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

mesh
103

The only information that I have for the subdomains per-se is the corresponding element-sets from the Abaqus mesh, each of which gets stored as a CellBlock in meshio. Ideally this should be analogous to a corresponding physical group in Gmsh. If I try to run

meshio-convert filename.inp filenameGmsh.msh

it gives out the following error:

meshio._exceptions.WriteError: Can only deal with one cell type for now

So I wanted to know if given mesh.points and mesh.cells can I generate subdomain-tags and export a “subdomains.xdmf” analogous to gmsh_physical_region.xdmf. Or let me ask in another way: what does “tetra_data” look like in your example? It is a simple list, with shape (numcells,) ?

104

Yes, the tetra_data is a simple list of shape (numcells,). To illustrate this, you could use this example:

import numpy as np
import pygmsh
import meshio
def create_mesh_gmsh():
    geom = pygmsh.built_in.Geometry()
    rect = geom.add_rectangle(0.0, 2.0, 0.0, 1.0, 0.0, lcar=0.2)
    geom.add_physical([rect.line_loop.lines[0], rect.line_loop.lines[2]], 1)
    geom.add_physical([rect.line_loop.lines[1]], 2)
    geom.add_physical([rect.line_loop.lines[3]], 3)
    geom.add_physical([rect.surface], 4)

    # Generate mesh
    mesh = pygmsh.generate_mesh(geom, dim=2, prune_z_0=True)
    cells = np.vstack(np.array([cells.data for cells in mesh.cells
                                if cells.type == "triangle"]))
    triangle_mesh = meshio.Mesh(points=mesh.points,
                                cells=[("triangle", cells)])

    facet_cells = np.vstack(np.array([cells.data for cells in mesh.cells
                                      if cells.type == "line"]))
    facet_data = mesh.cell_data_dict["gmsh:physical"]["line"]

    facet_mesh = meshio.Mesh(points=mesh.points,
                             cells=[("line", facet_cells)],
                             cell_data={"name_to_read": [facet_data]})

    # Write mesh
    meshio.xdmf.write("mesh.xdmf", triangle_mesh)
    meshio.xdmf.write("facet_mesh.xdmf", facet_mesh)

Here mesh.cells contains 5 CellBlock, four of them with facet data.

1 Like
105

Thanks! Will try it out.

106

Hi everyone,
I have followed all the suggestion you are giving using the last version of meshio. However, the code doesn’t recognize 3 of the 4 boundary conditions.
(*** Warning: Found no facets matching domain for boundary condition.)
I am trying to compute a mesh in 2D, following the code:

msh = meshio.read("presa.msh")

for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "line":
        line_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        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.write("mf.xdmf", line_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
File("Dolfin_circle_mesh.pvd").write(mesh)

mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("Dolfin_circle_facets.pvd").write(mf).

V = FunctionSpace(mesh, 'P', 1)

bc1 = DirichletBC(V, Constant(0.0), mf, 1)
bc2 = DirichletBC(V, Constant(10.0), mf, 2)
bc3 = DirichletBC(V, Constant(6.0), mf, 3)
bc4 = DirichletBC(V, Constant(3.0), mf, 4)

Additionally, I am using the next .geo file to generate the mesh file:

//#gmsh file
Point(6) = {0, 0, 0, 1};
Point(7) = {2, 0, 0, 1};
Point(8) = {0, 2, 0, 1};
Point(9) = {2, 2, 0, 1};
Line(10) = {6, 7};               
Line(20) = {7, 9};           
Line(30) = {9, 8};       
Line(40) = {8, 6};       
Line Loop(1) = {10, 20, 30, 40};

Plane Surface(1) = {1};
Physical Surface(0) = {1};

Physical Line(1) = {40};
Physical Line(2) = {20};
Physical Line(4) = {30};
Physical Line(3) = {10};

Could you please help me with the problem? I am not sure if it is because of the version of meshio, or I am doing something wrong here. Thank you so much.

107

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

108

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

109

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

110

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:

3 Likes
111

It works perfect! Thank you!!!

112

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

113

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.

114

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)
1 Like
115

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

117

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

118

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

1 Like
119

Thank you dokken for your help

120

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.

1 Like
121

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.

3 Likes
122

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!

123

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} ))
2 Likes