Hello everyone,
I am trying to calculate the reaction force through integration of stresses over a subdomain as done on the following link:
https://comet-fenics.readthedocs.io/en/latest/demo/tips_and_tricks/computing_reactions.html
I am attaching a minimal working code for 1D bar in tension. Left end of this bar is fixed and a displacement is applied on the right end. The program is providing the zero reaction force. The same method worked for me in dolfin but in dolfinx it is giving zero reaction force.
from dolfinx import fem, mesh # dolfinx library
from mpi4py import MPI
from ufl import replace, TestFunction, TrialFunction, nabla_div, Measure, grad, inner, dot, derivative, lhs, rhs
import dolfinx
import ufl
import numpy as np
from math import pi
from dolfinx.mesh import locate_entities, meshtags
from dolfinx.fem import FunctionSpace, Function, locate_dofs_geometrical, dirichletbc, Constant
from dolfinx.fem.petsc import LinearProblem
length = 200. # length of specimen
num_elem = 400 # number of elements in specimen
mesh = mesh.create_interval(MPI.COMM_WORLD,num_elem,[0., length])
def bound_cond(x):
tol = 1E-14
return x[0] <= 0. + tol
remain = locate_entities(mesh, mesh.topology.dim, bound_cond)
boundaries = meshtags(mesh,mesh.topology.dim-1,np.sort(remain),1)
ds_bound = Measure("ds", domain=mesh, subdomain_data=boundaries)
V_u = FunctionSpace(mesh, ('Lagrange', 1))
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)
u_R = Function(V_u)
disp = Function(V_u)
dx = Measure("dx",domain=mesh)
ds = Measure("ds",domain=mesh)
# Identify the boundaries to apply boundary conditions
def on_boundary_L(x): # Left edge to fix
return np.isclose(x[0],0.0)
dof_left_u = locate_dofs_geometrical(V_u, on_boundary_L)
bcl = dirichletbc(Constant(mesh,0.0), dof_left_u, V_u)
# Define Material Properties
E = Constant(mesh,30000.0) # Young's Modulus
nu = Constant(mesh,0.2) # Poisson ratio
def eps(u):
"""Strain tensor as a function of the displacement"""
return grad(u)
def sigma(u):
"""Stress tensor of the undamaged material as a function of the displacement"""
mu = E / (2.0*(1.0 + nu))
lmbda = (E * nu) / ((1.0 - 2.0*nu)*(1.0 + nu))
return E*(eps(u))
f_ext = Constant(mesh,0.0)
external_work = dot(f_ext, u) * dx
elastic_energy = 0.5*inner(sigma(u), eps(u))*dx
total_energy = elastic_energy - external_work
E_u = derivative(total_energy,u,v)
disp_a = Constant(mesh,0.00024)
def on_boundary(x):
return np.isclose(x[0],length)
dof_right_u = locate_dofs_geometrical(V_u, on_boundary)
bcr = dirichletbc(disp_a, dof_right_u, V_u)
bc_disp = [bcl, bcr]
E_du = replace(E_u,{u:du})
problem_u = LinearProblem(a = lhs(E_du), L = rhs(E_du), bcs = bc_disp,u = u, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "umfpack"})
problem_u.solve()
print(fem.assemble_scalar(dolfinx.fem.form(sigma(u)[0]*ds_bound(1)))) # reaction force
Please help me to find the solution of this issue.
Thank you
Manish