Hi, I am a fresher to fenics and plan to perform a cantilever example, considering gravity. according to F=Ku, I get the nodal force vector but the values of nodal force seems to be wrong because the value of nodal force at the node without external force not equals to zero. So i want to know how to obtain the nodal stress value? additionally, what should i do to get reaction based on F=Ku?
Thanks in advance
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
L = 25.
H = 5.
Nx = 5
Ny = 1
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)
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)
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)
vertex_values = u.compute_vertex_values()
A = assemble(a)
F = A*vertex_values
np.savetxt('./F.txt',F, fmt='%.2e')
np.savetxt('./vertex_values.txt',vertex_values, fmt='%.2e')