Problem reading mesh using meshio or interfacing msh with dolfinx

dolfinx
1

Hello everyone. I am trying to create simple square geometry using gmsh to test this approach. I used both the approaches and have issues.

  1. Using gmshio.model_to_mesh. On using this approach, I do not see any mesh when opened in Paraview although I have generated the mesh (picture below). What could be the issue here?
from mpi4py import MPI
from pathlib import Path
import dolfinx
from dolfinx.io import VTXWriter
from dolfinx.io import gmshio
import gmsh, meshio

#-------------------------------------------------------------------#

#-------------- gmshio---------------------------------#

# Geometry
gmsh.initialize() # Before using any functions in the Python API, Gmsh must be initialized
gmsh.model.add("box_test_1") # name of the model

# Create points
lc = 1e-2 # target mesh size close to the point
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(10.0, 0, 0, lc, 2)
gmsh.model.geo.addPoint(10.0, 10.0, 0, lc, 3)
gmsh.model.geo.addPoint(0, 10.0, 0, lc, 4)

# Create lines
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)

# Create Surface
gmsh.model.geo.addCurveLoop([4, 1, 2, 3], 1)
pl = gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize() # Before they can be meshed the CAD entities must be synchronized with the Gmsh model

# Add Physical Group
gmsh.model.addPhysicalGroup(1, [1], 5)
gmsh.model.addPhysicalGroup(1, [2], 6)
gmsh.model.addPhysicalGroup(1, [3], 7)
gmsh.model.addPhysicalGroup(1, [4], 8)

# Meshing options
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)
gmsh.model.mesh.generate(2)

gmsh.write("test_gm.msh")

# Interfacing with dolfinx
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=2)

with dolfinx.io.XDMFFile(domain.comm, "test_gmsh.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)

Screenshot from 2023-12-14 16-36-03
Screenshot from 2023-12-14 16-36-03932×920 2.92 KB

  1. Using meshio, I have the error given below. What could be the problem here?
from mpi4py import MPI
from pathlib import Path
import dolfinx
from dolfinx.io import VTXWriter
from dolfinx.io import gmshio
import gmsh, meshio

proc = MPI.COMM_WORLD.rank

# msh = meshio.read(Path("geometry")/"quad.msh")
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 proc == 0:
    # Read in mesh
    msh = meshio.read(Path("geometry")/"tri1.msh")

    # Create and save one file for the mesh, and one file for the facets
    triangle_mesh = create_mesh(msh, "triangle", 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 dolfinx.io.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 dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

Error on the cell_data line in def create_mesh:

Exception has occurred: ValueError
need at least one array to concatenate
  File "/home/igupta/studies/Basic_bi_copy/basic_bi.py", line 76, in create_mesh
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/igupta/studies/Basic_bi_copy/basic_bi.py", line 86, in <module>
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: need at least one array to concatenate

The tri1.msh file looks like this

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
4 4 1 0
1 0 0 0 0 
2 10 0 0 0 
3 10 10 0 0 
4 0 10 0 0 
1 0 0 0 10 0 0 1 1 2 1 -2 
2 10 0 0 10 10 0 1 2 2 2 -3 
3 0 10 0 10 10 0 1 3 2 3 -4 
4 0 0 0 0 10 0 1 4 2 4 -1 
1 0 0 0 10 10 0 0 4 4 1 2 3 
$EndEntities
$Nodes
8 40 1 40
0 1 0 1
1
0 0 0
0 2 0 1
2
10 0 0
0 3 0 1
3
10 10 0
0 4 0 1
4
0 10 0
1 1 0 9
5
6
7
8
9
10
11
12
13
0.9999999999991893 0 0
1.999999999996827 0 0
2.99999999999391 0 0
3.999999999991027 0 0
4.999999999992411 0 0
5.999999999993932 0 0
6.999999999995455 0 0
7.999999999996978 0 0
8.99999999999849 0 0
1 2 0 9
14
15
16
17
18
19
20
21
22
10 0.9999999999991893 0
10 1.999999999996827 0
10 2.99999999999391 0
10 3.999999999991027 0
10 4.999999999992411 0
10 5.999999999993932 0
10 6.999999999995455 0
10 7.999999999996978 0
10 8.99999999999849 0
1 3 0 9
23
24
25
26
27
28
29
30
31
8.999999999999998 10 0
7.999999999999996 10 0
6.999999999999994 10 0
5.999999999999991 10 0
4.999999999999989 10 0
3.999999999999988 10 0
2.999999999999985 10 0
1.999999999999982 10 0
0.9999999999999911 10 0
1 4 0 9
32
33
34
35
36
37
38
39
40
0 8.999999999999998 0
0 7.999999999999996 0
0 6.999999999999994 0
0 5.999999999999991 0
0 4.999999999999989 0
0 3.999999999999988 0
0 2.999999999999985 0
0 1.999999999999982 0
0 0.9999999999999911 0
$EndNodes
$Elements
4 40 1 40
1 1 1 10
1 1 5 
2 5 6 
3 6 7 
4 7 8 
5 8 9 
6 9 10 
7 10 11 
8 11 12 
9 12 13 
10 13 2 
1 2 1 10
11 2 14 
12 14 15 
13 15 16 
14 16 17 
15 17 18 
16 18 19 
17 19 20 
18 20 21 
19 21 22 
20 22 3 
1 3 1 10
21 3 23 
22 23 24 
23 24 25 
24 25 26 
25 26 27 
26 27 28 
27 28 29 
28 29 30 
29 30 31 
30 31 4 
1 4 1 10
31 4 32 
32 32 33 
33 33 34 
34 34 35 
35 35 36 
36 36 37 
37 37 38 
38 38 39 
39 39 40 
40 40 1 
$EndElements
2

You should add a physical group for the surface as well.

3

Thank you very much. That was the mistake :slight_smile: