[pradeepramesh]

Dear All,

How to access the Cell numbers in OpenFOAM framework?

I am working on a benchmark case, where i am interested to find the cell numbers. I have to apply "Influx boundary condition" at the right wall (Bold-Red-Line). Hence i m very keen to access the cell numbers of the described domain.

Please find the attached snapshot for a quick reference. Would be great to hear your views on this.

Thank you very much

Pradeep

[jhoepken]

Why do you need to access the cell ids (I think that this is what you are after), from the boundary condition?

[pradeepramesh]

I am trying to find a way to access the mesh information of the domain:
(Mesh - Cells - Faces - Edges - Points - Coordinates )

I am interested to assign source (Influx) on the cells at the boundary. Hence i am keen to know the Cell IDs.

Thank you Jens

[pradeepramesh]

One more aspect is to verify the cell IDs that are extracted from ParaView (SpreadSheet View) with the Cell IDs from the OpenFOAM framework.

[jhoepken]

So you need the coordinates of the vertices on the boundary?

[pradeepramesh]

Yes, coordinate of the vertices at the boundary.

As i mentioned previously, I am keen to apply a source term on a cell (ID xyz) in the domain.

Thank you Jens

[chrisb2244]

Firstly, boundary faces have the highest values for the faces. You can iterate through them using something like

labelList boundaryFaceList(mesh_.nFaces() - mesh_.nInternalFaces());
forAll(boundaryFaceList, faceI)
{
    // Do whatever. For example, to find the cell ID, use
    cellI = mesh_.faceOwner()[faceI + mesh_.nInternalFaces()];
}

This will only get you the list of all faces, however. If you know the physical boundary values of your mesh, you could compare using mesh_.cellCentres()[cellI] or similar, or if you don't know the size of your mesh (as in, for example, 0 left wall -> 10 cm on right wall), you can calculate the normals of the faces using something a little bit like:

Foam::vector Foam::hexRef4::calcSingleFaceNormal(const pointField& p, const label& faceI) const
{
    const labelList& f = mesh_.faces()[faceI];
    point fCentre = p[f[0]];
 
    vector sumN = vector::zero;
    scalar sumA = 0.0;
 
    label nPoints = f.size();
 
 
    for (label pi = 1; pi < nPoints; pi++)
    {
        fCentre += p[f[pi]];
    }
    fCentre /= nPoints;
 
    for (label pi = 0; pi < nPoints; pi++)
    {
        const point& nextPoint = p[f[(pi + 1) % nPoints]];
        vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
        scalar a = mag(n);
        sumN += n;
        sumA += a;
    }
    // This is to deal with zero-area faces. Mark very small faces
    // to be detected in e.g., processorPolyPatch.
    if (sumA < ROOTVSMALL)
    {
        return vector::zero;
    }
    else
    {
        return 0.5*sumN;
    }
}

The section above takes a pointField and a face index as arguments. faceI is pretty obvious, but for a pointField you want something like mesh_.points().


In all of these examples/sections, mesh_ is a reference to your polyMesh object. See http://www.openfoam.org/docs/cpp/ and search for 'polyMesh' (case sensitive). If it doesn't load a sensible list, refresh the page - that usually helps me.


Once you have the face normal, you can compare it to the direction you expect to be on your right wall. Should work fine. No idea how 'efficient' it is though, so consider that if you're conscious of such things.


(p.s. - you will need to change the class that that function belongs to when you place it where you need it - it is largely copied from an existing OF function, but if I recall correctly, it's from one of the applications, not a library, so copy-pasting and editing as needed into your own code, whilst causing duplication, is probably best. At least, it's what seemed easiest to me)

[chrisb2244]

For bonus points, ignore the long section I pasted and instead use

mesh_.faceAreas()[faceI + mesh_.nInternalFaces()];
// or
mesh_.faceAreas()[faceI];
// as appropriate

since that seems to have a lot of the same code. Can't remember why I didn't just use that, although I suspect it relates to being able to pass a different pointField. (For sub-sections, and new parts of the mesh being generated)

[pradeepramesh]

Thank you very much Christian