Issue in calculating reaction force by integration of stress over a subdomain

dolfinx
1

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

2

The dimension in remain does not match the dimension in your meshtag.

3

Thank you @dokken by reducing the dimension of remain it worked.
working code:

remain = locate_entities(mesh, mesh.topology.dim-1, bound_cond) 
boundaries = meshtags(mesh,mesh.topology.dim-1,np.sort(remain),1)

Thank you.

4

Hello. I tried to implement this example on a simple 3D-beam mesh. But I got following TypeError:
Can someone please tell me what I did wrong here?

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/nfs/home/zzhang/Dokumente/Zhang/DD_SOLVER/DD_SOLVER_I_BEAM/Elastoplastic/230710_plasticity_3D_beam_integratedForce.ipynb Cell 5 in 5
     52 problem = fem.petsc.LinearProblem(a, L, bcs=bcs, u = u, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
     53 u = problem.solve()
---> 54 print(fem.assemble_vector(dolfinx.fem.form(sigma(u)[0,0]*ds_bound)))

File ~/miniconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py:166, in form(form, dtype, form_compiler_params, jit_params)
    163         return list(map(lambda sub_form: _create_form(sub_form), form))
    164     return form
--> 166 return _create_form(form)

File ~/miniconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py:161, in form.._create_form(form)
    158 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    159 return form argument"""
    160 if isinstance(form, ufl.Form):
--> 161     return _form(form)
    162 elif isinstance(form, collections.abc.Iterable):
    163     return list(map(lambda sub_form: _create_form(sub_form), form))

File ~/miniconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py:155, in form.._form(form)
    149 # Subdomain markers (possibly None for some dimensions)
    150 subdomains = {_cpp.fem.IntegralType.cell: subdomains.get("cell"),
    151               _cpp.fem.IntegralType.exterior_facet: subdomains.get("exterior_facet"),
    152               _cpp.fem.IntegralType.interior_facet: subdomains.get("interior_facet"),
    153               _cpp.fem.IntegralType.vertex: subdomains.get("vertex")}
--> 155 return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)

File ~/miniconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py:51, in FormMetaClass.__init__(self, form, V, coeffs, constants, subdomains, mesh, ffi, code)
     49 self._code = code
     50 self._ufcx_form = form
---> 51 super().__init__(ffi.cast("uintptr_t", ffi.addressof(self._ufcx_form)),
     52                  V, coeffs, constants, subdomains, mesh)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.Form_float64(spaces: List[dolfinx::fem::FunctionSpace], integrals: Dict[dolfinx::fem::IntegralType, Tuple[List[Tuple[int, object]], dolfinx.cpp.mesh.MeshTags_int32]], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, mesh: dolfinx.cpp.mesh.Mesh = None)
    2. dolfinx.cpp.fem.Form_float64(form: int, spaces: List[dolfinx::fem::FunctionSpace], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], subdomains: Dict[dolfinx::fem::IntegralType, dolfinx.cpp.mesh.MeshTags_int32], mesh: dolfinx.cpp.mesh.Mesh)

Invoked with: , [], [], [], {: None, : array([  2,  12,  17,  35,  39,  44,  73,  77,  81,  86, 129, 133, 137,
       141, 146, 205, 209, 213, 217, 288, 292, 296, 376, 380, 466],
      dtype=int32), : None, : None},

Here is my code:

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.fem import Function, FunctionSpace, Constant
from ufl import TestFunction, dx, inner

L = 10.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[L,1,1]],[10,5,5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))
def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)

marked_facets = np.hstack([left_facets,right_facets])
marked_values = np.hstack([np.full_like(left_facets,1), np.full_like(right_facets,2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets],marked_values[sorted_facets])

from petsc4py.PETSc import ScalarType
# boundary condition on the left side: fixed
left_dofs = fem.locate_dofs_topological(V.sub(0),facet_tag.dim, facet_tag.find(1))
bc_l = fem.dirichletbc(ScalarType(0), left_dofs, V.sub(0))
right_dofs = fem.locate_dofs_topological(V.sub(0),facet_tag.dim, facet_tag.find(2))
bc_r = fem.dirichletbc(ScalarType(1.0),right_dofs,V.sub(0))
bcs =  [bc_l, bc_r]

from ufl import Measure
from dolfinx import fem, mesh          
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

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

def epsilon(u):
    return ufl.sym(ufl.grad(u)) 
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)


ds_bound = Measure("ds",domain=domain,subdomain_data=left_facets) 

u = Function(V)
du = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
u_R = Function(V)
disp = Function(V)
dx = ufl.Measure("dx",domain = domain)
ds = ufl.Measure("ds",domain = domain)

f = fem.Constant(domain, ScalarType((0, 0, 0)))
external_work = dot(f, u)* dx
elastic_energy = 0.5*inner(sigma(u), epsilon(u)) *dx
total_energy = elastic_energy - external_work

E_u = derivative(total_energy,u,v)
E_du = replace(E_u,{u:du})
a = lhs(E_du)
L = rhs(E_du)

# a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
# L = ufl.dot(f, v) * ufl.dx #+ ufl.dot(T, v) * ds

problem = fem.petsc.LinearProblem(a, L, bcs=bcs, u = u, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "umfpack"})
u = problem.solve()

print(fem.assemble_vector(dolfinx.fem.form(sigma(u)[0,0]*ds_bound)))

I would appreciate any help, tips, etc.

5

This is wrong.
It should be

ds_bound = Measure("ds",domain=domain,subdomain_data=facet_tag)

and

should be

print(fem.assemble_scalar(dolfinx.fem.form(sigma(u)[0,0]*ds_bound(1))))

or equivalently:

ds_bound = Measure("ds",domain=domain,subdomain_data=facet_tag, subdomain_id=1) 
print(fem.assemble_scalar(dolfinx.fem.form(sigma(u)[0,0]*ds_bound)))
1 Like
6

Hi Dokken,
Thank you very much for the quick answer. It works now.