mercredi 25 décembre 2013

Building hexagonal meshes with Gmsh

OpenFOAM's blockMesh is rather simple and efficient way of building meshes but it has certain annoying features:

Entities numbering. To construct a block, an edge or a boundary one has to know its vertices numbers, they are numbered automatically stating with 0. So while constructing a mesh one can, for example, use vim with +relativenumber and blockMeshDict opened in two splits (see screenshot). But when blockMesh reports about errors in the mesh description it also uses these automatic numbers of entities. Soon it becomes rather annoying to look for the block #4 or #17 to correct gradings or densities.

Mesh grading is a good thing but for example to make a cube mesh with a higher density near the walls, one has to cut the cube into 4 blocks and make consistent grading in each block.

Curved edges is something special. For a 2D case it's rather simple to make an utility to calculate coordinates of a point on the arc of a given radius passing through two vertices, in a 3D case it's less obvious.

And finally workflow vim -> blockMesh -> paraFoam is rather slow.

Gmsh usually used with the finite-element code GetDP so by default it generates tetrahedral meshes, like one below.
But it is possible to make hexagonal meshes with it. All GEO files from the sections below can be found in the repository.

Cube

To make a cube one has to create 8 vertices (with Gmsh's GUI, for example, but in fact it's simpler to create the points in an editor):

As the lengths of the cube sides are equal, densities of the 1D mesh on every side are equal (also I've used keyword Bump to grade the mesh near the corners). To use generated mesh in OpenFOAM one has to define physical groups also:

Tetrahedron

This mesh is less trivial than the previous one, the idea of a tetrahedron cut is taken from the article by David Eppstein. Let's assume tetrahedron with the following corners: (-5, 0, -5), (5, 0, -5) and (0, -5, 5), (0, 5, 5). To decompose the tetrahedron into 6-faced volumes one needs to know centers of the sides, centers of the faces and center of the tetrahedron:

After connection of the points with the lines, defining faces and volumes, we get the following geometry:
As the densities of the different parts of the mesh should be more or less the same one can add the same four lines as with the cube to have hexagonal mesh

and the final mesh will be like this
One of the characteristics of the mesh one should check before running a simulation is non-orthogonality (I omitted this step in case of cube as obviously for that case non-orthogonality is 0). Here is the report of the checkMesh utility for the generated mesh:

One can play with the points positions to get more orthogonal mesh but 44 is not so bad.

Cylinder

Trivial solution when one just take a cube and turn it into a cylinder by specifying certain sides as the arcs leads to the highly non-orthogonal meshes. Surely there're many ways to build mesh on a cylinder, below are several variants.

Version 1

Top view of the cylinder subdivision looks like
it is a square inside a circle. So we need 9 points (4 for the square, 4 on the circle and one for the center of the arcs) in each horizontal plane.

Alternatively we can start only with 9 points in one plane and then use extrude operation. After connection of the points with lines and arcs, defining planes and volumes cylinder geometry will look like
3 set of lines can have independent 1D mesh densities, so we define them

Version 2

The first idea that comes to mind in the quest for reducing non-orthogonality of the mesh is to use a square with the curved sides. Let them be arcs with a radius Rsq < Rcyl. Schematic top view of the cylinder division will look like
To have curved sides of the square we need to define additional points which will be the centers of the circle arcs

As usual we need to define all the lines connecting points, then planes, and finally volumes. Geometry will look like this
You can see those additional points which are the centers of the circle arcs of the "square" sides. Definition of the entities for transfinite meshing algorithm is almost the same as for the previous case. And here is the mesh
As the selection of the radius for the "square" arcs was more or less arbitrary we got a mesh with higher non-orthogonality:

Though the mesh is OK we'll make another mesh in the quest for the lower non-orthogonality.

Version 3

This time we'll convert the square into an octagon to gain control over additional points. Schematic top view of the cylinder subdivision will look like
So in addition to 18 points of the first subdivision scheme we get 16 more points

Spherical cap

Spherical cap meshing is a further development of the last variant of cylinder mesh. Schematic top and center-plane cut views of the geometry decomposition are shown below, top view is the same as for cylinder and center-plane cut is almost the same as a half of top view.
As the number of points in subdivision grows now it's easier to generate points programmatically than by hand. This is especially true for the points inside the cap and on "bottom" of the sphere. Still to connect points and to create the planes and the volumes one can use GUI (though it becomes less and less convenient). After these operations one get the following spherical cap decomposition into hexagonal volumes
In this case we have three sets of lines which can have different 1D mesh densities

And so after meshing spherical cap looks likes this
checkMesh report again indicates that we can improve our mesh, in general, it can be achieved by playing with points displacements in the top plane and especially by playing with the angles between octagonal cylindroid and points at the "bottom" of the cap.

And that's it. Basically dividing a shape into hexagonal (or pentagonal) volumes and using transfinite algorithm one can mesh anything with hexagons. Though when the number of entities in the geometry file increases usability of GUI goes down, it is still a little bit more convenient than blockMesh.

31 commentaires:

I just stumbled across your blog. Great post comparing blockMesh and gmsh.

I recently came across gmsh while going through the following blog http://lordvon64.blogspot.in/2012/03/hybrid-structuredunstructured-meshing.html?showComment=1431612709648#c4326128210914448246.

I have written a gmsh file to create a hybrid mesh for a thin flat plate. I have 9 surfaces in the 2D mesh. I'm not sure how to handle them during extrusion since they are creating various edges i.e more than the usual front,top,bottom,inlet,outlet and wing.

Can you help me get past this problem? Here is the github link to the geo file -----> https://github.com/pruthvi1991/gmsh/tree/master/pilot .

To capture the boundary layers I created a refined mesh near the wing surface. For the this I used the transfinite algorithm and I have 9 volumes in the mesh. These volumes intersect at internal faces. In OpenFOAM I don't know how to define these. When I left them undefined OpenFOAM named them as default faces.

defaultFaces{type patch;nFaces 284; startFace 562756;}

I don't know how to handle these internal faces. Do you know how to handle this problem?

Here are some pictures to make it clear. https://drive.google.com/open?id=0B0O_EM6xOrebfm5FZE1TYUhqZXNUVlF2d1RfUUdxWU5KenJlbmV4dUN3dXUySnUyTG42WEk&authuser=0

To capture the boundary layers I created a refined mesh near the wing surface. For the this I used the transfinite algorithm and I have 9 volumes in the mesh. These volumes intersect at internal faces. In OpenFOAM I don't know how to define these. When I left them undefined OpenFOAM named them as default faces.

defaultFaces{type patch;nFaces 284; startFace 562756;}

I don't know how to handle these internal faces. Do you know how to handle this problem?

Here are some pictures to make it clear. https://drive.google.com/open?id=0B0O_EM6xOrebfm5FZE1TYUhqZXNUVlF2d1RfUUdxWU5KenJlbmV4dUN3dXUySnUyTG42WEk&authuser=0

1. I have tried to convert msh (created from your geo-file), was not successful. First there was higher-order elements, then wrong oriented prisms and finally gmshToFoam aborted with sigSegv.2. I still do not understand your problem.3. I still do not understand what you are trying to achieve.4. You was not able to interest me enough to look for the errors in your geo-file.

ok thanks for the effort. I was able to convert without any such issues. I just updated the geo file by adding some physical surfaces and volumes but it shouldn't matter in the conversion.

"Inverting prism 147798" ------> Is this what you mean by wrong oriented prism?

I'm just trying to achieve a nice hybrid mesh for a pitching and plunging airfoil. If you have the time to look into this problem again, please do a git pull. There might have been some issues with the previous .geo file.

As mesh is generated using transfinite algorithm, 'Element size at point' really does not affect density of the mesh. Number of cells can be adjusted by modifying settings on lines 257-265 (https://bitbucket.org/mrklein/gmsh-meshing/src/db90a5048eb18c047eb98f2cf5b6c3eb283e45b2/cylinder-3.geo?at=default#cl-257). On lines 257-258 number of vertical cells is set (in fact number points on vertical lines, so it is number of cells + 1). Lines 259-260 control number of cells in radial direction in outer blocks. Finally lines 261-265 control subdivision of inner blocks and subdivision of outer block in azimuthal direction (and these values should be equal as lines are opposite of each other).

I do not know how to set exact locations of the points in transfinite algorithm. If you would like to have dense mesh in certain regions you can use "Using Progression" or "Using Bump" options. See http://www.geuz.org/gmsh/doc/texinfo/gmsh.html#Structured-grids.

You can use it like "gmsh file.geo -3 -o mesh.msh", or "gmsh - file.geo" (in this case you add "Mesh 3; Save "mesh.msh";" in the file.geo). See http://www.geuz.org/gmsh/doc/texinfo/gmsh.html#Non_002dinteractive-mode.

Hi Alexey, thanks for your tutorial! I have a gmsh script which creates a 2D airfoil cascade mesh. The boundary layer is meshed in hexagonal elements and everything else is tetrahedra. The only thing that I can't seem to get to work is mesh refinement of the trailing edge. Would it be possible for you to look at my script?

Thanks! Please see https://github.com/Quadrupole/airfoil2D. Basically my issue is now trying to get a good way to have natural refinement at the leading and trailing edges and a coarser mesh nearer the mid chord. I want a smooth transition between these two regions but the Bump and Progression methods of the Transfinite lines don't seem to offer this. It should be apparent from running the script. Let me know if there are any questions.

Actually I realised that I need Hyperbolic extrusion in the boundary layer - something like http://www.pointwise.com/DIY/Extruding-Structured-Domain-Out-of-Airfoil.shtml

Because I will always have the problem of the mesh becoming 'stretched' near the leading and trailing edge curvature. I am not sure if GMSH can handle that easily? I see that construct2d has this feature, but do you have any experience or knowledge of whether it can create a cascade geometry (where the there are two periodic boundaries above and below the airfoil)?

My experience with Construct2d is limited to: oh, it makes nice meshes, much better than ones I have tried to make with Gmsh ;) Since it was just for fun, I never went below default mesh generation settings.

The script is now updated to something that looks ok (at least for solving RANS). However I have some trouble importing to openFoam (the conversion gmshToFoam works fine) but when I run the simulation I get --> FOAM FATAL IO ERROR: size 12216 is not equal to the given value of 30962. I have added the openFoam case to the repo. I would greatly appreciate if you could have a look as I think the error is caused by improper number of cells being specified, possible related to the way I define my volume in GMSH.

Hi again Alexey - I have found that my mesh created as in the above example gives me incorrect results in OpenFoam and I believe it could be due to the mesh. When I tried the same case with a hyperbolic extruded boundary layer mesh I get the expected results. I am therefore wondering if you think it could be possible to mesh the blade in construct2D, export to GMSH and then "fill" in the region between the boundary layer hex grid and the periodic boundaries.

1. Generate airfoil mesh with Construct2D, it outputs mesh in p3d format. Since later you want to use it with OpenFOAM, general 3D mesh with 2 planes (see OOPT menu).2. Convert mesh to Gmsh format. Obviously, I used p3d2gmsh.py (https://github.com/mrklein/p3d2gmsh). For naca2315 Construct2D generated mesh with two ONE-TO-ONE boundaries, so they should be stitched later (names of the patches *-to-stitch-*).3. Generate mesh with periodic boundaries and a hole in the centre, where you put mesh generated with Construct2D (as a separate case).4. Use mergeMesh to merge meshes.5. Use stitchMesh to stitch extra-boundaries.

Unfortunately I have performed only steps 1 and 2 (and there are no problems), there could be certain difficulties during steps 3-5.

Hi Alexey,I have gone through your examples above and have found them really useful. Thank you very much for sharing them! I have got a problem when we specify the grades for bump since it doesn’t seem to be straightforward like progression. I have checked the source code and found the formula that does the calculation. However, when I did a simple test by providing values for the number of points, length as well as coef, the result doesn’t match what gmsh provided. I was wondering if you know how to calculate grades when using the bump in Gmsh. The formula for bump function in Gmsh is as follows; case 2: // Bump { double a; if(coef > 1.0) { a = -4. * sqrt(coef - 1.) * atan2(1., sqrt(coef - 1.)) / ((double)nbpt * length); } else { a = 2. * sqrt(1. - coef) * log(fabs((1. + 1. / sqrt(1. - coef)) / (1. - 1. / sqrt(1. - coef)))) / ((double)nbpt * length); } double b = -a * length * length / (4. * (coef - 1.)); val = d / (-a * SQU(t * length - (length) * 0.5) + b); }I would really appreciate any clarifications.

Thanks for the reply. The previous page that you sent is for progression in one side using the layer command. However, for bump, progression in both sides, I could not get the same node distributions using the following ad hoc function;

I would like to repeat: go to mailing list, pose question directly to guys, who develop Gmsh. Bump function, used in the software, is their invention, so either it is you, who reconstructs function using debugging, or you ask authors.