Group integrals with common integrand in IntegralData

.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
@@ -77,7 +79,7 @@ jobs:
       - name: Install Basix and FFCx
         run: |
           python3 -m pip install git+https://github.com/FEniCS/basix.git
-          python3 -m pip install git+https://github.com/FEniCS/ffcx.git
+          python3 -m pip install git+https://github.com/FEniCS/ffcx.git@dokken/group_integrals_v2
 
       - name: Clone DOLFINx
         uses: actions/checkout@v3

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

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
 
 

ufl/algorithms/domain_analysis.py

@@ -12,11 +12,13 @@
 import ufl
 from ufl.integral import Integral
 from ufl.form import Form
+from ufl.protocols import id_or_none
 from ufl.sorting import cmp_expr, sorted_expr
 from ufl.utils.sorting import canonicalize_metadata, sorted_by_key
 from ufl.algorithms.coordinate_derivative_helpers import (
     attach_coordinate_derivatives, strip_coordinate_derivatives)
 import numbers
+import typing
 
 
 class IntegralData(object):
@@ -134,14 +136,18 @@ def integral_subdomain_ids(integral):
         raise ValueError(f"Invalid domain id {did}.")
 
 
-def rearrange_integrals_by_single_subdomains(integrals, do_append_everywhere_integrals):
+def rearrange_integrals_by_single_subdomains(
+        integrals: typing.List[Integral],
+        do_append_everywhere_integrals: bool) -> typing.Dict[int, typing.List[Integral]]:
     """Rearrange integrals over multiple subdomains to single subdomain integrals.
 
     Input:
-        integrals: list(Integral)
+        integrals: List of integrals
+        do_append_everywhere_integrals: Boolean indicating if integrals defined on the whole domain should
+                                        just be restricted to the set of input subdomain ids.
 
     Output:
-        integrals: dict: subdomain_id -> list(Integral) (reconstructed with single subdomain_id)
+        integrals: The integrals reconstructed with single subdomain_id
     """
     # Split integrals into lists of everywhere and subdomain integrals
     everywhere_integrals = []
@@ -231,23 +237,37 @@ def build_integral_data(integrals):
     """
     itgs = defaultdict(list)
 
+    # --- Merge integral data that has the same integrals,
+    unique_integrals = defaultdict(tuple)
+    metadata_table = defaultdict(dict)
     for integral in integrals:
-        domain = integral.ufl_domain()
+        integrand = integral.integrand()
         integral_type = integral.integral_type()
+        ufl_domain = integral.ufl_domain()
+        metadata = integral.metadata()
+        meta_hash = hash(canonicalize_metadata(metadata))
         subdomain_id = integral.subdomain_id()
+        subdomain_data = id_or_none(integral.subdomain_data())
         if subdomain_id == "everywhere":
             raise ValueError("'everywhere' not a valid subdomain id.  Did you forget to call group_form_integrals?")
+        unique_integrals[(integral_type, ufl_domain, meta_hash, integrand, subdomain_data)] += (subdomain_id,)
+        metadata_table[(integral_type, ufl_domain, meta_hash, integrand, subdomain_data)] = metadata
+
+    for integral_data, subdomain_ids in unique_integrals.items():
+        (integral_type, ufl_domain, metadata, integrand, subdomain_data) = integral_data
+
+        integral = Integral(integrand, integral_type, ufl_domain, subdomain_ids,
+                            metadata_table[integral_data], subdomain_data)
         # Group for integral data (One integral data object for all
-        # integrals with same domain, itype, subdomain_id (but
-        # possibly different metadata).
-        itgs[(domain, integral_type, subdomain_id)].append(integral)
+        # integrals with same domain, itype, (but possibly different metadata).
+        itgs[(ufl_domain, integral_type, subdomain_ids)].append(integral)
 
     # Build list with canonical ordering, iteration over dicts
     # is not deterministic across python versions
     def keyfunc(item):
         (d, itype, sid), integrals = item
-        return (d._ufl_sort_key_(), itype, (type(sid).__name__, sid))
-
+        sid_int = tuple(-1 if i == "otherwise" else i for i in sid)
+        return (d._ufl_sort_key_(), itype, (type(sid).__name__, ), sid_int)
     integral_datas = []
     for (d, itype, sid), integrals in sorted(itgs.items(), key=keyfunc):
         integral_datas.append(IntegralData(d, itype, sid, integrals, {}))
@@ -262,8 +282,7 @@ def group_form_integrals(form, domains, do_append_everywhere_integrals=True):
     :returns: A new :class:`~.Form` with gathered integrands.
     """
     # Group integrals by domain and type
-    integrals_by_domain_and_type = \
-        group_integrals_by_domain_and_type(form.integrals(), domains)
+    integrals_by_domain_and_type = group_integrals_by_domain_and_type(form.integrals(), domains)
 
     integrals = []
     for domain in domains:
@@ -277,8 +296,8 @@ def group_form_integrals(form, domains, do_append_everywhere_integrals=True):
             #   f*dx((1,2)) + g*dx((2,3)) -> f*dx(1) + (f+g)*dx(2) + g*dx(3)
             # (note: before this call, 'everywhere' is a valid subdomain_id,
             # and after this call, 'otherwise' is a valid subdomain_id)
-            single_subdomain_integrals = \
-                rearrange_integrals_by_single_subdomains(ddt_integrals, do_append_everywhere_integrals)
+            single_subdomain_integrals = rearrange_integrals_by_single_subdomains(
+                ddt_integrals, do_append_everywhere_integrals)
 
             for subdomain_id, ss_integrals in sorted_by_key(single_subdomain_integrals):
 
@@ -307,8 +326,7 @@ def calc_hash(cd):
 
                     # Accumulate integrands of integrals that share the
                     # same compiler data
-                    integrands_and_cds = \
-                        accumulate_integrands_with_same_metadata(samecd_integrals[1])
+                    integrands_and_cds = accumulate_integrands_with_same_metadata(samecd_integrals[1])
 
                     for integrand, metadata in integrands_and_cds:
                         integral = Integral(integrand, integral_type, domain,

ufl/form.py

@@ -10,7 +10,9 @@
 # Modified by Anders Logg, 2009-2011.
 # Modified by Massimiliano Leoni, 2016.
 # Modified by Cecile Daversin-Catty, 2018.
+# Modified by Jørgen S. Dokken 2023.
 
+import numbers
 import warnings
 from collections import defaultdict
 from itertools import chain
@@ -51,10 +53,18 @@ def _sorted_integrals(integrals):
 
     all_integrals = []
 
+    def keyfunc(item):
+        if isinstance(item, numbers.Integral):
+            sid_int = item
+        else:
+            # As subdomain ids can be either int or tuples, we need to compare them
+            sid_int = tuple(-1 if i == "otherwise" else i for i in item)
+        return (type(item).__name__, sid_int)
+
     # Order integrals canonically to increase signature stability
     for d in sort_domains(integrals_dict):
         for it in sorted(integrals_dict[d]):  # str is sortable
-            for si in sorted(integrals_dict[d][it], key=lambda x: (type(x).__name__, x)):  # int/str are sortable
+            for si in sorted(integrals_dict[d][it], key=keyfunc):
                 unsorted_integrals = integrals_dict[d][it][si]
                 # TODO: At this point we could order integrals by
                 #       metadata and integrand, or even add the