Merge fix_dualspace

.github/workflows/fenicsx-tests.yml

@@ -41,6 +41,8 @@ jobs:
         with:
           path: ./ffcx
           repository: FEniCS/ffcx
+          ref: dokken/group_integrals_v2
+
       - name: Install FFCx
         run: |
           cd ffcx

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_cell.py

@@ -37,6 +37,20 @@ def test_hexahedron():
     assert cell.num_faces() == 6
 
 
+def test_pentatope():
+    cell = ufl.pentatope
+    assert cell.num_vertices() == 5
+    assert cell.num_edges() == 10
+    assert cell.num_faces() == 10
+
+
+def test_tesseract():
+    cell = ufl.tesseract
+    assert cell.num_vertices() == 16
+    assert cell.num_edges() == 32
+    assert cell.num_faces() == 24
+
+
 
 @pytest.mark.parametrize("cell", [ufl.interval])
 def test_cells_1d(cell):
@@ -53,12 +67,19 @@ def test_cells_2d(cell):
 
 
 @pytest.mark.parametrize("cell", [ufl.tetrahedron, ufl.hexahedron, ufl.prism, ufl.pyramid])
-def test_cells_2d(cell):
+def test_cells_3d(cell):
     assert cell.num_facets() == cell.num_faces()
     assert cell.num_ridges() == cell.num_edges()
     assert cell.num_peaks() == cell.num_vertices()
 
 
+@pytest.mark.parametrize("cell", [ufl.tesseract, ufl.pentatope])
+def test_cells_4d(cell):
+    assert cell.num_facets() == cell.num_sub_entities(3)
+    assert cell.num_ridges() == cell.num_faces()
+    assert cell.num_peaks() == cell.num_edges()
+
+
 def test_tensorproductcell():
     orig = ufl.TensorProductCell(ufl.interval, ufl.interval)
     cell = orig.reconstruct()

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_domains.py

@@ -4,6 +4,7 @@
 Tests of domain language and attaching domains to forms.
 """
 
+from mockobjects import MockMesh, MockMeshFunction
 import pytest
 
 from ufl import *
@@ -13,8 +14,6 @@
 all_cells = (interval, triangle, tetrahedron,
              quadrilateral, hexahedron)
 
-from mockobjects import MockMesh, MockMeshFunction
-
 
 def test_construct_domains_from_cells():
     for cell in all_cells:
@@ -198,7 +197,8 @@ def test_everywhere_integrals_with_backwards_compatibility():
 
     # Check some integral data
     assert ida.integral_type == "cell"
-    assert ida.subdomain_id == "otherwise"
+    assert(len(ida.subdomain_id) == 1)
+    assert ida.subdomain_id[0] == "otherwise"
     assert ida.metadata == {}
 
     # Integrands are not equal because of renumbering
@@ -208,6 +208,47 @@ def test_everywhere_integrals_with_backwards_compatibility():
     assert itg1.ufl_element() == itg2.ufl_element()
 
 
+def test_merge_sort_integral_data():
+    D = Mesh(triangle)
+
+    V = FunctionSpace(D, FiniteElement("CG", triangle, 1))
+
+    u = Coefficient(V)
+    c = Constant(D)
+    a = c * dS((2, 4)) + u * dx + u * ds + 2 * u * dx(3) + 2 * c * dS + 2 * u * dx((1, 4))
+    form_data = compute_form_data(a, do_append_everywhere_integrals=False).integral_data
+    assert(len(form_data) == 5)
+
+    # Check some integral data
+    assert form_data[0].integral_type == "cell"
+    assert(len(form_data[0].subdomain_id) == 1)
+    assert form_data[0].subdomain_id[0] == "otherwise"
+    assert form_data[0].metadata == {}
+
+    assert form_data[1].integral_type == "cell"
+    assert(len(form_data[1].subdomain_id) == 3)
+    assert form_data[1].subdomain_id[0] == 1
+    assert form_data[1].subdomain_id[1] == 3
+    assert form_data[1].subdomain_id[2] == 4
+    assert form_data[1].metadata == {}
+
+    assert form_data[2].integral_type == "exterior_facet"
+    assert(len(form_data[2].subdomain_id) == 1)
+    assert form_data[2].subdomain_id[0] == "otherwise"
+    assert form_data[2].metadata == {}
+
+    assert form_data[3].integral_type == "interior_facet"
+    assert(len(form_data[3].subdomain_id) == 1)
+    assert form_data[3].subdomain_id[0] == "otherwise"
+    assert form_data[3].metadata == {}
+
+    assert form_data[4].integral_type == "interior_facet"
+    assert(len(form_data[4].subdomain_id) == 2)
+    assert form_data[4].subdomain_id[0] == 2
+    assert form_data[4].subdomain_id[1] == 4
+    assert form_data[4].metadata == {}
+
+
 def xtest_mixed_elements_on_overlapping_regions():  # Old sketch, not working
 
     # Create domain and both disjoint and overlapping regions
@@ -366,10 +407,10 @@ def xtest_form_domain_model():  # Old sketch, not working
     # domains and regions to be part of their integrands
     dxb = dx('DB')   # Get Mesh by name
     dxbl = dx(Region(DB, (1, 4), 'DBL2'))
-              # Provide a region with different name but same subdomain ids as
-              # DBL
+    # Provide a region with different name but same subdomain ids as
+    # DBL
     dxbr = dx((1, 4))
-              # Assume unique Mesh and provide subdomain ids explicitly
+    # Assume unique Mesh and provide subdomain ids explicitly
 
     # Not checking measure objects in detail, as long as
     # they carry information to construct integrals below
@@ -466,8 +507,8 @@ def sub_elements_on_subdomains(W):
 
     # Disjunctified by UFL:
     alonly = dot(ul, vl) * dx(D1)
-                 # integral_1 knows that only subelement VL is active
+    # integral_1 knows that only subelement VL is active
     am = (dot(ul, vl) + ur * vr) * dx(D2)
-          # integral_2 knows that both subelements are active
+    # integral_2 knows that both subelements are active
     aronly = ur * vr * \
         dx(D3)  # integral_3 knows that only subelement VR is active

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/__init__.py

@@ -53,6 +53,8 @@
     - hexahedron
     - prism
     - pyramid
+    - pentatope
+    - tesseract
 
 * Domains::
 
@@ -183,7 +185,7 @@
     - sqrt
     - exp, ln, erf
     - cos, sin, tan
-    - acos, asin, atan, atan_2
+    - acos, asin, atan, atan2
     - cosh, sinh, tanh
     - bessel_J, bessel_Y, bessel_I, bessel_K
 
@@ -322,7 +324,7 @@
 from ufl.operators import (
     rank, shape, conj, real, imag, outer, inner, dot, cross, perp,
     det, inv, cofac, transpose, tr, diag, diag_vector, dev, skew, sym,
-    sqrt, exp, ln, erf, cos, sin, tan, acos, asin, atan, atan_2, cosh, sinh, tanh,
+    sqrt, exp, ln, erf, cos, sin, tan, acos, asin, atan, atan2, cosh, sinh, tanh,
     bessel_J, bessel_Y, bessel_I, bessel_K, eq, ne, le, ge, lt, gt, And, Or, Not,
     conditional, sign, max_value, min_value, variable, diff,
     Dx, grad, div, curl, rot, nabla_grad, nabla_div, Dn, exterior_derivative,
@@ -343,7 +345,7 @@
 
 # Predefined convenience objects
 from ufl.objects import (
-    vertex, interval, triangle, tetrahedron,
+    vertex, interval, triangle, tetrahedron, pentatope, tesseract,
     quadrilateral, hexahedron, prism, pyramid, facet,
     i, j, k, l, p, q, r, s,
     dx, ds, dS, dP,
@@ -387,7 +389,7 @@
     'transpose', 'tr', 'diag', 'diag_vector', 'dev', 'skew', 'sym',
     'sqrt', 'exp', 'ln', 'erf',
     'cos', 'sin', 'tan',
-    'acos', 'asin', 'atan', 'atan_2',
+    'acos', 'asin', 'atan', 'atan2',
     'cosh', 'sinh', 'tanh',
     'bessel_J', 'bessel_Y', 'bessel_I', 'bessel_K',
     'eq', 'ne', 'le', 'ge', 'lt', 'gt', 'And', 'Or', 'Not',
@@ -404,7 +406,7 @@
     'dc', 'dC', 'dO', 'dI', 'dX',
     'ds_b', 'ds_t', 'ds_tb', 'ds_v', 'dS_h', 'dS_v',
     'vertex', 'interval', 'triangle', 'tetrahedron',
-    'prism', 'pyramid',
+    'prism', 'pyramid', 'pentatope', 'tesseract',
     'quadrilateral', 'hexahedron', 'facet',
     'i', 'j', 'k', 'l', 'p', 'q', 'r', 's',
     'e', 'pi',

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

@@ -358,7 +358,7 @@ def atan(self, o, fp):
         f, = o.ufl_operands
         return fp / (1.0 + f**2)
 
-    def atan_2(self, o, fp, gp):
+    def atan2(self, o, fp, gp):
         f, g = o.ufl_operands
         return (g * fp - f * gp) / (f**2 + g**2)
 
@@ -850,9 +850,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/algorithms/compute_form_data.py

@@ -102,13 +102,12 @@ def _compute_max_subdomain_ids(integral_data):
     max_subdomain_ids = {}
     for itg_data in integral_data:
         it = itg_data.integral_type
-        si = itg_data.subdomain_id
-        if isinstance(si, int):
-            newmax = si + 1
-        else:
-            newmax = 0
-        prevmax = max_subdomain_ids.get(it, 0)
-        max_subdomain_ids[it] = max(prevmax, newmax)
+        for integral in itg_data.integrals:
+            # Convert string for default integral to -1
+            sids = (-1 if isinstance(si, str) else si for si in integral.subdomain_id())
+            newmax = max(sids) + 1
+            prevmax = max_subdomain_ids.get(it, 0)
+            max_subdomain_ids[it] = max(prevmax, newmax)
     return max_subdomain_ids