PEERS, BDMS or AFW implementation

I’m trying to implement weakly symmetric mixed linear elasticity in Fenics 2019.1. I tried the PEERS (Arnold), but also BDMS (Stenberg) and AFW (Arnold, Falk, Winther) elements in 2d. None of them seem to work. I’m aware of the discussions https://fenicsproject.org/qa/10359/enriched-spaces-for-peers-help-on-implementation/ and PEERS element implementation . None of them present the full implementation with boundary conditions, I guess the problem is there.

Consider below the MWE for a horizontal bar, clamped on the left, subjected to non-homogeneous traction on the right, and homogeneous traction on bottom and top.

Thanks,

import numpy as np
import dolfin as df
import matplotlib.pyplot as plt

# boundary condition based on the Mixed Poisson example
class BoundarySource(df.UserExpression):
    def __init__(self, mesh, g, **kwargs):
        self.mesh = mesh
        self.g = g
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        cell = df.Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        values[0] = self.g*n[0]
        values[1] = self.g*n[1]
    def value_shape(self):
        return (2,)

tx, ty =  0.0, 0.1 # Traction applied on the right end 
lamb, mu = 70., 45. # material parameters
Nx, Ny, Lx, Ly =  20, 5, 2.0, 0.5 # mesh parameters

# Mesh and mesh annotations
mesh = df.RectangleMesh(df.Point(0.0,0.0) , df.Point(Lx,Ly), Nx, Ny, 'left/right');
leftBnd = df.CompiledSubDomain('near(x[0], 0.0) && on_boundary')
rightBnd = df.CompiledSubDomain('near(x[0], Lx) && on_boundary', Lx=Lx)
topBnd = df.CompiledSubDomain('near(x[1], Ly) && on_boundary', Ly=Ly)
bottomBnd = df.CompiledSubDomain('near(x[1], 0.0) && on_boundary')
listBnds = [(0, bottomBnd), (1, rightBnd), (2, topBnd), (3, leftBnd)]
boundary_markers = df.MeshFunction("size_t", mesh, dim=1, value=0)
[bnd[1].mark(boundary_markers, bnd[0]) for bnd in listBnds]
dx = df.Measure('dx', domain=mesh)
ds = df.Measure('ds', domain=mesh, subdomain_data=boundary_markers)
normal = df.FacetNormal(mesh)
 
# For PEERS (see comments at the end for other spaces)
Ve = df.VectorElement("DG", mesh.ufl_cell(), 0)
Skew_e = df.FiniteElement("CG", mesh.ufl_cell(), 1)
Hdiv_e = df.FiniteElement("RT", mesh.ufl_cell(), 1) # naturally vectorial
Be = df.VectorElement("Bubble", mesh.ufl_cell(), 3) 
We = df.MixedElement(Hdiv_e, Hdiv_e, Ve, Skew_e,  Be)
Wh = df.FunctionSpace(mesh, We)

# sigma reconstruction : row-wise
sigma = lambda s1, s2, s0: ( df.as_tensor([[s1[0] , s1[1]], [s2[0] , s2[1]]]) + 
                            df.as_tensor([[s0[0].dx(1),-s0[0].dx(0)],[s0[1].dx(1),-s0[1].dx(0)]]) )

# asymetric measure : (asym(s), vec(W)) = (s, W), W is skew, vec give the axial vector (or just a scalar in 2d)
asym = lambda s: s[0,1] - s[1,0]


# Boundary Conditions
zero_traction = BoundarySource(mesh, df.Constant(0.), degree = 1)
traction_x = BoundarySource(mesh, df.Constant(tx), degree = 1)
traction_y = BoundarySource(mesh, df.Constant(ty), degree = 1)

bcs = []
bcs.append(df.DirichletBC(Wh.sub(0), zero_traction, boundary_markers, 0)) # stress row 0 , bottom
bcs.append(df.DirichletBC(Wh.sub(1), zero_traction, boundary_markers, 0)) # stress row 1 , bottom
bcs.append(df.DirichletBC(Wh.sub(0), zero_traction, boundary_markers, 2)) # stress row 0 , top
bcs.append(df.DirichletBC(Wh.sub(1), zero_traction, boundary_markers, 2)) # stress row 1 , top
bcs.append(df.DirichletBC(Wh.sub(0), traction_x, boundary_markers, 1)) # stress row 0, right
bcs.append(df.DirichletBC(Wh.sub(1), traction_y, boundary_markers, 1)) # stress row 1, right

# Variational formulation
w, dw = df.TrialFunction(Wh), df.TestFunction(Wh) # trial

(sig1, sig2, u, gamma, bsig) = df.split(w)
(tau1, tau2, v, theta, btau) = df.split(dw)

sig = sigma(sig1, sig2, bsig)
tau = sigma(tau1, tau2, btau)

# Note that A*sig = alpha*sig + beta*tr(sig)*Id 
alpha = 0.5/mu
beta = -lamb/(2*mu + 2*lamb)
a = lambda s, t : alpha*df.inner(s,t)*dx + beta*df.tr(s)*df.tr(t)*dx
b = lambda s, v, q: df.inner(df.div(s), v)*dx + df.inner(asym(s), q)*dx

A = a(sig,tau) + b(sig,v,theta) + b(tau,u,gamma)
F =  -df.inner(df.Constant((0.,0.)) , v)*dx  # no body forces
F+= df.inner(df.Constant((0.,0.)) , df.dot(tau, normal))*ds(3) # zero displacement on the left 

# # Solve
wh = df.Function(Wh)
df.solve(A == F, wh, bcs)

(sig1, sig2, u, gamma, bsig) = wh.split(deepcopy = True)

print(u.vector().get_local()[:].min())
print(u.vector().get_local()[:].max())

plt.figure(1)
df.plot(u, mode = 'displacement')
df.plot(mesh)


# For BDMS (Stenberg family)
# Ve = df.VectorElement("DG", mesh.ufl_cell(), 0)
# Skew_e = df.FiniteElement("DG", mesh.ufl_cell(), 1)
# Hdiv_e = df.FiniteElement("BDM", mesh.ufl_cell(), 1) # naturally vectorial
# Be = df.VectorElement("Bubble", mesh.ufl_cell(), 3)
# We = df.MixedElement(Hdiv_e, Hdiv_e, Ve, Skew_e,  Be)

# For AWF family (the bubble should be eliminated thereafter)
# k = 2
# Ve = df.VectorElement("DG", mesh.ufl_cell(), k-1)
# Skew_e = df.FiniteElement("DG", mesh.ufl_cell(), k-1)
# Hdiv_e = df.FiniteElement("BDM", mesh.ufl_cell(), k) # naturally vectorial
# We = df.MixedElement(Hdiv_e, Hdiv_e, Ve, Skew_e)