Hello,

I want to solve a linear elastic problem and I have problems to calculate the resulting stresses.

As far as I understand it, first the variational problem is solved to obtain the displacements. Then these displacements are used to calculate the stresses afterwards.

However, the calculation of the displacement is based on the function for the stresses, so one can assume that if the displacements look correct, the stresses must also be correct. Unfortunately, this is not the case, so I suspect that the error is in the subsequent calculation of the stresses.

To calculate the stresses, ‘project’ is used. Here I have already tried to change the degree (1,2,3) as well as the type of functions in the function space (P, CG, DG), but unfortunately this does not provide the correct stresses.

I would be very happy about a helpful comment.

Best regards

Felix

This is my code:

```
from fenics import *
import matplotlib.pyplot as plt
# Geometry and mesh
mesh = RectangleMesh(Point(0, 0), Point(1.0, 1.0), 10, 10)
# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 1)
T = TensorFunctionSpace(mesh, 'P', 1) # Space for stresses
u = TrialFunction(V)
v = TestFunction(V)
# Define material parameters
mu = Constant(1.666666667)
lambda_ = Constant(0.576923077)
# Helper function for deformation gradient
def epsilon(u):
return 0.5 * (grad(u) + grad(u).T)
# Function to compute the stresses
def stress(u):
return lambda_ * div(u) * Identity(2) + 2 * mu * epsilon(u)
# Define boundary conditions as functions
def left_boundary_condition():
return DirichletBC(V, Constant((0, 0)), "near(x[0], 0)")
def right_boundary_condition(disp_value):
return DirichletBC(V, Constant((disp_value, 0)), "near(x[0], 1.0)")
# Define weak form of the problem
F = inner(stress(u), epsilon(v)) * dx
a, L = lhs(F), rhs(F)
u_solution = Function(V) # Define the solution function
# Apply boundary conditions
bcs = [left_boundary_condition(), right_boundary_condition(0.1)]
# Solve the problem
solve(a == L, u_solution, bcs)
s = stress(u_solution)
sigma = project(s, T)
# Visualize the stresses
c = plot(sigma[0, 0]) # XX-component of the stress
plt.colorbar(c, orientation='vertical', label='Stress XX')
plt.title("Distribution of Stress XX-component")
plot(mesh, color='k', linewidth=0.5)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
# Visualize the displacements
c = plot(u_solution, mode='displacement')
plt.title("Final Displacement Distribution")
plot(mesh, color='k', linewidth=0.5)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
```