PyGmsh - Undesired Behavior with Transfinite Mesh

Question

I am trying to design a rather simple mesh with PyGmsh. I want to have two rectangles, both with transfinite (rectangular) meshes; they will share node points in the x-direction, but I want to be able to independently set the number of mesh points in the y-direction for each rectangle. See "Desired" mesh image below. I am able to produce the desired mesh using a .geo file - shown below, however, I am not able to produce the desired mesh using PyGmsh. Particularly with PyGmsh, without using recombine, the mesh produced is very close to the desired outcome, but I would like to recombine the triangles into rectangles. However, when the recombine command is used, the mesh warps at the interface of the two rectangles.

I believe the issue has to do with the fact that I cannot pair opposed transfinite curves in PyGmsh, which can be done directly in the .geo file (Transfinite Curve {1, 5} = 11 Using Progression 1;)

Is there a way to produce the desired mesh using PyGmsh? Looking through the documentation, I could not find a corresponding method to Transfinite Curve {1, 5} = 11 Using Progression 1;. It appears PyGmsh only allows individual lines to be set as Transfinite Curves, not pairs of lines.

import sys

# gmsh opens an external window which needs to be closed to finish the running code
import gmsh
factory = gmsh.model.geo    # alias to facilitate code writing
meshFact = gmsh.model.mesh  # meshing alias



gmsh.initialize()

# default mesh size (not necessary, since we are using transfinite curves
# and setting a certain number of points in all curves)
lc = 1.


# Geometry
height = 1.0
length = 2.0
thickness = 0.5

# points - clockwise
# Bottom Layer
p1 = factory.addPoint(0., 0., 0., lc)
p2 = factory.addPoint(0., height, 0., lc)
p3 = factory.addPoint(length, height, 0., lc)
p4 = factory.addPoint(length, 0., 0., lc)

# Top Layer
p5 = factory.addPoint(0., height + thickness, 0., lc)
p6 = factory.addPoint(length, height + thickness, 0., lc)


# lines - clockwise
l1 = factory.addLine(p1, p2)    # Inlet
l2 = factory.addLine(p2, p5)    # Top Layer West Wall
l3 = factory.addLine(p5, p6)    # Top Layer North Wall
l4 = factory.addLine(p6, p3)    # Top Layer East Wall
l5 = factory.addLine(p3, p4)    # Outlet
l6 = factory.addLine(p4, p1)    # bottom wall

# curve loops
cl1 = factory.addCurveLoop([l1, l2, l3, l4, l5, l6])

# surfaces
s1 = factory.addPlaneSurface([cl1])

gmsh.model.geo.mesh.setTransfiniteCurve(l3, 21)
gmsh.model.geo.mesh.setTransfiniteCurve(-l6, 21)

gmsh.model.geo.mesh.setTransfiniteCurve(l1, 11)
gmsh.model.geo.mesh.setTransfiniteCurve(-l5, 11)

gmsh.model.geo.mesh.setTransfiniteCurve(l2, 11)
gmsh.model.geo.mesh.setTransfiniteCurve(-l4, 11)

ts = gmsh.model.geo.mesh.setTransfiniteSurface(s1, "Left", [p1, p5, p6, p4])

factory.synchronize()

# mesh
# perform 2D mesh
meshFact.generate(2)

# Perhaps another method to do mesh recombination (?)
# gmsh.model.geo.mesh.setRecombine(2, tag=1, angle=45)      # dim = 2, tag = 1

# Toggling this line ON / OFF produces the attached figures
# meshFact.recombine()    # combines bisecting triangles - in order to get quadrangular 2D elements

gmsh.write("Rectangle_NoRecombine.msh")
gmsh.write("Rectangle_NoRecombine.geo_unrolled")

# Create graphical user interface
if 'close' not in sys.argv:
    gmsh.fltk.run()

# It finalizes the Gmsh API
gmsh.finalize()

Desired output can be achieved with the following .geo file

cl__1 = 1;

Point(1) = {0, 0, 0, cl__1};
Point(2) = {0, 1, 0, cl__1};
Point(3) = {2, 1, 0, cl__1};
Point(4) = {2, 0, 0, cl__1};
Point(5) = {0, 1.5, 0, cl__1};
Point(6) = {2, 1.5, 0, cl__1};

Line(1) = {1, 2};
Line(2) = {2, 5};
Line(3) = {5, 6};
Line(4) = {6, 3};
Line(5) = {3, 4};
Line(6) = {4, 1};

Curve Loop(1) = {1, 2, 3, 4, 5, 6};

Plane Surface(1) = {1};

Transfinite Curve {1, 5} = 11 Using Progression 1;
Transfinite Curve {2, 4} = 21 Using Progression 1;
Transfinite Curve {3, 6} = 21 Using Progression 1;

//Transfinite surface:
Transfinite Surface {1} = {1, 5, 6, 4};
Recombine Surface {1};

Images of the desired mesh versus the meshes produced via PyGmsh (red / green added to help with visualization).

Answer 1

The t11 tutorial found on github does outline some methods for using recombine. A section of that tutorial lists this comment section:

# Note that you could also apply the recombination algorithm and/or the
# subdivision step explicitly after meshing, as follows:
# gmsh.model.mesh.generate(2)# gmsh.model.mesh.recombine()
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
# gmsh.model.mesh.refine()

However, applying this method to the problem at hand does not produce the desired result. Modifying the code above to the following, does produce the desired result:

ts = gmsh.model.geo.mesh.setTransfiniteSurface(s1, "Left", [p1, p5, p6, p4])

factory.synchronize()

# mesh
# perform 2D mesh
# gmsh.model.geo.mesh.setRecombine(2, s1)       # This does not produce the desired result
gmsh.model.mesh.setRecombine(2, s1)             # dim = 2, tag = s1
# meshFact.recombine()                          # This method does not produce the desired result
meshFact.generate(2)
# meshFact.recombine()                          # This method

gmsh.write("Rectangle_NoRecombine.msh")
gmsh.write("Rectangle_NoRecombine.geo_unrolled")

This modification also changes the .geo_unrolled file, where the line Recombine Surface {1}; is appended to the end of the file now.

I am not sure if this is the designed behavior of the code, but this does produce the desired mesh.

Answer 2

If you're fine to use another package for your problem you could use FElupe (spoiler: I'm the author). You could then convert the FElupe mesh to a meshio-mesh object and write it to a gmsh file. More details about FElupe meshes can be found in the docs. FElupe is available on PyPI: pip install felupe[all].

import felupe as fem

height = 1.0
length = 2.0
thickness = 0.5

rectangles = [
    fem.Rectangle(a=(0, 0), b=(length, height), n=(21, 11)),
    fem.Rectangle(a=(0, height), b=(length, height + thickness), n=(21, 11)),
]

container = fem.MeshContainer(rectangles, merge=True)
mesh = fem.mesh.stack(container.meshes)
mesh.as_meshio().write("mesh.msh")
mesh

A summary of the metadata of the mesh:

<felupe Mesh object>
  Number of points: 441
  Number of cells:
    quad: 400

The mesh may be visualized with PyVista / matplotlib.

ax = mesh.imshow()