Hello.

I am doing a simulation of elastic bodies and ran into some strange results.

This is my MRE:

I have an elastic bar (10x10x100) and am compressing it with a set force on the long axis.

Now I wrote code to apply my force not on the whole area (10x10) on the top but only a square in the middle of it with a side length of l.

Now when I set l < 2.82 I get a maximum displacement of 0.0.

I think that is an error. If I apply an absurd force of 1e90N to a 2x2 square on a 10x10 square, there should be some visible displacement but fenics says theres nothing.

Is that a mistake in my code and can I fix that somehow?

```
import fenics as fx
import mshr
#define the mesh in m
domain = mshr.Box(fx.Point(0, 0, 0), fx.Point(10, 10, 100))
mesh = mshr.generate_mesh(domain, 64)
# unscaled variables for an elastic body
E = 300
G = 103
mu = G
lambda_ = G*(E-2*G) / (3*G - E)
# this follows linearElasticity.ipynb of the fenics tutorial closely
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 = 2.0 #side length of the applied square
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 = 1e90/(l*l) #absurd 1e90 force
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.solve(a == L, u, bc, solver_parameters={'linear_solver':'mumps'})
V = fx.FunctionSpace(mesh, 'P', 1)
u_magnitude = fx.sqrt(fx.dot(u, u))
u_magnitude = fx.project(u_magnitude, V)
print('min/max u:',
u_magnitude.vector().min(),
u_magnitude.vector().max())
```

This prints me 0.0, 0.0.

Setting l=10 the maximum displacement prints the theoretically expected 0.003m (when appliying 1N Force) so the code works there.

Another problem is visible when setting l=5 (black square) and applying a force of 1e3. The bar is bend not only straight down but also to the side and the side seems to be random when retrying with different values for l.

Thanks in advance!