Issue of computing boundary force from sigma field

dolfinx
1

Based on the tutorial case for a linear elasticity problem: Implementation — FEniCSx tutorial

the demo is modified a little, that is, ignore the body force f, and add a x-direction force on the outside face (x[0] == 1).

The computation is done well, however, when I tried to compute the boundary force back from the displacement field for validation, the force i got is much smaller, why? What’s wrong with my code? Thanks in advance for any help!

the code is listed below:

  • this part is the same as the tutorial demo, until setting L
import dolfinx as dfx

# Scaled variable
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

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

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [20, 6, 6], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

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

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)

def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


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


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
  • this part is my custom code for solving
def traction_boundary(x):
    return np.isclose(x[0], 1)

fdim = 2
boundary_facets_load = dfx.mesh.locate_entities_boundary(domain, fdim, traction_boundary)


Tface = fem.Constant(domain, default_scalar_type((-0.013, 0, 0)))


facet_marker = dfx.mesh.meshtags(domain, fdim, boundary_facets_load, 1)
ds2 = ufl.Measure("ds", domain=domain, subdomain_data=facet_marker)

L = ufl.dot(Tface, v) * ds2

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
  • Then re-compute boundary force

V2 = fem.functionspace(domain, ("CG", 1, (domain.geometry.dim, domain.geometry.dim)))

sigma_expr = sigma(uh)

expr_sigma = fem.Expression(sigma_expr, V2.element.interpolation_points())

sigma_u = fem.Function(V2)

sigma_u.interpolate(expr_sigma)

n = ufl.FacetNormal(domain)

t = ufl.dot(sigma_u, n)

force_x_expr = t[0] * ds2(1)

total_force_x = fem.assemble_scalar(fem.form(force_x_expr))

total_force_x 
>>> -0.0006366250943262476

However, I set the Tface as (-0.013, 0, 0), but I got the total_force_x of -0.0006366. Anything wrong with my code for computing force from uh?

2

Please note that this is not the correct space for putting sigma in, as it will be a DG-1 quantity.
Secondly, note that you do not need the interpolation at all, i.e.

uh = problem.solve()

n = ufl.FacetNormal(domain)
t = ufl.dot(sigma(uh), n)

force_x_expr = t[0] * ds2(1)

total_force_x = fem.assemble_scalar(fem.form(force_x_expr))

print(total_force_x)
print(fem.assemble_scalar(fem.form(Tface[0] * ds2(1))))

which with refinement looks quite close.