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)