How to extract stress/strain at integration point and force at nodes(like F=K*u)?

1

Hi, Is it possible to get strain/stress at all integration points and force at all nodes (like F=K*u). the type of element is triangular element with 1 degree.

from dolfin import *
import matplotlib.pyplot as plt
L = 25.
H = 5.
Nx = 50
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

plot(mesh)
plt.show()
dim  = mesh.topology().dim()

def eps(v):return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
   lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
rho_g = 1e-3
f = Constant((0, -rho_g))

V = VectorFunctionSpace(mesh, 'Lagrange', 1,dim)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):return near(x[0], 0.)
bc = DirichletBC(V, Constant((0.,0.)), left)

u = Function(V, name="Displacement")
solve(a == l, u, bc)

plot(1e3*u, mode="displacement")
print("Maximal deflection:", -u(L,H/2.)[1])
print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3))
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigma(u), Vsig))
print("Stress at (0,H):", sig(0, H))