Issue reading the mesh tags imported from gmsh

General
, , , ,
1

Hello all,

I know topic (reading the mesh from gmsh) has been posted and answered several times, even the documentation gives great clarity in explaining the subdomains. I heave tried with the different geometries before and worked fine.

But I am facing the wierd issue now of malloc in reading the mesh tags from the mesh imported from gmsh for this particular case although it writes the tags smoothly which I verified.

Here is MWE

import dolfinx
from mpi4py import MPI
import meshio
import numpy as np
from dolfinx.io import XDMFFile

comm = MPI.COMM_WORLD
rank = comm.rank

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
    return out_mesh

if rank == 0:
    domain = meshio.read("dp.msh")
    # Create and save one file for the mesh, and one file for the facets
    triangle_mesh = create_mesh(domain, "triangle", prune_z=True)
    line_mesh = create_mesh(domain, "line", prune_z=True)
    meshio.write("mesh.xdmf", triangle_mesh)
    meshio.write("mt.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

I don’t see an option of uploading the file here to upload my .msh file, in case you need to try on ur end, please let me know, I will just paste it here in comments.

Getting error while reading the tags

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see FAQ — PETSc 3.21.1 documentation and FAQ — PETSc 3.21.1 documentation
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0

Thanks!

2

Please try using dolfinx.io.gmshio as you see in Defining subdomains for different materials — FEniCSx tutorial

Regardless of this suggestion, it’s impossible for us to comment any further without accessing the mesh file.

3

This is usually related to the mesh generation not doing the right thing, i.e. you have tagged some surfaces in your model that are not included in the actual mesh. This sometimes happens when using gmsh. Please provide a minimal script that shows how you created your msh file.

4

Thanks @francesco-ballarin and @dokken for your kind response.

After just restarting the gmsh and deleting the commented lines in .geo file to generate the .msh file, its working now.

5

Hello @dokken and @francesco-ballarin,

I am getitng petsc malloc error again while reading the mesh tags, while it can write the xdmf file correctly as I can view the tags on paraview.

MWE

import dolfinx
from dolfinx import default_scalar_type
from mpi4py import MPI
import meshio
from petsc4py import PETSc
import numpy as np
from dolfinx.io import XDMFFile, gmshio

dtype = PETSc.ScalarType  # type: ignore
# dtype = default_scalar_type  # type: ignore
comm = MPI.COMM_WORLD
rank = comm.rank

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells},
                           cell_data={"name_to_read": [cell_data.astype(np.int32)]})
    return out_mesh

if rank == 0:
    # Read in mesh
    msh = meshio.read("Tshape.msh")
    # Create and save one file for the mesh, and one file for the facets
    triangle_mesh = create_mesh(msh, "quad", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write("mesh.xdmf", triangle_mesh)
    meshio.write("mt.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()


with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

I have uploaded the .geo and .msh files, can you please check from your end, if it results in same problem, if yes, can you please let me know how to correct this.

https://drive.google.com/drive/folders/1D2GfsxqVKjls8pIKSCSWkLEW4IwqMOYT?usp=sharing

Thanks.

6

I would like to ask you to put the geo-file here, as google drive links tend to die over time.

7

Inspecting your mesh, i see that there are more points in your mesh than those that are part of the cells:

(Pdb++) triangle_mesh.points.shape
(64010, 2)
(Pdb++) len(np.unique(triangle_mesh.cells_dict["quad"].reshape(-1)))
63593

which is problematic.

This is probably due to the definition of Transfinite surfaces after making Physical Surface and Physical Curve.

8

Thanks @dokken for response.
Once I added the physical groups to all the surfaces, the issue was lifted and was able to read the mesh and mesh tags generated through Gmsh.

Regards.