Hello,

I would like to simulate how an elastic sphere gets compressed when squished against two surfaces.

The goal is to get to know the force (and deformation work) needed to speeze the sphere into a certain shape

Here an illustration of what I mean:

In practice I would probably only push from the top and have a 0 deformation dirichlet condition at the bottom. That should yield the same result using the doubled force at the top.

I know how to compress an elastic box but dont get how I can do the same to a sphere as the loadSurface should change with force as the contact area between the top surface at the sphere grows.

```
import fenics as fx
import mshr
domain = mshr.Box(fx.Point(0, 0, 0), fx.Point(10, 10, 100))
mesh = mshr.generate_mesh(domain, 64)
# Elastic variables
E = 300 #young modulus
nu = 0.5 #poisson ratio
G = E / (2*(1+nu))+3 #+3 to not be just on the border of physicality
mu = G # mu and lambda_ are the factors in the naviers equation
lambda_ = G*(E-2*G) / (3*G - E)
V = fx.VectorFunctionSpace(mesh, 'Lagrange', 1)
tol = 1e-6
def boundary(x, on_boundary):
return on_boundary and (fx.near(x[2], 0, tol))
bc = fx.DirichletBC(V, fx.Constant((0, 0, 0)), boundary)
# Define strain and stress
def epsilon(u):
return 0.5*(fx.nabla_grad(u) + fx.nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
#return lambda_*fx.nabla_div(u)*fx.Identity(d) + 2*mu*epsilon(u)
return lambda_* fx.div(u) * fx.Identity(d) + 2*mu*epsilon(u)
# Define variational problem
u = fx.TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = fx.TestFunction(V)
# define where force applies
l = 10
class LoadSurface(fx.SubDomain):
def inside(self, x, on_boundary):
tol = 1e-6
return on_boundary and fx.near(x[2], 100, tol) and \
(x[0] >= 5 - l/2 - tol and x[0] <= 5 + l/2 + tol) and \
(x[1] >= 5 - l/2 - tol and x[1] <= 5 + l/2 + tol)
bnd_markers = fx.MeshFunction('size_t', mesh, d-1, 0)
loadSurface = LoadSurface()
loadSurface.mark(bnd_markers, 1)
ds = fx.Measure('ds', domain=mesh, subdomain_data=bnd_markers)
force = 1e3/(l*l)
T = fx.Constant((0, 0, -force))
a = fx.inner(sigma(u), epsilon(v))*fx.dx
L = fx.dot(T, v)*ds(1)
# Compute solution
u = fx.Function(V)
fx.parameters['ghost_mode'] = 'shared_facets'
fx.solve(a == L, u, bc, solver_parameters={'linear_solver':'mumps'})
```

Is there a way to simulate that with a changing load surface or are the better methods?