Build a torus with gmsh - Error when loading mesh into DOLFINx

mesh
1

Hi all,

I am trying to build and create a volume discretization of a torus using gmsh. While I can successfully create the mesh with gmsh API, I have issues when I try to load the mesh directly into DOLFINx, using the approach suggested in this tutorial.

Please note that, although I am aware that a torus can be built by making use of native OCC functions, this is just a toy problem that I need to extend to more complicated geometries. That is why I prefer setting up an arbitrary geometry from the beginning.

However, when I invoke the function dolfinx.io.gmshio.read_from_msh, I get an out of bound error.

import gmsh
import numpy as np
from mpi4py import MPI
import os

_OUTPUT_PATH = "/root/shared/out"

r = 0.1   # Radius of torus cross section [m]
a = 1     # Radius of torus [m]
h = a     # Axial distance between the two tori [m]
gdim = 3  # Geometric dimension of the mesh
mesh_density = 20E-3  # Mesh density [m]

rank = MPI.COMM_WORLD.rank

gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 0)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:

    # Define torus geometry
    top_cross_section = gmsh.model.occ.addCircle(a, 0, 0.5*h, r, zAxis = [0, 1, 0])
    top_cross_section_loop = gmsh.model.occ.addCurveLoop([top_cross_section])
    top_cross_section_surface = gmsh.model.occ.addPlaneSurface([top_cross_section_loop])
    gmsh.model.occ.synchronize()
    
    gmsh.model.occ.revolve([(2, top_cross_section_surface)], 0,0,0, 0,0,1, np.pi)
    gmsh.model.occ.revolve([(2, top_cross_section_surface)], 0,0,0, 0,0,1, -np.pi)

    # glue volumes, set mesh size and sync to model
    gmsh.model.occ.fragment(gmsh.model.occ.getEntities(3), [])
    gmsh.model.occ.mesh.setSize(gmsh.model.occ.getEntities(0), mesh_density)  # Only entities of dimension 0 (points) are handled by OCC setSize.

    gmsh.model.occ.synchronize()

    # Generate mesh
    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.option.setNumber("Mesh.Algorithm3D", 4)
    gmsh.model.mesh.generate(gdim)
    
    meshfile = os.path.join(_OUTPUT_PATH, 'torus.msh')
    gmsh.write(meshfile)

from dolfinx.io.gmshio import model_to_mesh
mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

The error I receive is the following:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [2], in <cell line: 2>()
      1 from dolfinx.io.gmshio import model_to_mesh
----> 2 mesh, ct, ft = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
      3 gmsh.finalize()

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/gmshio.py:205, in model_to_mesh(model, comm, rank, gdim)
    202 perm_sort = np.argsort(cell_dimensions)
    204 # Broadcast cell type data and geometric dimension
--> 205 cell_id = cell_information[perm_sort[-1]]["id"]
    206 tdim = cell_information[perm_sort[-1]]["dim"]
    207 num_nodes = cell_information[perm_sort[-1]]["num_nodes"]

IndexError: index -1 is out of bounds for axis 0 with size 0

The mesh looks correctly generated, in the sense the mesh is successful and I can open and inspect it with gmsh (after saving to file for inspection).

I am using the above snippet in the docker container dokken92/dolfinx_custom:15072022, running DOLFINx version: 0.4.2.0.

Can someone please help me figure out what I am doing wrong?

Thanks!

2

You need to add physical groups to your mesh, as the function read_from_msh states

Reads a mesh from a msh-file and returns the distributed DOLFINx
mesh and cell and facet markers associated with physical groups
in the msh file.

1 Like