Gmsh: generate 3D meshes with spherical caps

mesh
1

Hi all,

I’m using gmsh. I need to mesh a 3D shape formed by a hemisphere centred at the origin of coordinates with radius R_{\beta} and a spherical cap centred at c_{\alpha} with radius R_{\alpha}, as schematised in the following figure:

Copy of diagram

There doesn’t seem to be any built-in method for creating 3D spherical caps. My idea is then to generate the 3D volume as a solid of revolution from its 2D profile, which I could define quite easily:

import gmsh
model = gmsh.model; geom = gmsh.model.geo
from mpi4py import MPI
mesh_comm = MPI.COMM_WORLD

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model.add("3D model")

Rbeta, Ralpha, calpha = 0.6976318239461305, 0.8720397799296024, -0.5232238679565371
ms = 0.1

p = [
    geom.addPoint(0, 0, Ralpha+calpha, ms),
    geom.addPoint(Rbeta, 0, 0, ms),
    geom.addPoint(0, 0, -Rbeta, ms),
    geom.addPoint(0, 0, 0, ms),
    geom.addPoint(0, 0, calpha, ms)
]
l = [
    geom.addCircleArc(p[0], p[4], p[1]),
    geom.addCircleArc(p[1], p[3], p[2])
]

Is it possible to do something like that? I’ve read that you could achieve something like this with the extrude function, but I can’t find any examples of such rotations and I don’t know how to proceed.

Any suggestions would be appreciated, as well as other ways of approaching the problem. Thanks in advance!

2

Please note that this is not a forum about gmsh, but about FEniCS. Some users may have enough experience with both gmsh and FEniCS to help, but this sort of questions are better asked to gmsh developers.

1 Like
3

You should follow gmsh tutorials. To perform revolve, you should use addPlaneSurface to generate surface from arcs. Then you need to get the geometrical entity tag of that surface to revolve it. You can check tutorial 3 for more information.

2 Likes
4

Following @Ekrem_Ekici’s reply, I have been able (I think) to revolve the surface. However, now I’m having trouble with translating the mesh to dolfinx with meshio. Below are the code and the error message. Any ideas what is wrong?

from dolfinx.io import gmshio
import gmsh
model = gmsh.model; geom = gmsh.model.geo
from mpi4py import MPI
mesh_comm = MPI.COMM_WORLD

Rbeta, Ralpha, calpha = 0.6976318239461305, 0.8720397799296024, -0.5232238679565371
ms = 0.1
gdim = 2

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model.add("3D acorn")

p = [
    geom.addPoint(0, 0, Ralpha+calpha, ms),
    geom.addPoint(Rbeta, 0, 0, ms),
    geom.addPoint(0, 0, -Rbeta, ms),
    geom.addPoint(0, 0, 0, ms),
    geom.addPoint(0, 0, calpha, ms)
]
l = [
    geom.addCircleArc(p[0], p[4], p[1]),
    geom.addCircleArc(p[1], p[3], p[2]),
    geom.addLine(p[2], p[0])
]

cl = geom.addCurveLoop([l[j] for j in range(len(l))])
s = geom.addPlaneSurface([cl], 10) 

# Rotations are specified by an axis point (0, 0, 0), the axis direction (0, 0, 1), and a
# rotation angle (Pi/2). We input the dimension and tag of the surface to rotate [(2, 10)]
# and ask for 7 subdivisions.
revsol = geom.revolve([(2, 10)], 0, 0, 0, 0, 0, 1, np.pi/2, [7])

geom.synchronize()
model.add_physical_group(dim=gdim+1, tags=[revsol[1][1]]) # tag of the volume revsol

geom.synchronize()
model.mesh.generate(gdim+1)

domain = gmshio.model_to_mesh(model, mesh_comm, 0, gdim=gdim+1)

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Input In [64], in <cell line: 43>()
     40 geom.synchronize()
     41 model.mesh.generate(gdim+1)
---> 43 domain = gmshio.model_to_mesh(model, mesh_comm, 0, gdim=gdim+1)
     44 plotMesh(domain)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/gmshio.py:188, in model_to_mesh(model, comm, rank, gdim)
    185 x = extract_geometry(model)
    187 # Get mesh topology for each element
--> 188 topologies = extract_topology_and_markers(model)
    190 # Extract Gmsh cell id, dimension of cell and number of
    191 # nodes to cell for each
    192 num_cell_types = len(topologies.keys())

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/gmshio.py:99, in extract_topology_and_markers(model, name)
     92 for entity in entities:
     93     # Get cell type, list of cells with given tag and
     94     # topology of tagged cells
     95     # NOTE: Assumes that each entity only have one
     96     # cell-type, i.e. facets of prisms and pyramid meshes
     97     # are not supported
     98     (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity)
---> 99     assert len(entity_types) == 1
    101     # Determine number of local nodes per element to create the
    102     # topology of the elements
    103     properties = model.mesh.getElementProperties(entity_types[0])

AssertionError:
5

Why are you not using the Open Cascade interface? A hemisphere can be created in three lines using the Boolean operations.

        gmsh.model.occ.addSphere(0, 0, 0, 1.0, 1) # The sphere
        gmsh.model.occ.addBox(0, -1.0, 0, 1.0, 2.0, 1.0, 2) # Box for truncating sphere
        gmsh.model.occ.intersect([(3,1)], [(3,2)], 3, removeObject=True, removeTool=True)
6

Judging from the NOTE in the line before the assert, it probably means that your mesh contains more than one element type (e.g., tet and hex), which isn’t supported (yet) by dolfinx.

7

You should extract return the output of this line as a tuple, see here.