# How to import xdmf mesh generated by gmsh into fenics

@dokken

Maybe you can help me out.
I have created simple square geometry in Gmsh and converted the .msh file into xdmf. When i import it in fenics using your code I’m getting different results compare to the simple command "UnitSquareMesh(30,30). I’m trying tonsolve the simple linear elastic problem. Maybe you can help me out. How i can import the external 2d mesh into fenics?

is it the correct code to import 2d mesh?

mesh = Mesh()

with XDMFFile(“lavender.xdmf”) as infile:

mvc = MeshValueCollection(“size_t”, mesh, 1)

with XDMFFile(“bound.xdmf”) as infile:

mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

1 Like

Without an example illustrating the difference you are experiencing, I cant help you.

When I use Fenics built in mesh UnitSquareMesh(30,30) I’m getting dispacment 0.3.

When i create geometry of 30x30 in gmsh and import the mesh in Fenics the resulting displacement is 8.3e-9.

I’m using clamped boundary conditions at the bottom and displacement at top.

Maybe I’m importing the mesh wrongly. Maybe you can help me out how to import the 2d mesh in Fenics with some descriptions?

Please provide the code that you are running (and the mesh, or the scripts that generated the mesh).

There are plenty of topics covering this, for instance

and many more

here you go


to generate the square mesh:

//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {30, 0, 0, 1.0};
//+
Point(3) = {30, 30, 0, 1.0};
//+
Point(4) = {0, 30, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {3, 4, 1, 2};
//+
Plane Surface(1) = {1};
//+
Physical Curve(5) = {1};
//+
Physical Curve(6) = {2};
//+
Physical Curve(7) = {3};
//+
Physical Curve(8) = {4};
//+
Physical Surface(9) = {1};


code to generate the mesh in xdmf

import meshio
import numpy
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)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
return out_mesh

line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
meshio.write(“bound.xdmf”, line_mesh)

triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)
meshio.write(“lavender.xdmf”, triangle_mesh)


Fenics code:

from dolfin import *
from fenics import *
from ufl import nabla_grad
from ufl import nabla_div
import numpy as np

mesh=Mesh()
with XDMFFile(‘lavender.xdmf’) as infile:
mvc = MeshValueCollection(“size_t”, mesh, 2)
with XDMFFile(‘bound.xdmf’) as infile:
sub = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# Strain function

def epsilon(u):

# Stress function
def sigma(u):
return lmbda_*div(u)*Identity(2) + 2*mu*epsilon(u)

# Define material properties

E = Constant(100000)
nu = Constant(0.3)
mu = E/(2*(1+nu))
lmbda_ = E*nu/((1+nu)*(1-2*nu))
P2 = VectorElement("Lagrange", mesh.ufl_cell(), flag_quad)
TH = P2
W = FunctionSpace(mesh, TH)

# Define traction on the boundary and body forces
T = Constant((0.0, 0.0))
B = Constant((0.0, 0.0))

u = Function(W)
du = TrialFunction(W)
v = TestFunction(W)

# Boundary Conditions
def clamped_boundary(x, on_boundary):
return on_boundary and x[1] < tol

bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary)
tol=1E-5
top = CompiledSubDomain('on_boundary && near(x[1], 1, tol)', tol=1E-5)
applied_disp = 0.5
bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top)

bcs = [bc1,bc2]

a = inner(sigma(du), epsilon(v))*dx
l = dot(B, v)*dx + inner(T, v)*ds(1)
A_ass, L_ass = assemble_system(a, l, bcs)
solve(A_ass, u.vector(), L_ass)
plot(u, title="Displacement", mode="displacement")
to_print = True

file_results = XDMFFile("+-elasticity_reet.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

file_results.write(u,0.)


Please edit your posts and make sure They use 3x encapsulation surrounding codes or files, i.e.


//+
Point(….)



and

python
from dolfin import *
#….



here you go


to generate the square mesh:

//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {30, 0, 0, 1.0};
//+
Point(3) = {30, 30, 0, 1.0};
//+
Point(4) = {0, 30, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {3, 4, 1, 2};
//+
Plane Surface(1) = {1};
//+
Physical Curve(5) = {1};
//+
Physical Curve(6) = {2};
//+
Physical Curve(7) = {3};
//+
Physical Curve(8) = {4};
//+
Physical Surface(9) = {1};


code to generate the mesh in xdmf

import meshio
import numpy
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)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
return out_mesh

line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
meshio.write(“bound.xdmf”, line_mesh)

triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)
meshio.write(“lavender.xdmf”, triangle_mesh)


Fenics code:

from dolfin import *
from fenics import *
from ufl import nabla_grad
from ufl import nabla_div
import numpy as np

mesh=Mesh()
with XDMFFile(‘lavender.xdmf’) as infile:
mvc = MeshValueCollection(“size_t”, mesh, 2)
with XDMFFile(‘bound.xdmf’) as infile:
sub = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# Strain function

def epsilon(u):

# Stress function
def sigma(u):
return lmbda_*div(u)*Identity(2) + 2*mu*epsilon(u)

# Define material properties

E = Constant(100000)
nu = Constant(0.3)
mu = E/(2*(1+nu))
lmbda_ = E*nu/((1+nu)*(1-2*nu))
P2 = VectorElement("Lagrange", mesh.ufl_cell(), flag_quad)
TH = P2
W = FunctionSpace(mesh, TH)

# Define traction on the boundary and body forces
T = Constant((0.0, 0.0))
B = Constant((0.0, 0.0))

u = Function(W)
du = TrialFunction(W)
v = TestFunction(W)

# Boundary Conditions
def clamped_boundary(x, on_boundary):
return on_boundary and x[1] < tol

bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary)
tol=1E-5
top = CompiledSubDomain('on_boundary && near(x[1], 1, tol)', tol=1E-5)
applied_disp = 0.5
bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top)

bcs = [bc1,bc2]

a = inner(sigma(du), epsilon(v))*dx
l = dot(B, v)*dx + inner(T, v)*ds(1)
A_ass, L_ass = assemble_system(a, l, bcs)
solve(A_ass, u.vector(), L_ass)
plot(u, title="Displacement", mode="displacement")
to_print = True

file_results = XDMFFile("+-elasticity_reet.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

file_results.write(u,0.)


It is quite obvious as you have hardcoded the top boundary condition, which assumes that your mesh is a unit square, while what you actually have in your mesh is a [0,30]\times[0,30] mesh.
Thus using

top = CompiledSubDomain('on_boundary && near(x[1], 30, tol)', tol=1E-5)


Here is a cleaned up version of your code:


import matplotlib.pyplot as plt
from dolfin import *
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)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("bound.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("lavender.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("lavender.xdmf") as infile:
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("bound.xdmf") as infile:
sub = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# Strain function

def epsilon(u):

# Stress function

def sigma(u):
return lmbda_*div(u)*Identity(2) + 2*mu*epsilon(u)

# Define material properties

E = Constant(100000)
nu = Constant(0.3)
mu = E/(2*(1+nu))
lmbda_ = E*nu/((1+nu)*(1-2*nu))
P2 = VectorElement("Lagrange", mesh.ufl_cell(), flag_quad)
TH = P2
W = FunctionSpace(mesh, TH)

# Define traction on the boundary and body forces
T = Constant((0.0, 0.0))
B = Constant((0.0, 0.0))

u = Function(W)
du = TrialFunction(W)
v = TestFunction(W)

# Boundary Conditions

def clamped_boundary(x, on_boundary):
return on_boundary and x[1] < tol

bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary)
tol = 1E-5
top = CompiledSubDomain('on_boundary && near(x[1], 30, tol)', tol=1E-5)
applied_disp = 0.5
bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top)

bcs = [bc1, bc2]

a = inner(sigma(du), epsilon(v))*dx
l = dot(B, v)*dx + inner(T, v)*ds(1)
A_ass, L_ass = assemble_system(a, l, bcs)
solve(A_ass, u.vector(), L_ass)

plot(u, title="Displacement", mode="displacement")
plt.savefig("u.png")

file_results = XDMFFile("+-elasticity_reet.xdmf")
file_results.write(u, 0.)
file_results.close()


Note that I’ve made sure that the indentation of all functions are correct, which wasn’t the case in your code

Also note that you could use the marked boundaries from the gmsh mesh (“bound.xdmf”) to set boundary conditions.

@dokken thanks for your kind help. I really appreciate that.

A little modification. I have this image domain in which i have 2 different materials the circles in domain have different material. how can I give different values of E to different domains in fenics? The applied boundary conditions are same as above:

circles have E = 1e6
and othe part has E= 2e3

def clamped_boundary(x, on_boundary):
return on_boundary and x[1] < tol

bc1 = DirichletBC(W, Constant((0, 0)), clamped_boundary)
tol = 1E-5
top = CompiledSubDomain('on_boundary && near(x[1], 30, tol)', tol=1E-5)
applied_disp = 0.5
bc2 = DirichletBC(W.sub(1), Constant((applied_disp/2.0)), top)

bcs = [bc1, bc2]


Everything is same as above i have attached the geometry below:

One question do i have to change this xdmf file?


import matplotlib.pyplot as plt
from dolfin import *
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)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("bound.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("lavender.xdmf", triangle_mesh)



As long as you have Added physical tags in Gmsh, They are now saved in the xdmf files.

Let’s suppose the image given above is meshed using GMSH. The circles in images are extracted using OpenCv Contours detection. And I add physical tags to only boundary than whats the procedure to detect subdomains. In simple I use Gmsh for meshing the Image which has different domains how will fenics detect the subdoamins?

You can loop over the read in facets, find the connecting cells on each side and mark them.

However, I would strongly advice you to marker the surfaces/subdomains (cell volumes) in Gmsh/your mesh generator, as it is alot more expensive to do recursive searches through neighbouring cells if you do not have specific geometrical information.

It would be amazing if you can elaborate a little bit how i can proceed in GmSh for this. I have marked the boundary only and gave the whole image an only Physical surface. As circles are in the image I can’t give them physical tag. How can we proceed in case of Images ?

I do not understand what you mean by this being an image. GMSH is a mesh generator, that has a graphical user interface where you can manually tag each subdomain.

You could also do this with Python coding, see
https://jsdokken.com/dolfinx-tutorial/chapter3/em.html#meshing-a-complex-structure-with-subdomains

What I want to say is that I want to mesh a Image. The above Image is just a dummy. Actual Image doesn’t have circle. It has quite random geometry which is quite complex to define. In this way what’s the procedure for physical tags and import the mesh to Fenics with white background has different material and black has different material. Inhave attached the Image below:

You could either:
1, Use a graphical user interface of GMSH to give each object a physical surface tag.
2. Write a recursive script looping through
a. For each interface of an object find each connected cell.
b. Check each connected cells loop through neighbours (going through each of their respective facets, making sure you never cross the interface again), marking each side with a single number.

If I use Option 1 can I give physical tag to random geometries too?

How to proceed in this case? Any tutorial?

2nd option is quite complex.

As i do not know what meshing tool you use to go from an image to a mesh, it is hard to give you any specific tips.

If you use Gmsh for this, there might be a tutorial on their webpage:
https://gmsh.info/

The second option doesn’t take many lines of Python to program. See for instance: