# FEniCS: both normal and shear stress boundary conditions for elasticity? [closed]

I would like to have both the normal (xx) component and shear (xy) component of a 2D (stress) tensor defined on a boundary (y=const, for instance) for an equation which is of the type

$$\nabla \cdot \sigma = 0$$

The normal value of $\sigma_{xx}=P$ at $x=0,1$ can be accounted for via a Neumann boundary condition, in the term obtained by integration by parts in weak form:

$$\int_\Omega \sigma : \nabla \tau = \int_{\partial\Omega} P \cdot \tau \, d s .$$

$\tau$ is a vector test function for the trial function $\vec{u}$, which physically is the displacement. However, I am having difficulty fixing the shear stress at the boundaries $x=0,1$ of the unit square, as a simple example of my goal. How can one best accomplish this in FEniCS?

My question is similar to this one, but as far as I can tell, they are distinct.

A minimal running code in FEniCS (in python) which gives incorrect results for $\sigma$ is below. Thank you!

from dolfin import *

# Mesh Params ################
E, nu = 1.0, 0.4 #material params
mesh = UnitSquareMesh(100, 100)

print 'Creating Lagrangian function spaces...'
Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
Vf = FunctionSpace(mesh, "Lagrange", 1)
Vt = TensorFunctionSpace(mesh, "Lagrange", 1)

# P is stress vector of stress values normal and tangent to boundary
P = Expression(('(x[0]-xc)','(x[0]-xc)'),xc=0.5)

class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-7
return on_boundary and abs(x[0]) < tol

class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-7
return on_boundary and abs(x[0] - 1) < tol

left_boundary = LeftBoundary()
right_boundary = RightBoundary()

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
right_boundary.mark(boundary_parts, 1)
left_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]

##############################
# Define variational problem #
##############################
print 'Defining test and trial functions...'
v = TestFunction(Vf)
tau = TestFunction(Vv)
phi = TrialFunction(Vf)
u = TrialFunction(Vv)

# functions for variational problem
d = u.geometric_dimension()
# PLANE STRESS
sigma = E/(1+nu)* (epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)

# Define normal and tangent directions on mesh
norm = FacetNormal(mesh)
tang = as_vector([norm[1], -norm[0]])

tau_n = dot(tau,norm)
tau_t = dot(tau,tang)
P_n = dot(P,norm)
P_t = dot(P,tang)

#################
# SOLVE
#################
Lu = dot(P_n,tau_n)*ds(1) + dot(P_n,tau_n)*ds(2) + dot(P_t,tau_t)*ds(1) + dot(P_t,tau_t)*ds(2)
u = Function(Vv)
solve(au == Lu, u)

## View the results