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

1

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.

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

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

And we obtain the force vector as

[-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!

1 Like
2

Use

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

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

1 Like
3

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:

# 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.

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.

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.

4

You can simply integrate the tractions over surface:

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.

2 Likes