I am trying to make a mesh consisting of 3 quarters of an annulus completed with a distorted quarter annulus. I have a working FEniCSx code with a working mesh consisting of a full annulus with the 2 Physical Surface regions. The working mesh can be built with the following file.geo:
// Gmsh project created on Sat Mar 8 22:37:56 2025
SetFactory("OpenCASCADE");
//+
Mesh.CharacteristicLengthMax = 0.0009;
Mesh.CharacteristicLengthMin = 0.0001;
Disk(1) = {0, 0, 0, 0.05, 0.05};
//+
Disk(2) = {0, 0, 0, 0.035, 0.035};
BooleanDifference{ Surface{1}; Delete;}{ Surface{2}; Delete;}
// Surface 3 is the full annulus. I need to split it into 2.
Point(4) = {0, 0.05, 0};
//+
Point(5) = {0, 0.035, 0};
//+
// Origin.
Point(6) = {0, 0, 0};
Line(4) = {2,3};
Line(5) = {4,5};
Circle(6) = {5, 6, 2};
Circle(7) = {4, 6, 3};
Curve Loop(8) = {5,6,4,7};
Plane Surface(9) = {8};
// 3-quarter surface
BooleanDifference{ Surface{1}; Delete;}{ Surface{9};}
BooleanFragments{ Surface{1}; Surface{9}; Delete; }{ }
Physical Surface("three_quarter_material", 20) = {1};
Physical Surface("quarter_material", 21) = {9};
Physical Curve("theta_0", 22) = {4};
Physical Curve("theta_pi_over_2", 23) = {5};
Physical Curve("r_inner_quarter", 24) = {6};
Physical Curve("r_outer_quarter", 25) = {7};
Physical Curve("r_inner_three_quarter", 26) = {8};
Physical Curve("r_outer_three_quarter", 27) = {9};
I can examine the annulus with gmsh and see that all Physical Curves and surfaces are fine.
Here is my Python code that produces the modified quarter-annulus msh file:
import gmsh
import numpy as np
from math import pi, sqrt
# -----------------------------------------------------------------------------
# User parameters
# -----------------------------------------------------------------------------
r_o, r_i = 0.05, 0.035
A_i, B, n_pt = 0.008, 0.003, 30
gmsh.initialize()
gmsh.model.add("annulus_full")
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.0009)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.0001)
# -----------------------------------------------------------------------------
# 1) Build 3/4‐annulus in OCC
# -----------------------------------------------------------------------------
c1 = gmsh.model.occ.addDisk(0, 0, 0, r_o, r_o)
c2 = gmsh.model.occ.addDisk(0, 0, 0, r_i, r_i)
ann, _ = gmsh.model.occ.cut([(2, c1)], [(2, c2)])
ann_tag = ann[0][1]
L = 2*r_o
p1 = gmsh.model.occ.addPoint(0, 0, 0)
p2 = gmsh.model.occ.addPoint(0, L, 0)
p3 = gmsh.model.occ.addPoint(L, L, 0)
p4 = gmsh.model.occ.addPoint(L, 0, 0)
l1 = gmsh.model.occ.addLine(p1, p2) # θ=π/2 radial
l2 = gmsh.model.occ.addLine(p2, p3)
l3 = gmsh.model.occ.addLine(p3, p4)
l4 = gmsh.model.occ.addLine(p4, p1) # θ=0 radial
wloop = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4])
wedge = gmsh.model.occ.addPlaneSurface([wloop])
three, _ = gmsh.model.occ.cut([(2, ann_tag)], [(2, wedge)])
three_tag = three[0][1]
gmsh.model.occ.synchronize()
# -----------------------------------------------------------------------------
# 2) Build quarter patch in OCC
# -----------------------------------------------------------------------------
k = r_i*A_i + (pi/8)*A_i**2
a,b,c= pi/8, r_o, -k
disc = b*b - 4*a*c
if disc < 0:
raise RuntimeError("No real solution for A_i")
A_o = (-b + sqrt(disc))/(2*a)
θ = np.linspace(0, 0.5*pi, n_pt)
r_out = r_o + A_o*np.sin(2*θ) + B*np.sin(4*θ)
r_in = r_i + A_i*np.sin(2*θ) + B*np.sin(4*θ)
x_o,y_o = r_out*np.cos(θ), r_out*np.sin(θ)
x_i,y_i = r_in *np.cos(θ), r_in *np.sin(θ)
ms = (r_out.max() - r_in.min())/5.0
pts_o = [gmsh.model.occ.addPoint(x,y,0,ms) for x,y in zip(x_o,y_o)]
pts_i = [gmsh.model.occ.addPoint(x,y,0,ms) for x,y in zip(x_i,y_i)]
curves_q = []
for i in range(n_pt-1):
curves_q.append(gmsh.model.occ.addLine(pts_o[i], pts_o[i+1])) # outer arc
curves_q.append(gmsh.model.occ.addLine(pts_o[-1], pts_i[-1])) # θ=π/2
for i in range(n_pt-1,0,-1):
curves_q.append(gmsh.model.occ.addLine(pts_i[i], pts_i[i-1])) # inner arc
curves_q.append(gmsh.model.occ.addLine(pts_i[0], pts_o[0])) # θ=0
loop_q = gmsh.model.occ.addCurveLoop(curves_q)
quarter = gmsh.model.occ.addPlaneSurface([loop_q])
gmsh.model.occ.synchronize()
outer_q = curves_q[:n_pt-1]
thetaPi2_q = curves_q[n_pt-1]
inner_q = curves_q[n_pt:-1] # <-- this slice was fixed earlier
theta0_q = curves_q[-1]
# -----------------------------------------------------------------------------
# 3) Fragment copy of 3/4 with quarter patch
# -----------------------------------------------------------------------------
cp = gmsh.model.occ.copy([(2, three_tag)])
copy_tag = next(t for d,t in cp if d==2)
gmsh.model.occ.synchronize()
out, mapEnt = gmsh.model.occ.fragment(
[(2, copy_tag)], [(2, quarter)],
removeObject=False, removeTool=False
)
gmsh.model.occ.synchronize()
# Tag surfaces
surfs = [t for d,t in out if d==2]
surf3, surf1 = surfs
pg = gmsh.model.addPhysicalGroup(2, [surf3], tag=20)
gmsh.model.setPhysicalName(2, pg, "three_quarter_material")
pg = gmsh.model.addPhysicalGroup(2, [surf1], tag=21)
gmsh.model.setPhysicalName(2, pg, "quarter_material")
# -----------------------------------------------------------------------------
# 4) Radial curve tags via mapEnt
# -----------------------------------------------------------------------------
theta0_new, thetaPi2_new = [], []
for (dNew, tNew), elders in zip(out, mapEnt):
if dNew!=1: continue
old1D = [t for (d,t) in elders if d==1]
if l4 in old1D: theta0_new.append(tNew)
elif l1 in old1D: thetaPi2_new.append(tNew)
theta0_all = [theta0_q] + theta0_new
thetaPi2_all = [thetaPi2_q] + thetaPi2_new
# -----------------------------------------------------------------------------
# 5) Classify the 3/4 arcs, *first* dropping any radials by angle
# -----------------------------------------------------------------------------
# ensure 1D mesh so every curve has nodes
gmsh.model.mesh.generate(1)
# grab all boundary‐curve tags of the 3/4‐piece
b3all = gmsh.model.getBoundary([(2, surf3)],
oriented=False, recursive=False)
curves3 = [t for (d,t) in b3all if d==1]
# now drop the θ≈0 and θ≈π/2 radials by checking the mean angle of their nodes
tol = 1e-3
filtered = []
for c in curves3:
_, coords, _ = gmsh.model.mesh.getNodes(1, c)
xs, ys = coords[0::3], coords[1::3]
theta_mean = np.mean(np.arctan2(ys, xs))
# normalize to [0,2π)
if theta_mean < 0:
theta_mean += 2*np.pi
# if it's within tol of 0 or π/2, skip it
if abs(theta_mean) < tol or abs(theta_mean - 0.5*np.pi) < tol:
continue
filtered.append(c)
# now only the two arc‐curves remain
curves3 = filtered
rinner3, router3 = [], []
for c in curves3:
_, coords, _ = gmsh.model.mesh.getNodes(1, c)
xs, ys = coords[0::3], coords[1::3]
ravg = np.mean(np.hypot(xs, ys))
if abs(ravg - r_i) < abs(ravg - r_o):
rinner3.append(c)
else:
router3.append(c)
# -----------------------------------------------------------------------------
# 6) Create the six Physical Curve groups
# -----------------------------------------------------------------------------
groups = {
22: ("theta_0", theta0_all),
23: ("theta_pi_over_2", thetaPi2_all),
24: ("r_inner_quarter", inner_q),
25: ("r_outer_quarter", outer_q),
26: ("r_inner_three_quarter", rinner3),
27: ("r_outer_three_quarter", router3),
}
for tag,(name,clist) in groups.items():
if not clist:
raise RuntimeError(f"No curves found for {name}")
pg = gmsh.model.addPhysicalGroup(1, clist, tag=tag)
gmsh.model.setPhysicalName(1, pg, name)
# -----------------------------------------------------------------------------
# 7) Mesh 2D & write
# -----------------------------------------------------------------------------
gmsh.model.mesh.generate(1)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.partition(0) # collapse everything into partition 0
gmsh.write("annulus_full.msh")
gmsh.finalize()
The produced mesh looks totally fine with Gmsh GUI. All the physical curves and surfaces seem to be OK.
However, my working FEnicsx code now yields:
Info : Reading 'meshes/annulus_full.msh'...
Info : 153 entities
Info : 6782 nodes
Info : 13508 elements
Info : Done reading 'meshes/annulus_full.msh'
Invalid rank, error stack:
internal_Issend(117): MPI_Issend(buf=0x55c7aed13b21, count=1, MPI_BYTE, 1, 1, comm=0x84000005, request=0x55c7ae7e87c4) failed
internal_Issend(78).: Invalid rank has value 1 but must be nonnegative and less than 1
Abort(656902) on node 0 (rank 0 in comm 416): application called MPI_Abort(comm=0x84000005, 656902) - process 0
Which makes me think something is wrong with my mesh, but it’s very hard to debug, I don’t know what line/feature is faulty.
A MWE (untested):
import gmsh
from dolfinx.io import gmshio, VTXWriter
from mpi4py import MPI
from ufl import Measure, SpatialCoordinate, dx, FacetNormal
class TheMesh:
def __init__(self, mesh_file="meshes/quarter_annulus.msh") -> None:
self.mesh, self.cell_markers, self.facet_markers = gmshio.read_from_msh(mesh_file, MPI.COMM_WORLD, gdim=2)
self.x = SpatialCoordinate(self.mesh)
self.dx = Measure("dx", domain=self.mesh, subdomain_data=self.cell_markers)
self.ds = Measure("ds", domain=self.mesh, subdomain_data=self.facet_markers)
self.dS = Measure("dS", domain=self.mesh, subdomain_data=self.facet_markers)
self.n = FacetNormal(self.mesh)
the_mesh=TheMesh(mesh_file="meshes/annulus_full.msh")
Any pointer would be very much welcome.