Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

Quentin · November 21, 2023, 11:11am

Thanks for the link.

but it seems to me that this post and the link that you shared mostly covers geometry that are generated from scratch with gmsh, whereas I deal with a CAD in a .step format. Is anyone has already dealt with such a thing ?

Thanks,

Quentin

dokken · November 21, 2023, 11:46am

The whole point of the reference is that it contains the updated syntax for reading in a msh file and converting it to xdmf, replacing

Quentin:

mesh = meshio.read("SLICE.msh")

line_cells = []
for cell in mesh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "line":
        if len(line_cells) == 0:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])

line_data = []
for key in mesh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = mesh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, mesh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = mesh.cell_data_dict["gmsh:physical"][key]

triangle_mesh = meshio.Mesh(points=mesh.points, cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=mesh.points,
                           cells=[("line", line_cells)],
                           cell_data={"name_to_read":[line_data]})
meshio.write("mesh.xdmf", triangle_mesh)

meshio.xdmf.write("mf.xdmf", line_mesh)

with

mesh_from_file = meshio.read("mesh.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]})
    return out_mesh
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

or similar

Quentin · November 21, 2023, 1:32pm

Dear dokken,

thank you for you reply.

Even though the facet mesh is generated sucessfully whit the comand above, for some reason, If I use

tetrahedron_mesh = create_mesh(mesh, "tetrahedron", prune_z=False)
meshio.write("SLICE.xdmf", tetrahedron_mesh)

I get the following error

Traceback (most recent call last):
  File "/mnt/c/Users/quent/Documents/Codes/Python/Mesh_GMSH2MeshIO/XDMF_OUTPUT_BOUNDARIES_TAGS/Conversion_GMSH2MeshIO.py", line 151, in <module>
    tetrahedron_mesh = create_mesh(MAILLAGE, "tetrahedron", prune_z=False)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/mnt/c/Users/quent/Documents/Codes/Python/Mesh_GMSH2MeshIO/XDMF_OUTPUT_BOUNDARIES_TAGS/Conversion_GMSH2MeshIO.py", line 141, in create_mesh
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/quentin_digimind/anaconda3/envs/MESHING/lib/python3.12/site-packages/meshio/_mesh.py", line 227, in get_cell_data
    return np.concatenate(
           ^^^^^^^^^^^^^^^
ValueError: need at least one array to concatenate

even though I added physical groups for volumes as

volumes = gmsh.model.getEntities(dim=3)
for volume in volumes:
    gmsh.model.addPhysicalGroup(3, [volume[1]] ,tag=volume[1] )
    gmsh.model.setPhysicalName(3, volume[1], str(volume[1]))

whereas

meshio.write("SLICE.xdmf", meshio.Mesh(points=mesh.points,cells=mesh.cells))

provides a mesh that seems OK.

Just to let you know .

best regards,

Quentin

dokken · November 21, 2023, 2:55pm

As I do not have your step file, I cannot reproduce your error.

Note that it is not adviced to call
gmsh.model.occ.synchronize()
after adding physical groups, as it tends to remove them,

Quentin · November 22, 2023, 5:11pm

I finally fixe this issue ( was related to an issue with volume physical groups), so I am able to write the mesh as described in the post, but when I try to import the mesh into FeniCS I get the error

mesh = Mesh()
with XDMFFile("SLICE.xdmf") as infile:
    infile.read(mesh)
    
mvc = MeshValueCollection("size_t", mesh, 2) 

with XDMFFile("facet_SLICE.xdmf") as infile:
    infile.read(mvc, "name_to_read")
    
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

I get the error

Traceback (most recent call last):
  File "/mnt/c/Users/quent/Documents/Codes/FENICS/CaseSetup/HeatTransfer/SliceCalculation_FirstGo/OneSliceCalculus.py", line 25, in <module>
    mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
         ^^^^^^^^
AttributeError: module 'mshr.cpp' has no attribute 'mesh'

inspired from other scripts, I tried to declare something like

mf = MeshFunction("size_t",mesh, 2)

but I cannot find a way to get the mvc values into mf …

dokken · November 22, 2023, 5:47pm

Please include the full code, with imports, as this is likely due to the use of multiple wildcard imports within the same script.

Quentin · November 22, 2023, 5:49pm

here it is

import numpy as np
from fenics import *
from dolfin import *
from mshr import *


tol=1E-8

# read a mesh properly written
mesh = Mesh()
with XDMFFile("SLICE.xdmf") as infile:
    infile.read(mesh)
    
mvc = MeshValueCollection("size_t", mesh, 1) 

with XDMFFile("facet_SLICE.xdmf") as infile:
    infile.read(mvc, "name_to_read")
    
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
```

Quentin · November 22, 2023, 5:52pm

indeed only importing fenics solve the problem thanks

dokken · November 22, 2023, 6:40pm

This should be reduced to
from dolfin import *, as fenics and dolfin is the same package.

NataliaKor · November 30, 2023, 11:23am

Dear Dokken. When I try to open this tutorial it shows that “File not found”. Could you send the link again?

dokken · November 30, 2023, 11:32am

Not sure what link you are referring to

NataliaKor · December 1, 2023, 5:24am

https://jsdokken.com/converted_files/tutorial_pygmsh.html This link.

dokken · December 1, 2023, 5:51am

Try Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken

NataliaKor · December 4, 2023, 3:04am

Dear @dokken how we can merge two meshes together for getting one XMl file?

Devin · December 8, 2023, 1:22am

Hi dokken,

I tried this but unfortunately it didn’t work for my case. My python script and the error message can be found below. What mistakes did I make?

from future import print_function
from fenics import *
from dolfin import *
from ufl import nabla_div
from mpi4py import MPI
import numpy as np
import meshio

Scaled variables

E = 2.5e11
nv = 0.27

lambda_ = E * nv / ((1 + nv) * (1 - 2 * nv))
mu = 0.5 * E / (1 + nv)

Create mesh and define function space

mesh = BoxMesh(Point(-1, -1, -1), Point(1, 1, 1), 10, 10, 10)

msh = meshio.read(“S1_V8.msh”)
meshio.write(“mesh.xdmf”, meshio.Mesh(points=msh.points, cells={“tetra”: msh.cells[“tetra”]}))
meshio.write(“mf.xdmf”, meshio.Mesh(points=msh.points, cells={“triangle”: msh.cells[“triangle”]},
cell_data={“triangle”: {“name_to_read”: msh.cell_data[“triangle”][“gmsh:physical”]}}))

meshio.write(“cf.xdmf”, meshio.Mesh(
points=msh.points, cells={“tetra”: msh.cells[“tetra”]},
cell_data={“tetra”: {“name_to_read”: msh.cell_data[“tetra”][“gmsh:physical”]}}))

mesh = Mesh()
with XDMFFile(“mesh.xdmf”) as infile:
infile.read(mesh)
mvc = MeshValueCollection(“size_t”, mesh, 2)
with XDMFFile(“mf.xdmf”) as infile:
infile.read(mvc, “name_to_read”)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

meshio.write(“S1_V8.xml”,msh)
mesh = Mesh(“mesh.xdmf”)

V = VectorFunctionSpace(mesh, ‘P’, 1)

dokken · December 8, 2023, 8:50am

As referred to in the previous post:

The syntax has updated:

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]})
    return out_mesh

Devin · December 8, 2023, 6:45pm

I updated my codes and now I got another issue saying "KeyError:‘gmsh:physical’ ". My msh file is generated from gmsh based on a stl file. My python scrip for gmsh is pasted below.

import gmsh
import sys

gmsh.initialize()
gmsh.merge(‘S1_V8.stl’)

s = gmsh.model.getEntities(2)
l = gmsh.model.geo.addSurfaceLoop([e[1] for e in s])
V0 = gmsh.model.geo.addVolume([l])
gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(3, [V0], 1)

gmsh.model.mesh.generate(3)

#obtain the model dimension
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(-1, -1)

gmsh.write(“S1_V8_tests.msh”)
gmsh.finalize()

dokken · December 8, 2023, 7:42pm

What read command did you now use?
Could you print mesh_from_file and mesh_from_file.cell_data

Devin · December 8, 2023, 8:46pm

I modified something and also tried to print those as you suggested. The read codes and the outputs are:

dokken · December 9, 2023, 6:59pm

The good news is that there is cell data in your grid, marked with “gmsh:physical”, so the error you had in:

does not make sense, and as you can see from your new code, you get way further, as you get to the write stage. Please remove the “tetra_mesh.xdmf” and “tetra_mesh.h5” file and try to run the script again.

Topic		Replies	Views
Extracting mesh regions for gmsh mesh mesh	24	5244	December 5, 2023
Transitioning from Triangle -> GMSH -> XDMF I/O	8	939	April 16, 2023
Need help converting GMSH to FEniCS installation	61	6620	September 28, 2022
Create mesh with gmsh, convert it and import it in dolfin, problems mesh	5	10983	April 1, 2019
How to import mixed triangle and quad meshes from gmsh to fenics mesh	22	4880	June 20, 2020

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

Scaled variables

Create mesh and define function space

Related topics