Hi all!
I just implemented the exact “Orthotropic linear elasticity” example from Orthotropic linear elasticity — Numerical tours of continuum mechanics using FEniCS master documentation. After solving the equation, I want to check whether my traction boundary condition on the Top boundary.
With the following extension (where T0 and T1 are the components of the traction force T which is 0, 1e-03 per default in this example):
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigma(u), Vsig))
import numpy as np
for right_boundary_point in [Point(x[0], x[1]) for x in mesh.coordinates() if np.isclose(x[0], L)]:
sigma_tensor = np.array(sig(right_boundary_point)).reshape(2, 2)
print(f"sigma_tensor[:, 1]: {sigma_tensor[:, 1]}")
print(f"T: [{T0}, {T1}]")
I checked the values of the stress tensor on the boundary, multiplied by (0, 1), which is the outer normal. The results are not the exact values of T, but something like [-3.83081901e-07, 1.00894364e-03] for example.
Now in my application I have a lot larger values and T is around [0, 1e6]. With this value I get the following results (just examplewise for one point on the boundary): [-8.29929781e+02, 1.02096271e+06], which in the first component doesn’t match T at all.
Can someone tell me the reason? Is it just a numerical rounding error (maybe related to the mesh or something) or am I understanding the boundary condition wrong?
Best regards and thank you in advance!