2D Mesh Tutorial using GMSH
This tutorial was created to show how to generate a 2D mesh for OpenFOAM using the GMSH Open Source Mesh Generator. OpenFOAM by default only works with 3D mesh elements, so some special steps need to be applied to create a 2D mesh. This is not meant to be a tutorial on GMSH or OpenFOAM, but just some useful steps to get the two tools to work together. Note that two ways of geometry generation are shown. One uses the GMSH GUI to create the geometry and requires some modifications of the resulting geometry (.geo) file. The other one directly creates the geometry file. Meshing might be done with or without GUI in both ways.
Contents
1 Defining the Geometry and Mesh in GMSH
1.1 Graphical User Interface (GUI) and *.geo-modification
1. Open GMSH and create a new file.
2. In a single plane (2D), create the geometry by first creating all points, then combining the points into lines, and then the lines into a surface. These can all be done from the Geometry Menu under Elementary entities->Add->New.
3. Extrude the final surface into 3D by selecting Elementary entities->Extrude->Translate under the Geometry menu. Select Surface, then click the surface in the viewer. This will extend the surface into 3D space. The distance and direction of the extension is defined under the Contextual Geometry Definitions window that appears, under the Translate tab.
4. Now that the shape is 3D, the boundaries can be defined. To define a boundary, go to Physical Groups->Add->Surface from Geometry and add each surface that will be a boundary. Remember the order of the boundaries, as the GMSH GUI will label them with numbers that later need to be changed to a text name for OpenFOAM. For example, if the shape is a square, select the bottom, left, top, right, front, and back faces and individually add them as a Physical Group (If your boundary consists of multiple faces, you can select them together and then add the Physical Group). Later the names will be changed from a number to “bottom”, “left”, “top”, etc... These names will be used in OpenFOAM.
5. The volume of the geometry must also be named. Select Physical Groups->Add->Volume from Geometry and select the volume.
6. The shape should be saved as a .geo file automatically, but if not, save it now. Then open the .geo file in a text editor. Rename the Physical Surfaces to the names that will be used in OpenFOAM as the boundary names. You can also assign internal surfaces to physical surfaces to make internal walls in the geometry. For example (the numbers may be different):
Physical Surface("front") = {28}; Physical Surface("back") = {6}; Physical Surface("bottom") = {27}; Physical Surface("left") = {15}; Physical Surface("top") = {19}; Physical Surface("right") = {23};
Do the same for the volume:
Physical Volume("internal") = {1};
7. In the .geo file, the Extrude statement must be modified to have a single layer, and then should be recombined. Do this by adding the Layers{1} and Recombine commands to the Extrude command. This ensures that the mesh blocks are created only in 2D and then extended to 3D, and not divided in the 3rd axis. For example:
Extrude {0, 0, 10} { Surface{6}; Layers{1}; Recombine; }
8. Save the .geo file in the text editor and then reload in GMSH using Reload under the Geometry menu. There should be no errors produced by GMSH.
Below is an example of a 2D mesh .geo file of a simple square:
Point(1) = {-100, 100, 0, 1e+22}; Point(2) = {100, 100, 0, 1e+22}; Point(3) = {100, -100, 0, 1e+22}; Point(4) = {-100, -100, 0, 1e+22}; Line(1) = {1, 2}; Line(2) = {2, 3}; Line(3) = {3, 4}; Line(4) = {4, 1}; Line Loop(6) = {4, 1, 2, 3}; Plane Surface(6) = {6}; Physical Volume("internal") = {1}; Extrude {0, 0, 10} { Surface{6}; Layers{1}; Recombine; } Physical Surface("front") = {28}; Physical Surface("back") = {6}; Physical Surface("bottom") = {27}; Physical Surface("left") = {15}; Physical Surface("top") = {19}; Physical Surface("right") = {23};
9. Create the mesh in GMSH by selecting 3D from the Mesh menu. You can also do it from the command line:
gmsh -3 <fileName>.geo
which will automatically produce a fileName.msh file in the same directory. If everything is set up correctly (via the Extrude modifications), the sides of the 3D mesh should only have straight lines connecting the front face to the back face. This is required for OpenFOAM because for 2D analysis the mesh blocks should only change in two dimensions. OpenFOAM will produce an error if the mesh does not meet this requirement.
The size of the mesh elements can be controlled by setting the Element Size Factor under Mesh in the Options window (select Tools->Options from the main menu). The smaller the number, the smaller the mesh elements. For the changes in element size factor to take place, just go back to Geometry, Reload, then back to Mesh and click 3D.
10. Save the mesh by selecting ‘Save’ from the Mesh menu. This will create a file using the same name as the .geo file, but ending in .msh. This will be used by the OpenFOAM utility gmshToFoam that is shipped with OpenFOAM.
1.2 Creating *.geo from scratch
Instead of creating the geometry via GUI, you could also create the *.geo file from scratch in any editor.
1. In an editor of your choice, create a new file with a .geo file extension.
2. Define the 2D-geometry using the syntax below.
Point(<pointNr>) = {<x>, <y>, <z>, <cellSize>};
Line(<lineNr>) = {<startPointNr>, <endPointNr>};
Line Loop(<lineLoopNr>) = {(+/-)<lineNr1>, (+/-)<lineNr2>, ...};
Plane Surface(<surfaceNr>) = {<lineLoopNr>};
Be aware that the Line Loop function needs sorted input. The end- and start-points of adjacent lines have to match and the direction (e.g. clockwise) must be consistent. The direction of lines may be reversed with a minus.
Note that the order of appearance of lines in Line Loop will be used again after extrusion.
3. Extrude the surface using
surfaceVector[] = Extrude {<x>, <y>, <z>} {
Surface{<surfaceNr>};
Layers{1}; // create only one layer of elements in the direction of extrusion
Recombine;}; // recombine triangular mesh to quadrangular mesh
where 'surfaceVector[]' is a vector containing
- surfaceVector[0] = surfaceNumber of opposing Plane
- surfaceVector[1] = volumeNumber of extruded Volume
- surfaceVector[...] = surfaceNumber of surrounding surfaces in the same order as in Line Loop
after extrusion.
4. Define physicals for later use as boundaries using the syntax below
Physical Surface("<boundaryName>") = {<surfaceNr>}
Physical Surface("<boundaryName>") = surfaceVector[<integerAsShownAbove>];
Physical Volume("<interiorName>") = surfaceVector[1];
5. Define variables for your parametric mesh
squareSide = 200; //m
meshThickness = squareSide / 10;
gridsize = squareSide / 20;
6. The *.geo file should now look similar to this example.
// Inputs
squareSide = 200; //m
meshThickness = squareSide / 10;
gridsize = squareSide / 20;
// All numbering counterclockwise from bottom-left corner
Point(1) = {-squareSide/2, -squareSide/2, 0, gridsize};
Point(2) = {squareSide/2, -squareSide/2, 0, gridsize};
Point(3) = {squareSide/2, squareSide/2, 0, gridsize};
Point(4) = {-squareSide/2, squareSide/2, 0, gridsize};
Line(1) = {1, 2}; // bottom line
Line(2) = {2, 3}; // right line
Line(3) = {3, 4}; // top line
Line(4) = {4, 1}; // left line
Line Loop(5) = {1, 2, 3, 4};
// the order of lines in Line Loop is used again in surfaceVector[]
Plane Surface(6) = {5};
surfaceVector[] = Extrude {0, 0, meshThickness} {
Surface{6};
Layers{1};
Recombine;
};
/* surfaceVector contains in the following order:
[0] - front surface (opposed to source surface)
[1] - extruded volume
[2] - bottom surface (belonging to 1st line in "Line Loop (6)")
[3] - right surface (belonging to 2nd line in "Line Loop (6)")
[4] - top surface (belonging to 3rd line in "Line Loop (6)")
[5] - left surface (belonging to 4th line in "Line Loop (6)") */
Physical Surface("front") = surfaceVector[0];
Physical Volume("internal") = surfaceVector[1];
Physical Surface("bottom") = surfaceVector[2];
Physical Surface("right") = surfaceVector[3];
Physical Surface("top") = surfaceVector[4];
Physical Surface("left") = surfaceVector[5];
Physical Surface("back") = {6}; // from Plane Surface (6) ...
/* On my PC, the file must end with a free line to avoid errors which might come from different control characters in UNIX, Mac and Windows! Just to be save, insert one line below this comment!*/
7. Mesh the *.geo-File using either the GUI or a command line as shown in the previous chapter.
1.3 Optional: Structured Mesh
There are several ways to obtain a structured mesh in GMSH. One easy way for simple geometries is adding
Transfinite Surface{<surfaceNr>}={<edgePointsNr1>, <edgePointsNr2>, ...};
// forces later meshing to contain structured triangles
// e.g. Transfinite Surface{6} = {1,2,3,4};
Recombine Surface{<surfaceNr>};
//combine triangles to quadrangles
// e.g.Recombine Surface{6};
just before the extrusion command in the *.geo-file.
Example:
// Inputs
squareSide = 200; //m
meshThickness = squareSide / 10;
gridsize = squareSide / 10;
// Geometry
Point(1) = {-squareSide/2, -squareSide/2, 0, gridsize};
Point(2) = {squareSide/2, -squareSide/2, 0, gridsize};
Point(3) = {squareSide/2, squareSide/2, 0, gridsize};
Point(4) = {-squareSide/2, squareSide/2, 0, gridsize};
Line(1) = {1, 2}; // bottom line
Line(2) = {2, 3}; // right line
Line(3) = {3, 4}; // top line
Line(4) = {4, 1}; // left line
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};
//Transfinite surface:
Transfinite Surface {6};
Recombine Surface {6};
surfaceVector[] = Extrude {0, 0, meshThickness} {
Surface{6};
Layers{1};
Recombine;
};
surfaceVector contains in the following order:
[0] - front surface (opposed to source surface)
[1] - extruded volume
[2] - bottom surface (belonging to 1st line in "Line Loop (6)")
[3] - right surface (belonging to 2nd line in "Line Loop (6)")
[4] - top surface (belonging to 3rd line in "Line Loop (6)")
[5] - left surface (belonging to 4th line in "Line Loop (6)")
Physical Surface("front") = surfaceVector[0];
Physical Volume("internal") = surfaceVector[1];
Physical Surface("bottom") = surfaceVector[2];
Physical Surface("right") = surfaceVector[3];
Physical Surface("top") = surfaceVector[4];
Physical Surface("left") = surfaceVector[5];
Physical Surface("back") = {6};
Result:
2 Converting the Mesh to OpenFOAM
1. An OpenFOAM case directory with a controlDict file must already be created. Copy the .msh file into the case directory.
2. Open a command console and go to the case directory, and enter gmshToFoam <file.msh>. For example:
gmshToFoam square.msh
This will create the mesh used by OpenFOAM under the constant/polyMesh directory. If there are any errors, ensure that a proper case directory is set up (refer to OpenFOAM documentation, or use a pre-shipped case directory such as icoFoam or simpleFoam). If the errors indicate no 3D elements exist, or there is another mesh problem, make sure all steps were followed in the GMSH section of this tutorial.
3. Check that the mesh conforms to OpenFOAM standards by running checkMesh from the console.
4. To confirm that the boundaries made it to the OpenFOAM mesh, open the boundary file in the case directory under constant/polyMesh. Each boundary is defined here. For the faces that are along the 3rd dimension, set the type to ‘empty’. You can also use the Contrib modifyPatches to automate that. For example:
back { type empty; nFaces 2858; startFace 4217; } front { type empty; nFaces 2858; startFace 7075; }
This will indicate to OpenFOAM that this is a 2D case. Ensure that all other boundaries as defined in the GMSH .geo file as physical surfaces are present. There could also be a defaultFaces boundary that should have nFaces = 0 (it is unclear what will happen if nFaces is not 0):
defaultFaces { type patch; nFaces 0; startFace 10073; }
Alternatively, you can edit the constant/polyMesh/boundary dictionary with an awk command on the shell, or in a shell script.
$ awk -v RS="\0" -v ORS="" '{ gsub(/(patch name 1)\n {\n type patch/, "(patch name 1)\n {\n type (geometric patch type)"); print }' constant/polyMesh/boundary > tmp $ mv tmp constant/polyMesh/boundary $ awk -v RS="\0" -v ORS="" '{ gsub(/(patch name 2)\n {\n type patch/, "(patch name 2)\n {\n type (geometric patch type)"); print }' constant/polyMesh/boundary > tmp $ mv tmp constant/polyMesh/boundary ... $ awk -v RS="\0" -v ORS="" '{ gsub(/(patch name n)\n {\n type patch/, "(patch name n)\n {\n type (geometric patch type)"); print }' constant/polyMesh/boundary > tmp $ mv tmp constant/polyMesh/boundary
Notice that the number of spaces or tabs (\t sign) matter for gsub(input, output) to find what has to be changes. You could also name all your boundaries of the same kind with a suffix (e.g. surf1_wall, surf2_wall, surf3_empty...), and awk will substitute all parts where this pattern (e.g. _wall, _empty, _symmetry...) is found.
5. To generate an internal wall you need to remove the zero sized boundary from "constant/polyMesh/boundary" and add two zero sized patches at the end of the file named internalWallMaster and internalWallSlave and update the total number of patches (you can use Contrib modifyPatches to automate that) then use splitMesh to cut the mesh at the internalWall. If you have more than one internal wall use splitMeshWithSets. More details can be found in importing a fluent mesh with internal walls
6. The 2D mesh is now ready to be used by OpenFOAM. Be sure that the boundaries names in the initial conditions files (for example p, U, etc...) match those if the boundary file. For the empty faces (3rd dimension), set the boundary type to ‘empty’ in these files.
7. If all went well, OpenFOAM should be able to use the new mesh for a 2D case. If things don’t go well... OpenFOAM will let you know.
3 Mesh Examples
3.1 2D plane-symmetric stenosis
/*
Profile of the axisymmetric stenosis following a cosine
function dependent on the axial coordinate x [1].
f(x) = R * (1 - f0/2 *(1 + Cos(2.0*Pi * (x-x0)/L) )
-- x0 maximum stenosis position.
-- L stenosis length.
-- R vessel radius
-- f0 obstruction fraction, ranging 0--1.
References:
[1] Varghese, S. S., Frankel, S. H., Fischer, P. F., 'Direct numerical simulation of steotic flows. Part 1. Steady flow', Journal of Fluid Mechanics, vol. 582, pp. 253 - 280.
*/
ls = 5;
Xi = 100; // um
Xo = 100; // um
L = 100.0; // um
x0 = Xi + L/2.0;
R = 50.0; // um
f0 = 0.5; // 0--1
Z = 5;
Point(1) = {0, 0, 0, ls};
Point(2) = {Xi, 0, 0, ls};
Point(3) = {Xi, R, 0, ls};
Point(4) = {0, R, 0, ls};
Point(5) = {Xi + L, 0, 0, ls};
Point(6) = {Xi + L + Xo, 0, 0, ls};
Point(7) = {Xi + L + Xo, R, 0, ls};
Point(8) = {Xi + L, R, 0, ls};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {2, 5};
pList[0] = 3; // First point label
nPoints = 21; // Number of discretization points (top-right point of the inlet region)
For i In {1 : nPoints}
x = Xi + L*i/(nPoints + 1);
pList[i] = newp;
Point(pList[i]) = {x,
( R * (1 - f0/2 *(1 + Cos(2.0*Pi * (x-x0)/L) ) )),
0,
ls};
EndFor
pList[nPoints+1] = 8; // Last point label (top-left point of the outlet region)
Spline(newl) = pList[];
Transfinite Line {9, 10} = Ceil(L/ls) Using Progression 1;
Transfinite Line {4, -2, 8, -6} = Ceil(R/ls) Using Progression 1.1;
Transfinite Line {1, 3} = Ceil(Xi/ls) Using Progression 1;
Transfinite Line {5, 7} = Ceil(Xo/ls) Using Progression 1;
Line Loop(11) = {4, 1, 2, 3};
Plane Surface(12) = {11};
Line Loop(13) = {2, 10, 8, -9};
Plane Surface(14) = {13};
Line Loop(15) = {8, 5, 6, 7};
Plane Surface(16) = {15};
Transfinite Surface {14,12,16};
Recombine Surface {14,12,16};
Extrude {0,0,Z} {
Surface{14,12,16}; Layers{1}; Recombine;
}
Coherence;
Physical Surface("symmetryLine") = {51, 37, 73};
Physical Surface("frontAndBack") = {60, 38, 82, 16, 14, 12};
Physical Surface("wall") = {59, 29, 81};
Physical Surface("inlet") = {47};
Physical Surface("outlet") = {77};
Physical Volume("volume") = {2, 1, 3};
