Merge branch 'main' into pbrubeck/merge-upstream

.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_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()

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/core/base_form_operator.py

@@ -16,6 +16,7 @@
 from collections import OrderedDict
 
 from ufl.argument import Argument, Coargument
+from ufl.coefficient import BaseCoefficient
 from ufl.constantvalue import as_ufl
 from ufl.core.operator import Operator
 from ufl.core.ufl_type import ufl_type
@@ -134,20 +135,19 @@ def count(self):
     def ufl_shape(self):
         """Return the UFL shape of the coefficient.produced by the operator."""
         arg, *_ = self.argument_slots()
-        if isinstance(arg, BaseForm):
+        if not isinstance(arg, BaseCoefficient) and isinstance(arg, (BaseForm, Coargument)):
             arg, *_ = arg.arguments()
         return arg._ufl_shape
 
     def ufl_function_space(self):
         """Return the function space associated to the operator.
 
-        I.e. return the dual of the base form operator's Coargument.
+        I.e. return the dual of the base form operator's Coargument space.
         """
         arg, *_ = self.argument_slots()
-        if isinstance(arg, BaseForm):
+        if not isinstance(arg, BaseCoefficient) and isinstance(arg, (BaseForm, Coargument)):
             arg, *_ = arg.arguments()
-            return arg._ufl_function_space
-        return arg._ufl_function_space.dual()
+        return arg.ufl_function_space()
 
     def _ufl_expr_reconstruct_(
         self, *operands, function_space=None, derivatives=None, argument_slots=None

ufl/core/interpolate.py

@@ -8,14 +8,12 @@
 #
 # Modified by Nacime Bouziani, 2021-2022
 
-from ufl.action import Action
 from ufl.argument import Argument, Coargument
-from ufl.coefficient import Cofunction
 from ufl.constantvalue import as_ufl
 from ufl.core.base_form_operator import BaseFormOperator
 from ufl.core.ufl_type import ufl_type
 from ufl.duals import is_dual
-from ufl.form import BaseForm, Form
+from ufl.form import BaseForm
 from ufl.functionspace import AbstractFunctionSpace
 
 
@@ -35,16 +33,15 @@ def __init__(self, expr, v):
             v: the FunctionSpace to interpolate into or the Coargument
                 defined on the dual of the FunctionSpace to interpolate into.
         """
-        # This check could be more rigorous.
-        dual_args = (Coargument, Cofunction, Form, Action, BaseFormOperator)
+        dual_args = (Coargument, BaseForm)
 
         if isinstance(v, AbstractFunctionSpace):
             if is_dual(v):
                 raise ValueError("Expecting a primal function space.")
             v = Argument(v.dual(), 0)
         elif not isinstance(v, dual_args):
             raise ValueError(
-                "Expecting the second argument to be FunctionSpace, FiniteElement or dual."
+                "Expecting the second argument to be FunctionSpace, Coargument, or BaseForm."
             )
 
         expr = as_ufl(expr)
@@ -54,11 +51,9 @@ def __init__(self, expr, v):
         # Reversed order convention
         argument_slots = (v, expr)
         # Get the primal space (V** = V)
-        if isinstance(v, BaseForm):
-            arg, *_ = v.arguments()
-            function_space = arg.ufl_function_space()
-        else:
-            function_space = v.ufl_function_space().dual()
+        arg, *_ = v.arguments()
+        function_space = arg.ufl_function_space()
+
         # Set the operand as `expr` for DAG traversal purpose.
         operand = expr
         BaseFormOperator.__init__(

ufl/tensors.py

@@ -74,7 +74,6 @@ def sub(e, *indices):
                     return sub(e0, 0) if j == () else sub(e0, 0)[(*j, slice(None))]
             except ValueError:
                 pass
-
         # Simplify [v[0,:], v[1,:], ..., v[k,:]] -> v
         if (
             all(
@@ -85,7 +84,10 @@ def sub(e, *indices):
             and all(sub(e, 0, 0) == sub(e0, 0, 0) for e in expressions[1:])
         ):
             indices = [sub(e, 0, 1).indices() for e in expressions]
-            if all(i[0] == k for k, i in enumerate(indices)):
+            if all(
+                i[0] == k and all(isinstance(subindex, Index) for subindex in i[1:])
+                for k, i in enumerate(indices)
+            ):
                 return sub(e0, 0, 0)
 
         # Construct a new instance to be initialised