On a mission to remove the groups hack

.github/workflows/test.yml

@@ -34,6 +34,7 @@ jobs:
       PYOP2_SPMD_STRICT: 1
       EXTRA_PYTEST_ARGS: --splitting-algorithm least_duration --timeout=600 --timeout-method=thread -o faulthandler_timeout=660 --durations-path=./fuse-repo/test/test_durations.json --durations=50
       PYTEST_MPI_MAX_NPROCS: 8
+      FIREDRAKE_RUN_SPLIT_TESTS_TIMEOUT: 2400
     steps:
       - name: Fix HOME
         # For unknown reasons GitHub actions overwrite HOME to /github/home

fuse/cells.py

@@ -440,22 +440,6 @@ def get_topology(self, renumber=False):
             return self._topology
         return self._topology_unrelabelled
 
-    def get_renumbered_topology(self):
-        structure = [generation for generation in nx.topological_generations(self.graph())]
-        structure.reverse()
-
-        min_ids = [min(dimension) for dimension in structure]
-        vertices = self.ordered_vertices()
-        relabelled_verts = {vertices[i]: i for i in range(len(vertices))}
-
-        self._topology = {}
-        for i in range(len(structure)):
-            dimension = structure[i]
-            self._topology[i] = {}
-            for node in dimension:
-                self._topology[i][node - min_ids[i]] = tuple([relabelled_verts[vert] for vert in self.get_node(node).ordered_vertices()])
-        return self._topology
-
     def get_sub_entities(self):
         min_ids = self.get_starter_ids()
         sub_entities = {d: {e.id - min_ids[d]: [] for e in self.d_entities(d)} for d in range(self.get_spatial_dimension() + 1)}
@@ -540,7 +524,6 @@ def ordered_vertices(self, get_class=False):
             flatten = itertools.chain.from_iterable(full_list)
             verts = list(dict.fromkeys(flatten))
             if self.oriented:
-                # make sure this is necessary
                 return self.oriented.permute(verts)
             return verts
 

fuse/groups.py

@@ -29,6 +29,30 @@ def perm_list_to_matrix(identity, perm):
     return res
 
 
+def is_normal_subgroup(subgroup, group):
+    members = {h.perm for h in subgroup.members()}
+    return all(((~x) * h * x).perm in members
+               for x in group.members() for h in subgroup.members())
+
+
+def sub_entity_cone_offset(cell):
+    """Offset between how a cell lists its sub entities and the convention that sub
+    entity k is the one opposite vertex k.
+
+    An orientation value is computed by comparing a cell's cone of sub entities, while
+    numeric_rep labels group members by vertex order. The two agree only under the
+    opposite-vertex convention, so this offset relates them.
+    """
+    verts = cell.ordered_vertices()
+    opposite = []
+    for connection in cell.connections:
+        missing = [v for v in verts if v not in set(connection.ordered_vertices())]
+        if len(missing) != 1:
+            raise ValueError("Cone offset is only defined for a simplex")
+        opposite.append(verts.index(missing[0]))
+    return Permutation(opposite)
+
+
 class GroupMemberRep(object):
 
     def __init__(self, perm, M, group):
@@ -117,40 +141,29 @@ def matrix_form_subgroup(self, group):
         if group.size() == 1:
             # Trivial case
             return np.array([1])
-        if group.size() == 6 and self.group.size() == 6:
-            # horrible hack for S3
-            members = [m.numeric_rep() for m in group.members()]
-            permuted_members = [((m)*(~self)).numeric_rep() for m in group.members()]
-            mapping = {4: 4, 3: 0, 0: 3}
-            if (~self).numeric_rep() in mapping.keys():
-                n = self.group.get_member_by_val(mapping[(~self).numeric_rep()])
-                permuted_members = [((m)*(~n)).numeric_rep() for m in group.members()]
-            mat = perm_list_to_matrix(members, permuted_members)
-        elif group.size() == self.group.size():
+        if group.size() == self.group.size():
+            # A free orbit of the full symmetry group
+            # Uses sub_entity_cone_offset to reconcile FUSE facet numbering with FIAT.
+            w = self.group.get_member(sub_entity_cone_offset(self.group.cell))
+            oriented = w * (~self) * w
             members = [m.numeric_rep() for m in group.members()]
-            permuted_members = [(m*(~self)).numeric_rep() for m in group.members()]
+            permuted_members = [(m*oriented).numeric_rep() for m in group.members()]
             mat = perm_list_to_matrix(members, permuted_members)
         elif group.size() == self.perm.size:
-            if self.perm.size == 3:
+            if is_normal_subgroup(group, self.group):
+                # Products leave the orbit group, so they are projected back through the
+                # coset section. That section commutes with the translation only when the
+                # subgroup is normal
                 cosets = self.group.cosets_by_submember(group)
                 members = [cosets[m.array_form].numeric_rep() for m in group.members()]
                 permuted_members = [cosets[(m*(~self)).array_form].numeric_rep() for m in group.members()]
-                mapping = {4: 3, 3: 4, 0: 0}
-                # mapping = {4: 4, 3: 0, 0: 3}
-                if (~self).numeric_rep() in mapping.keys():
-                    n = self.group.get_member_by_val(mapping[(~self).numeric_rep()])
-                    permuted_members = [cosets[(m*(~n)).array_form].numeric_rep() for m in group.members()]
                 mat = perm_list_to_matrix(members, permuted_members)
             else:
                 mat = np.array(PermutationMatrix(self.perm).as_explicit()).astype(np.float64)
         elif group.size() < self.group.size():
             members = [m.numeric_rep() for m in group.members()]
             permuted_members = [(m*(~self)).numeric_rep() for m in group.members()]
             mat = perm_list_to_matrix(members, permuted_members)
-            # cosets = self.group.cosets_by_submember(group)
-            # members = [cosets[m.array_form].numeric_rep() for m in group.members()]
-            # permuted_members = [cosets[(m*(~self)).array_form].numeric_rep() for m in group.members()]
-            # mat = perm_list_to_matrix(members, permuted_members)
         else:
             raise NotImplementedError("Complex subgroups where group size is not the same as perm size are not supported")
         return mat
@@ -217,7 +230,9 @@ def conjugacy_class(self, g):
         return conj_class
 
     def cosets(self, subset):
-        # Divides current group by given subset
+        # Divides current group by given subset into left cosets gH.
+        # The g*h order is load bearing: cosets_by_submember relies on it to
+        # recover the right factor h, and reversing it is not a relabelling.
         # can be modified to allow members of given subset not to exist in group self
         seen = self.members().copy()
         cosets = []
@@ -235,6 +250,7 @@ def cosets(self, subset):
         return cosets
 
     def cosets_by_submember(self, subset):
+        # Maps each member x = g*h of self to the right factor h in subset.
         cosets = self.cosets(subset)
         cosets_by_submember = {}
         for i, m in enumerate(subset.members()):

test/test_cells.py

@@ -220,6 +220,70 @@ def test_tensor_connectivity():
         assert all(connectivity[i] == t for i, t in topology.items())
 
 
+@pytest.mark.parametrize(["cell"], [(polygon(3),), (make_tetrahedron(),)])
+def test_fiat_topology_vertex_frame(cell):
+    # The vertex coordinates handed to FIAT and the vertex labels in the topology have
+    # to be in the same frame. get_topology(renumber=False) labels vertices by id, so
+    # verts must be in id order too. Mixing the frames silently attaches entities to
+    # the wrong vertices rather than raising.
+    fiat_cell = cell.to_fiat()
+    verts = np.array(fiat_cell.get_vertices(), dtype=float)
+    topology = fiat_cell.get_topology()
+
+    assert all(key == label for key, (label,) in topology[0].items())
+
+    min_ids = cell.get_starter_ids()
+    for dim in range(1, cell.get_spatial_dimension() + 1):
+        for entity in cell.d_entities(dim):
+            labels = topology[dim][entity.id - min_ids[dim]]
+            expected = np.array([cell.get_node(v, return_coords=True)
+                                 for v in entity.ordered_vertices()], dtype=float)
+            assert np.allclose(np.sort(verts[list(labels)], axis=0),
+                               np.sort(expected, axis=0))
+
+
+def test_vertex_frames_differ():
+    # ordered_vertices() is a connectivity traversal, vertices() is id sorted. They are
+    # not the same on the tetrahedron, so anything consuming one of them is committed to
+    # that frame. The group representations use ordered_vertices, to_fiat uses vertices.
+    cell = make_tetrahedron()
+    ids = [v.id for v in cell.vertices()]
+    ordered = cell.ordered_vertices()
+
+    assert sorted(ordered) == ids
+    assert ordered != ids
+    assert [ids.index(v) for v in ordered] == [2, 0, 1, 3]
+
+    assert np.allclose(cell.to_fiat().get_vertices(), cell.vertices(return_coords=True))
+
+
+@pytest.mark.xfail(strict=True,
+                   reason="construct_subelement builds sub entities with renumber=True "
+                          "(ordered vertex labels) while the cell uses renumber=False (id "
+                          "labels) and verts never follows renumber, so get_entity_transform "
+                          "maps a sub entity's vertices onto a permutation of the parent's "
+                          "view of them. The vertex set is right, the correspondence is not.")
+@pytest.mark.parametrize(["cell"], [(polygon(3),), (make_tetrahedron(),)])
+def test_entity_transform_preserves_vertex_order(cell):
+    # get_entity_transform maps a sub entity's reference vertices into the parent. It has
+    # to land on the parent's own view of that entity in the same order, otherwise the two
+    # disagree about which vertex is which.
+    fiat_cell = cell.to_fiat()
+    topology = fiat_cell.get_topology()
+    disagree = []
+    for dim in range(1, cell.get_spatial_dimension()):
+        for entity in range(len(topology[dim])):
+            sub = fiat_cell.construct_subelement(dim, entity)
+            transform = fiat_cell.get_entity_transform(dim, entity)
+            mapped = np.array([transform(v) for v in sub.get_vertices()], dtype=float)
+            parent = np.array(fiat_cell.get_vertices_of_subcomplex(topology[dim][entity]),
+                              dtype=float)
+            assert np.allclose(np.sort(mapped, axis=0), np.sort(parent, axis=0))
+            if not np.allclose(mapped, parent):
+                disagree.append((dim, entity))
+    assert disagree == []
+
+
 @pytest.mark.parametrize(["cell"], [(ufc_triangle(),), (polygon(3),), (make_tetrahedron(), ), (make_tetrahedron(), )])
 def test_new_connectivity(cell):
     cell = cell.to_fiat()

test/test_convert_to_fiat.py

@@ -668,7 +668,7 @@ def test_projection_convergence_3d(elem_gen, elem_code, deg, conv_rate):
                                                               (construct_tet_ned, "N1curl", 1, 0.8),
                                                               (construct_tet_rt2, "RT", 2, 1.8),
                                                               (construct_tet_ned2, "N1curl", 2, 1.8),
-                                                              (lambda cell: periodic_table(1, 3, 1, 3), "N2curl", 3, 3.8)])
+                                                              pytest.param(lambda cell: periodic_table(1, 3, 1, 3), "N2curl", 3, 3.8, id="N1-3"),])
 def test_const_vec(elem_gen, elem_code, deg, conv_rate):
     cell = make_tetrahedron()
     elem = elem_gen(cell)
@@ -1122,6 +1122,60 @@ def expr(mesh):
     assert all([res < max_err for res in errors])
 
 
+def _one_form_norm_spread(ufl_elem, is_vector, mesh_factory, perms):
+    # For a transformation that is a signed permutation, the norm of the assembled
+    # vector should be fixed.
+    norms = []
+    for g in perms:
+        mesh = mesh_factory(g)
+        V = FunctionSpace(mesh, ufl_elem)
+        v = TestFunction(V)
+        x = SpatialCoordinate(mesh)
+        f = as_vector((2, 3, 5)) if is_vector else x[0] + 2*x[1] + 3*x[2]
+        b = assemble(inner(f, v)*dx)
+        norms.append(float(np.linalg.norm(np.asarray(b.dat.data_ro).reshape(-1))))
+        print(g, norms[-1])
+    norms = np.array(norms)
+    return norms, norms.max() - norms.min()
+
+
+_TET_ONE_FORM_PERMS = [
+    Permutation([0, 1, 2, 3]),
+    Permutation([0, 2, 3, 1]),
+    Permutation([0, 3, 1, 2]),
+    Permutation([0, 1, 3, 2]),
+    Permutation([0, 3, 2, 1]),
+    Permutation([0, 2, 1, 3]),
+]
+
+
+@pytest.mark.parametrize("elem_gen", [
+    pytest.param(construct_tet_cg4, id="CG-4"),
+    pytest.param(construct_tet_rt2, id="RT-2"),
+    pytest.param(construct_tet_ned_2nd_kind_2, id="N2curl-2"),
+    # CG-6 is the cheapest element reaching the six member branch of
+    # matrix_form_subgroup, which the three above do not exercise.
+    pytest.param(construct_tet_cg6, id="CG-6"),
+])
+def test_two_tet_one_form_orientation_invariance(elem_gen):
+    # construct_tet_ned2 (1st-kind Nedelec deg 2) and construct_tet_ned_2nd_kind_3
+    # are deliberately excluded: their face orientation matrices are not signed
+    # permutations, so the norm is not expected to be fixed
+    #
+    # Blind spot worth knowing: TwoTetMesh shares only face 0, whatever the permutation,
+    # so this test and test_two_tet_projection cannot see a face dependent orientation
+    # error. They also only check that two cells agree, which fixes the orientation
+    # family up to a constant. test_const_vec[N1-3] interpolates a constant on
+    # UnitCubeMesh and does catch both.
+    from firedrake.utility_meshes import TwoTetMesh
+    elem = elem_gen()
+    ufl_elem = elem.to_ufl()
+    is_vector = len(elem.get_value_shape()) > 0
+    _, spread = _one_form_norm_spread(
+        ufl_elem, is_vector, lambda g: TwoTetMesh(perm=g, use_fuse=True), _TET_ONE_FORM_PERMS)
+    assert spread < 1e-10
+
+
 @pytest.mark.parametrize("elem_gen,elem_code,deg",
                          [(construct_tet_cg4, "CG", 4), (construct_tet_rt2, "RT", 2), (construct_tet_ned2, "N1curl", 2), (construct_tet_bdm2, "BDM", 2),
                           ])

test/test_groups.py

@@ -3,6 +3,9 @@
 from fuse.groups import perm_matrix_to_perm_array
 from sympy.combinatorics import Permutation
 from test_convert_to_fiat import create_dg1
+from test_3d_examples_docs import (construct_tet_cg4, construct_tet_cg6, construct_tet_rt2,
+                                   construct_tet_ned_2nd_kind_2, construct_tet_ned_2nd_kind_3)
+from FIAT.orientation_utils import make_entity_permutations_simplex
 import numpy as np
 
 
@@ -117,3 +120,108 @@ def test_perm_mat_conversion():
         mat_form = g.matrix_form()
         array_form = perm_matrix_to_perm_array(mat_form)
         assert np.allclose(g.perm.array_form, array_form)
+
+
+def test_coset_convention():
+    # Cosets are the left cosets gH, and cosets_by_submember recovers the right
+    # factor h of x = g*h. matrix_form_subgroup depends on this order, so it is
+    # pinned here: reversing the product is not a relabelling.
+    cell = polygon(3)
+    cS3 = S3.add_cell(cell)
+    cC3 = C3.add_cell(cell)
+
+    sub_members = cC3.members()
+    assert sub_members[0].perm.is_Identity
+
+    cosets = cS3.cosets(cC3)
+    by_submember = cS3.cosets_by_submember(cC3)
+
+    seen = []
+    for coset in cosets:
+        assert len(coset) == len(sub_members)
+        for i, h in enumerate(sub_members):
+            assert coset[i] == coset[0] * h
+            assert by_submember[coset[i].array_form] == h
+        seen += [m.array_form for m in coset]
+    assert sorted(seen) == sorted([m.array_form for m in cS3.members()])
+
+    assert {k: v.array_form for k, v in by_submember.items()} == {
+        (0, 1, 2): (0, 1, 2), (0, 2, 1): (0, 1, 2),
+        (1, 2, 0): (1, 2, 0), (1, 0, 2): (1, 2, 0),
+        (2, 0, 1): (2, 0, 1), (2, 1, 0): (2, 0, 1)}
+
+
+def test_face_matrices_match_fiat():
+    # The orientation numbering FUSE produces for a triangular face has to agree
+    # with the one FIAT indexes its entity permutations by.
+    elem = construct_tet_cg4()
+    elem.to_ufl()
+    fiat_perms = make_entity_permutations_simplex(2, 2)
+
+    for face, dof_ids in elem.entity_ids[2].items():
+        for o, mat in elem.matrices[2][face].items():
+            block = np.asarray(mat)[np.ix_(dof_ids, dof_ids)]
+            assert perm_matrix_to_perm_array(block) == fiat_perms[o]
+
+
+SIX_MEMBER_XFAIL = mark.xfail(
+    strict=True,
+    reason="A free orbit of the full symmetry group legitimately violates this. The six "
+           "member branch of matrix_form_subgroup sends orientation 3, not 0, to the "
+           "identity permutation, so M(0) is not the identity and the family is not a "
+           "homomorphism. A replacement satisfying both was tried (6dacc09) and broke "
+           "test_const_vec[N1-3]: facet matching only fixes the family up to a left "
+           "factor, and that replacement picked a different, face dependent one.")
+
+
+def entity_orientation_matrices(elem, dim, entity=0):
+    dof_ids = elem.entity_ids[dim][entity]
+    return {o: np.asarray(mat)[np.ix_(dof_ids, dof_ids)]
+            for o, mat in elem.matrices[dim][entity].items()}
+
+
+@mark.parametrize("elem_gen", [
+    param(construct_tet_cg4, id="CG-4"),
+    param(construct_tet_rt2, id="RT-2"),
+    param(construct_tet_ned_2nd_kind_2, id="N2curl-2"),
+    param(construct_tet_cg6, id="CG-6", marks=SIX_MEMBER_XFAIL),
+    param(construct_tet_ned_2nd_kind_3, id="N2curl-3", marks=SIX_MEMBER_XFAIL),
+])
+def test_orientation_matrix_is_representation(elem_gen):
+    # Entities whose DOFs come from a coset orbit carry a representation of their
+    # symmetry group: the identity orientation does not move DOFs, and composing two
+    # orientations composes their matrices. A free orbit of the full symmetry group
+    # does not - see SIX_MEMBER_XFAIL.
+    elem = elem_gen()
+    elem.to_ufl()
+
+    for dim in range(1, elem.cell.dim()):
+        if not elem.entity_ids[dim][0]:
+            continue
+        mats = entity_orientation_matrices(elem, dim)
+        assert np.allclose(mats[0], np.eye(mats[0].shape[0]))
+
+        members = elem.cell.d_entities(dim)[0].group.members()
+        for a in members:
+            for b in members:
+                composed = mats[(a * b).numeric_rep()]
+                assert np.allclose(composed, mats[b.numeric_rep()] @ mats[a.numeric_rep()])
+
+
+def test_s3():
+    # Golden pin on the six member branch. These matrices are the ones proven correct by
+    # the full suite, so any change to that branch has to reproduce them exactly. Runs in
+    # milliseconds, unlike the assembly tests that actually discriminate.
+    from fuse.groups import perm_list_to_matrix
+
+    s3 = S3.add_cell(polygon(3))
+    for g in s3.members():
+        members = [m.numeric_rep() for m in s3.members()]
+        permuted_members = [((m)*(~g)).numeric_rep() for m in s3.members()]
+        mapping = {4: 4, 3: 0, 0: 3}
+        if (~g).numeric_rep() in mapping.keys():
+            n = g.group.get_member_by_val(mapping[(~g).numeric_rep()])
+            permuted_members = [((m)*(~n)).numeric_rep() for m in s3.members()]
+        expected = perm_list_to_matrix(members, permuted_members)
+
+        assert np.allclose(g.matrix_form_subgroup(s3), expected), g