Dear All,
I have been writing a python script to produce a 3D mesh for a pipe with a number of alternating bends and straight sections. The unstructured meshing algorithms built into gmsh produce only triangles and tetrahedra. I wish to create a pipe mesh that has hexahedral elements regularly spaced: through the thickness, around the circumference, and along the length with a mesh cross-section like that below:
{see Fig1 in comments}
I concluded that the best way to do this, using Gmsh, would be using transfinite curves/surfaces. The issue I ran into was that in order to create transfinite surfaces for structured meshing, it is required that surfaces have a minimum of 3 corner tags. An annulus obviously does not have this. This was preventing me from forming the desired mesh for a pipe with bends.
My workaround for this issue is shown in the MWE below. It involves splitting the desired pipe in half so that all surfaces have 4 corners so that they can be structurally meshed as transfinite surfaces forming the desired 3D mesh of hexahedra. The goal was to then combine the two halves again, as the meshes on the surfaces that join them are identical but I can’t work out how to fuse them without destroying the structured mesh. The .msh file produced by this MWE is shown below:
{Fig2}
Although this MWE creates a .msh file that looks really nice, the two volumes are separate and accordingly this mesh cannot be used for FE using FEniCS. This is shown to be true in the screenshot included below of an attempt to excite waves in this pipe mesh, the MWE has clearly formed two separate volumes.
{see Fig3 in comments}
I tried using “gmsh.model.occ.fuse([(3, 1)], [(3, 2)])” to get around this issue but then it turns the annulus ends of the pipe back into a single surface which cannot be meshed using transfinite curves because it has only 2 corners so I am back to where I started.
I have spent a long time trying to solve this problem in a large number of different ways but it always goes back to the underlying issue; that I am unable to mesh an annulus in a structured way to produce the desired mesh cross-section without introducing extra boundaries.
The piece of code below is the MWE that produces a pipe model and mesh with two straight sections and a “bend”. This mesh is hexahedral and has mesh elements of the desired shape but is two separate volumes.
import gmsh
import numpy as np
import sys
import meshio
# Initialisation
file_suffix = "fenics_forum_MWE"
r1 = 0.05
r2 = 0.04
Nthickness = 3
Ncirc = 128
Nlen = 300
gmsh.initialize()
gmsh.model.add(f"pipe_creation_{file_suffix}")
# Make spline
p = []
for m in range(10):
gmsh.model.occ.addPoint(0, 0, m/10, 1, tag=100+m)
p.append(100+m)
for m in range(50):
gmsh.model.occ.addPoint(0, 0.1*np.cos(2*np.pi*m/50)-0.1, 1+(m/50), 1, tag=200+m)
p.append(200+m)
for m in range(10):
gmsh.model.occ.addPoint(0, 0, 2+(m/10), 1, tag=300+m)
p.append(300+m)
spline_tag = gmsh.model.occ.addSpline(p, 900)
gmsh.model.occ.addWire([900], 900)
#####################################################################################################################################################################
# We define the shapes we would like to extrude along the spline (two half annulus):
#####################################################################################################################################################################
point1_tag = gmsh.model.occ.addPoint(r1, 0, 0)
point2_tag = gmsh.model.occ.addPoint(r2, 0, 0)
point3_tag = gmsh.model.occ.addPoint(-r1, 0, 0)
point4_tag = gmsh.model.occ.addPoint(-r2, 0, 0)
line12_tag = gmsh.model.occ.addLine(point1_tag, point2_tag)
line43_tag = gmsh.model.occ.addLine(point4_tag, point3_tag)
outer_pipe_tag_0pi = gmsh.model.occ.addCircle(0, 0, 0, r1, angle1=0, angle2=np.pi, tag=1000)
inner_pipe_tag_0pi = gmsh.model.occ.addCircle(0, 0, 0, r2, angle1=0, angle2=np.pi)
curveloop1_tag = gmsh.model.occ.addCurveLoop([line12_tag, inner_pipe_tag_0pi, line43_tag, outer_pipe_tag_0pi], tag=1050)
pipe_cs_1 = gmsh.model.occ.addPlaneSurface([curveloop1_tag], tag=1)
gmsh.model.occ.addPipe([(2, pipe_cs_1)], spline_tag, 'DiscreteTrihedron')
gmsh.model.occ.remove([(2, pipe_cs_1)])
gmsh.model.occ.synchronize()
outer_pipe_tag_pi2pi = gmsh.model.occ.addCircle(0, 0, 0, r1, angle1=np.pi, angle2=2*np.pi, tag=2000)
inner_pipe_tag_pi2pi = gmsh.model.occ.addCircle(0, 0, 0, r2, angle1=np.pi, angle2=2*np.pi)
curveloop2_tag = gmsh.model.occ.addCurveLoop([line12_tag, inner_pipe_tag_pi2pi, line43_tag, outer_pipe_tag_pi2pi], tag=2050)
pipe_cs_2 = gmsh.model.occ.addPlaneSurface([curveloop2_tag], tag=11)
gmsh.model.occ.addPipe([(2, pipe_cs_2)], spline_tag, 'DiscreteTrihedron')
gmsh.model.occ.remove([(2, pipe_cs_2)])
gmsh.model.occ.synchronize()
#####################################################################################################################################################################
# delete all non-necessary entities
entities = gmsh.model.getEntities()
print(entities)
entities_to_remove = []
for e in entities:
if e[0] in range(3):
entities_to_remove.append(e)
gmsh.model.removeEntities(entities_to_remove)
gmsh.model.occ.remove(entities_to_remove)
gmsh.model.occ.synchronize()
# make the structured mesh
for c in [(1, 1005), (1, 1007), (1, 1011), (1, 1015), (1, 2005), (1, 2006), (1, 2011), (1, 2015)]:
gmsh.model.mesh.setTransfiniteCurve(c[1], Nthickness+1)
for c in [(1, 1006), (1, 1008), (1, 1016), (1, 1013), (1, 2007), (1, 2008), (1, 2016), (1, 2013)]:
gmsh.model.mesh.setTransfiniteCurve(c[1], (Ncirc//2)+1)
for c in [(1, 1009), (1, 1010), (1, 1012), (1, 1014), (1, 2009), (1, 2010), (1, 2012), (1, 2014)]:
gmsh.model.mesh.setTransfiniteCurve(c[1], Nlen)
for s in [(2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7)]:
gmsh.model.mesh.setTransfiniteSurface(s[1])
gmsh.model.mesh.setRecombine(s[0], s[1])
gmsh.model.mesh.setSmoothing(s[0], s[1], 100)
gmsh.model.mesh.setTransfiniteVolume(1)
for s in [(2, 12), (2, 13), (2, 14), (2, 15), (2, 16), (2, 17)]:
gmsh.model.mesh.setTransfiniteSurface(s[1])
gmsh.model.mesh.setRecombine(s[0], s[1])
gmsh.model.mesh.setSmoothing(s[0], s[1], 100)
gmsh.model.mesh.setTransfiniteVolume(2)
gmsh.model.addPhysicalGroup(2, [2, 3, 4, 5, 6, 7, 12, 13, 14, 15, 16, 17], 1, "pipe surface") # 2, 9, 7, 14 are ends. 4, 6, 11, 13 are curved surfaces
gmsh.model.addPhysicalGroup(3, [1, 2], 2, "pipe volume")
gmsh.model.mesh.generate(3)
gmsh.write(f"pipe_creation_{file_suffix}.msh")
# Launch the GUI to see the results:
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.finalize()
msh = meshio.read(f"pipe_creation_{file_suffix}.msh")
print(msh.cell_data_dict)
quad_data = msh.cell_data_dict["gmsh:physical"]["quad"]
hexahedron_data = msh.cell_data_dict["gmsh:physical"]["hexahedron"]
quad_cells = msh.cells_dict["quad"]
hexahedron_cells = msh.cells_dict["hexahedron"]
hexahedron_mesh = meshio.Mesh(points=msh.points, cells={"hexahedron": hexahedron_cells}, cell_data={"name_to_read":[hexahedron_data]})
quad_mesh =meshio.Mesh(points=msh.points, cells=[("quad", quad_cells)], cell_data={"name_to_read":[quad_data]})
meshio.write(f"hex_{file_suffix}.xdmf", hexahedron_mesh)
meshio.write(f"quad_{file_suffix}.xdmf", quad_mesh)
Please advise on how I can produce a 3D mesh of a bent pipe, as produced using an extrusion of a cross-section along a spline, where the mesh is formed of regular/structured hexahedral elements as shown in the images included!
Thank you very much!! 
