[Gmsh] gmsh or gmshToFoam problem?

[Madeleine Vincent]

gmsh friends -

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::polyMesh(... 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 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/gmsh/attachments/20120619/cd5b22cc/attachment.html>