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

mesh
165

Dear,

I believe I followed all the procedures written in the posts above, but I am not able to generate the .xdmf file

note: I made the mesh in Gmsh version 4.8.3 on macOS.

Follow the code below and also the file .geo:

import meshio
msh = meshio.read("two_cantilever.msh")
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]
for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                         cell_data={"Volume":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"Boundary":[triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

See the error:

 File "/Users/xyz/Desktop/Topopt_final/3dmodel/gmsh_to_xdmf.py", line 66, in <module>
    tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},

  File "/Users/xyz/anaconda3/envs/fenics2/lib/python3.8/site-packages/meshio/_mesh.py", line 77, in __init__
    raise ValueError(

ValueError: Incompatible cell data. Cell block 0 has length 214, but corresponding cell data Volume item has length 66420.

two_cantilever.geo:

// Gmsh project created on Wed Apr 28 09:20:18 2021
SetFactory("OpenCASCADE");
//+
Box(1) = {0, 0, 0, 3, 0.01, 1};
//+
Box(2) = {-2, -0.03, 0, 3, 0.01, 1};
//+
Cylinder(3) = {0.5, -0.04, 0.25, 0, 0.05, 0, 0.05, 2*Pi};
//+
Cylinder(4) = {0.5, -0.04, 0.75, 0, 0.05, 0, 0.05, 2*Pi};
//+
Physical Surface(1) = {2};
//+
Physical Surface(2) = {7};
//+
Physical Surface(3) = {13};
//+
Physical Surface(4) = {16};
//+
Physical Volume(31) = {1};
//+
Physical Volume(32) = {2};
//+
Physical Volume(33) = {4};
//+
Physical Volume(34) = {3};

For my code I need to apply boundary considerations on the extreme faces of the plates and on the outer faces of the cylinders.

Edit .: The .msh file was larger than the size allowed here. Follow the file link:

can anybody help me? @dokken

166

You are not concatenating data stemming from different cell blocks.
You have not followed the instructions in:

Need help converting GMSH to FEniCS - #8 by dokken

To repost this again:

import meshio


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("two_cantilever.msh")
facet_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("facet_mesh.xdmf", facet_mesh)

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

Note that you have not used boolean operations on your geometry, resulting in another mesh that what you would like. See: Using subdomains in diffusion equation - #6 by dokken

1 Like
167

Dear Dokken,

I really appreciate your help.

Well, I made an example with a simpler geometry, a single cantilever without the “holes” (I will study more about gmsh) that follows below. The two_cantilever here didn’t work! Did you ever run there?

The new code works, however in the result it seems that both the left side BCdirichlet restriction and the right side load insertion are not being considered.

bc = DirichletBC (V, Constant ((0.0, 0.0, 0.0)), mf, 14)
ds = Measure ("ds", domain = mesh, subdomain_data = mf, subdomain_id = 13)

see the code:

from fenics import *
import meshio

L = 3.0
B = 0.01
H = 1.0

# Load
F = 30
ds_size = 0.1

# Define 3D geometry
mesh = Mesh()
with XDMFFile("untitled.xdmf") as infile:
    infile.read(mesh)

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

mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)

dx = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=15)

# Define function space and base functions
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)

# Boundary Condition
"""
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], -2.0)
left = Left()
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), left)
"""
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), mf, 14)

# Load
class LoadSurface(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], L) and \
               between(x[1], ((B-ds_size)/2.0, (B+ds_size)/2.0)) and \
               between(x[2], ((H-ds_size)/2.0, (H+ds_size)/2.0))
#edited               
loadSurface = LoadSurface()
mf.set_all(0)
loadSurface.mark(mf, 1)
ds = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=13)
t = Expression(('q_x', 'q_y', 'q_z'), degree = 3, q_x = 0.0, q_y = 0.0, q_z = -F/ds_size)

## Constitutive relation

def eps(v):
    return sym(grad(v))

E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
    dim = v.geometric_dimension()
    return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)

## Variational formulation

## Resolution

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(t, u_)*ds(1)

u = Function(V, name="Displacement")
solve(a == L, u, bc)

file_results = XDMFFile("u_new_3D.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)

.geo File:

//+
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 3, 0.01, 1};
//+
Physical Surface(13) = {2};
//+
Physical Surface(14) = {1};
//+
Physical Volume(15) = {1};

.msh File: (Using the code you provided above)

168

There are several things in your code that strikes me as odd.
First, why are you loading the mesh from another file than the volume tag, i.e.

should be

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

Secondly, you are overwriting the facet tags you are reading in

Thus you are not setting the load on any subdomain, as you removed the tags 13 by using set_all. Before posting more code, please carefully go through your code line for line, and check that you get the expected result for each operation.

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