Hello all,

I have a very simple structure with three cylinders. This is the geo file:

```
//+
SetFactory("OpenCASCADE");
//+
Cylinder(1) = {0.0, 0, 0, 0, 0, 30, 0.5, 2*Pi};
//+
Cylinder(2) = {0.0, 1.0, 0, 0, 0, 30, 0.5, 2*Pi};
//+
Cylinder(3) = {0.8, 0, 0, 0, 0, 30, 0.3, 2*Pi};
//+
Physical Volume(10) = {2};
//+
Physical Volume(11) = {1};
//+
Physical Volume(12) = {3};
//+
Physical Surface(13) = {1};
//+
Physical Surface(14) = {2};
//+
Physical Surface(15) = {3};
//+
Physical Surface(16) = {4};
//+
Physical Surface(17) = {5};
//+
Physical Surface(18) = {6};
//+
Physical Surface(19) = {7};
//+
Physical Surface(20) = {8};
//+
Physical Surface(21) = {9};
```

I want to solve a simple elastic problem by fixing the top base of all cylinders (14, 17, 20) and apply a load along x direction in the bottom base of the smallest cylinder (21).

I use dolphin-converter to get the xml files.

This is the dolfin code:

```
from dolfin import *
from ufl import nabla_div
import numpy as np
import time
# Define the material properties of steel
E = 210e3 # Young's modulus
nu = 0.28 # Poisson's ratio
lame_mu = E/(2*(1+nu)) # Lame parameters
lame_lambda = (E*nu)/((1+nu)*(1-2*nu))
#############################################################
# Load mesh
mesh = Mesh('three_circle.geo.xml')
Volume = MeshFunction('size_t', mesh, "three_circle.geo_physical_region.xml")
bnd_mesh = MeshFunction('size_t', mesh, "three_circle.geo_facet_region.xml")
xdmf = XDMFFile(mesh.mpi_comm(),"Mesh.xdmf")
xdmf.write(mesh)
xdmf.write(Volume)
xdmf = XDMFFile(mesh.mpi_comm(),"boundaries.xdmf")
xdmf.write(bnd_mesh)
xdmf.close()
mesh= Mesh()
xdmf = XDMFFile(mesh.mpi_comm(), "Mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("Mesh.xdmf") as infile:
infile.read(mvc, "f")
Volume = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc2 = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("boundaries.xdmf") as infile:
infile.read(mvc2, "f")
bnd_mesh = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = dx(metadata={'quadrature_degree': 2},
subdomain_data=Volume)
ds = ds(metadata={'quadrature_degree': 2},
subdomain_data=bnd_mesh)
##########################################################################
# Define function space
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
Area = 0.269959
# 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 lame_lambda*nabla_div(u)*Identity(d) + 2*lame_mu*epsilon(u)
# Dirichelet BCs
bc1 = DirichletBC(V, Constant((0, 0, 0)), bnd_mesh, 14)
bc2 = DirichletBC(V, Constant((0, 0, 0)), bnd_mesh, 17)
bc3 = DirichletBC(V, Constant((0, 0, 0)), bnd_mesh, 20)
bcs = [bc1, bc2, bc3]
# Boundary domain
ds = Measure('ds', subdomain_data=bnd_mesh)
load_force = 5
g_B= Constant((-load_force/Area, 0.0, 0.0))
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(v, g_B)*ds(21)
start_time = time.time()
# Compute solution
u = Function(V)
solve(a == L, u, bcs)
end_time = time.time()
elapsed_time = end_time - start_time
print("Elapsed time:", elapsed_time, "seconds")
# Compute von Mises stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'Lagrange', 1)
von_Mises = project(von_Mises, V)
u.rename('Displacemente field vector (mm)', '')
von_Mises.rename('von Mises stress (MPa)', '')
cylinder_steel_file = XDMFFile('3_cylinders.xdmf')
cylinder_steel_file.parameters["flush_output"] = True
cylinder_steel_file.parameters["functions_share_mesh"] = True
cylinder_steel_file.write(u, 0.0)
cylinder_steel_file.write(von_Mises, 0.0)
```

This code is working but I can see that the smallest cylinder is penetrating one of the other cylinders. What would be the best approach to avoid this? I have done my research and found that I could include a contact penalty. Is this feasible for problems with many physical volumes (ex: 10 cylinders touching each other).