[Madeleine P. Vincent]

All -

I have a question regarding gmsh and the gmshToFoam utility within OpenFOAM, one of which seems to be not acting as I had presumed.

I have a simple mesh problem which I've distilled down to the it's basic form.
I am running a gmsh script (.geo) to create a 3D rectangular flow field with a circular cylinder in the middle.
When I run gmsh, it meshes succesfully. However, when I run gmshToFoam (conversion to the OpenFOAM polyMesh format), I end up with a mesh that is not usable.

First of all, I get a warning saying:


[warning]
--> FOAM Warning :
From function polyMesh:olyMesh(... construct from shapes...)
in file meshes/polyMesh/polyMeshFromShapeMesh.C at line 619
Found 9022 undefined faces in mesh; adding to default patch.
[/warning]



Invoking the "checkMesh" utility results in the following errors:


[error]
Checking geometry...
Overall domain bounding box (0 0 -0.75) (9 6 0.75)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (1.15651e-17 5.07766e-19 9.032e-18) OK.
***High aspect ratio cells found, Max aspect ratio: 7.58943e+198, number of cells 3
<<Writing 3 cells with high aspect ratio to set highAspectRatioCells
Minumum face area = 6.81219e-06. Maximum face area = 0.177405. Face area magnitudes OK.
Min volume = 1.33333e-300. Max volume = 0.0208743. Total volume = 80.8949. Cell volumes OK.
Mesh non-orthogonality Max: 89.2144 average: 30.4407
*Number of severely non-orthogonal faces: 1046.
Non-orthogonality check OK.
<<Writing 1046 non-orthogonal faces to set nonOrthoFaces
***Error in face pyramids: 8 faces are incorrectly oriented.
<<Writing 8 faces with incorrect orientation to set wrongOrientedFaces
Max skewness = 2.23263 OK.
Coupled point location match (average 0) OK.

Failed 2 mesh checks.
[/error]



The contents of the gmsh script (.geo file) is as follows:


[script]
// Everything is parameterized on the cylinder radius
rad = 0.15;

//Characteristic length - used for mesh sizing
lc = 4*rad;

//Plane dimesions
planeLength = 60*rad;
planeWidth = 40*rad;

//Cyclinder constants
cir_x = 15*rad; //The absolute X distance of the center
cir_y = planeWidth / 2; //The absolute Y distance of the center

//The bottom plane number
bottomPlaneNum = 16;

//The extrusion thickness
extThickness = rad*10;
zBottom = -(extThickness/2.0);
zTop = (extThickness/2.0);
z = zBottom; //The z-height of this plane.
// NOTE - For 2D, the plane is centered at 0, so
// the z value used will be -(h/2)

//Rectangular plane (Counter-clockwise, starting in lower left corner)
Point(1) = {0, 0, z, lc};
Point(2) = {planeLength, 0, z, lc};
Point(3) = {planeLength, planeWidth, z, lc};
Point(4) = {0, planeWidth, z, lc};

//Circle
Point(5) = {cir_x, cir_y, z, lc}; // Center point
Point(6) = {cir_x, cir_y + rad, z, lc}; // 12 o'clock
Point(7) = {cir_x + rad, cir_y, z, lc}; // 3
Point(8) = {cir_x, cir_y - rad, z, lc}; // 6
Point(9) = {cir_x - rad, cir_y, z, lc}; // 9

//Lines
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Circle(5) = {6, 5, 7};
Circle(6) = {7, 5, 8};
Circle(7) = {8, 5, 9};
Circle(8) = {9, 5, 6};
Line Loop(9) = {1, 2, 3, 4};
Line Loop(10) = {5, 6, 7, 8};

//Surface
// Don't add any line loop(s) specific to meshing
Plane Surface(bottomPlaneNum) = {9, 10};

//Extrude the bottom surface in the z direction a total of extThickness units.
// This extrusion is centered on 0 (from -h/2 to +h/2)
newSurfs[] = Extrude {0,0,extThickness} {
Surface{bottomPlaneNum};
};

//Name the surfaces
Physical Surface("topAndBottom") = {bottomPlaneNum, newSurfs[0]};
Physical Surface("nearAndFar") = {newSurfs[2], newSurfs[4]};
Physical Surface("outlet") = {newSurfs[3]};
Physical Surface("inlet") = {newSurfs[5]};
Physical Surface("cylinder") = {newSurfs[9], newSurfs[8], newSurfs[7], newSurfs[6]};

Physical Volume(1) = {newSurfs[1]};

//Meshing
fineness = 15;

Field[4] = Cylinder;
Field[4].Radius = 1.5*rad;
Field[4].VIn = lc / fineness;
Field[4].VOut = lc / 2;
Field[4].XAxis = 0.0;
Field[4].XCenter = cir_x;
Field[4].YAxis = 0.0;
Field[4].YCenter = cir_y;
Field[4].ZAxis = 1.0;
Field[4].ZCenter = 0.0;

Field[6] = Box;
Field[6].VIn = lc / (fineness / 3);
Field[6].VOut = lc;
Field[6].XMin = cir_x - 3*rad;
Field[6].XMax = cir_x + 10*rad;
Field[6].YMin = cir_y - 3*rad;
Field[6].YMax = cir_y + 3*rad;
Field[6].ZMin = zBottom -0.000001;
Field[6].ZMax = zTop +0.000001;

// Use minimum of all the fields as the background field
Field[100] = Min;
Field[100].FieldsList = {4, 6};
Background Field = 100;

[/script]



I have no problem with the 2D case -- I use Layers{1} and Recombine within the extrusion command. This works well.

I have also seen that using many layers within the 3D case will also produce a mesh that passes the checkMesh utility. But this doesn't provide a good 3D mesh, as the domain is segregated throughout in these distinct layers, and it doesn't take into account the finer meshing (in the Z direction) around the cylinder.

So I'm wondering if this is an issue in the meshing of gmsh, or the gmshToFoam utility? And what remedy might be available.

Thank you very much for any help you are able to provide.

Madeleine P. Vincent