When developing your own code in OpenFOAM, it is almost certain that access to mesh information will be required in order to evaluate various parameters. In this post, I have provided a short description on how to access some of this information which could come in handy while programming. This post is comprised of the following topics:


Mesh connectivities

Cells

const labelListList& cellPoints = mesh.cellPoints();  // Cell to node
const labelListList& cellEdges = mesh.cellEdges();    // Cell to edge
const cellList& cells = mesh.cells();                 // Cell to face
const labelListList& cellCells = mesh.cellCells();    // Cell to cell

Faces

const faceList& faces = mesh.faces();                   // Face to node
const labelListList& faceEdges = mesh.faceEdges();      // Face to edge
const labelList& faceOwner = mesh.faceOwner();          // Face to owner cell
const labelList& faceNeighbour = mesh.faceNeighbour();  // Face to neighbour cell 

Nodes

const labelListList& pointPoints = mesh.pointPoints();  // Node to node 
const labelListList& pointEdges = mesh.pointEdges();    // Point to edge 
const labelListList& pointFaces = mesh.pointFaces();    // Point to face 
const labelListList& pointCells = mesh.pointCells();    // Point to cell 

Edges

const edgeList& edges = mesh.edges();                   // Edge to node 
const labelListList& edgeFaces = mesh.edgeFaces();      // Edge to face 
const labelListList& edgeCells = mesh.edgeCells();      // Edge to cell

Mesh coordinates

Cell center, face center and node

const volVectorField& C = mesh.C();         // Cell center coordinates
const surfaceVectorField& Cf = mesh.Cf(); // Face center coordinates
const pointField& points = mesh.points(); // Node coordinates

Edge center

Since OpenFOAM is based on the conventional cell-centred Finite Volume Method, edge coordinates are not computed by the fvMesh class. However, if necessary they can be obtained very easily by averaging the 2 node coordinates associated to each edge.

// Include vectorList definition to store edge coordinates
#include "vectorList.H"

// Create edge list which contains edge to node connectivity
const edgeList& edges = mesh.edges();
   
// Create a vectorList to store edge centre coordinates 
// This list has size of edgeList 'edges' and is initialized to zero 
vectorList Ce(edges.size(), vector::zero);             

// Loop over all edges
forAll(edges, edge)
{
   const label& own = edges[edge][0];     // Index of edge owner node
   const label& nei = edges[edge][1];     // Index of edge neighbour node
  
   // Calculate edge center coordinates by averaging node coordinates
   Ce[edge] = 0.5*(points[own] + points[nei]);
}

Boundary mesh data

// Store boundary mesh information
const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();

Distinguish between internal and boundary faces

// Loop over all faces (internal and boundary)
forAll(faces, face)
{
  if (mesh.isInternalFace(face)) // Internal face is found
  {
    // Do your calculations for internal faces i.e.
// U[face] = vector::zero;
} else // Boundary face is found { const label& patch = boundaryMesh.whichPatch(face); // Boundary patch index const label& facei = boundaryMesh[patch].whichFace(face); // Local boundary face index
// Do your calculations for boundary faces i.e.
// U.boundaryField()[patch][facei] = vector::zero;

}
}

Boundary patches

// Loop over boundary patches
forAll(mesh.boundary(), patch)
{
  const word& patchName = mesh.boundary()[patch].name();            // Boundary patch name
    
  // Loop over all faces of boundary patch
  forAll(mesh.boundary()[patch], facei)
  {
    const label& bCell = boundaryMesh[patch].faceCells()[facei];    // Boundary cell index
    const label& face = boundaryMesh[patch].start() + facei;        // Face index
        
// Do your calculations e.g. // U.boundaryField()[patch][facei] = vector::zero;
}
// Loop over all nodes of boundary patch forAll(boundaryMesh[patch].meshPoints(), pointi) {
const label& point = boundaryMesh[patch].meshPoints()[pointi]; // Node index
// Do your calculations e.g. // U[point] = vector::one;
}
}

Other useful mesh parameters

// Cell volumes
const scalarField& V = mesh.V();

// Face area normal vectors const surfaceVectorField& Sf = mesh.Sf();
// Face areas const surfaceScalarField& magSf = mesh.magSf();
// Face normals const surfaceVectorField& N = Sf/magSf;
// Total number of cells const label& nCells = mesh.nCells();
// Total number of nodes const label& nPoints = mesh.nPoints();
// Total number of internal faces const label& nInternalFaces = mesh.nInternalFaces();
// Total number of internal nodes const label& nInternalPoints = mesh.nInternalPoints();

OpenFOAM classes

For more information the following classes in OpenFOAM Foundation release version 7  can be referred:

Similar Posts

24 Comments

    1. Tobi, the aim of this post is to provide a brief description of the mesh information that I believe one is most likely to use while programming in OpenFOAM. I know that everything is there in the code but it is not always easy for novice OpenFOAM programmers to dig the relevant information. Nevertheless, based on your suggestion, I have added the link to the relevant class.

  1. Hello,
    This seems to be useful for my problem. I need to extract the face normal of a patch in Openfoam. Where can I put the codes you provided above in my folder case? Thanks

  2. Thanks for the concise summary. I was programming a solver that needs to check neighbor cells for a certain criteria, and found that cellCells doesn’t work in parallel if you try to address across processor boundaries, that is, the cells across decomposition boundaries aren’t included in cellCells() output. Hopefully using some of the other stuff in here will get this to work! 😀

    1. You are welcome and apologies for the delayed response! You are right because cellCells() is only aware of the cells who belong to that particular processor. If you want to access the complete mesh information then you could refer to globalMeshData class. It might have something useful for your needs.

  3. Thanks a lot for this!
    I am developing a solver which will generate boundary condition for solving Laplacian using an already calculated field. I calculated the boundary condition but I am struggling on how to use it for solving the Laplacian. Can you suggest something?

  4. Thanks a lot for your post. This was very helpful.
    Could you please help me calculate the face areas of cells in a boundary.
    scalar localArea = mesh.magSf().boundaryFieldRef()[patchi][facei];
    I have tried as above. But, am unsuccessful. [patchi] in my boundary patch.

  5. Hello,

    I found this incredibly helpful!. Thank you for putting the effort and clarifying not only the syntax but the detailed implementation. I am struggling with the idea of cell labels. It seems that all fields are saved at the center of any given cell, in the case of finite volume methods, this makes sense. However, in the PolyMesh folder there is no file such as cells. However, there is owner and neighbor information. How can I construct a connectivity diagram myself just to improve my understanding of creating meshes?.

    Thank you!

    1. Hi, I am glad that you found the post helpful. OpenFOAM is based on the cell-centred Finite Volume Method where the primary unknown variables are certainly stored at the cell centres. Then there are additional variables such as fluxes and boundary values which are generally stored at the face centres. In some cases you might have other variables which could be stored at the nodes.

      When it comes to mesh generation the information you generally require are the nodal coordinates and element connectivities. Therefore, cell centre coordinates are derived from the nodal coordinates by OpenFOAM and not part of the ‘polyMesh’ folder.

      Now the question is what are these files named ‘owner’ and ‘neighbour’? These essentially give you the connectivity of faces to cells. So every internal face will have an owner cell and a neighbour cell (generally speaking), however, boundary faces will only have one cell (owner cell) since it is at the boundary. So to cut it short, the size of owner file will correspond to the total number of faces in the mesh and each face will have an owner cell associated to it. On the other hand, the size of neighbour file will be equal to the internal faces of the mesh and each internal face will have a neighbour cell. I hope that it is clear now!

  6. Hi Jibran,

    First of all – thank you very much for this post, this helped me A LOT since I am just getting started with programming in OpenFOAM. I have some questions that I hope you could help me clarify. But first you should know – I am building my own functionObject to be able to output some desired data in a file in the end of the simulation and I used the command foamNewFunctionObject to get it going. So to my questions:

    1. I had to change all your “mesh.” statements to “mesh_” why is that? (I saw someone using “mesh_” so I tried it out and it worked)
    2. If I want to find a collection of all extensions that the boundaryMesh() has (e.g. .faceCells(), .start() etc,) where can I find that?
    3. You are writing U.boundaryField()[patch][facei] … I was trying to write that in the console because what I was expecting was that I would get an output showing me the velocity value for each face for each time step. However it gave me an error saying that “U was not declared in this scope”. How would I change it to make it work?

    I am super thankful for any help on this.

    Kind regards,
    David

  7. Hi Haidar,

    Thanks for the post. I have a few questions for you.

    1. Look at the line below from your post:
      const edgeList& edges = mesh.edges();
      What is the first “edges”? Is it only a name you are giving to this data type? Can one use something else. I see this pattern in all of the data you created in the post.
      I do understand that the “edges” in the end of the line is a member function.
    2. I am currently studying one code in OpenFoam to better understand how I can write my own application.
      In one section of this sample application, there is this line:
      // Get the fvPatch
      const fvPatch& patch = U.boundaryField()[patchID].patch();

      I know it is trying to pull some mesh info out for a specific patch to use later, but my question is that what is the type of the data which is output from the above line? and how I can view that data for that particular patch?
  8. Hi Haidar,

    How to loop over cells in parallel run.
    I want to give ‘if’ condition for some cells. Please help me.

    Regards,
    Bibin

    1. Hi Bibin,
      OpenFOAM code is executed in exactly the same manner on all processors unless you explicitly specify it to do so. Therefore, if you loop over cells in you code, it will loop over cells in all processors separately. Now depending on what you want to implement in the if condition, additional stuff might be necessary (e.g. MPI reduce operations etc.).

      1. Hi Haider,

        I loop over all the cells inside my time loop and replace the field (e.g. T) in the cell of interest with the value I want. It works fine with the static mesh. But it doesn’t work with the adaptive mesh. Could you please tell me why and how can I make it work for the adaptive mesh? I have written the piece of code that works with the static mesh below.

        const pointField& ctrs = mesh.C();
        forAll(ctrs, cellI)
        {
        scalar xc= magSqr(x0.value() – ctrs[cellI]);
        if (xc <= radius)
        {
        T[cellI] = 2000;
        }

        }

  9. Thank you very much for your post above. It has been very useful for me. In my experiments with OpenFoam, I have a pointVectorField defined. Here in the following piece of code I am trying to access the boundary field of the defined pointVectorField. But I keep getting compilation errors. It would be very helpful for me to know how to access and modify the boundary field.
    auto& mappedPatchField =
    mappedField.boundaryFieldRef()[patchId_];
    mappedPatchField[i_cell] = ....

    here mappedField is the pointVectorField and I want to access the field on the nodes of patch with patchId_ and modify it. When I try to compile it I get the following error message
    no match for ‘operator[]’ (operand types are ‘Foam::pointPatchField<Foam::Vector<double> >’ and ‘Foam::label {aka int}’)

    Thank you very much in advance.

  10. Hi Haidar,

    I’ve been trying to apply a nonuniform profile for my dynamic mesh with velocityComponentLaplacian but for some reason I can’t access the mesh. When I tried to construct the mesh point coordinates I got ‘mesh’ not declared in this scope, please let me know if I missed something and thanks for your time.

    Best regards,
    Howard

Leave a Reply

Your email address will not be published. Required fields are marked *