Refactor permutation code

.github/workflows/test.yml

@@ -28,6 +28,8 @@ jobs:
           git fetch
           git checkout indiamai/fuse-quads
           git pull
+          pip install pybind11 Cython
+          make
       - name: Install checkedout fuse
         run: |
           python3 -m pip install --break-system-packages -e '.[dev]'

fuse/cells.py

@@ -10,7 +10,7 @@
 from matplotlib.patches import FancyArrowPatch
 from mpl_toolkits.mplot3d import proj3d
 from sympy.combinatorics.named_groups import SymmetricGroup
-from fuse.utils import sympy_to_numpy, fold_reduce, numpy_to_str_tuple
+from fuse.utils import sympy_to_numpy, fold_reduce, numpy_to_str_tuple, orientation_value
 from FIAT.reference_element import Simplex, TensorProductCell as FiatTensorProductCell, Hypercube
 from ufl.cell import Cell, TensorProductCell
 
@@ -546,9 +546,8 @@ def permute_entities(self, g, d):
     def basis_vectors(self, return_coords=True, entity=None):
         if not entity:
             entity = self
-        entity_levels = [sorted(generation) for generation in nx.topological_generations(entity.G)]
         self_levels = [sorted(generation) for generation in nx.topological_generations(self.G)]
-        vertices = entity_levels[entity.graph_dim()]
+        vertices = entity.ordered_vertices()
         if self.dimension == 0:
             # return [[]
             raise ValueError("Dimension 0 entities cannot have Basis Vectors")
@@ -592,15 +591,13 @@ def plot(self, show=True, plain=False, ax=None, filename=None):
         xs = np.linspace(-1, 1, 20)
         if ax is None:
             ax = plt.gca()
-        print("plotting")
         if self.dimension == 1:
             # line plot in 1D case
             nodes = self.d_entities(0, get_class=False)
             points = []
             for node in nodes:
                 attach = self.attachment(top_level_node, node)
                 points.extend(attach())
-            print(np.array(points))
             plt.plot(np.array(points), np.zeros_like(points), color="black")
 
         for i in range(self.dimension - 1, -1, -1):
@@ -613,7 +610,6 @@ def plot(self, show=True, plain=False, ax=None, filename=None):
                     if len(plotted) < 2:
                         plotted = (plotted[0], 0)
                     vert_coords += [plotted]
-                    print(np.array(plotted))
                     if not plain:
                         plt.plot(plotted[0], plotted[1], 'bo')
                         plt.annotate(node, (plotted[0], plotted[1]))
@@ -856,13 +852,18 @@ def __init__(self, cell, name=None):
     def cellname(self):
         return self.name
 
-    def construct_subelement(self, dimension):
+    def construct_subelement(self, dimension, e_id=0, o=None):
         """Constructs the reference element of a cell
         specified by subelement dimension.
 
         :arg dimension: subentity dimension (integer)
+        :arg e_id: subentity id, default 0, (integer)
+        :arg o: orientation of subentity, default None (GroupMemberRep)
         """
-        return self.fe_cell.d_entities(dimension)[0].to_fiat()
+        if o:
+            return self.fe_cell.d_entities(dimension)[e_id].orient(o).to_fiat()
+        else:
+            return self.fe_cell.d_entities(dimension)[e_id].to_fiat()
 
     def get_facet_element(self):
         dimension = self.get_spatial_dimension()
@@ -1027,3 +1028,33 @@ def constructCellComplex(name):
         return TensorProductPoint(*components).to_ufl(name)
     else:
         raise TypeError("Cell complex construction undefined for {}".format(str(name)))
+
+
+def compare_topologies(base, new):
+    """Compute orientations of sub entities against a base topology
+
+       base topology is assumed to follow the FIAT numbering convention,
+       ie edges are ordered from lower to higher
+       """
+    if list(base.keys()) != list(new.keys()):
+        raise ValueError("Topologies of different sizes cannot be compared")
+    orientations = []
+    for dim in sorted(base.keys()):
+        for entity in sorted(base[dim].keys()):
+            assert entity in new[dim].keys()
+            base_array = list(base[dim][entity])
+            new_array = list(new[dim][entity])
+            if sorted(base_array) == sorted(new_array):
+                orientations += [orientation_value(base_array, new_array)]
+            else:
+                # numbering does not match - renumber
+                # base is treated as ordered list
+                # new is renumbered by sorting the values
+                base_numbering = {base_array[i]: i for i in range(len(base_array))}
+                sorted_new = sorted(new_array)
+                new_numbering = {sorted_new[i]: i for i in range(len(new_array))}
+                base_renumbered = [base_numbering[b] for b in base_array]
+                new_renumbered = [new_numbering[n] for n in new_array]
+                orientations += [orientation_value(base_renumbered, new_renumbered)]
+
+    return orientations

fuse/dof.py

@@ -13,6 +13,7 @@ class Pairing():
 
     def __init__(self):
         self.entity = None
+        self.orientation = None
 
     def _to_dict(self):
         o_dict = {"entity": self.entity}
@@ -36,11 +37,22 @@ def __call__(self, kernel, v, cell):
 
     def convert_to_fiat(self, ref_el, dof, interpolant_deg):
         pt = dof.eval(FuseFunction(lambda *x: x))
+        # pt1 = dof.tabulate([[1]])
         return PointEvaluation(ref_el, pt)
+        # return PointEvaluation(ref_el, tuple(pt1[0]))
 
     def add_entity(self, entity):
         res = DeltaPairing()
         res.entity = entity
+        if self.orientation:
+            res = res.permute(self.orientation)
+        return res
+
+    def permute(self, g):
+        res = DeltaPairing()
+        if self.entity:
+            res.entity = self.entity.orient(g)
+        res.orientation = g
         return res
 
     def __repr__(self):
@@ -56,58 +68,71 @@ def _from_dict(obj_dict):
 
 
 class L2Pairing(Pairing):
-    """ need to think about the abstraction level here -
-    are we wanting to define them as quadrature now? or defer this?
-    """
+
     def __init__(self):
         super(L2Pairing, self).__init__()
 
     def __call__(self, kernel, v, cell):
-        # print(self.entity)
+        # TODO get degree of v
         # if cell == self.entity:
+        #     ref_el = self.entity.to_fiat()
         #     # print("evaluating", kernel, v, "on", self.entity)
-        #     quadrature = create_quadrature(self.entity.to_fiat(), 5)
+        #     Q = create_quadrature(self.entity.to_fiat(), 5)
         #     # need quadrature here too - therefore need the information from the triple.
         # else:
         #     ref_el = cell.to_fiat()
-        #     print(cell)
         #     ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()]
-        #     entity = ref_el.construct_subelement(self.entity.dim())
+        #     entity_ref = ref_el.construct_subelement(self.entity.dim())
+        #     entity = ref_el.construct_subelement(self.entity.dim(), ent_id, self.orientation)
         #     Q_ref = create_quadrature(entity, 5)
-        #     quadrature = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref)
-        quadrature = create_quadrature(self.entity.to_fiat(), 5)
+        #     Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref, self.orientation)
+        Q = create_quadrature(self.entity.to_fiat(), 5)
 
         def kernel_dot(x):
             return np.dot(kernel(*x), v(*x))
 
-        return quadrature.integrate(kernel_dot)
+        return Q.integrate(kernel_dot)
 
     def tabulate(self):
         pass
 
     def add_entity(self, entity):
         res = L2Pairing()
         res.entity = entity
+        if self.orientation:
+            res = res.permute(self.orientation)
+        return res
+
+    def permute(self, g):
+        if g.perm.is_Identity:
+            return self
+
+        res = L2Pairing()
+        if self.entity:
+            res.entity = self.entity.orient(g)
+        res.orientation = g
         return res
 
     def convert_to_fiat(self, ref_el, dof, interpolant_degree):
         total_deg = interpolant_degree + dof.kernel.degree()
         ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()]
-        entity = ref_el.construct_subelement(self.entity.dim())
+        # entity_ref = ref_el.construct_subelement(self.entity.dim())
+        entity = ref_el.construct_subelement(self.entity.dim(), ent_id, self.orientation)
         Q_ref = create_quadrature(entity, total_deg)
-        Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref)
+        # pts_ref, wts_ref = Q_ref.get_points(), Q_ref.get_weights()
+
+        # pts, wts, J = map_quadrature(pts_ref, wts_ref, Q_ref.ref_el, entity_ref, jacobian=True)
+        Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref, self.orientation)
         Jdet = Q.jacobian_determinant()
+        # Jdet = pseudo_determinant(J)
         qpts, _ = Q.get_points(), Q.get_weights()
-        print(qpts)
-        print(dof.tabulate(qpts))
+
         f_at_qpts = dof.tabulate(qpts).T / Jdet
-        print(len(Q.pts))
-        print(f_at_qpts.shape)
         functional = FrobeniusIntegralMoment(ref_el, Q, f_at_qpts)
         return functional
 
     def __repr__(self):
-        return "integral_{}({{kernel}} * {{fn}}) dx)".format(str(self.entity))
+        return "integral_{}({{kernel}} * {{fn}}) dx) ".format(str(self.entity))
 
     def dict_id(self):
         return "L2Inner"
@@ -154,7 +179,9 @@ def permute(self, g):
     def __call__(self, *args):
         return self.pt
 
-    def tabulate(self, Qpts):
+    def tabulate(self, Qpts, attachment=None):
+        if attachment:
+            return np.array([attachment(*self.pt) for _ in Qpts]).astype(np.float64)
         return np.array([self.pt for _ in Qpts]).astype(np.float64)
 
     def _to_dict(self):
@@ -186,16 +213,20 @@ def degree(self):
         return self.fn.as_poly().total_degree()
 
     def permute(self, g):
-        new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))})
-        return PolynomialKernel(new_fn, symbols=self.syms)
+        return self
+        # new_fn = self.fn.subs({self.syms[i]: g(self.syms)[i] for i in range(len(self.syms))})
+        # return PolynomialKernel(new_fn, symbols=self.syms)
 
     def __call__(self, *args):
         res = sympy_to_numpy(self.fn, self.syms, args[:len(self.syms)])
         if not hasattr(res, '__iter__'):
             return [res]
         return res
 
-    def tabulate(self, Qpts):
+    def tabulate(self, Qpts, attachment=None):
+        # TODO do we need to attach qpts
+        # if attachment:
+        #     return np.array([self(*attachment(*pt)) for pt in Qpts]).astype(np.float64)
         return np.array([self(*pt) for pt in Qpts]).astype(np.float64)
 
     def _to_dict(self):
@@ -232,7 +263,7 @@ def __init__(self, pairing, kernel, entity=None, attachment=None, target_space=N
 
     def __call__(self, g):
         new_generation = self.generation.copy()
-        return DOF(self.pairing, self.kernel.permute(g), self.trace_entity, self.attachment, self.target_space, g, self.immersed, new_generation, self.sub_id, self.cell)
+        return DOF(self.pairing.permute(g), self.kernel.permute(g), self.trace_entity, self.attachment, self.target_space, g, self.immersed, new_generation, self.sub_id, self.cell)
 
     def eval(self, fn, pullback=True):
         return self.pairing(self.kernel, fn, self.cell)
@@ -256,7 +287,6 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non
 
     def convert_to_fiat(self, ref_el, interpolant_degree):
         return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree)
-        raise NotImplementedError("Fiat conversion only implemented for Point eval")
 
     def __repr__(self, fn="v"):
         return str(self.pairing).format(fn=fn, kernel=self.kernel)
@@ -294,16 +324,16 @@ def eval(self, fn, pullback=True):
 
     def tabulate(self, Qpts):
         immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g)
-        res = self.kernel.tabulate(Qpts)
+        res = self.kernel.tabulate(Qpts, self.attachment)
         return immersion*res
 
     def __call__(self, g):
-        permuted = self.cell.permute_entities(g, self.trace_entity.dim())
         index_trace = self.cell.d_entities_ids(self.trace_entity.dim()).index(self.trace_entity.id)
-        new_trace_entity = self.cell.get_node(permuted[index_trace][0]).orient(permuted[index_trace][1])
-
-        return ImmersedDOF(self.pairing, self.kernel.permute(permuted[index_trace][1]), new_trace_entity,
-                           self.attachment, self.target_space, g, self.triple, self.generation, self.sub_id, self.cell)
+        permuted_e, permuted_g = self.cell.permute_entities(g, self.trace_entity.dim())[index_trace]
+        new_trace_entity = self.cell.get_node(permuted_e).orient(permuted_g)
+        new_attach = lambda *x: g(self.attachment(*x))
+        return ImmersedDOF(self.pairing.permute(permuted_g), self.kernel.permute(permuted_g), new_trace_entity,
+                           new_attach, self.target_space, g, self.triple, self.generation, self.sub_id, self.cell)
 
     def __repr__(self):
         fn = "tr_{1}_{0}(v)".format(str(self.trace_entity), str(self.target_space))

fuse/groups.py

@@ -1,10 +1,10 @@
 import fuse.cells as cells
+from fuse.utils import orientation_value
 from sympy.combinatorics import PermutationGroup, Permutation
 from sympy.combinatorics.named_groups import SymmetricGroup, DihedralGroup, CyclicGroup, AlternatingGroup
 from sympy.matrices.expressions import PermutationMatrix
 import numpy as np
 import sympy as sp
-import math
 
 
 def perm_matrix_to_perm_array(p_mat):
@@ -57,12 +57,7 @@ def compute_perm(self, base_val=None):
     def numeric_rep(self):
         identity = self.group.identity.perm.array_form
         m_array = self.perm.array_form
-        val = 0
-        for i in range(len(identity)):
-            loc = m_array.index(identity[i])
-            m_array.remove(identity[i])
-            val += loc * math.factorial(len(identity) - i - 1)
-        return val
+        return orientation_value(identity, m_array)
 
     def __eq__(self, x):
         assert isinstance(x, GroupMemberRep)