I would like to triangulate a domain composed of a box cut by two closed ellipsoids, see the following figure (blue part).
In python, when I set the coeff
below 15, I got the following error:
Error : Impossible to mesh periodic surface 2
terminate called after throwing an instance of 'std::runtime_error'
what(): Impossible to mesh periodic surface 2
After I make it bigger than 16, everything will be fine. Why did this happen? How can I avoid this error?
Another question is how can I catch this error in python?
Setting the option General.AbortOnError
didn't work.
import gmsh
import numpy as np
gmsh.initialize()
gmsh.model.add("test")
alpha, kappa = 1.2, 1
h = 1/64
a_upper = (0.5, 0.5, 0.5)
a_lower = (0.5, 0.5, 0.5)
epsilon = 1e-5
sep = epsilon/2
s = (kappa*h)**(1/(2 - alpha))
d = a_upper[2] + a_lower[2] + epsilon \
- a_upper[2]*np.sqrt(1 - s**2/a_upper[0]**2) - a_lower[2]*np.sqrt(1 - s**2/a_lower[0]**2)
def make_ellipsoid(cx, cy, cz, ax, ay, az):
ellipsoid = gmsh.model.occ.add_sphere(cx, cy, cz, 1)
gmsh.model.occ.dilate([(3, ellipsoid)], cx, cy, cz, ax, ay, az)
return ellipsoid
c_upper = make_ellipsoid(0, 0, a_upper[2] + sep, *a_upper)
c_lower = make_ellipsoid(0, 0, - a_lower[2] - sep, *a_lower)
d = a_upper[2] + a_lower[2] + epsilon \
- a_upper[2]*np.sqrt(1 - s**2/a_upper[0]**2) - a_lower[2]*np.sqrt(1 - s**2/a_lower[0]**2)
# inner = gmsh.model.occ.add_cylinder(0, 0, -a_lower[2], 0, 0, a_lower[2]+a_upper[2], s)
inner = gmsh.model.occ.add_cylinder(0, 0, -d, 0, 0, 2*d, s)
omega = gmsh.model.occ.cut([[3, inner]], [[3, c_upper], [3, c_lower]], removeTool=True)
gmsh.model.occ.synchronize()
coeff = 15 # 16 is ok
def meshSizeCallback(dim, tag, x, y, z, lc):
return max((2*np.sqrt(x**2 + y**2))**alpha*h, coeff*epsilon)
gmsh.model.mesh.setSizeCallback(meshSizeCallback)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
# 0: no, 1: abort meshing, 2: throw an exception unless in interactive mode,
# 3: throw an exception always, 4: exit
gmsh.option.setNumber("General.AbortOnError", 3)
try:
ret = gmsh.model.mesh.generate()
except:
print('Error in mesh generation!')