# Weak form of the resultant force on a surface in FEniCS-X

**URL:** <https://fenicsproject.discourse.group/t/weak-form-of-the-resultant-force-on-a-surface-in-fenics-x/6242>\
**Category:** Uncategorized\
**Created:** [July 16, 2021, 10:02am UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-resultant-force-on-a-surface-in-fenics-x/6242 "2021-07-16T10:02:09Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Yingqi\_Jia](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/yingqi_jia/32/2316_2.png) [@Yingqi\_Jia](https://fenicsproject.discourse.group/u/Yingqi_Jia)\
**Post date:** [July 16, 2021, 10:02am UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-resultant-force-on-a-surface-in-fenics-x/6242/1 "2021-07-16T10:02:09Z")

</div>

Hi all,

I have a small question about the weak form of the resultant force. For example, we have a brick element with the bottom surface fixed. We add the prescribed vertical displacement on the top surface.  
And this 3D elasticity problem can be easily solved as follows.

```auto
import numpy as np
import dolfinx
from dolfinx.cpp.mesh import CellType
from mpi4py import MPI
from petsc4py import PETSc
import ufl

mesh = dolfinx.BoxMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([1, 1, 1])], \
    [1, 1, 1], CellType.hexahedron)
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
fdim = mesh.topology.dim - 1
u = dolfinx.Function(V) # trial function (displacment field)
v = ufl.TestFunction(V) # test unction
mu = 1
lmbda = 1

def BottomSurface(x):
    return np.isclose(x[2], 0)
def TopSurface(x):
    return np.isclose(x[2], 1)

u1 = dolfinx.Function(V)
with u1.vector.localForm() as u1_loc:
    u1_loc.set(0.0)
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, BottomSurface)
bc1 = dolfinx.DirichletBC(u1, dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets))

disp = dolfinx.Function(V)
Vz = V.sub(2).collapse()
u2 = dolfinx.Function(Vz)
with u2.vector.localForm() as u2_loc:
    u2_loc.set(1.0)
boundary_dofs_2 = dolfinx.fem.locate_dofs_geometrical((V.sub(2), Vz), TopSurface)
bc2 = dolfinx.DirichletBC(u2, boundary_dofs_2, V.sub(2))

boundaries = [(1, lambda x: np.isclose(x[2], 0)), # bottom surface
              (2, lambda x: np.isclose(x[2], 1))] # top surface

# loop through all the boundary conditions to identify the facets
facet_indices, facet_markers = [], []
for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
facet_tag = dolfinx.MeshTags(mesh, fdim, np.array(np.hstack(facet_indices), dtype=np.int32), \
    np.array(np.hstack(facet_markers), dtype=np.int32))
metadata = {"quadrature_degree": 4}

P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * ufl.Identity(len(u))

ds = ufl.Measure('ds', subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", metadata=metadata)
R = ufl.inner(ufl.grad(v), P)*dx

problem = dolfinx.fem.NonlinearProblem(R, u, [bc1, bc2])
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
num_its, converged = solver.solve(u)

```

In addition, we can compute the _matrix form_ of the internal force as

```auto
F = dolfinx.fem.assemble_vector(R)
F.assemble()
F = F.getArray()
print(F.round(3))

```

And we obtain the force vector as

```auto
[-0.147 -0.147 -0.662 0.147 -0.147 -0.662 -0.147 0.147 -0.662 0.147
  0.147 -0.662 0. -0. 0.662 0. 0. 0.662 0. 0.
  0.662 0. -0. 0.662]

```

In this case, the resultant reaction force in the _z-direction_ on the _top surface_ should be 0.662\*4=2.648. So how can we construct the corresponding weak form such that `dolfinx.fem.assemble_scalar(W)=2.648`?

Thanks!

---

<div class="post-metadata">

**Author:** ![garth](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/garth/32/5_2.png) [@garth](https://fenicsproject.discourse.group/u/garth)\
**Post date:** [July 17, 2021, 12:43pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-resultant-force-on-a-surface-in-fenics-x/6242/2 "2021-07-17T12:43:28Z")

</div>

Use

```auto
np.sum(F[boundary_dofs_2[0]])

```

Using `dolfinx.fem.assemble_scalar` won’t support this ‘variational’ (conservative) computation of the reaction forces.

---

<div class="post-metadata">

**Author:** ![Yingqi\_Jia](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/yingqi_jia/32/2316_2.png) [@Yingqi\_Jia](https://fenicsproject.discourse.group/u/Yingqi_Jia)\
**Post date:** [July 18, 2021, 7:44am UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-resultant-force-on-a-surface-in-fenics-x/6242/3 "2021-07-18T07:44:14Z")

</div>

Hi @garth, the code `np.sum(F[boundary_dofs_2[0]])` does give us the **value** of the reaction force, but I still want to find the weak form. Because I want to take some derivatives with this form and it is very convenient to do so in FEniCSX.

I have an idea to construct this weak form by calculating the integral of P\_{33} over the top surface:

```auto
# direct calculation
print(np.sum(F[boundary_dofs_2[0]]))

# weak form
n_vec = dolfinx.Constant(mesh, (0,0,1))
L = ufl.inner(P, ufl.outer(n_vec, n_vec))*ds(2)
F = dolfinx.fem.assemble_scalar(L)
print(F)

```

But the result is somehow incorrect.

```auto
2.647058823529411
2.294117647058822

```

An interesting thing is that we can obtain the correct reaction force from the average of integrals over the top and bottom surfaces.

```auto
n_vec = dolfinx.Constant(mesh, (0,0,1))
L1 = ufl.inner(P, ufl.outer(n_vec, n_vec))*ds(1)
L2 = ufl.inner(P, ufl.outer(n_vec, n_vec))*ds(2)
F = dolfinx.fem.assemble_scalar(L1+L2)
print(F/2)

```

So I am not if we can follow this idea to establish the weak form of reaction force.

---

<div class="post-metadata">

**Author:** ![garth](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/garth/32/5_2.png) [@garth](https://fenicsproject.discourse.group/u/garth)\
**Post date:** [July 19, 2021, 7:53am UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-resultant-force-on-a-surface-in-fenics-x/6242/4 "2021-07-19T07:53:25Z")

</div>

You can simply integrate the tractions over surface:

```auto
t = ufl.inner(P * n, n1)*ds(1)
R = dolfinx.fem.assemble_scalar(t)

```

where `P` is the stress tensor, `n` is the unit outward normal vector and `n1` is the direction that you want to resolve the traction into (above snippet as not been tested). But, this computation is not conservative, i.e. the reactions and applied forces will not sum to exactly zero. This is a characteristic of Galerkin methods with continuous bases, and may or may not be an issue for your application.