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
- Mesh coordinates
- Boundary mesh data
- Other useful mesh parameters
- OpenFOAM classes for further reference
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:
- PrimitiveMesh class lists all functions
- FvMesh class contains mesh information for Finite Volume discretisation
- PolyMesh class
- PolyBoundaryMesh class
Nice but you could directly refer to doxygen and thats it 🙂
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.
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
Hi Emma. Generally you would code this in your solver file but it depends on what you want to evaluate.
Are there terminal commands that are equivalent to check? This would prevent from having to compile.
Hello, some of this information is already stored in the ‘polymesh’ folder of the ‘constant’ directory. There might be some OpenFOAM utilities as well.
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! 😀
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.
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?
You are welcome and sorry for the delayed response. Not sure what exactly your problem is. Perhaps you can send me a detailed email to discuss this.
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.
You are welcome! Try this without Ref.
scalar localArea = mesh.magSf().boundaryField()[patchi][facei];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!
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!
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:
I am super thankful for any help on this.
Kind regards,
David
Hi David,
You are welcome. I will try to answer your questions here:
Hope this helps!
Hi Haidar,
Thanks for the post. I have a few questions for you.
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.
In one section of this sample application, there is this line:
// Get the fvPatchconst 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?
Hello Mohammad. Apologies for the late response.
Hi Haidar,
How to loop over cells in parallel run.
I want to give ‘if’ condition for some cells. Please help me.
Regards,
Bibin
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.).
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;
}
}
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.
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
Actually please ignore my question, I realised I forgot to construct the polyMesh.
Best regards,
Howard
Hello,
I am building my own finite volume solver but I want to use the mesh data from OpenFOAM.
I make the .unv mesh in Salome and then convert it to an Foam mesh using the command ideasUnvToFoam.
Now, I would like to extract the cell-to-node connectivity matrix. How do I do that?
Where do I need to input the following statement?
const labelListList& cellPoints = mesh.cellPoints(); // Cell to node
I tried in controlDict but it does not work.
Best,
Alessandro
Hi Alessandro,
You need to use them in your custom OpenFOAM solver not in the input OpenFOAM files!
Thanks for the reply.
But I do not want to solve the equations, I only need to do use the above commands as a preprocessing step. How do I write the subroutine in such a way that mesh data is extracted without running the simulation?
Alessandro
You can write a short OpenFOAM utility where you create your mesh and dump the required connectivity to a file.
Then you can use it in your own FVM solver (I am presuming it is not in OpenFOAM).
Have a look at some OpenFOAM utilities here: https://www.openfoam.com/documentation/user-guide/standard-utilities.php
Hi Haidar
Thanks for the post.
I want to access volume of cells in volScalarField but in openFoam there is two ways to access to the mesh volume and none of them are volScalarField:
const DimensionedField& V() const;
or
const scalarField& cellVolumes() const;
Hi Mojtaba,
Glad that the post helped you. May I ask why would you want to store volume in a volScalarField?
In OpenFOAM volFields have an internal field and a boundary field. Since cell volumes are only associated to the cells themselves they are not associated to a volField.
Depending on your need, you could always work around and create a volScalarField whose internalField() is filled with cellVolumes().
Thank you very much for your reply.
I want to add a source term in my continuity equation. In this source term, there is a constant called “specific area” which has a value of non-zero on the boundary and zero value in other locations. My problem is 2D, and the value of the specific area is defined as the ratio of the area of the face, connected to the boundary on a cell boundary, to the volume of the boundary cell. I was wondering how could I program this definition to calculate the specific area, and which type of specific area should I add to the continuity equation?
Thank you very much in advance!
I think what you need is mentioned in the Boundary patches part of this post.
You can loop over boundary patches and then faces of those patches.
You will have access to the global face index of this patch face and also the boundary cell index.
The code snippet is repeated below for convenience.
// 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;
} // End loop over faces of boundary patch
} // End loop over boundary patches
The post is great, man!
I have a boundary patch with a polygon shape. But with the fvPatch instance, I could not find a way to obtain the edges of this polygon. Could you please kindly share your thoughts on this issue? Thanks a lot!
Hi,
Thanks for the post!
A question, should the the surface normal vector be a reference?
Shouldn’t it be:
const surfaceVectorField nf = Sf/magSf;
Hi,
Thank you for your post !! I have a query to ask you.
My domain is divided into N elements (I= 1 to N). In my work, I need to duplicate the node as (IB, IC) for a purpose. Let’s say,
for(I = 1 to N)
{
if (T[I] == 100)
IB = I;
}
for(I=IB+1 to N)
{
if (T[I] < 100)
IC = I;
}
(1) How to declare the for loop specifically in OpenFOAM?
(2) How to address IB and IC ? Like if its scalar variable we can declare "scalar IB" but that's not gonna work here.
Kindly someone share the ideas.
Thank you
Hello!
My question is regarding ‘multiRegion’ in OpenFOAM. Suppose I have two regions – solid and liquid.
The two regions share a common boundary, known as the coupled patch. From being in the solid region or vice versa, through which command can I find out the patchID of the coupled patch?