How to create a grid mesh using pygmsh?

[amine-aboufirass]

I am using the following code to generate a simple grid mesh, also using fipy:

from fipy import CellVariable, Gmsh2D
import pygmsh, pyvista

geom = pygmsh.built_in.Geometry()
p1 = geom.add_point((0.,0.,0.), lcar = 1.0)
p2 = geom.add_point((2.,0.,0.), lcar = 1.0)
p3 = geom.add_point((0.,2.,0.), lcar = 1.0)
p4 = geom.add_point((2.,2.,0.), lcar = 1.0)
l6 = geom.add_line(p1, p2)
l7 = geom.add_line(p2, p4)
l8 = geom.add_line(p4, p3)
l9 = geom.add_line(p3, p1)

ll10 = geom.add_line_loop((l6, l7, l8, l9))
s11 = geom.add_surface(ll10)

geom_string = geom.get_code()
mesh = Gmsh2D(geom_string)

phi = CellVariable(name = "solution variable", mesh = mesh, value = 0.)

ugrid= pyvista.UnstructuredGrid(mesh.VTKCellDataSet._vtk_obj)
plotter = pyvista.Plotter()
plotter.set_background('white')
plotter.add_mesh(ugrid, style='wireframe', color='black')
plotter.add_bounding_box(color='red')
plotter.show_grid(color="red")
plotter.view_xy()
plotter.show()

However the resulting plot shows triangular elements and not grid elements. Is there anything I can do at the pygmsh level in order to enforce grid elements, or is this something that can only be solved at the fipy level. In other words where does the decision for which kind of elements occur?

[amine-aboufirass]

Nvm, I think I found where the mesh is created, it is line 218 in gmshMesh.py from the fipy source code. It seems to call the following command:

gmsh geo_file.geo -2 -nopopup -format msh2 -o msh_file.msh

However I'm still not seeing any indication of element type. Does gmsh only do triangular elements? Thanks.

[nschloe]

What do you mean by "grid elements" ?

[amine-aboufirass]

rectangular elements...

[nschloe]

Well, by default gmsh does triangular cells. Just googling "gmsh quad" gives you many indications of how to create quad meshes, e.g., https://github.com/cycheung/gmsh/blob/master/tutorial/t11.geo.

[amine-aboufirass]

Using this video and not without some seriously annoying glitches on the part of the gmesh GUI (it is SO difficult to select entities!) I was able to generate this which will result in a rectangular mesh:

// Gmsh project created on Wed Feb 12 14:36:38 2020
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {2, 0, 0, 1.0};
//+
Point(3) = {0, 1, 0, 1.0};
//+
Point(4) = {2, 1, 0, 1.0};
//+
Line(1) = {2, 4};
//+
Line(2) = {4, 3};
//+
Line(3) = {3, 1};
//+
Line(4) = {1, 2};
//+
Curve Loop(1) = {3, 4, 1, 2};
//+
Plane Surface(1) = {1};
//+
Transfinite Curve {3, 1} = 10 Using Progression 1;
//+
Transfinite Curve {2, 4} = 20 Using Progression 1;
//+
Transfinite Surface {1};
//+
Recombine Surface {1};

Then I tried to adapt my code accordingly:

from fipy import CellVariable, Gmsh2D
import pygmsh, pyvista

geom = pygmsh.built_in.Geometry()
p1 = geom.add_point((0.,0.,0.), lcar = 1.0)
p2 = geom.add_point((2.,0.,0.), lcar = 1.0)
p3 = geom.add_point((0.,2.,0.), lcar = 1.0)
p4 = geom.add_point((2.,2.,0.), lcar = 1.0)
l6 = geom.add_line(p1, p2)
l7 = geom.add_line(p2, p4)
l8 = geom.add_line(p4, p3)
l9 = geom.add_line(p3, p1)
ll10 = geom.add_line_loop((l6, l7, l8, l9))
s11 = geom.add_surface(ll10)
geom.set_transfinite_lines([l6, l7, l8, l9], 6)
geom.set_transfinite_surface(s11)

geom_string = geom.get_code()
mesh = Gmsh2D(geom_string)

phi = CellVariable(name = "solution variable", mesh = mesh, value = 0.)

ugrid= pyvista.UnstructuredGrid(mesh.VTKCellDataSet._vtk_obj)
plotter = pyvista.Plotter()
plotter.set_background('white')
plotter.add_mesh(ugrid, style='wireframe', color='black')
plotter.add_bounding_box(color='red')
plotter.show_grid(color="red")
plotter.view_xy()
plotter.show()

This was partially successful however I still need to "Recombine Surface" as shown in the .geo script. Anything in pygmsh that can allow me to do this? A quick perusal of available methods in the geom object did not reveal anything regarding Recombining Surfaces... Thanks.

[nschloe]

There's always add_raw_code().

[amine-aboufirass]

Right. Is there a more pythonic way to do it rather than adding the raw code?

I see this pull request which might be interesting/relevant: https://github.com/nschloe/pygmsh/pull/120

[nschloe]

Right now, no. Feel free to PR.

[amine-aboufirass]

I could try. Would it be very difficult to implement? General guidelines/instructions?

[nschloe]

It should perhaps go into https://github.com/nschloe/pygmsh/blob/master/pygmsh/built_in/geometry.py, recombine(surface), then assert if surface is of SurfaceBase type and add the code with its ID.

[gdmcbain]

Another way to create a grid mesh is to use pygmsh.Geometry.extrude, passing it recombine=True; see

https://github.com/nschloe/pygmsh/blob/f667b73e599fbbf6c67bbcc1150f11e50000b533/test/test_hex.py#L14-L22

[gdmcbain]

Besides invoking Gmsh's Recombine with add_raw_code as in

https://github.com/nschloe/pygmsh/blob/f667b73e599fbbf6c67bbcc1150f11e50000b533/test/test_quads.py#L10

there's another little catch that you might like to know about, if you expect a mesh consisting purely of quadrilaterals with no admixture of triangles. This is Gmsh's Mesh.RecombinationAlgorithm. This can give rise to a nasty heisenbug; see, e.g., kinnala/scikit-fem#275. Some values of Mesh.RecombinationAlgorithm only try to use quadrilaterals but don't guarantee it. Worse, some values are only available in more recent versions of Gmsh and Gmsh doesn't complain if it's passed an unknown value but just silently degrades to using one of the algorithms it knows which might not be one that guarantees no triangles! To be robust against various versions of Gmsh, I'm using

    geom.add_raw_code('Mesh.RecombinationAlgorithm=2;\n')

but if you know that you have a more recent one, I believe 3 is better.

I don't think this applies to extrusions with recombine=True; I think they always produce quadrilaterals (or indeed hexahedra in three dimensions).

[amine-aboufirass]

@nschloe I have just submitted a pull request for this. It seems to work nicely, though I do see that most checks are failing. What needs to be done about those?

@gdmcbain I am working with simple 2D geometries for the moment, so no extrusions. I'm sure your comment will come in useful for more complex situations. I just haven't faced these yet though.

[amine-aboufirass]

@nschloe I fixed the formatting problems but the pytest tests are still failing. This is strange because when I run python -m pytest in my local machine things seem to run fine... If anyone could take a quick look that would be great. Thanks.

[nschloe]

Fixed via https://github.com/nschloe/pygmsh/pull/305.