Merge branch 'main' into pbrubeck/fix/call-form

test/test_apply_algebra_lowering.py

@@ -55,7 +55,7 @@ def test_determinant2(A2):
 
 def test_determinant3(A3):
     assert determinant_expr(A3) == (A3[0, 0]*(A3[1, 1]*A3[2, 2] - A3[1, 2]*A3[2, 1])
-                                    + A3[0, 1]*(A3[1, 2]*A3[2, 0] - A3[1, 0]*A3[2, 2])
+                                    + (A3[1, 0]*A3[2, 2] - A3[1, 2]*A3[2, 0])*(-A3[0, 1])
                                     + A3[0, 2]*(A3[1, 0]*A3[2, 1] - A3[1, 1]*A3[2, 0]))
 
 

test/test_derivative.py

@@ -483,6 +483,30 @@ def test_coefficient_derivatives(self):
     self.assertEqual(replace(actual, fd.function_replace_map), expected)
 
 
+def test_vector_coefficient_scalar_derivatives(self):
+    V = FiniteElement("Lagrange", triangle, 1)
+    VV = VectorElement("Lagrange", triangle, 1)
+
+    dv = TestFunction(V)
+
+    df = Coefficient(VV, count=0)
+    g = Coefficient(VV, count=1)
+    f = Coefficient(VV, count=2)
+    u = Coefficient(V, count=3)
+    cd = {f: df}
+
+    integrand = inner(f, g)
+
+    i0, i1, i2, i3, i4 = [Index(count=c) for c in range(5)]
+    expected = as_tensor(df[i1]*dv, (i1,))[i0]*g[i0]
+
+    F = integrand*dx
+    J = derivative(F, u, dv, cd)
+    fd = compute_form_data(J)
+    actual = fd.preprocessed_form.integrals()[0].integrand()
+    assert (actual*dx).signature() == (expected*dx).signature()
+
+
 def test_vector_coefficient_derivatives(self):
     V = VectorElement("Lagrange", triangle, 1)
     VV = TensorElement("Lagrange", triangle, 1)

test/test_tensoralgebra.py

@@ -105,6 +105,27 @@ def test_cross(self):
     self.assertEqualValues(C, D)
 
 
+def test_perp(self):
+    # Test perp is generally doing the correct thing
+    u = as_vector([3, 1])
+    v = perp(u)
+    w = as_vector([-1, 3])
+    self.assertEqualValues(v, w)
+
+    # Test that a perp does the correct thing to Zero
+    u = zero(2)
+    v = perp(u)
+    self.assertEqualValues(u, v)
+
+    # Test that perp throws an error if the wrong thing is provided
+    u = as_vector([3, 1, -1])  # 3D vector instead of 2D
+    with pytest.raises(ValueError):
+        v = perp(u)
+    u = as_matrix([[1, 3], [0, 4]])  # Matrix instead of vector
+    with pytest.raises(ValueError):
+        v = perp(u)
+
+
 def xtest_dev(self, A):
     C = dev(A)
     D = 0*C  # FIXME: Add expected value here

ufl/algorithms/apply_algebra_lowering.py

@@ -55,6 +55,9 @@ def c(i, j):
             return Product(a[i], b[j]) - Product(a[j], b[i])
         return as_vector((c(1, 2), c(2, 0), c(0, 1)))
 
+    def perp(self, o, a):
+        return as_vector([-a[1], a[0]])
+
     def dot(self, o, a, b):
         ai = indices(len(a.ufl_shape) - 1)
         bi = indices(len(b.ufl_shape) - 1)

ufl/algorithms/apply_derivatives.py

@@ -835,9 +835,9 @@ def coefficient(self, o):
             dosum = Zero(o.ufl_shape)
             for do, v in zip(dos, self._v):
                 so, oi = as_scalar(do)
-                rv = len(v.ufl_shape)
-                oi1 = oi[:-rv]
-                oi2 = oi[-rv:]
+                rv = len(oi) - len(v.ufl_shape)
+                oi1 = oi[:rv]
+                oi2 = oi[rv:]
                 prod = so * v[oi2]
                 if oi1:
                     dosum += as_tensor(prod, oi1)

ufl/compound_expressions.py

@@ -93,6 +93,8 @@ def determinant_expr(A):
             return determinant_expr_2x2(A)
         elif sh[0] == 3:
             return determinant_expr_3x3(A)
+        else:
+            return determinant_expr_nxn(A)
     else:
         return pseudo_determinant_expr(A)
 
@@ -116,15 +118,21 @@ def determinant_expr_3x3(A):
     return codeterminant_expr_nxn(A, [0, 1, 2], [0, 1, 2])
 
 
+def determinant_expr_nxn(A):
+    nrow, ncol = A.ufl_shape
+    assert nrow == ncol
+    return codeterminant_expr_nxn(A, list(range(nrow)), list(range(ncol)))
+
+
 def codeterminant_expr_nxn(A, rows, cols):
     if len(rows) == 2:
         return _det_2x2(A, rows[0], rows[1], cols[0], cols[1])
     codet = 0.0
     r = rows[0]
     subrows = rows[1:]
     for i, c in enumerate(cols):
-        subcols = cols[i + 1:] + cols[:i]
-        codet += A[r, c] * codeterminant_expr_nxn(A, subrows, subcols)
+        subcols = cols[:i] + cols[i + 1:]
+        codet += (-1)**i * A[r, c] * codeterminant_expr_nxn(A, subrows, subcols)
     return codet
 
 

ufl/finiteelement/mixedelement.py

@@ -311,8 +311,10 @@ def __init__(self, family, cell=None, degree=None, dim=None,
         # Initialize element data
         MixedElement.__init__(self, sub_elements, value_shape=value_shape,
                               reference_value_shape=reference_value_shape)
-        FiniteElementBase.__init__(self, sub_element.family(), cell, sub_element.degree(), quad_scheme,
-                                   value_shape, reference_value_shape)
+
+        FiniteElementBase.__init__(self, sub_element.family(), sub_element.cell(), sub_element.degree(),
+                                   sub_element.quadrature_scheme(), value_shape, reference_value_shape)
+
         self._sub_element = sub_element
 
         if variant is None:

ufl/finiteelement/restrictedelement.py

@@ -39,7 +39,10 @@ def __repr__(self):
         return f"RestrictedElement({repr(self._element)}, {repr(self._restriction_domain)})"
 
     def sobolev_space(self):
-        return L2
+        if self._restriction_domain == "interior":
+            return L2
+        else:
+            return self._element.sobolev_space()
 
     def is_cellwise_constant(self):
         """Return whether the basis functions of this element is spatially

ufl/operators.py

@@ -20,7 +20,7 @@
 from ufl.constantvalue import Zero, RealValue, ComplexValue, as_ufl
 from ufl.differentiation import VariableDerivative, Grad, Div, Curl, NablaGrad, NablaDiv
 from ufl.tensoralgebra import (
-    Transposed, Inner, Outer, Dot, Cross,
+    Transposed, Inner, Outer, Dot, Cross, Perp,
     Determinant, Inverse, Cofactor, Trace, Deviatoric, Skew, Sym)
 from ufl.coefficient import Coefficient
 from ufl.variable import Variable
@@ -175,7 +175,7 @@ def perp(v):
     v = as_ufl(v)
     if v.ufl_shape != (2,):
         raise ValueError("Expecting a 2D vector expression.")
-    return as_vector((-v[1], v[0]))
+    return Perp(v)
 
 
 def cross(a, b):

ufl/tensoralgebra.py

@@ -217,6 +217,34 @@ def __str__(self):
                             parstr(self.ufl_operands[1], self))
 
 
+@ufl_type(is_index_free=True, num_ops=1)
+class Perp(CompoundTensorOperator):
+    __slots__ = ()
+
+    def __new__(cls, A):
+        sh = A.ufl_shape
+
+        # Checks
+        if not len(sh) == 1:
+            raise ValueError(f"Perp requires arguments of rank 1, got {ufl_err_str(A)}")
+        if not sh[0] == 2:
+            raise ValueError(f"Perp can only work on 2D vectors, got {ufl_err_str(A)}")
+
+        # Simplification
+        if isinstance(A, Zero):
+            return Zero(sh, A.ufl_free_indices, A.ufl_index_dimensions)
+
+        return CompoundTensorOperator.__new__(cls)
+
+    def __init__(self, A):
+        CompoundTensorOperator.__init__(self, (A,))
+
+    ufl_shape = (2,)
+
+    def __str__(self):
+        return "perp(%s)" % self.ufl_operands[0]
+
+
 @ufl_type(num_ops=2)
 class Cross(CompoundTensorOperator):
     __slots__ = ("ufl_free_indices", "ufl_index_dimensions")