Symmetry of boundary terms in nonlinear mechanics [bug?]

mathematics
1

Hi,

I am working on nonlinear mechanics, where there is pressure acting on a boundary. This yields a term looking like

\int_\Gamma p\, J \mathbf{F}^{-T} n \cdot v\,ds

with the standard notation (J=\text{det}(\mathbf d), \mathbf F = \mathbf I + \nabla \mathbf d). This term is not symmetric (at least I concluded so after some lengthy computations), and I have also verified this by doing some computations in another FEM library. Thing is, UFL yields a symmetric jacobian no matter what, so I was wondering if this has something to do with how UFL computes things. Here is an MWE:

import numpy as np
from dolfin import *

# Create mesh and define function space
mesh = UnitCubeMesh(10, 10, 10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
left = CompiledSubDomain("near(x[0], side) && on_boundary", side=0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side=1.0)
c_bc = Constant(("0.0", "0.0", "0.0"))
r_bc = Expression(("scale*0.0",
               "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
               "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
              scale=0.5, y0=0.5, z0=0.5, theta=pi/3, degree=4)
bcl = DirichletBC(V, c_bc, left)
bcr = DirichletBC(V, r_bc, right)
bcs = [bcl, bcr]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v = TestFunction(V)             # Test function
u = Function(V)                 # Displacement from previous iteration

# Kinematics
I = Identity(3)    # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F
E = 0.5 * (C - I)
J = det(F)

# Elasticity parameters
E_elast, nu = 10.0, 0.3
mu, lmbda = Constant(E_elast/(2*(1 + nu))), Constant(E_elast*nu/((1 + nu)*(1 - 2*nu)))
P = mu * F + (-mu + lmbda * ln(J)) * inv(F.T)
Res = inner(P, grad(v)) * dx - dot(Constant(1e0) * J * inv(F.T) * FacetNormal(mesh), v) * ds
Jac = derivative(Res, u, du)
solve(Res == 0, u, bcs, J=Jac)

# Test symmetry
Jnpy = assemble(Jac).array()
ten = np.abs(Jnpy-Jnpy.T) < 1e1 * np.finfo(float).eps
print("First Piola Symmetric: ", np.all(ten))

Which of course in the end gives ‘True’. Any ideas on why this is happening?

Thanks in advance,
Nicolas

2

Hi,
I think that this is due to the fact that you perform the integration on the whole boundary. Hence your loading work is only \int_\Omega p n \cdot u ds = p \int_\Omega div(u) dx for a constant pressure which is conservative. Just change the pressure by some spatially varying expression, or integrate only on a subpart of the boundary and the term should be symmetric.

3

Thank you! Just added this and it is indeed not symmetric:

top = CompiledSubDomain("near(x[2], side) && on_boundary", side=1.0)
markers = MeshFunction("size_t", mesh, 1)
markers.set_all(0)
TOP = 7
top.mark(markers, TOP)
ds = Measure('ds', domain=mesh, subdomain_data=markers,
         metadata={'optimize': True})
Res = inner(P, grad(v)) * dx - dot(Constant(1e0) * J *
                               inv(F.T) * FacetNormal(mesh), v) * ds(TOP)

Marked as solved.