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?

msh = meshio.read("./contact_beams.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)
# 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
load = -20
F = load
g = -9.81
# Create mesh and define function space
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
File("contact_beams_mesh.pvd").write(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("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):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(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)
```