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

Try naming all of your Lines as Physical Line, and not Physical Curve.

Hi Sebastian. Yes, I already did that modification, but I have the same result. I think that my problem is not there. Thank you

What are the boundary conditions you are using? Are they all Dirichlet type?

The problem is that meshio defines multiple cell block for facets. The following works:

import meshio
import numpy as np
msh = meshio.read("juanfe360.msh")

line_cells = []
for cell in msh.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 msh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

triangle_mesh = meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=msh.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)

I’ve explained the same case for 3D-2D in the post above:

2 Likes

It works perfect! Thank you!!!

I am currently trying to covert a simple a .msh to XDMF mesh with the following code from the tutorials above:

msh = meshio.read("test_box.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif cell.type == "tetra":
        tetra_cells = cell.data
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                 cell_data={"name_to_read":[tetra_data]})
triangle_mesh = meshio.Mesh(points=msh.points,
                            cells=[("triangle", triangle_cells)],
                            cell_data={"name_to_read": [triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)
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)

cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

but when it tries to read in the file in infield.read(mesh) I get a ValueError: vector. I cannot figure out why, as the same code works on my colleagues computer. I have attached the beginning of my mesh file! Any help would be appreciated.

mesh file:
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
358
1 0 0 1
2 0 0 0
3 0 1 1
4 0 1 0
5 1 0 1
6 1 0 0
7 1 1 1
8 1 1 0
9 0 0 0.1666666666666662
10 0 0 0.3333333333333324
11 0 0 0.4999999999999986
12 0 0 0.6666666666666646
13 0 0 0.8333333333333321
14 0 0.1666666666666662 1
15 0 0.3333333333333324 1
16 0 0.4999999999999986 1
17 0 0.6666666666666646 1
18 0 0.8333333333333321 1
19 0 1 0.1666666666666662
20 0 1 0.3333333333333324
21 0 1 0.4999999999999986
22 0 1 0.6666666666666646
23 0 1 0.8333333333333321
24 0 0.1666666666666662 0
25 0 0.3333333333333324 0
26 0 0.4999999999999986 0
27 0 0.6666666666666646 0
28 0 0.8333333333333321 0
29 1 0 0.1666666666666662
30 1 0 0.3333333333333324
31 1 0 0.4999999999999986
32 1 0 0.6666666666666646
33 1 0 0.8333333333333321
34 1 0.1666666666666662 1
35 1 0.3333333333333324 1
36 1 0.4999999999999986 1
37 1 0.6666666666666646 1
38 1 0.8333333333333321 1
39 1 1 0.1666666666666662
40 1 1 0.3333333333333324
41 1 1 0.4999999999999986
42 1 1 0.6666666666666646
43 1 1 0.8333333333333321
44 1 0.1666666666666662 0
45 1 0.3333333333333324 0
46 1 0.4999999999999986 0
47 1 0.6666666666666646 0

If it works on your colleagues computers, it is probably a version issue. So which versions of
Meshio, FEniCS and gmsh are you using? And what versions are they using?

Additionally, you need to supply the whole mesh file (make the coarsest possible mesh) for anyone to be able to debug it.

Hi dokken,

With this code I got a 2D mesh embedded in a 3D mesh. Replacing “msh.points” by
“points=msh.points[:, :2]” in the lines declaring triangle_mesh, line_mesh I obtained a true 2D mesh when reading it from Fenics. I am including the full script.

Thanks for your help

import meshio
import numpy as np
msh = meshio.read("juanfe360.msh")

line_cells = []
for cell in msh.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 msh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

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

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

I am not sure what you are looking for of guidance here, as you are not asking a question.

Hello dokken
I have an error as following : KeyError: ‘tetra’
when I change “tetra” by triangle , it is working.

That means that there are no Tetrahedron cells in your mesh. Make sure to use physical markers on cells and facets

1 Like

Thank you dokken for your help

Hi,
How can I set boundary conditions from mf file?
For example:
bc_bottom = DirichletBC(V, Constant(0.0), mf, X), where X is the tag of boundary (that name in Gmsh is Bottom_boundary) in mf file.

You should not use string tags in Gmsh, as They cannot be read in to a mesh function. Use integers to tag the domains, and the syntax you proposed above will work.

1 Like

Hi,
I am having a problem in converting a 1D mesh generated with gmsh to dolfin.
(I am using gmsh 4.5.5, meshio 4.2.0 and dolfin 2019.1.0)
I use the following code to convert the .msh file to XDMF

import numpy as np
import meshio
#read msh
msh = meshio.read('example1D.msh')
fn = 'example1D'
dim_keep = 2
#get cells
line_cells = []
vertex_cells = []
for cell in msh.cells:
    if cell.type == 'line':
        if len(line_cells) == 0:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])

    if cell.type == 'vertex':
        if len(vertex_cells) == 0:
            vertex_cells = cell.data
        else:
            vertex_cells = np.vstack([vertex_cells, cell.data])
#get cells data
line_data = None
vertex_data = None
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "vertex":
        vertex_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "line":
        line_data = msh.cell_data_dict["gmsh:physical"][key]

line_mesh = meshio.Mesh(points = np.atleast_2d(msh.points[:,:dim_keep]),
                        cells = [("line", line_cells)],
                        cell_data = {"name_to_read":[line_data]})

meshio.write(fn + '_1D_mesh.xdmf', line_mesh)
if not (vertex_data is None):
    vertex_mesh = meshio.Mesh(points = msh.points[:,:dim_keep], 
                              cells = [("vertex", vertex_cells)],
                              cell_data = {"name_to_read":[vertex_data]})
    meshio.write(fn + '_0D_mesh.xdmf', vertex_mesh)

note that I am using the parameter dim_keep to keep only the first dimensions.
I then import the mesh and the markers using

import dolfin as dol
mesh = dol.Mesh()
with dol.XDMFFile(fn + '_1D_mesh.xdmf') as infile:
    infile.read(mesh)
mvc = dol.MeshValueCollection('size_t', mesh= mesh, dim = 1) 
with dol.XDMFFile(fn + '_1D_mesh.xdmf') as infile:
    infile.read(mvc, 'name_to_read')
marker_1D = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = dol.MeshValueCollection('size_t', mesh= mesh, dim = 0) 
with dol.XDMFFile(fn + '_0D_mesh.xdmf') as infile:
    infile.read(mvc, 'name_to_read')
marker_0D = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)

This works as long as I use dim_keep =2,3 but if I use dim_keep =1 (I want to have a truly 1D mesh, not embedded in 2 or 3D space) I get the following error when trying to load the mesh (infile.read(mesh)):

*** Error:   Unable to determine geometric dimension.
*** Reason:  GeometryType "X" in XDMF file is unknown or unsupported.
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown

It is not a big deal since i can always embed my 1D mesh in a 2d space but I don’t understand what is happening.

For the sake of completeness my .geo file reads

SetFactory("OpenCASCADE");
Point(1) = {0,0,0,0.1};
Point(2) = {1,0,0,0.2};
Line(1) ={1,2};
Physical Line('domain',1) = {1};
Physical Point('left',1) = {1};
Physical Point('right', 2) = {2};

Thanks for helping!

There was a change in the structure of cells to a list of cell blocks.
This is my script (meshio 4.0.16, gmsh 4.4.1) to convert volume, boundary and physical groups from gmsh to vtu. For xdmf there is not much too change.

#Read gmsh file (2D) and write vtu-files (or other formats) for domain, 
#boundaries and possibly subdomains according to physical groups in gmsh. 
#All mesh entities must belong to a physical group!

#Tested with Gmsh 4.4.1  (file format seems to change quite often)

#Example file and initial version adopted from
#https://github.com/iitrabhi/dolfinx/tree/iitrabhi/mvc-xdmf/python/demo


import meshio    # currently 4.0.16
from numpy import array

mesh = meshio.read('poisson_subdomain.msh')
# points: 2Darray( point, xyz)
# cells: list( cellblocks( type='line'/'triangle', data=2Darray(element, points) ))
# cell_data: dict('gmsh:physical', 2Darray(geo_tag, ph_tag))
# field_data: dict('top'/'bottom'/PH_NAME, 1Darray (ph_tag, geo_type)   
points, cells = mesh.points, mesh.cells
cell_data, field_data = mesh.cell_data, mesh.field_data


# some variable declarations
type_index=0    # to access type ('line' or 'triangle') of a cellblock
data_index=1    # to access data (nodes of elements) of a cellblock
gmshdict={1:'line', 2: 'triangle'}      # gmsh convention
gmsh_string='gmsh:physical'


# A mesh consists of cellblocks, 
# now we go through them at first for domain and boundary, 
# which are easily recognizable by celltype
domain_cells=[]         # start with empty lists
domain_cell_data=[]
boundary_cells=[]    
boundary_cell_data=[]

for cellblock_data, cellblock in zip(physical_cell_data, cells):

        if cellblock[type_index]=='line':
                boundary_cells.append(cellblock)
                boundary_cell_data.append(cellblock_data)

        if cellblock[type_index]=='triangle':
                domain_cells.append(cellblock)
                domain_cell_data.append(cellblock_data)

# write results to file
if len(domain_cell_data):  # only if there are data to write
        meshio.write(output_basename+"_domain.vtu",
        meshio.Mesh(points=points, cells=domain_cells, cell_data={domain_data_string: domain_cell_data}))

if len(boundary_cell_data): # only if there are data to write
        meshio.write(output_basename+"_boundary.vtu",
        meshio.Mesh(points=points, cells=boundary_cells, cell_data={boundary_data_string: boundary_cell_data}) )


# now we want to extract subdomains given by physical groups in gmsh
# so we need an additional loop 
for name, data in field_data.items():
# name=user-defined name of physical group, data=[physical_id, geometry_type]
        selected_cells=[]
        selected_cell_data=[]
        ph_id=data[0]   # selection by physical id
        cell_type=gmshdict[data[1]]     # 1='line' or 2='triangle'

        for cellblock_data, cellblock in zip(physical_cell_data, cells):

                # access matching data by an index field
                selected_data=array(cellblock[data_index][cellblock_data==ph_id]) #
                if len(selected_data):  # append only nonzero data
                        selected_cells.append(meshio.CellBlock(cell_type, selected_data))
                        selected_cell_data.append(cellblock_data)

        if len(selected_cell_data): # write only nonzero data
                outputfilename=output_basename+"_physical_group_"+name+".vtu"
                meshio.write(outputfilename,
                meshio.Mesh(points=points, cells=selected_cells,
                cell_data={selection_data_string: selected_cell_data} ))
1 Like

I am afraid, the first for-loop picks only a subset of triangle/tetra, if there are multiple cell blocks.
Instead of for cell in msh.cells, there should be for key in mesh.cells_dict and so on, as it is done for the cell_data. Also a mesh.prune() might be useful to get rid of orphaned points, since in the current example code, all points are saved for both submeshes.

Hi everyone, hi dokken,

I have a geo file, then a 2D msh file created with gmsh, I converted it in xdmf with

meshio-convert whistle.msh whistle.xdmf

I had the same initial problem
“*** -------------------------------------------------------------------------
*** Error: Unable to recognise cell type.
*** Reason: Unknown value “mixed”.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0”

So i read this post and i tried different routines to convert in xdmf properly, but each one failed.
When I tried

meshio.write(“whistle.xdmf”, meshio.Mesh(points=msh.points, cells={“line”: msh.cells[“line”]}))

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”]}}))

I get the error

File “<tmp 1>”, line 9, in
meshio.write(“whistle.xdmf”, meshio.Mesh(points=msh.points, cells={“line”: msh.cells[“line”]}))
TypeError: list indices must be integers or slices, not str

I tried with “tetra” "triangle and “line”, each time same result.

I also tried with

for cell in msh.cells:
if cell.type == “triangle”:
triangle_cells = cell.data
elif cell.type == “tetra”:
tetra_cells = cell.data
for key in msh.cell_data_dict[“gmsh:physical”].keys():
if key == “triangle”:
triangle_data = msh.cell_data_dict[“gmsh:physical”][key]
elif key == “tetra”:
tetra_data = msh.cell_data_dict[“gmsh:physical”][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={“tetra”: tetra_cells},
cell_data={“name_to_read”:[tetra_data]})
triangle_mesh = meshio.Mesh(points=msh.points,
cells=[(“triangle”, triangle_cells)],
cell_data={“name_to_read”: [triangle_data]})
meshio.write(“mesh.xdmf”, tetra_mesh)

I get this mistake :

Traceback (most recent call last):
File “<tmp 1>”, line 13, in
for key in msh.cell_data_dict[“gmsh:physical”].keys():
KeyError: ‘gmsh:physical’

I am not clearly understanding what i am doing, i’m very new to gmsh, meshio and dolfin. Thanks for your answer.

Here is my geo file :

//Dimensions (mm)
h = 7.69;
D=25.0;
delta = 3.0;
e = 0.5;
L1 = 50;
H1 = 50;
L2 = 50;
H2 = 50;
epsilon = 0.5;

// Points
Point(1) = {-L1, H1/2, 0};
Point(2) = {0, H1/2, 0};
Point(3) = {0, delta/2, 0};
Point(4) = {e, delta/2+epsilon, 0};
Point(5) = {e, D/2, 0};
Point(6) = {h+e, D/2, 0};
Point(7) = {h+e, delta/2, 0};
Point(8) = {h+2e, delta/2+epsilon, 0};
Point(9) = {h+2
e, H2/2, 0};
Point(10) = {h+2*e+L2, H2/2, 0};

Point(11) = {-L1, -H1/2, 0};
Point(12) = {0, -H1/2, 0};
Point(13) = {0, -delta/2, 0};
Point(14) = {e, -delta/2-epsilon, 0};
Point(15) = {e, -D/2, 0};
Point(16) = {h+e, -D/2, 0};
Point(17) = {h+e, -delta/2, 0};
Point(18) = {h+2e, -delta/2-epsilon, 0};
Point(19) = {h+2
e, -H2/2, 0};
Point(20) = {h+2*e+L2, -H2/2, 0};

// Lines
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,5};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,9};
Line(9) = {9,10};
Line(10) = {10,20};
Line(11) = {20,19};
Line(12) = {19,18};
Line(13) = {18,17};
Line(14) = {17,16};
Line(15) = {16,15};
Line(16) = {15,14};
Line(17) = {14,13};
Line(18) = {13,12};
Line(19) = {12,11};
Line(20) = {11,1};

// Loops
Line Loop(1) = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};

// Surface
Plane Surface(1) = {1};
Recombine Surface{1};
Transfinite Line{2:8} = 15;
Transfinite Line{12:18} = 15;

First, you need to use physical markers on your lines and surfaces

See for instance this tutorial: https://jsdokken.com/converted_files/tutorial_pygmsh.html
In particular have a look at part 3.
What you should also do is inspect the objects that you are trying to access.
I.e. print msh.cells[“line”] etc.

1 Like

Thanks a lot for your explanation and furthermore for the tutorial. It worked really well and was very clear