Meshing fails for small domain

Netgen
1

I am trying to mesh a geometry that consists of a cube with a sphere and cylinders cut out. I want to use the physical sizes of the domain, which means there are some small features I want to resolve. The attached script works fine if the geometry is, for example, in a unit cube, but crashes or gets stuck if letting the size of the domain be smaller than 1/1000. (1/1000 works, but 1/1100 fails. The size I want to use is 1/2800).
All features that I want to include are scaled with the side length of the cube, hence it is only a matter of scaling the size of the gridded domain.

I have tried several things, which makes it fail in different ways.

  1. If I do not specify a maxh, the script hangs and gives a segementation fault when I stop it.
NETGEN-6.2.2403
Developed by Joachim Schoeberl at
2010-xxxx Vienna University of Technology
2006-2010 RWTH Aachen University
1996-2006 Johannes Kepler University Linz
Including OpenCascade geometry kernel
optfile ./ng.opt does not exist - using default values
togl-version : 2
OCC module loaded
loading ngsolve library
NGSolve-6.2.2403
Using Lapack
Including sparse direct solver UMFPACK
Running parallel using 16 thread(s)
 Start Findpoints
 Analyze spec points
 Find edges
 Start Findpoints
 Analyze spec points
 Find edges
 Start Findpoints
 Analyze spec points
 Find edges
 Surface 1 / 16
Thank you for using NGSolve
Segmentation fault (core dumped)
  1. If specifying maxh (relative to the cube side length) it gives an error and then hangs.
NETGEN-6.2.2403
Developed by Joachim Schoeberl at
2010-xxxx Vienna University of Technology
2006-2010 RWTH Aachen University
1996-2006 Johannes Kepler University Linz
Including OpenCascade geometry kernel
optfile ./ng.opt does not exist - using default values
togl-version : 2
OCC module loaded
loading ngsolve library
NGSolve-6.2.2403
Using Lapack
Including sparse direct solver UMFPACK
Running parallel using 16 thread(s)
 Start Findpoints
 Analyze spec points
 Find edges
 Start Findpoints
 Analyze spec points
 Find edges
 Start Findpoints
 Analyze spec points
 Find edges
 Surface 1 / 16
 ERROR: Problem in Surface mesh generation
Traceback (most recent call last):
  File "<string>", line 38, in <module>
netgen.libngpy._meshing.NgException: Problem in Surface mesh generation

Do I have to respect anything in terms of mesh size? Or are there any constraints or limitations for gridding small features in netgen?

I found a similar post in the old forum, but it did not explain the solution nor solve my problem.

smallmesh.py (1.1 KB)

2

Hi,
Can I suggest you to move to the newer, more modern opencascade interface? This is more robust here than the csg kernel (and also more powerful):

from netgen.occ import *

sx = 1/1100
sy = sx
sz = sx

# size of sylinder and sphere inside cell
cyl_radius = 0.14*sx
sphere_radius = 0.364*sx

max_mesh_size = sx*0.1

#solid_domain = OrthoBrick(Pnt(0,0,0), Pnt(sx,sy,sz))
# that gives a cube
solid_domain = Box((0,0,0), (sx,sy,sz))

sphere = Sphere((sx/2,sy/2,sz/2), sphere_radius)
cylinders = Cylinder((sx/2,sy/2,0.0), Z, cyl_radius, sz)
cylinders = cylinders + Cylinder((0.0, sy/2, sz/2), X, cyl_radius, sx)
cylinders = cylinders + Cylinder((sx/2, 0.0, sz/2), Y, cyl_radius, sy)

fluid_domain = sphere + cylinders
fluid_domain = solid_domain * fluid_domain

solid_domain = solid_domain - fluid_domain

geo = OCCGeometry(Glue([fluid_domain, solid_domain]))
mesh = Mesh(geo.GenerateMesh(maxh=max_mesh_size))
mesh.Curve(3)

Draw(mesh)

best, Christopher

3

Thanks a lot for point for pointing me to the opencascade interface (I was not aware of it) and the revised code. This works much better!

4

One further trick here: Sometimes I run into troubles still with very small geometries (it might be that we do not use tolerances in OCC always 100% correct in our interface or occ is not always 100% stable there - not sure). But what helps then is to design the geometry in O(1) size and then scale mesh and geometry to size after meshing:

mesh = geo.GenerateMesh(...)
mesh.Scale(1e-9)
geo = OCCGeometry(geo.shape.Scale((0,0,0), 1e-9))
mesh.SetGeometry(geo)
mesh = Mesh(mesh)
mesh.Curve(5)

this usually works then.

best