Compute force on a surface and object weight in post-processing

variational formulation
1

Hi,

I have been following the tutorial of clamped 3D-beam bending by its own weight; the codes seem to work. Now I like to play with it, by computing the external clamping force on its left end. Also I like to compute the weight of the whole beam, since we know the rho (density) and g . The clamping force on the left end should be equal to the weight of the beam. My thought is to sum all inner(sigma, n) (sigma is stress, n is normal) at left surface, but how do I know sigma at a surface? For the weight, my thought is to add (rhogV) for all elements, but how do we know the volume of each element, and how to go through all elements?

I found a similar question in this following thread, but am not sure if I get its point. Could anyone point me the right direction? Thank you.
https://fenicsproject.org/qa/12464/integrating-forces-on-the-boundary/

My 3d-beam bending codes, almost identical to the tutorial.

from fenics import *
from ufl import nabla_div, nabla_grad
#import matplotlib.pyplot as plt

# scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# create mesh and define function space
#Boxmesh(x0,y0,z0, x,y,z,nx,ny,nz) create tetrahegonal mesh between two points with number of cells
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

# DirichletBC(FunctionSpace, value, boundary)
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# define strain and stress
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# define vairational problem
u = TrialFunction(V)
d = u.geometric_dimension()  #space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))   # gravity
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# compute solution
u = Function(V)
solve(a == L, u, bc)

# plot solution
#plot(u, title='Displacement', mode='displacement')  #Implemtn name error

# plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # dviatoric stress
von_Mises = sqrt(3./2*inner(s,s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
#plot(von_Mises, title='Stress Intensity')

# compute magnitue of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print('min/max u: ', u_magnitude.vector().get_local().min(), u_magnitude.vector().get_local().max())

# compute strain energy density
strain_energy_density = 0.5*inner(sigma(u), epsilon(u))
SED= project(strain_energy_density, V)

File('ft06_elasticity_u.pvd') << u
File('ft06_elasticity_vonMises.pvd') << von_Mises
File('ft06_elasticity_magnitude.pvd') << u_magnitude