# How to define different materials / importet 3D geometry (gmsh)

hi

i have a problem with defining different materials for two different volumes.
i generated a mesh (with gmsh) which includes 2 beams. beam 1 is clamped on the left side (physical surface 1) beam 2 is clamped on the right side (physical surface 2), and they do not touch. I defined the physical volume 1 for the left beam and the physical volume 2 for the right beam. how can i use this to create two different materials for the beams?

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)],
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)

# Scaled variables
E1 = Constant(210* 1e9)
E2 = Constant(70* 1e9)
nu1 = Constant(0.3)
nu2 = Constant(0.3)
mu1 = E1/2/(1+nu1)
mu2 = E2/2/(1+nu2)
lambda1_ = E1*nu1/(1+nu1)/(1-2*nu1)
lambda2_ = E2*nu2/(1+nu2)/(1-2*nu2)
rho1 = 7800
rho2 = 2700
g = -9.81

# Create mesh and define function space
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
File("contact_beams_mesh.pvd").write(mesh)

mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("mf.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("beam_facets.pvd").write(mf)

V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

bcl = DirichletBC(V, Constant((0, 0, 0)), mf, 1)
bcr = DirichletBC(V, Constant((0, 0, 0)), mf, 2)
bcs = [bcl, bcr]

# Define strain and stress
def epsilon(u):

def sigma(u):
return lambda1_*nabla_div(u)*Identity(d) + 2*mu1*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, rho1*g+F),)
#T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx # + dot(T, v)*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs)``````

Following Dolfinx discontinous expression, this code can be modified for the syntax in dolfin as follows:

``````import dolfin

mesh = dolfin.UnitCubeMesh(9, 9, 9)
V = dolfin.FunctionSpace(mesh, "DG", 0)
v = dolfin.Function(V)
x = V.tabulate_dof_coordinates()
mf = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
class Obstacle(dolfin.SubDomain):
def inside(self, x, on_boundary):
return (dolfin.between(x[1], (0.2, 0.7)) and
dolfin.between(x[0], (0.3, 0.7)))
Obstacle().mark(mf, 2)
dolfin.XDMFFile("CF.xdmf").write(mf)
for i in range(x.shape[0]):
if mf.array()[i] == 2:
v.vector().vec().setValueLocal(i, 2)
else:
v.vector().vec().setValueLocal(i, 1)

dolfin.XDMFFile(dolfin.MPI.comm_world, "DG.xdmf").write_checkpoint(v, "v", 0)
``````
2 Likes

I tried to implement your code, but there is a thing I don’t understand:
Don’t I also have to know beforehand, where the inclusions are?
I tried it like this:

``````> class Obstacle(SubDomain):
>     def inside(self, x, on_boundary):
>         return (between(x[1], (0, np.max(Mesh_Koord))) and
>                 between(x[0], (0, np.max(Mesh_Koord))))
>
> Obstacle().mark(subdomains, 2)
``````

But then everything is subdomain(2).

I used this to create the mesh function marking each cell, as I do not have data from an XDMF or HDF5 file to base my geometry on. As you have loaded your mesh function from file, you should only consider.
where I in my case have the subdomain 2 that has the material parameter 2, and the rest of the domain having the material parameter 1.

Thank you very much, I think that this is what I wanted!

Hello dokken,

I am trying something similar using dolfinx but encounter a small problem. I used gmsh to write the subdomains in the mesh using

for key in mesh.cell_data_dict[“gmsh:physical”].keys():
if key == “triangle”:
triangle_data = mesh.cell_data_dict[“gmsh:physical”][key]
elif key == “tetra”:
tetra_data = mesh.cell_data_dict[“gmsh:physical”][key]

tetra_mesh = meshio.Mesh(points=mesh.points, cells={“tetra”: tetra_cells},
cell_data={“volume”:[tetra_data]})
triangle_mesh =meshio.Mesh(points=mesh.points,
cells=[(“triangle”, triangle_cells)],
cell_data={“boundary”:[triangle_data]})

Afterwards I read the meashtags as follows

with XDMFFile(MPI.COMM_WORLD,“mesh.xdmf”, “r”) as infile:
volumes = infile.read_meshtags(mesh, name = “Grid”)

Up to this everything works fine, however while writing my material property using something like

Its getting complicated, the function V and x.shape[0] seem to refer the Points of the mesh but the meshtags I read in seem to refer the Cells. Thus this approach fails, so I am wondering, is there any way to read/write the meshtags in relation to the points? Or am I missing the point completely and this approach is simply not the right one?

Kind regards
Philipp

Please supply a minimal working example, as it is not clear to me what you would like to do.
The post above, similar to Dolfinx discontinous expression shows how to write discontinuous cell data to a dolfin/dolfinx function. This means that you should use a DG 0 function space to represent the data.
If you are trying to write for anything else than cells, the function space and mesh function has to be adapted.

Sorry seems like I thought it would be kind of trivial

Well okay at first I generate my generate a mesh using pygmsh using the following Code
import numpy as np
import pygmsh as pg
import meshio

``````geom = pg.opencascade.Geometry()
iso = geom.add_rectangle([0, 0, 0.0], 10e-3, 2.5e-3, char_length = .1e-3)
air = geom.add_rectangle([0.0, 0,  5e-3], 10e-3, 2.5e-3, char_length = .1e-3)
axis = [0,0,5e-3]
isok = geom.extrude(iso, axis)
axis_air = [0,0,1e-3]
airk = geom.extrude(air, axis_air)
sens = geom.add_rectangle([2.1e-3, 0 , 4.2e-3],0.8e-3, 0.8e-3, char_length = .1e-3)
axis_sens = [0,0,0.8e-3]
sensk = geom.extrude(sens,axis_sens)
geom.add_raw_code('BooleanDifference{ Volume{1}; Delete; }{ Volume{3}; }')
geom.add_raw_code('Physical Surface(3) = {39, 36};') #outside
geom.add_raw_code('Physical Surface(5) = {35, 40};') #inflow

mesh = pg.generate_mesh(geom)

tetra_cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "tetra"]))

for key in mesh.cell_data_dict["gmsh:physical"].keys():
if key == "tetra":
tetra_data = mesh.cell_data_dict["gmsh:physical"][key]

tetra_mesh = meshio.Mesh(points=mesh.points, cells={"tetra": tetra_cells},
cell_data={"volume":[tetra_data]})
meshio.write("mesh.xdmf", tetra_mesh)
``````

And then I´ll try to read the mesh and write the conductivities to the three different volumes

``````from dolfinx import Mesh, geometry, Function, DirichletBC, Constant, solve, FunctionSpace, VectorFunctionSpace
from dolfinx.io import XDMFFile, VTKFile
import meshio
from ufl import TrialFunction, TestFunction, dot, grad, dx, inner, Measure
from mpi4py import MPI

mesh = Mesh
with XDMFFile(MPI.COMM_WORLD,"mesh.xdmf", "r") as infile:
vols = infile.read_meshtags(mesh, name = "Grid")
element =("CG", 1)
F = FunctionSpace(mesh, element)
x = F.tabulate_dof_coordinates()
cond = Function(F)
for i in range(x.shape[0]):
if vols.values[i] == 6:
cond.vector.setValueLocal(i,1)
elif vols.values[i] == 8:
cond.vector.setValueLocal(i,2)
else:
cond.vector.setValueLocal(i,3)
``````

Afterwards I want to solve a heat equation with this system, unfortunately everything has an cond of 1.
I removed quite a lot of the Code, I hope this is enough for you to see the mistake i made.

This should be DG-0 as the mesh tags for tetrahedrons are cell based, as I said earlier.

Okay thanks, that helps quite a lot. Guess I´ll have to modify the saved mesh too, as I am getting an exception when changing the element type

``UFLException: Order "0" invalid for "Lagrange" finite element.``

It should not be Lagrange, it should be “DG”, Which are discontinuous lagrangian elements (piecewise constants per element).

Oh sorry got a typo in there. Now its working, thanks a lot

Okay,

The calculations are running so far. Unfortunately I did encounter a new problem.
When using the following code

``````from dolfinx import Mesh, geometry, Function, DirichletBC, Constant, solve, FunctionSpace, VectorFunctionSpace
from dolfinx.io import XDMFFile, VTKFile
from dolfinx.fem import assemble
from ufl import TrialFunction, TestFunction, dot, grad, dx, ds, inner, Measure, nabla_grad
from mpi4py import MPI
import meshio
import pygmsh as pg

lam = [0.195, 3, 26e-3]
pth = 5.6
rho = 1.18
zair = 1e-3
esize =[2e-4,.5e-4,.5e-4]
>#Meshing
iso = geom.add_rectangle([0, 0, 0.0], 5e-3, 2.5e-3, char_length = esize[0])
air = geom.add_rectangle([0.0, 0,  5e-3], 5e-3, 2.5e-3, char_length = esize[1])
axis = [0,0,5e-3]
isok = geom.extrude(iso, axis)
axis_air = [0,0,zair]
airk = geom.extrude(air, axis_air)
sens = geom.add_rectangle([2.1e-3, 0 , 4.2e-3],0.8e-3, 0.8e-3, char_length = esize[2])
axis_sens = [0,0,0.8e-3]
sensk = geom.extrude(sens,axis_sens)
geom.add_raw_code('BooleanDifference{ Volume{1}; Delete; }{ Volume{3}; }')
geom.add_raw_code('Physical Surface(6) = {33};') #out Mat
geom.add_raw_code('Physical Surface(8) = {36};') #outside Mat

mesh = pg.generate_mesh(geom, verbose=False)
tetra_cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "tetra"]))

triangle_cells = np.vstack(np.array([cells.data for cells in mesh.cells
if cells.type == "triangle"]))

for key in mesh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = mesh.cell_data_dict["gmsh:physical"][key]

tetra_mesh = meshio.Mesh(points=mesh.points, cells={"tetra": tetra_cells},
cell_data={"volume":[tetra_data]})
triangle_mesh =meshio.Mesh(points=mesh.points,
cells=[("triangle", triangle_cells)],
cell_data={"boundary":[triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

>#calculating
ptherm = pth/(1000*1.6e-3*0.8e-3*0.8e-3)
#Mesh lesen
mesh = Mesh
with XDMFFile(MPI.COMM_WORLD,"mesh.xdmf", "r") as infile:
vols = infile.read_meshtags(mesh, name = "Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "mf.xdmf", "r") as xdmf_infile:
element =("DG", 0)
F = FunctionSpace(mesh, element)
element =("CG", 2)
V = FunctionSpace(mesh, element)
#Definde boundary condition
ind = [1,2,3,4,5,6,7,8]
dom =[]
#Ordnet Knoten den Flächen der RB zu
for element in ind:
dom.append(mt.indices[mt.values == element])
u_D = Function(V)
u_D.vector.set(0)
bc = []
#RB schreiben
for element in dom:
bc.append(dolfinx.DirichletBC(u_D, dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, element)))
body = []
ind =[9,10,11]
#Knoten und deren Koord holen
x = F.tabulate_dof_coordinates()
cond = Function(F)
vel = Function(F)
f = Function(F)
for i in range(x.shape[0]):
if vols.values[i] == 9:
cond.vector.setValueLocal(i,lam[0])
f.vector.setValueLocal(i,0)
elif vols.values[i] == 11:
cond.vector.setValueLocal(i,lam[1])
f.vector.setValueLocal(i,ptherm)
else:
cond.vector.setValueLocal(i,lam[2])
f.vector.setValueLocal(i,0)
#Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
L = f*v*dx
#Compute solution
u = Function(V)
solve(a==L, u, bc)
#save
m = u*dx2
t = assemble.assemble_scalar(m)
overtemp = t/(.8e-3*0.8e-3*0.8e-3)
vtkfile = VTKFile('wuek_automata.pvd')
vtkfile.write(u)
``````

Sorry the code is a bit messy, i chopped huge parts away to simplify it.

Okay now the deal, I need a quite fine mesh, at least for esize[1] and esize[2], something like 0.1e-4. Pc can handle it no problem, but with such small element sizes the calcutation somehow “crashes”.

Crashing means, the calculated temperature is something like 5e-4 instead of lets say 1.5K

I tried several things and if I use a constant thermal conductivity for all bodys this problem does not occur. The problem seems to be a combination of small elements and the huge difference in thermal conductivity betwen “sensor” and “air”.

Is there any way to avoid, whatever is happening there?
Regards

Philipp

Please use ``` for appropriate code formatting (indentation and such). This is the markdown syntax, which makes it possible to run your code.
Note that you should be able to scale your problem such that your mesh size does not have to be of order 10^{-4}.

Oh sorry for false fomatting.

I do have scaled the problem allready, the code is unfortunately at the workstation, but the problem occurs anyway. The only difference is an factor of 1000 in the esize. (from meter to mm)