# How to import xdmf mesh generated by gmsh into fenics

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

I have cell neighbours of each element. How can I differentiate the two and define different material.

As I told you

So start at a marked interface, find the two cells connected to that facet.
Give each of the cells a unique number (say cell0 is marked with 5 cell1 is marked with 7).
in sudo-code this is what you need to do for each cell:

def mark_neigbours(cell, tag, interface_markers, volume_markers, visited_cells):
if c in visited_cells:
return
c_facets = facets(cell)
for facet in c_facets:
if facet in (interface markers):
pass
else
connected_cells=cells(facet)
for c in connected_cells:
if c!= cell
volume_markers[c] = tag
visited_cells.append(c)
mark_neigbours(c, tag, interface_markers, visited_cells)


As I said, this is a recursive algorithm, which is going to be quite complex (and I might have made mistakes in the pseudo code). You should really let the mesh generator handle the tagging of cells.

after marking neighbours how will define the different material for each cell?

That is the whole point of volume_markers, which should be a MeshFunction, which you can modify the array of, with the appropriate tags, such that you can compute integrals over subdomains.