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

mesh
#99

@dokken I had a quick question regarding the usage of meshio to write meshes for subdomains that can be loaded into dolfin. As of now I have 2 numpy arrays each containing the cells for a subdomain. Something like:

import meshio, os
meshFolder = os.getcwd()
msh = meshio.read("abaqusMesh.inp")
#  msh.points -- all the nodes
#  msh.cells[0] -- cells corresponding to one of the subdomain
#  msh.cells[1] -- cells corresponding to the other
entireMesh = meshio.Mesh(msh.points, np.vstack((msh.cells[0], msh.cells[1])))
particlesMesh = meshio.Mesh(msh.points, msh.cells[1])

#  subdomainMesh = meshio.Mesh(msh.points, ...) # what could be done here?
 
meshio.write(os.path.join(meshFolder, "entireMesh.xdmf"), entireMesh)
meshio.write(os.path.join(meshFolder, "particles.xdmf"), particlesMesh) # -- only the second subdomain
#  and consequently this...
#  meshio.write(os.path.join(meshFolder, "subdomains.xdmf"), subdomainMesh)

I want to be able to assemble forms on subdomains like

from dolfin import *
msh = Mesh()
with XDMFFile("entireMesh.xdmf") as sdf:
    sdf.read(msh)
subdomains = MeshFunction("size_t", msh, "subdomains.xdmf")
dx = Measure('dx')(subdomain_data=subdomains)

Could you guide me on this or point to appropriate resources?

#100

So, first of all, you only need the entire mesh.
For the entire mesh, you should be able to save the mesh with the corresponding subdomains as shown in the post above and define the appropriate measures: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

#101

I did see the above comment. I think that is based off of a premise that the mesh is read from a Gmsh.-generated one. For instance,

import meshio
msh = meshio.read("filename.inp") 
for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data

wouldn’t work in my case (as I have two CellBlocks each with type “tetra” and pointing to a separate subdomain). Also msh.cell_data_dict is an empty dictionary in my case. That is why in order to write the entire mesh I have to do

completeMesh = meshio.Mesh(msh.points, cells = {"tetra":np.vstack((msh.cells[0].data, msh.cells[1].data)) })
meshio.write("entireMesh.xdmf", completeMesh)
#102

But then you do not save any subdomain information, so how should you Get any data to Measure. The code above is easily generalized to multiple cell blocks, by stacking inside the for loop.

#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.

Marking surfaces of a 3D mesh (mshr) to apply b.c
Converting gmsh to xdmf
Marking surfaces of a 3D mesh (mshr) to apply b.c
#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?

Code from Poisson subdomain problem tutorial not working
#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:

2 Likes
Single 2D Truss mesh
Meshio Convert missing files
Converting 3D mesh to xdmf/h5
Creating subdomains in mshr with an stl file
Importing 3D .msh file into Fenics as mesh
Navier Stokes equiation 3D
Why to update from dolfin-convert to meshio?
Converting mesh from XML format to hd5 format
#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)
#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