diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py
index 320b3e32..0e6e741e 100644
--- a/tsfc/finatinterface.py
+++ b/tsfc/finatinterface.py
@@ -61,6 +61,7 @@
     "Nedelec 2nd kind H(curl)": finat.NedelecSecondKind,
     "Raviart-Thomas": finat.RaviartThomas,
     "Regge": finat.Regge,
+    "Gopalakrishnan-Lederer-Schoberl": finat.GopalakrishnanLedererSchoberl,
     "BDMCE": finat.BrezziDouglasMariniCubeEdge,
     "BDMCF": finat.BrezziDouglasMariniCubeFace,
     # These require special treatment below
diff --git a/tsfc/ufl_utils.py b/tsfc/ufl_utils.py
index 0cec4dbf..bc5069b6 100644
--- a/tsfc/ufl_utils.py
+++ b/tsfc/ufl_utils.py
@@ -403,6 +403,10 @@ def apply_mapping(expression, element, domain):
 
       G(X) = det(J)^2 K g(x) K^T  i.e. G_il(X)=(detJ)^2 K_ij g_jk K_lk
 
+    'covariant contravariant piola' mapping for g:
+
+      G(X) = det(J) J^T g(x) K^T     i.e. G_il(X) = det(J) J_ji g_jk(x) K_lk
+
     If 'contravariant piola' or 'covariant piola' (or their double
     variants) are applied to a matrix-valued function, the appropriate
     mappings are applied row-by-row.
@@ -442,6 +446,13 @@ def apply_mapping(expression, element, domain):
         *k, i, j, m, n = indices(len(expression.ufl_shape) + 2)
         kmn = (*k, m, n)
         rexpression = as_tensor(detJ**2 * K[i, m] * expression[kmn] * K[j, n], (*k, i, j))
+    elif mapping == "covariant contravariant piola":
+        J = Jacobian(mesh)
+        K = JacobianInverse(mesh)
+        detJ = JacobianDeterminant(mesh)
+        *k, i, j, m, n = indices(len(expression.ufl_shape) + 2)
+        kmn = (*k, m, n)
+        rexpression = as_tensor(detJ * J[m, i] * expression[kmn] * K[j, n], (*k, i, j))
     elif mapping == "symmetries":
         # This tells us how to get from the pieces of the reference
         # space expression to the physical space one.