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

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

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

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.

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

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

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)

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

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.

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)

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

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)

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

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.

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.



You need to make your code reproducible for anyone to be able to help you. Please supply a minimal script for generating the mesh.

Sorry, Sir,

Here is the .goe file(since the .msh file is too big in size, I am sending .geo file)

// Gmsh project created by Jørgen S. Dokken 2020
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1, 1, .285};
Box(2) = {0.4, 0.4, 0.285, .2, .2, .03};
Box(3) = {0., 0, 0.315, 1, 1, .01};
Box(4) = {0.0, 0.0, 0.325, 1, 1, 1};
v() = BooleanFragments{ Volume{1}; Delete;}{ Volume{2,3,4}; Delete; };
Physical Volume("air", 1) ={#v()};
Physical Volume("hBN", 2) ={#v()-1};
Physical Volume("cavity", 3) ={#v()-2};
Physical Volume("substrate", 4) ={#v()-3};
Physical Surface("backgate", 5) = {29};
Physical Surface("patterned_gate", 6) = {18, 23};

Try editing the opacity values of your ParaView colormap for the volume plot. With the default values, i.e.
default
The volume marked “1” is shown as completely transparent. If you edit your colormap to make all values opaque, i.e.
opaque
you will see that all volumes are present.

1 Like

To add to @conpierce8 excellent comment, I would suggest changing the colormap to Viridis (matplotlib), as it is a much better color map than Cool to Warm, see: Introduction to the viridis color maps

1 Like

Thank you so much. it is working now.

1 Like

Hello. I have tried this code but I get the following error:
ModuleNotFoundError: No module named ‘gmsh’

I have installed gmsh in fenicsproject as well as in my system, but the error still persists.