I would like to be able to calculate and visualize the traction vector field for a given surface in a problem with hyperelasticity and Dirichlet boundary conditions. Below is sample code for an example problem in which a displacement is enforced on one side of a cube while the other is held in place.

```
import dolfin as dlf
import matplotlib.pyplot as plt
import numpy as np
def main():
# Optimization options for the form compiler
dlf.parameters['form_compiler']['cpp_optimize'] = True
dlf.parameters['form_compiler']['representation'] = 'uflacs'
# Define system
E = 10.0
nu = 0.3
mu = dlf.Constant(E/(2*(1 + nu)))
lmbda = dlf.Constant(E*nu/((1 + nu)*(1 - 2*nu)))
L = 0.1
W = 0.1
mesh = dlf.BoxMesh(dlf.Point(0, 0, 0), dlf.Point(L, W, W), 2, 2, 2)
V = dlf.VectorFunctionSpace(mesh, 'CG', 1)
# Define boundaries
facets = dlf.MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
facets.set_all(0)
left = dlf.CompiledSubDomain('near(x[0], side) && on_boundary', side=0.0)
left.mark(facets, 1)
right = dlf.CompiledSubDomain('near(x[0], side) && on_boundary', side=L)
right.mark(facets, 2)
# Create boundary conditions
cl = dlf.Constant((0.0, 0.0, 0.0))
cr = dlf.Expression(('gamma', '0.0', '0.0'), gamma=0.1, degree=1)
bcl = dlf.DirichletBC(V, cl, left)
bcr = dlf.DirichletBC(V, cr, right)
bcs = [bcl, bcr]
# Solve for displacement vector field
u = dlf.Function(V)
du = dlf.TrialFunction(V)
v = dlf.TestFunction(V)
dim = len(u)
F = dlf.Identity(dim) + dlf.grad(u)
body_force = dlf.Constant((0.0, 0.0, 0.0))
traction_force = dlf.Constant((0.0, 0.0, 0.0))
E = green_lagrange_strain(u)
energy_density = lmbda/2*dlf.tr(E)**2 + mu*dlf.tr(dlf.dot(E, E))
pot = (energy_density*dlf.dx - dlf.dot(body_force, u)*dlf.dx -
dlf.dot(traction_force, u)*dlf.ds)
F = dlf.derivative(pot, u, v)
J = dlf.derivative(F, u, du)
dlf.solve(F == 0, u, bcs, J=J)
# Save solution in VTK format
out_file = dlf.File('outs/saint_venant-kirchhoff_displacement.pvd')
out_file << u
# Calculate second Piola-Kirchhoff pseudo traction vector field
N = dlf.FacetNormal(mesh)
T = dlf.dot(second_piola_kirchhoff_stress(u, lmbda, mu), N)
def green_lagrange_strain(u):
linear_part = dlf.grad(u) + dlf.grad(u).T
nonlinear_part = dlf.dot(dlf.grad(u), dlf.grad(u).T)
return 1/2*(linear_part + nonlinear_part)
def second_piola_kirchhoff_stress(u, lmbda, mu):
E = green_lagrange_strain(u)
return 2*mu*E + lmbda*dlf.tr(E)*dlf.Identity(len(u))
if __name__ == '__main__':
main()
```

I seem to be able to integrate over the surface to calculate the mean normal force vector,

```
Tnn_c = dlf.dot(T, N)
Tnn = Tnn_c*N
ds = dlf.ds(domain=mesh, subdomain_data=facets)
left_normal_force = dlf.assemble(Tnn_c*ds(1))
right_normal_force = dlf.assemble(Tnn_c*ds(2))
```

However, attempts to output the traction vector so that it can be visualized with Paraview fail. If I naively try to project the traction vector field,

```
W = dlf.VectorFunctionSpace(mesh, 'DG', 1)
dlf.project(T, W)
```

I end up with the error,

```
ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal.
```

If I instead try to manually project and integrate over surfaces only,

```
W = dlf.VectorFunctionSpace(mesh, 'DG', 1)
w = dlf.TrialFunction(W)
v = dlf.TestFunction(W)
a = dlf.dot(w, v)*dlf.ds
L = dlf.dot(T, v)*dlf.ds
w = dlf.Function(W)
dlf.solve(a == L, w)
out_file = dlf.File('outs/saint_venant-kirchhoff_traction.pvd')
out_file << w
```

I get the following error

```
*** Error: Unable to compute matrix factorisation.
*** Reason: The provided data did not satisfy the prerequisites.
*** Where: This error was encountered inside EigenLUSolver.cpp.
```

I wondered whether I should be trying to use the boundary mesh, so I ran

```
bmesh = dlf.cpp.mesh.BoundaryMesh(mesh, 'exterior')
N = dlf.FacetNormal(bmesh)
T = dlf.dot(second_piola_kirchhoff_stress(u, lmbda, mu), N)
```

but ended with the error

```
TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'
```

The closest posts on the forums that I could find are here, which involves a different variational form but seems to be specific to a Neumann boundary condition, and here, which has much more going on and many things that are not clear to me and don’t seem to be documented anywhere.

I am using dolfin version 2019.1.0, which I installed directly (i.e. I am not running it with docker).

Any input would be greatly appreciated.