[Voulet]

Hi there. I would like to get the faceIDs of some cells whose neighbour cell and owner cell have a specific scalar value (volscalarfield)

• I loop throw my cell (with the help of my 'box' voldscalarfield).

• if my cell scalar has a certain value:
  • i get the 6 faces of my cell
  • for all these faces:
    • if the adjacant cell has a scalar of a certain value :
      • i append the faceid to a list

My issue is that the faceid i'm receving are made of 0, 1, 2, 3, 4, 5 : ie they are the relative faces of my cell.
How can i get the absolute faceId in order to use them after ?


Thanks a lot.


Here is my code :

    DynamicList<label> Omega2FaceList;         // Facelist

    // loop over the cell of my mesh
    forAll(box, celli) {
        // if my scalar is equal to one
        if (box[celli] = 1.0){
            // get facesid of the cell
            const cell& faces = mesh.cells()[celli];
            // Info << "cell faces : " << faces << endl;

            forAll(faces,facei) {
                const scalar& sOwner = box[mesh.owner()[facei]];
                const scalar& sNeighbour = box[mesh.neighbour()[facei]];
                // Info << "cell owner value : " << box[mesh.owner()[facei]] << endl;
                // Info << "cell neighbour value: " << box[mesh.neighbour()[facei]] << endl;
                if ((sOwner == 1.0) && (sNeighbour == 0.0)) {
                    Omega2FaceList.append(facei);
                    Info << "facei: " << facei << endl;

                }
            }
        }
    }

[Voulet]

No one to help me ?

[dlahaye]

Not sure here.

However, loops over faces (instead of over cells) in OpenFoam are more intuitive (given the face-based structure of OpenFoam).

Each face has one owner and one neighbor cell.

Thus you could loop over all faces and check if the field value in the owner and neighbor cell has target value. If so, store value of the given face.

The following was adapted from Listing 6.10 in Chapter 6 of the book by Moukalled e.a.

DynamicList<label> Omega2FaceList; // Facelist

const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();

forAll (own,facei)
{
if ( (box[own[facei]]==1.) && (box[nei[facei]]==1.) )
{
Omega2FaceList.append(facei);
}
}

Hope this gives ideas .

[olesen]

Originally Posted by Voulet

    DynamicList<label> Omega2FaceList;         // Facelist

    // loop over the cell of my mesh
    forAll(box, celli) {
        // if my scalar is equal to one
        if (box[celli] = 1.0){
            // get facesid of the cell
            const cell& faces = mesh.cells()[celli];
            // Info << "cell faces : " << faces << endl;

            forAll(faces,facei) {
                const scalar& sOwner = box[mesh.owner()[facei]];
                const scalar& sNeighbour = box[mesh.neighbour()[facei]];
                // Info << "cell owner value : " << box[mesh.owner()[facei]] << endl;
                // Info << "cell neighbour value: " << box[mesh.neighbour()[facei]] << endl;
                if ((sOwner == 1.0) && (sNeighbour == 0.0)) {
                    Omega2FaceList.append(facei);
                    Info << "facei: " << facei << endl;

                }
            }
        }
     }

I have no idea what it is about the forAll macro that misleads so many people. It is simply a for-loop over the size of some list-like item. You are reusing the facei, which is an index, as if it is the face ID!


Without questioning if there is a better way, here is a cleaner version of your code. Please examine the changes to see your error.

DynamicList<label> Omega2FaceList;         // Facelist

// loop over the cell of my mesh
forAll(box, celli)
{
    if (!equal(box[celli], 1))
     {
        continue;
    }

     for (const label facei : mesh.cells()[celli])
    {
        if 
        (
            equal(box[mesh.owner()[facei]], 1)
         &&
            (
                !mesh.isInternalFace(facei)
              || equal(box[mesh.neighbour()[facei]], 0)
         )
         {
             Omega2FaceList.append(facei);
         }
      }
    }
}

I've just made an assumption about what you want to do on boundaries.


However, since you are essentially just looking for faces to select, it makes much more sense to be checking them directly.

DynamicList<label> Omega2FaceList;         // Facelist
for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
{    const label own = mesh.owner()[facei];
      const label nei = mesh.neighbour()[facei];

    if (equal(box[own], 1) && equal(box[nei], 0))
    {
        Omega2FaceList.append(facei);
    }
}

for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
 {
    const label own = mesh.owner()[facei];


   if (equal(box[own], 1))
     {
        Omega2FaceList.append(facei);
    }
}

BTW: in both cases I have used the equal() method when checking values, to account for possible rounding.



Depending on the density of faces to be selected and what you want to do with the information afterwards, you may find a bitSet a convenient container. For example,

bitSet selectFaces(mesh.nFaces());


 for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
{
     const label own = mesh.owner()[facei];
     const label nei = mesh.neighbour()[facei];

    if (equal(box[own], 1) && equal(box[nei], 0))
    {
        selectFaces.set(facei);
    }
}


... later on


Info<< "found " << selectFaces.count() << " faces" << nl;


if (selectFaces.test(1000)) Info<< "face 1000 is selected" << nl;


for (const label facei : selectFaces) ...




// Pick the other faces

selectFaces.flip();

[Voulet]

Thanks a lot. Your answer help me to understand more deeply the openfoam mesh structuration.

 DynamicList<label> Omega2FaceList;         // Facelist
for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
{    const label own = mesh.owner()[facei];
      const label nei = mesh.neighbour()[facei];

    if (equal(box[own], 1) && equal(box[nei], 0))
    {
        Omega2FaceList.append(facei);
    }
}

for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
 {
    const label own = mesh.owner()[facei];


   if (equal(box[own], 1))
     {
        Omega2FaceList.append(facei);
    }
}

I don't understand the second loop. I guess that is intended to get the cells on the boundary whose owner cell have a box == 1 value. That's what i'm looking for. I do undertand that the faces indexing start with internal faces and then swith to boundary faces (that's why for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)).


However i get a list out of bound. It seems that the mesh.owner.()size() is the same than the mesh.nInternalFaces() and that the boundary cells are not present in the mesh.owner. How is it possible ? Looking at my constant/polymesh/owner, the size of the list should is however equal to mesh.nFaces (nearly twice as mesh.nInternalFaces).



Looking at the documentation, line 278 https://cpp.openfoam.org/v8/fvMesh_8H_source.html that the owner list only concern the internal faces.


Am I wrongly construcing my mesh ? I only #include fvMesh.H to do it.

[Voulet]

Yes, it seems that the mesh.owner() list is only okay for the internal faces. So how can i know to which cell is connected a face whose id correspond to a boundary face ?