How to use gmesh to partition a sphere with high-quality tetrahedra

General
1

Hi everyone !
I am trying to dissect a concentric spherical ring, which is a large ball digging out a small ball.In order to make the grid more uniform, I added some extra points. But it didn’t work in the final mesh generation.

import gmsh
from mpi4py import MPI
from dolfinx.io import XDMFFile, gmshio
import numpy as np
gmsh.initialize()     

model = gmsh.model()     
model.add("bigSphere-smallSphere")
model.setCurrent("bigSphere-smallSphere")

bigsphere_dim_tags = model.occ.addSphere(0, 0, 0, 1)
smallsphere_dim_tags=model.occ.addSphere(0, 0, 0, 0.35)
model_dim_tags = model.occ.cut([(3, bigsphere_dim_tags)], [(3, smallsphere_dim_tags)])   
rr=np.arange(0.35, 1.01, 0.01)
d1=np.arange(0, np.pi+0.1, 0.1)       # [0，pi]
d2=np.arange(0, 2*np.pi+0.1, 0.01)    #  [0,2pi]

for i in rr:
    for j in d1:
        for k in d2:
            point_tags=model.occ.addPoint(i*np.cos(j)*np.cos(k),i*np.cos(j)*np.cos(k),i*np.sin(j))

lc=0.1        
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)

model.occ.synchronize()  


boundary = model.getBoundary(model_dim_tags[0], oriented=True)   
boundary_ids = [b[1] for b in boundary]
model.addPhysicalGroup(2, boundary_ids, tag=1)     
model.setPhysicalName(2, 1, "Sphere surface")

volume_entities = [model[1] for model in model.getEntities(3)]  
model.addPhysicalGroup(3, volume_entities, tag=2)
model.setPhysicalName(3, 2, "Sphere volume")
model.mesh.generate(3)

mesh_comm = MPI.COMM_WORLD
model_rank = 0

msh, mt, ft = gmshio.model_to_mesh(model, mesh_comm, model_rank)
msh.name = "ball_d1"
mt.name = f"{msh.name}_cells"
ft.name = f"{msh.name}_facets"

with XDMFFile(msh.comm, "out_gmsh/mesh.xdmf", "w") as file:
    file.write_mesh(msh)
    msh.topology.create_connectivity(2, 3)
    file.write_meshtags(mt, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
    file.write_meshtags(ft, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=3)
ft.name = "Facet markers"
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=3)

image
image700×700 235 KB

2

You have to explain us what you mean by didn’t work.

If you mean that the mesh you get in paraview (below your code) is indeed very odd, just be aware that what you are showing has a visualization artifact. When you do the clip in paraview, go to the properties (typically, bottom left of the screen), scroll down till you find the option “Crinkle clip”, and enable it. This ensures that full tethraedra are displayed and not clipped by the plane.

3

Thanks francesco-ballarin.When I choose the option “Crinkle clip”.The grid is not as strange as it looked before.

image
image918×775 215 KB

image
image960×723 227 KB

But it didn’t work in the final mesh generation
It is mean that I noticed which the points I added did not appear in the mesh generation results.
This is reflected in the code, which means that commenting out the following lines does not affect the results.

rr=np.arange(0.35, 1.01, 0.01)
d1=np.arange(0, np.pi+0.1, 0.1)       # [0，pi]
d2=np.arange(0, 2*np.pi+0.1, 0.01)    #  [0,2pi]

for i in rr:
    for j in d1:
        for k in d2:
            point_tags=model.occ.addPoint(i*np.cos(j)*np.cos(k),i*np.cos(j)*np.cos(k),i*np.sin(j))

This is the result of commenting out these few lines of code

image
image819×166 17.2 KB

This is the result of not commenting on these few lines of code
image
image765×132 13.1 KB

The multiple nodes are the ones I added, but they do not appear in my visualization results.

4

See: [Gmsh] How to add a Point to a Surface with the gmsh Python API?