# Two materials and Mechanism

Dear,

I have two doubts:

Dear,

I have two doubts:

1st question - According to the example below, I have two volumes and I would like to consider two different types of mtarial in each volume. The model considers the individual volumes in gmsh using BooleanFragments. However, how do I consider each volume to be a material in the code below? I’m sorry, I did some research on discursse but I still can’t make this identification of two materials in the code.

2nd question - Is it possible in fenics to consider this model as a mechanism? Consider the contact between the two volumes? or consider some kind of rotation between them?

gsmh code:

``````//+
Box(1) = {0, 0, 0, 3, 1, 1};
//+
Cylinder(2) = {1, -0.2, 0.5, 0, 2.2, 0, 0.3, 2*Pi};
//+
BooleanFragments{ Volume{1}; Delete; }{ Volume{2}; Delete; }
//+
Physical Surface(20) = {11};
//+
Physical Surface(21) = {6};
//+
Physical Volume(22) = {1};
//+
Physical Volume(23) = {4, 2, 3};
``````

Example code:

``````from fenics import *
import meshio

F = 2e2

# Define 3D geometry
mesh = Mesh()
with XDMFFile("mecanismo_mesh.xdmf") as infile:
mvc2 = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mecanismo_facet.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
mvc3 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mecanismo_mesh.xdmf") as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc3)

dx1 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=22)
dx2 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=23)

# Define function space and base functions
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

# Boundary Condition
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), mf, 20)

BCs = [bc]

ds = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=21)
t = Expression(('q_x', 'q_y', 'q_z'), degree = 1, q_x = 0.0, q_y = 0.0, q_z = -F)

## Constitutive relation

def eps(v):

E = Constant(210e3)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
dim = v.geometric_dimension()
return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)

## Variational formulation

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(t, u_)*ds

A = assemble(a)
B = assemble(L)

[bc.apply(A, B) for bc in BCs]

u = Function(V, name="Displacement")
solve(A ,u.vector() ,B)

#file_results
file_results = XDMFFile("u_new_3D.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
``````

Dear @dokken,

It would be like the code below:

My question was on line 32, do I call dx1? Is correct?

I did a test by inverting the materials and the displacement result was the same as if it were all the material with the greater Em. I see that there should be a larger displacement in the load region and not equal to the displacement of a single material.

``````from fenics import *
import meshio

F = 2e2

# Define 3D geometry
mesh = Mesh()
with XDMFFile("mecanismo_mesh.xdmf") as infile:
mvc2 = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mecanismo_facet.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
mvc3 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mecanismo_mesh.xdmf") as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc3)

dx1 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=22) # box
dx2 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=23) # cilinder

# Define Material
Em = [210e3, 100]
#Em = [100, 210e3] # inverting

Vm = FunctionSpace(mesh, "DG", 0)
vm = Function(Vm)
xm = Vm.tabulate_dof_coordinates()

for i in range(xm.shape):
if cf.array()[i] == dx1:
vm.vector().vec().setValueLocal(i, Em)
else:
vm.vector().vec().setValueLocal(i, Em)

# Define function space and base functions
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

# Boundary Condition
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), mf, 20)

BCs = [bc]

ds = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=21)
t = Expression(('q_x', 'q_y', 'q_z'), degree = 1, q_x = 0.0, q_y = 0.0, q_z = -F)

## Constitutive relation

def eps(v):

#E = Constant(210e3)
E = vm
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
dim = v.geometric_dimension()
return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)

## Variational formulation

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(t, u_)*ds

A = assemble(a)
B = assemble(L)

[bc.apply(A, B) for bc in BCs]

u = Function(V, name="Displacement")
solve(A ,u.vector() ,B)

#file_results
file_results = XDMFFile("u_two_material.xdmf")
#file_results = XDMFFile("u_two_material_inverting.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
``````

Its not very clear what you are asking.
You are not using `dx1` and `dx2` in the correct way.

Note that your example is far from minimal (as you could have stopped the whole example after generating `vm`, and visualize it in Paraview to see if it is correct.

Here is a complete example, that also includes the mesh conversion.

``````from fenics import *
import meshio
from pathlib import Path
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

mesh_file = Path("mecanismo_mesh.msh")
xdmf_file = mesh_file.with_suffix(".xdmf")
facet_file = xdmf_file.with_stem("mecanismo_facet")

mesh3D = create_mesh(msh, "tetra")
mesh2D = create_mesh(msh, "triangle")
meshio.write(xdmf_file, mesh3D)
meshio.write(facet_file, mesh2D)

F = 2e2

# Define 3D geometry
mesh = Mesh()
with XDMFFile(str(xdmf_file)) as infile:
mvc2 = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(str(facet_file)) as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
mvc3 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(str(xdmf_file)) as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc3)

# Define Material
Em = [210e3, 100]
#Em = [100, 210e3] # inverting

Vm = FunctionSpace(mesh, "DG", 0)
vm = Function(Vm)

for cell in cells(mesh):
i = cell.index()
if cf.array()[i] == 22:
vm.vector().vec().setValueLocal(i, Em)
elif cf.array()[i]==23:
vm.vector().vec().setValueLocal(i, Em)

# Define function space and base functions
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

# Boundary Condition
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), mf, 20)

BCs = [bc]

ds = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=21)
t = Expression(('q_x', 'q_y', 'q_z'), degree = 1, q_x = 0.0, q_y = 0.0, q_z = -F)

## Constitutive relation

def eps(v):

#E = Constant(210e3)
E = vm
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
dim = v.geometric_dimension()
return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)

## Variational formulation

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(t, u_)*ds

A = assemble(a)
B = assemble(L)

[bc.apply(A, B) for bc in BCs]

u = Function(V, name="Displacement")
solve(A ,u.vector() ,B)

#file_results
file_results = XDMFFile("u_two_material.xdmf")
#file_results = XDMFFile("u_two_material_inverting.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.close()
``````

Choosing the inverse values gives a significantly different solution.

For any further inquiries, please make an effort in isolating what you have problems with (in this case, the definition of `E` through `vm`).

Dear @Dokken,

Thanks a lot for your help, it works now.

What about considering this model as a mechanism? Create a form of contact between the volumes, effectively “separating” the meshes, or even a way to reduce the degrees of freedom between the nodes that are shared between the two volumes, releasing movement in relation to some specific axis, something of the like, like some structural software does. Is it possible to do this with the Fenics?

@bleyerj or someone who works with mechanism could also help?

Dear @Dokken,

As I said before the code is working. However, I went to test in parallel mode and there was the following error:

``````Traceback (most recent call last):
File "OTE_Two_cantilever_static_3D_two_material.py", line 181, in <module>
rho_n = interpolate(Constant(float(Vol)), W) #Control of density
File "/opt/anaconda3/envs/fenics2/lib/python3.8/site-packages/fenics_adjoint/interpolation.py", line 11, in interpolate
output = backend.interpolate(*args, **kwargs)
File "/opt/anaconda3/envs/fenics2/lib/python3.8/site-packages/dolfin/fem/interpolation.py", line 71, in interpolate
Pv.interpolate(v._cpp_object)
File "/opt/anaconda3/envs/fenics2/lib/python3.8/site-packages/dolfin/function/function.py", line 365, in interpolate
self._cpp_object.interpolate(u)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to interpolate function into function space.
*** Reason:  Wrong size of vector.
*** Where:   This error was encountered inside FunctionSpace.cpp.
*** Process: 5
***
*** DOLFIN version: 2019.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpiexec detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

Process name: [[52824,1],6]
Exit code:    1
--------------------------------------------------------------------------
``````

According to the error, could you tell me what happens in this case?

What code is throwing this error?
I cant help you by just seeing an error message and not a minimal code that reproduces the error.