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

mesh
169

One little follow-up question, somewhat unrelated to the main thread…
I’ve used mesh refinement in gmsh recently and then “Recombine 2D” more by accident and it resulted in a mesh of mixed cell types (triangles and quads).
Are there any considerations to prefer one over the other or even mix them?

170

Dolfin does not support Mixed type meshes. (It is in the long term plans of dolfinx).

1 Like
171

Hello everyone,

I was working with the code given here earlier by @dokken


import meshio
msh = meshio.read("cylinder.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})
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)

And am receiving the following error code.

Traceback (most recent call last):
  File "convert.py", line 25, in <module>
    cells=[("triangle", triangle_cells)],
NameError: name 'triangle_cells' is not defined

File for the cylinder.msh was created on GMSH with version 4.8.4

The “msh” files was to large for this post but you can find it at the following link.

Not sure how the triangles are not defined as the msh is using triangles? What am I missing here?

172

Your msh file does not contain any physical markers (i.e Physical Line/Volume/Surface), thus you cannot read in cell_data.

I would suggest using the following code for reading meshes (from Accessing and marking imported boundaries - #8 by dokken). Note that here I have defined a physical surface and physical curve.

import meshio


msh = meshio.read("cylinder.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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh


tetra_mesh = create_mesh(msh, "tetra")
triangle_mesh = create_mesh(msh, "triangle")
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)
1 Like
173

@dokken

I appreciate your help and I think that got me one step closer, unfortunately I have hit another issue.

I am using the code you provided from above and have updated my design to have surfaces and volumes and more.

I did take a look at the original post that you directed me to and that one does not cover this error I am now getting.

Traceback (most recent call last):
  File "convert.py", line 11, in <module>
    msh = meshio.read("cylinderv5.msh")
  File "/usr/local/lib/python3.6/dist-packages/meshio/_helpers.py", line 69, in read
    return reader_map[file_format](filename)
  File "/usr/local/lib/python3.6/dist-packages/meshio/gmsh/main.py", line 19, in read
    mesh = read_buffer(f)
  File "/usr/local/lib/python3.6/dist-packages/meshio/gmsh/main.py", line 48, in read_buffer
    return reader.read_buffer(f, is_ascii, data_size)
  File "/usr/local/lib/python3.6/dist-packages/meshio/gmsh/_gmsh41.py", line 73, in read_buffer
    points, point_tags, point_entities = _read_nodes(f, is_ascii, data_size)
  File "/usr/local/lib/python3.6/dist-packages/meshio/gmsh/_gmsh41.py", line 174, in _read_nodes
    raise ReadError("parametric nodes not implemented")
meshio._exceptions.ReadError: parametric nodes not implemented

This error does not seem well documented on the web and the few things I have found refer to using abaqus mesh as the mesher and not GMSH.

Any idea where this error is coming from. I do realize that the mesh is really unrefined but was less worried about resolutions and more worried about converting it.

The new mesh file can be found here

Also here is the geo files if that is more helpful

Thank you much,

1 Like
174

As far as I can tell you would like to have a 3D mesh (tetrahedrons). Then you also need to mark the volumes:

// Gmsh project created on Tue May 18 16:34:55 2021
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 00, 00, 0.6, 0, 2*Pi};
//+
Circle(2) = {0, 00, 00, 0.4, 0, 2*Pi};
//+
Curve Loop(1) = {1};
//+
Curve Loop(2) = {2};
//+
Plane Surface(1) = {1, 2};
//+
Extrude {0, 0, 10} {
  Curve{2}; Curve{1}; Layers {1}; Recombine;
}
//+
Curve Loop(5) = {6};
//+
Curve Loop(6) = {4};
//+
Plane Surface(4) = {5, 6};
//+
Physical Surface("Brime", 7) = {2};
//+
Physical Surface("AirSide", 8) = {3};
//+
Physical Surface("Top", 9) = {1};
//+
Physical Surface("Bottom", 10) = {4};
//+
Surface Loop(1) = {2, 4, 3, 1};
//+
Volume(1) = {1};
//+
Physical Volume(11) = {1};

This file, converted to msh with gmsh 4.6.0,
reads nicely with:

import meshio


msh = meshio.read("cylinderv5.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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh


tetra_mesh = create_mesh(msh, "tetra")
triangle_mesh = create_mesh(msh, "triangle")
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)
2 Likes
175 
Hello everyone, How I can use GMSH file on FEniCS (1D mesh) ? 
my problem was divided into two regions, each region has a different equation, how I can divide the two regions in fenics ?. 

I am using the follows code, but unfortunately I got an error.

 ##########################################################
################### Python code #############################
###########################################################
from fenics import *
import meshio

msh = meshio.read('QW.msh')
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"line": msh.cells["line"]}, cell_data={"line": {"name_to_read":
                        msh.cell_data["line"]["gmsh:physical"]}}))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 1)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mf)


print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))

######################################
################################# Error #############################################
############################################
Traceback (most recent call last):
File “/home/jaouane/Desktop/works/gooog/gmsh.py”, line 17, in
print(assemble(Constant(1)*dx))
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dolfin/fem/assembling.py”, line 198, in assemble
dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dolfin/fem/assembling.py”, line 56, in _create_dolfin_form
return Form(form,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dolfin/fem/form.py”, line 43, in init
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dolfin/jit/jit.py”, line 47, in mpi_jit
return local_jit(*args, **kwargs)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dolfin/jit/jit.py”, line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/ffc/jitcompiler.py”, line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dijitso/jit.py”, line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/ffc/jitcompiler.py”, line 76, in jit_generate
dep_module_name = jit(dep, parameters, indirect=True)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/ffc/jitcompiler.py”, line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dijitso/jit.py”, line 177, in jit
build_shared_library(signature, header, source, dependencies,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dijitso/build.py”, line 153, in build_shared_library
status, output = get_status_output(cmd)
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/site-packages/dijitso/system.py”, line 40, in _get_status_output_subprocess
pipe = subprocess.Popen(cmd, shell=False, cwd=cwd, env=env,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/subprocess.py”, line 951, in init
self._execute_child(args, executable, preexec_fn, close_fds,
File “/home/jaouane/anaconda3/envs/spyder/lib/python3.9/subprocess.py”, line 1821, in _execute_child
raise child_exception_type(errno_num, err_msg, err_filename)
FileNotFoundError: [Errno 2] No such file or directory: ‘c++’

Process finished with exit code 1
176

What happens if you try to use a built in mesh (i.e. IntervalMesh)? Do you obtain the same error?

Here is a minimal example of reading in a 1D mesh that works with dolfin:
File test.geo

//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Physical Curve(1) = {1};

from fenics import *
import meshio


def create_mesh_1D(mesh, cell_type):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh

mesh_out = create_mesh_1D(meshio.read("test.msh"), "line")
meshio.write("mesh.xdmf", mesh_out)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 1)
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mf)

print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
178

Thank you for your response, I try it with the above code, unfortunately I got same error.
I am going to try it by using IntervalMesh.

179

Then it Seems like you have an installation issue. I suspected this from the beginning as your error is in the assembly, not reading in the mesh.

1 Like
180

Thank you very much Dokken, You are right , I have an installation issue, now it’s working fine.

181

Thank you. This works perfectly.
If anyone is interested in adapting this for a 2D mesh, only minor modifications are necessary, and the code snippet is as follows:

import gmsh
import pygmsh
import meshio
import numpy
from dolfin import *
msh = meshio.read("mesh2d.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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mf.xdmf", line_mesh)

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


mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(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("dolfinmesh.pvd").write(mesh)
File("dolfincellfunc.pvd").write(mf)
1 Like
182

Is the above code will work for 3d domain.?

183

Simply set prune_z=False in the create_mesh function, as long as replacing line with triangle and triangle with tetrahedron. The MeshValueCollection should have dimension of mesh.topology().dim() if you read in cell markers, and mesh.topology().dim()-1 if you read in facet markers.

184

Thank you for your reply,
Is it I modified correctly as you said.?

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)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh

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

triangle_mesh = create_mesh(msh, "triangle", False)
tetra_mesh = create_mesh(msh, "tetra", False)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)
#print(msh.cell_data_dict)
from dolfin import *
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["optimize"] = True

mesh=Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds = Measure("ds", domain=mesh, subdomain_data=mf)

dx=Measure("dx", domain=mesh, subdomain_data=cf)
185

You can easily verify this yourself by looking at the XDMF file in Paraview, or integrate over the subdomains (for instance the volumes or surface areas, using

dx_C = Measure("dx", domain=mesh, subdomain_data=cf)
print(assemble(1*dx_C(value_of_domain))

and similarly for ds

The code you have supplied is also not properly formatted (missing indentation).

186

Thank you for your help.
Here is the code with proper indentation.

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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh



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


triangle_mesh = create_mesh(msh, "triangle", False)
tetra_mesh = create_mesh(msh, "tetra", False)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh) 
#print(msh.cell_data_dict)
from dolfin import *
parameters["allow_extrapolation"] = True 
parameters["form_compiler"]["optimize"] = True 

mesh=Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)


mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")   
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc) 
  
ds = Measure("ds", domain=mesh, subdomain_data=mf)

dx=Measure("dx", domain=mesh, subdomain_data=cf)
1 Like
187

Yes, I have verified with paraview by reading xdmf file. when I have seen the file in wireframe or surface it shows correctly. But, if click for the volume it will not show anything

188

There are two different files to consider here

  1. mesh.xdmf
  2. mf.xdmf

The mesh.xdmf file contains the cells (tetrahedrons) and corresponding markers, while mf.xdmf will show you the facets of the mesh, colored with the corresponding markers.

189

Yes, I considered mesh.xdmf only.
In the following figures, you can see that, the file is perfectly matching with the gmsh file when I will see it as a wireframe or surface but not in the case of volume.

surface
surface1144×795 16 KB

wireframe
wireframe1144×795 70.6 KB

volume
volume1144×795 18.9 KB

1 Like