In the meantime, I figured out the solution to my problem.
What I suggested in the question was that differentially, the following equality must hold
$$
\int_{\partial \Omega^t} \delta \mathbf{u} \cdot (\boldsymbol{\sigma}\mathbf{n})\;\mathrm{d}A
= \int_{\partial \Omega^t} \delta \mathbf{u} \cdot \boldsymbol{\gamma}_\text{known} \;\mathrm{d}A
\tag{1}
$$
where $\boldsymbol{\gamma}_\text{known}$ is the traction on $\partial\Omega^t$, incorporated as a Neumann boundary condition.
However, this does not work very well when you are working with finite elements. The reason is that when you solve the problem
$$
\int_{\Omega} \delta \boldsymbol{\varepsilon}:\boldsymbol{\sigma}\;\mathrm{d}V
- \int_{\partial \Omega^{t}} \delta \mathbf{u} \cdot \boldsymbol{\gamma}_\text{known} \;\mathrm{d}A = 0
\tag{2}
$$
you obtain the solution vector $\mathbf{u}$ that satisfies (1).
The fact that the weak solution $\mathbf{u}$ satisfies (2) does not imply it will also satisfy (1), because the solution is a discretized one, and solving (1) would require
- the solution (2) for $\mathbf{u}$ by projecting it to the Gauss integration points,
- and re-projection of the stresses (evaluated at the nodes), back to the Gauss points.
The projections are where this error is introduced, which causes the slight difference from the expected result.
If we desire that our integrated value and given traction match, we need not look at any expression other than (2), just because (2) is the only equality that our discretized solution $\mathbf{u}$ satisfies. Since $\partial\Omega = \partial\Omega^t \cup \partial\Omega^u$, and (2) is satisfied for the given traction $\boldsymbol{\gamma}_\text{known}$, we must subtract the vectors assembled for (2) to obtain the unknown forces, i.e.
$$
F_\text{external, unknown} = F_\text{internal} - F_\text{external, known},.
$$
Then, we can use the degree of freedom map to sum the forces in the horizontal and vertical directions respectively. Here is a code that works for this problem
from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Sub domain for clamp at left end
def left(x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary
# Sub domain for rotation at right end
def right(x, on_boundary):
return x[0] > 1 - DOLFIN_EPS and on_boundary
# Material parameters
kappa = 1.0
mu = 1.0
Du = TrialFunction(V)
du = TestFunction(V)
u = Function(V)
# The force on the right boundary
f = Constant((1.0, 1.0))
I = Identity(len(u))
n = FacetNormal(mesh)
# Create mesh function over the cell facets
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(left).mark(boundary_subdomains, 1)
AutoSubDomain(right).mark(boundary_subdomains, 2)
dss = ds(subdomain_data=boundary_subdomains)
zero = Constant((0.0, 0.0))
bcs = DirichletBC(V, zero, left)
eps = sym(grad(u))
deps = sym(grad(du))
sigma = kappa * tr(eps) * I + 2 * mu * dev(eps)
f_int = inner(deps, sigma) * dx
f_ext = inner(f, du) * dss(2)
F = f_int - f_ext
J = derivative(F, u, Du)
solve(F == 0, u, bcs, J=J)
f_ext_known = assemble(f_ext)
f_ext_unknown = assemble(f_int) - f_ext_known
x_dofs = V.sub(0).dofmap().dofs()
y_dofs = V.sub(1).dofmap().dofs()
Fx = 0
for i in x_dofs:
Fx += f_ext_unknown[i]
Fy = 0
for i in y_dofs:
Fy += f_ext_unknown[i]
# The following must be equal to -1 * f
print("Horizontal reaction force: %f"%(Fx))
print("Vertical reaction force: %f"%(Fy))
Note that if you a displacement controlled problem, namely only Dirichlet boundary conditions and no Neumann boundary conditions, you should use a different approach since you do not have an external force term in the problem.
