Merge branch 'main' into dokken/update-3D-1D

.github/workflows/fenicsx-tests.yml

@@ -53,6 +53,8 @@ jobs:
           pip install .[ci]
       - name: Run FFCx unit tests
         run: python3 -m pytest -n auto ffcx/test
+      - name: Run FFCx demos
+        run: python3 -m pytest -n auto ffcx/demo/test_demos.py
 
   dolfinx-tests:
     name: Run DOLFINx tests

test/test_check_arities.py

@@ -9,7 +9,9 @@
     TestFunction,
     TrialFunction,
     adjoint,
+    as_tensor,
     cofac,
+    conditional,
     conj,
     derivative,
     ds,
@@ -84,3 +86,41 @@ def test_product_arity():
     with pytest.raises(ArityMismatch):
         L = inner(v, v) * dx
         compute_form_data(L, complex_mode=False)
+
+
+def test_zero_simplify_arity():
+    """
+    Test that adding verious zero-like expressions to a form is simplified,
+    such that one can compute form data for the integral.
+    """
+    cell = tetrahedron
+    D = Mesh(FiniteElement("Lagrange", cell, 1, (3,), identity_pullback, H1))
+    V = FunctionSpace(D, FiniteElement("Lagrange", cell, 2, (), identity_pullback, H1))
+    v = TestFunction(V)
+    u = Coefficient(V)
+
+    nonzero = 1
+    with pytest.raises(ArityMismatch):
+        F = inner(u, v + nonzero) * dx
+        compute_form_data(F)
+    z = Coefficient(V)
+
+    # Add a Zero-component (rank-0) of a tensor to a rank-1 tensor
+    zero = as_tensor([0, z])[0]
+    F = inner(u, v + zero) * dx
+    fd = compute_form_data(F)
+    assert fd.num_coefficients == 1
+
+    # Add a conditional that should have been simplified to zero (rank-0)
+    # to a rank-1 tensor
+    zero = conditional(z < 0, 0, 0)
+    F = inner(u, v + zero) * dx
+    fd = compute_form_data(F)
+    assert fd.num_coefficients == 1
+
+    # Check that nested zero conditionals are simplifed to zero (rank-0)
+    # and can be added to a rank-1 tensor
+    zero = conditional(z < 0, 0, conditional(z == 0, 0, 0))
+    F = inner(u, v + zero) * dx
+    fd = compute_form_data(F)
+    assert fd.num_coefficients == 1

test/test_evaluate.py

@@ -3,6 +3,8 @@
 
 import math
 
+import numpy as np
+
 from ufl import (
     Argument,
     Coefficient,
@@ -33,12 +35,28 @@
     tr,
     triangle,
 )
-from ufl.constantvalue import as_ufl
+from ufl.constantvalue import ConstantValue, as_ufl
 from ufl.finiteelement import FiniteElement
 from ufl.pullback import identity_pullback
 from ufl.sobolevspace import H1
 
 
+class CustomConstant(ConstantValue):
+    def __init__(self, value):
+        super().__init__()
+        self._value = value
+
+    @property
+    def ufl_shape(self):
+        return ()
+
+    def evaluate(self, x, mapping, component, index_values):
+        return self._value
+
+    def __repr__(self):
+        return f"CustomConstant({self._value})"
+
+
 def testScalars():
     s = as_ufl(123)
     e = s((5, 7))
@@ -132,6 +150,21 @@ def testAlgebra():
     assert e == v
 
 
+def testConstant():
+    """Test that constant division doesn't discard the complex type in the case the value is
+    a numpy complex type, not a native python complex type.
+    """
+    _a = np.complex128(1 + 1j)
+    _b = np.complex128(-3 + 2j)
+    a = CustomConstant(_a)
+    b = CustomConstant(_b)
+    expr = a / b
+    e = expr(())
+
+    expected = complex(_a) / complex(_b)
+    assert e == expected
+
+
 def testIndexSum():
     cell = triangle
     domain = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))

test/test_interpolate.py

@@ -10,6 +10,7 @@
     Adjoint,
     Argument,
     Coefficient,
+    Cofunction,
     FunctionSpace,
     Mesh,
     TestFunction,
@@ -69,6 +70,20 @@ def test_symbolic(V1, V2):
     assert Iu.ufl_operands == (u,)
 
 
+def test_symbolic_adjoint(V1, V2):
+    # Set dual of V2
+    V2_dual = V2.dual()
+
+    u = Argument(V1, 1)
+    vstar = Cofunction(V2_dual)
+    Iu = Interpolate(u, vstar)
+
+    assert Iu.ufl_function_space() == V2_dual
+    assert Iu.argument_slots() == (vstar, u)
+    assert Iu.arguments() == (u,)
+    assert Iu.ufl_operands == (u,)
+
+
 def test_action_adjoint(V1, V2):
     # Set dual of V2
     V2_dual = V2.dual()

test/test_simplify.py

@@ -32,11 +32,13 @@
     triangle,
 )
 from ufl.algorithms import compute_form_data
-from ufl.core.multiindex import FixedIndex, MultiIndex
+from ufl.constantvalue import Zero
+from ufl.core.multiindex import FixedIndex, Index, MultiIndex, indices
 from ufl.finiteelement import FiniteElement
 from ufl.indexed import Indexed
 from ufl.pullback import identity_pullback
 from ufl.sobolevspace import H1
+from ufl.tensors import ComponentTensor, ListTensor
 
 
 def xtest_zero_times_argument(self):
@@ -193,3 +195,76 @@ def test_nested_indexed(self):
     multiindex = MultiIndex((FixedIndex(0),))
     assert Indexed(expr, multiindex) is expr[0]
     assert Indexed(expr, multiindex) is comps[1]
+
+
+def test_repeated_indexing(self):
+    # Test that an Indexed with repeated indices does not contract indices
+    shape = (2, 2)
+    element = FiniteElement("Lagrange", triangle, 1, shape, identity_pullback, H1)
+    domain = Mesh(FiniteElement("Lagrange", triangle, 1, (2,), identity_pullback, H1))
+    space = FunctionSpace(domain, element)
+    x = Coefficient(space)
+    C = as_tensor([x, x])
+
+    fi = FixedIndex(0)
+    i = Index()
+    ii = MultiIndex((fi, i, i))
+    expr = Indexed(C, ii)
+    assert i.count() in expr.ufl_free_indices
+    assert isinstance(expr, Indexed)
+    B, jj = expr.ufl_operands
+    assert B is x
+    assert tuple(jj) == tuple(ii[1:])
+
+
+def test_untangle_indexed_component_tensor(self):
+    shape = (2, 2, 2, 2)
+    element = FiniteElement("Lagrange", triangle, 1, shape, identity_pullback, H1)
+    domain = Mesh(FiniteElement("Lagrange", triangle, 1, (2,), identity_pullback, H1))
+    space = FunctionSpace(domain, element)
+    C = Coefficient(space)
+
+    r = len(shape)
+    kk = indices(r)
+
+    # Untangle as_tensor(C[kk], kk) -> C
+    B = as_tensor(Indexed(C, MultiIndex(kk)), kk)
+    assert B is C
+
+    # Untangle as_tensor(C[kk], jj)[ii] -> C[ll]
+    jj = kk[2:]
+    A = as_tensor(Indexed(C, MultiIndex(kk)), jj)
+    assert A is not C
+
+    ii = kk
+    expr = Indexed(A, MultiIndex(ii))
+    assert isinstance(expr, Indexed)
+    B, ll = expr.ufl_operands
+    assert B is C
+
+    rep = dict(zip(jj, ii))
+    expected = tuple(rep.get(k, k) for k in kk)
+    assert tuple(ll) == expected
+
+
+def test_simplify_indexed(self):
+    element = FiniteElement("Lagrange", triangle, 1, (3,), identity_pullback, H1)
+    domain = Mesh(FiniteElement("Lagrange", triangle, 1, (2,), identity_pullback, H1))
+    space = FunctionSpace(domain, element)
+    u = Coefficient(space)
+    z = Zero(())
+    i = Index()
+    j = Index()
+    # ListTensor
+    lt = ListTensor(z, z, u[1])
+    assert Indexed(lt, MultiIndex((FixedIndex(2),))) == u[1]
+    # ListTensor -- nested
+    l0 = ListTensor(z, u[1], z)
+    l1 = ListTensor(z, z, u[2])
+    l2 = ListTensor(u[0], z, z)
+    ll = ListTensor(l0, l1, l2)
+    assert Indexed(ll, MultiIndex((FixedIndex(1), FixedIndex(2)))) == u[2]
+    assert Indexed(ll, MultiIndex((FixedIndex(2), i))) == l2[i]
+    # ComponentTensor + ListTensor
+    c = ComponentTensor(Indexed(ll, MultiIndex((i, j))), MultiIndex((j, i)))
+    assert Indexed(c, MultiIndex((FixedIndex(1), FixedIndex(2)))) == l2[1]

ufl/algebra.py

@@ -253,12 +253,7 @@ def evaluate(self, x, mapping, component, index_values):
         a, b = self.ufl_operands
         a = a.evaluate(x, mapping, component, index_values)
         b = b.evaluate(x, mapping, component, index_values)
-        # Avoiding integer division by casting to float
-        try:
-            e = float(a) / float(b)
-        except TypeError:
-            e = complex(a) / complex(b)
-        return e
+        return a / b
 
     def __str__(self):
         """Format as a string."""

ufl/algorithms/apply_derivatives.py

@@ -18,7 +18,6 @@
 from ufl.checks import is_cellwise_constant
 from ufl.classes import (
     Coefficient,
-    ComponentTensor,
     Conj,
     ConstantValue,
     ExprList,
@@ -219,28 +218,12 @@ def variable(self, o, df, unused_l):
 
     # --- Indexing and component handling
 
-    def indexed(self, o, Ap, ii):  # TODO: (Partially) duplicated in nesting rules
+    def indexed(self, o, Ap, ii):
         """Differentiate an indexed."""
         # Propagate zeros
         if isinstance(Ap, Zero):
             return self.independent_operator(o)
 
-        # Untangle as_tensor(C[kk], jj)[ii] -> C[ll] to simplify
-        # resulting expression
-        if isinstance(Ap, ComponentTensor):
-            B, jj = Ap.ufl_operands
-            if isinstance(B, Indexed):
-                C, kk = B.ufl_operands
-                kk = list(kk)
-                if all(j in kk for j in jj):
-                    rep = dict(zip(jj, ii))
-                    Cind = [rep.get(k, k) for k in kk]
-                    expr = Indexed(C, MultiIndex(tuple(Cind)))
-                    assert expr.ufl_free_indices == o.ufl_free_indices
-                    assert expr.ufl_shape == o.ufl_shape
-                    return expr
-
-        # Otherwise a more generic approach
         r = len(Ap.ufl_shape) - len(ii)
         if r:
             kk = indices(r)
@@ -1450,29 +1433,12 @@ def base_form_coordinate_derivative(self, o, f, dummy_w, dummy_v, dummy_cd):
             o_[3],
         )
 
-    def indexed(self, o, Ap, ii):  # TODO: (Partially) duplicated in generic rules
+    def indexed(self, o, Ap, ii):
         """Apply to an indexed."""
         # Reuse if untouched
         if Ap is o.ufl_operands[0]:
             return o
 
-        # Untangle as_tensor(C[kk], jj)[ii] -> C[ll] to simplify
-        # resulting expression
-        if isinstance(Ap, ComponentTensor):
-            B, jj = Ap.ufl_operands
-            if isinstance(B, Indexed):
-                C, kk = B.ufl_operands
-
-                kk = list(kk)
-                if all(j in kk for j in jj):
-                    rep = dict(zip(jj, ii))
-                    Cind = [rep.get(k, k) for k in kk]
-                    expr = Indexed(C, MultiIndex(tuple(Cind)))
-                    assert expr.ufl_free_indices == o.ufl_free_indices
-                    assert expr.ufl_shape == o.ufl_shape
-                    return expr
-
-        # Otherwise a more generic approach
         r = len(Ap.ufl_shape) - len(ii)
         if r:
             kk = indices(r)

ufl/algorithms/check_arities.py

@@ -57,7 +57,8 @@ def sum(self, o, a, b):
         """Apply to sum."""
         if a != b:
             raise ArityMismatch(
-                f"Adding expressions with non-matching form arguments {_afmt(a)} vs {_afmt(b)}."
+                f"Adding expressions with non-matching form arguments "
+                f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
             )
         return a
 
@@ -86,7 +87,7 @@ def product(self, o, a, b):
             if len(c) != len(a) + len(b) or len(c) != len({x[0] for x in c}):
                 raise ArityMismatch(
                     "Multiplying expressions with overlapping form arguments "
-                    f"{_afmt(a)} vs {_afmt(b)}."
+                    f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
                 )
             # It's fine for argument parts to overlap
             return c
@@ -138,7 +139,7 @@ def variable(self, o, f, a):
     def conditional(self, o, c, a, b):
         """Apply to conditional."""
         if c:
-            raise ArityMismatch(f"Condition cannot depend on form arguments ({_afmt(a)}).")
+            raise ArityMismatch("Condition cannot depend on form arguments.")
         if a and isinstance(o.ufl_operands[2], Zero):
             # Allow conditional(c, arg, 0)
             return a
@@ -153,7 +154,7 @@ def conditional(self, o, c, a, b):
             # conditional(c, test, nonzeroconstant)
             raise ArityMismatch(
                 "Conditional subexpressions with non-matching form arguments "
-                f"{_afmt(a)} vs {_afmt(b)}."
+                f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
             )
 
     def linear_indexed_type(self, o, a, i):

ufl/algorithms/expand_indices.py

@@ -137,7 +137,7 @@ def index_sum(self, x):
 
         # TODO: For the list tensor purging algorithm, do something like:
         # if index not in self._to_expand:
-        #     return self.expr(x, *[self.visit(o) for o in x.ufl_operands])
+        #     return self.expr(x, *map(self.visit, x.ufl_operands))
 
         for value in range(x.dimension()):
             self._index2value.push(index, value)