Hi everyone, 

I was trying to use a mixed-formulation for an incompressible hyperelasticity problem. See the following MWE which uses a Taylor-Hood approximation -- P2 for displacement (u) and P1 for pressure (p): 

mesh = UnitCubeMesh(5,5,5)
Vh = VectorFunctionSpace(mesh, "CG", 2)
Qh = FunctionSpace(mesh, "CG", 1)
Wh = Vho * Qh
w = Function(Wh)
dw = TrialFunction(Wh)
v = TestFunction(Wh)
u, p = split(w)
mu, lmbda = Constant(1.), Constant(1.e3)
Js = (lmbda + p + sqrt((lmbda + p)**2. + 4.*lmbda*mu ))/(2.*lmbda)
F = Identity(mesh.geometric_dimension()) + grad(u)
J = det(F)
I1 = inner(F, F)
Pi = p * (J - Js) + mu/2 * (I1 - 3) - mu * ln(Js) + lmbda/2 * (Js - 1)**2
Pi *= dx
a = derivative(Pi, w, v)
Jac = derivative(a, w, dw)
bcs = [    
DirichletBC(Wh.sub(0).sub(0), Constant(0.), 1),
DirichletBC(Wh.sub(0).sub(1), Constant(0.), 3),    
DirichletBC(Wh.sub(0).sub(2), Constant(0.), 5),    
DirichletBC(Wh.sub(0).sub(0), Constant(1.), 2)
]
problem = NonlinearVariationalProblem(a, w, bcs, J=Jac)
solver = NonlinearVariationalSolver(problem)
solver.solve()
u, p = w.split()
File("displacement.pvd").write(u)


Everything works fine as expected. Now, I wanted to try another pair of approximations that is inf-sup stable, namely the one in the celebrated paper by Crouzeix-Raviart (1973). See (4.36) in http://archive.numdam.org/article/M2AN_1973__7_3_33_0.pdf. 

This amounts to approximating `u` by a "P2 + Bubble" and `p` by a "DG-1" space. This is also given in (3.2.3) in https://link.springer.com/content/pdf/10.1007%2F978-3-642-87047-7_5.pdf 
and (2.38) in https://link.springer.com/content/pdf/10.1007%2F978-3-642-61623-5_2.pdf

Now, in 3D as far as I can tell this amounds to approximating `u` with "P2 + an internal bubble + FacetBubble". Would anyone know how this could be achieved within firedrake? Because as far as I can see, firedrake supports `FacetBubble`. However, when I try something like: 

B3 = FiniteElement("Bubble", "tetrahedron", 3)
FB3 = FiniteElement("FacetBubble", "tetrahedron", 3)
P2 = FiniteElement("CG", "tetrahedron", 2)
Vh = FunctionSpace(mesh, P2 + B3 + FB3)

I get the following error

  File "/usr/lib/python3.8/functools.py", line 875, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 189, in convert_enrichedelement
    elements, deps = zip(*[_create_element(elem, **kwargs)
  File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 189, in <listcomp>
    elements, deps = zip(*[_create_element(elem, **kwargs)
  File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 310, in _create_element
    finat_element, deps = convert(ufl_element, **kwargs)
  File "/usr/lib/python3.8/functools.py", line 875, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/home/firedrake/firedrake/src/tsfc/tsfc/finatinterface.py", line 177, in convert_finiteelement
    return lmbda(cell, element.degree()), set()
  File "/home/firedrake/firedrake/src/FInAT/finat/fiat_elements.py", line 318, in __init__
    super(Bubble, self).__init__(FIAT.Bubble(cell, degree))
  File "/home/firedrake/firedrake/src/fiat/FIAT/bubble.py", line 33, in __init__
    super(Bubble, self).__init__(ref_el, degree, codim=0)
  File "/home/firedrake/firedrake/src/fiat/FIAT/bubble.py", line 24, in __init__
    raise RuntimeError('Bubble element of degree %d and codimension %d has no dofs' % (degree, codim))
RuntimeError: Bubble element of degree 3 and codimension 0 has no dofs

Any pointers would be appreciated
--
Bhavesh Shrimali
Ph.D. Student, Civil and Environmental Engineering,
University of Illinois Urbana-Champaign
Champaign, IL, United States

Email : bshrima2@illinois.edu