Add HCurlDiv space and covariant-contravariant Piola mapping

ufl/__init__.py

@@ -74,6 +74,7 @@
     -HCurl
     -HEin
     -HDivDiv
+    -HCurlDiv
 
 
 * Pull backs::
@@ -84,6 +85,7 @@
     -l2_piola
     -double_contravariant_piola
     -double_covariant_piola
+    -covariant_contravariant_piola
 
 * Function spaces::
 
@@ -418,7 +420,7 @@
     identity_pullback,
     l2_piola,
 )
-from ufl.sobolevspace import H1, H2, L2, HCurl, HDiv, HDivDiv, HEin, HInf
+from ufl.sobolevspace import H1, H2, L2, HCurl, HCurlDiv, HDiv, HDivDiv, HEin, HInf
 from ufl.split_functions import split
 from ufl.tensors import (
     as_matrix,
@@ -448,12 +450,14 @@
     "HInf",
     "HEin",
     "HDivDiv",
+    "HCurlDiv",
     "identity_pullback",
     "l2_piola",
     "contravariant_piola",
     "covariant_piola",
     "double_contravariant_piola",
     "double_covariant_piola",
+    "covariant_contravariant_piola",
     "l2_piola",
     "MixedPullback",
     "SymmetricPullback",

ufl/pullback.py

@@ -31,6 +31,7 @@
     "L2Piola",
     "DoubleContravariantPiola",
     "DoubleCovariantPiola",
+    "CovariantContravariantPiola",
     "MixedPullback",
     "SymmetricPullback",
     "PhysicalPullback",
@@ -328,6 +329,51 @@ def physical_value_shape(self, element, domain) -> typing.Tuple[int, ...]:
         return (gdim, gdim)
 
 
+class CovariantContravariantPiola(AbstractPullback):
+    """The covariant contravariant Piola pull back."""
+
+    def __repr__(self) -> str:
+        """Return a representation of the object."""
+        return "CovariantContravariantPiola()"
+
+    @property
+    def is_identity(self) -> bool:
+        """Is this pull back the identity (or the identity applied to mutliple components)."""
+        return False
+
+    def apply(self, expr):
+        """Apply the pull back.
+
+        Args:
+            expr: A function on a physical cell
+
+        Returns: The function pulled back to the reference cell
+        """
+        from ufl.classes import Jacobian, JacobianDeterminant, JacobianInverse
+
+        domain = extract_unique_domain(expr)
+        J = Jacobian(domain)
+        detJ = JacobianDeterminant(J)
+        K = JacobianInverse(domain)
+        # Apply transform "row-wise" to TensorElement(PiolaMapped, ...)
+        *k, i, j, m, n = indices(len(expr.ufl_shape) + 2)
+        kmn = (*k, m, n)
+        return as_tensor((1.0 / detJ) * K[m, i] * expr[kmn] * J[j, n], (*k, i, j))
+
+    def physical_value_shape(self, element, domain) -> typing.Tuple[int, ...]:
+        """Get the physical value shape when this pull back is applied to an element.
+
+        Args:
+            element: The element that the pull back is applied to
+            domain: The domain
+
+        Returns:
+            The value shape when the pull back is applied to the given element
+        """
+        gdim = domain.geometric_dimension()
+        return (gdim, gdim)
+
+
 class MixedPullback(AbstractPullback):
     """Pull back for a mixed element."""
 
@@ -589,6 +635,7 @@ def physical_value_shape(self, element, domain) -> typing.Tuple[int, ...]:
 l2_piola = L2Piola()
 double_covariant_piola = DoubleCovariantPiola()
 double_contravariant_piola = DoubleContravariantPiola()
+covariant_contravariant_piola = CovariantContravariantPiola()
 physical_pullback = PhysicalPullback()
 custom_pullback = CustomPullback()
 undefined_pullback = UndefinedPullback()

ufl/sobolevspace.py

@@ -53,6 +53,7 @@ def __init__(self, name, parents=None):
             "HCurl": 0,
             "HEin": 0,
             "HDivDiv": 0,
+            "HCurlDiv": 0,
             "DirectionalH": 0,
         }[self.name]
 
@@ -156,7 +157,7 @@ def __lt__(self, other):
 
         if other in [HDiv, HCurl]:
             return all(self._orders[i] >= 1 for i in self._spatial_indices)
-        elif other.name in ["HDivDiv", "HEin"]:
+        elif other.name in ["HDivDiv", "HEin", "HCurlDiv"]:
             # Don't know how these spaces compare
             return NotImplementedError(f"Don't know how to compare with {other.name}")
         else:
@@ -175,3 +176,4 @@ def __str__(self):
 HInf = SobolevSpace("HInf", [H2, H1, HDiv, HCurl, L2])
 HEin = SobolevSpace("HEin", [L2])
 HDivDiv = SobolevSpace("HDivDiv", [L2])
+HCurlDiv = SobolevSpace("HCurlDiv", [L2])