First Class Quads

fuse/__init__.py

@@ -1,5 +1,5 @@
 
-from fuse.cells import Point, Edge, polygon, make_tetrahedron, constructCellComplex
+from fuse.cells import Point, Edge, polygon, line, make_tetrahedron, constructCellComplex, TensorProductPoint
 from fuse.groups import S1, S2, S3, D4, Z3, Z4, C3, C4, S4, A4, tri_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group
 from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, PolynomialKernel
 from fuse.triples import ElementTriple, DOFGenerator, immerse

fuse/cells.py

@@ -148,6 +148,13 @@ def compute_scaled_verts(d, n):
         raise ValueError("Dimension {} not supported".format(d))
 
 
+def line():
+    """
+    Constructs the default 1D interval
+    """
+    return Point(1, [Point(0), Point(0)], vertex_num=2)
+
+
 def polygon(n):
     """
     Constructs the 2D default cell with n sides/vertices
@@ -336,6 +343,7 @@ def compute_cell_group(self):
         """
         verts = self.ordered_vertices()
         v_coords = [self.get_node(v, return_coords=True) for v in verts]
+
         n = len(verts)
         max_group = SymmetricGroup(n)
         edges = [edge.ordered_vertices() for edge in self.edges()]
@@ -407,6 +415,15 @@ def get_starter_ids(self):
         min_ids = [min(dimension) for dimension in structure]
         return min_ids
 
+    def local_id(self, node):
+        structure = [sorted(generation) for generation in nx.topological_generations(self.G)]
+        structure.reverse()
+        min_id = self.get_starter_ids()
+        for d in range(len(structure)):
+            if node.id in structure[d]:
+                return node.id - min_id[d]
+        raise ValueError("Node not found in cell")
+
     def graph_dim(self):
         if self.oriented:
             dim = self.dimension + 1
@@ -449,6 +466,9 @@ def ordered_vertices(self, get_class=False):
                 return self.oriented.permute(verts)
             return verts
 
+    def ordered_vertex_coords(self):
+        return [self.get_node(o, return_coords=True) for o in self.ordered_vertices()]
+
     def d_entities_ids(self, d):
         return self.d_entities(d, get_class=False)
 
@@ -471,7 +491,7 @@ def get_node(self, node, return_coords=False):
         if return_coords:
             top_level_node = self.d_entities_ids(self.graph_dim())[0]
             if self.dimension == 0:
-                return [()]
+                return ()
             return self.attachment(top_level_node, node)()
         return self.G.nodes.data("point_class")[node]
 
@@ -550,7 +570,6 @@ def basis_vectors(self, return_coords=True, entity=None):
         self_levels = [sorted(generation) for generation in nx.topological_generations(self.G)]
         vertices = entity_levels[entity.graph_dim()]
         if self.dimension == 0:
-            # return [[]
             raise ValueError("Dimension 0 entities cannot have Basis Vectors")
         top_level_node = self_levels[0][0]
         v_0 = vertices[0]
@@ -715,8 +734,8 @@ def copy(self):
     def to_fiat(self, name=None):
         if len(self.vertices()) == self.dimension + 1:
             return CellComplexToFiatSimplex(self, name)
-        if len(self.vertices()) == 2 ** self.dimension:
-            return CellComplexToFiatHypercube(self, name)
+        # if len(self.vertices()) == 2 ** self.dimension:
+        #     return CellComplexToFiatHypercube(self, name)
         raise NotImplementedError("Custom shape elements/ First class quads are not yet supported")
 
     def to_ufl(self, name=None):
@@ -736,6 +755,13 @@ def dict_id(self):
     def _from_dict(o_dict):
         return Point(o_dict["dim"], o_dict["edges"], oriented=o_dict["oriented"], cell_id=o_dict["id"])
 
+    def equivalent(self, other):
+        if self.dimension != other.dimension:
+            return False
+        if set(self.ordered_vertex_coords()) != set(other.ordered_vertex_coords()):
+            return False
+        return self.get_topology() == other.get_topology()
+
 
 class Edge():
     """
@@ -759,7 +785,7 @@ def __call__(self, *x):
             if hasattr(self.attachment, '__iter__'):
                 res = []
                 for attach_comp in self.attachment:
-                    if len(attach_comp.atoms(sp.Symbol)) == len(x):
+                    if len(attach_comp.atoms(sp.Symbol)) <= len(x):
                         res.append(sympy_to_numpy(attach_comp, syms, x))
                     else:
                         res.append(attach_comp.subs({syms[i]: x[i] for i in range(len(x))}))
@@ -794,17 +820,17 @@ def _from_dict(o_dict):
 
 class TensorProductPoint():
 
-    def __init__(self, A, B, flat=False):
+    def __init__(self, A, B):
         self.A = A
         self.B = B
         self.dimension = self.A.dimension + self.B.dimension
-        self.flat = flat
+        self.flat = False
 
     def get_spatial_dimension(self):
         return self.dimension
 
-    def dimension(self):
-        return tuple(self.A.dimension, self.B.dimension)
+    def dim(self):
+        return (self.A.dimension, self.B.dimension)
 
     def d_entities(self, d, get_class=True):
         return self.A.d_entities(d, get_class) + self.B.d_entities(d, get_class)
@@ -819,17 +845,91 @@ def vertices(self, get_class=True, return_coords=False):
         return verts
 
     def to_ufl(self, name=None):
-        if self.flat:
-            return CellComplexToUFL(self, "quadrilateral")
         return TensorProductCell(self.A.to_ufl(), self.B.to_ufl())
 
     def to_fiat(self, name=None):
-        if self.flat:
-            return CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name))
         return CellComplexToFiatTensorProduct(self, name)
 
     def flatten(self):
-        return TensorProductPoint(self.A, self.B, True)
+        assert self.A.equivalent(self.B)
+        return FlattenedPoint(self.A, self.B)
+
+
+class FlattenedPoint(Point, TensorProductPoint):
+
+    def __init__(self, A, B):
+        self.A = A
+        self.B = B
+        self.dimension = self.A.dimension + self.B.dimension
+        self.flat = True
+        fuse_edges = self.construct_fuse_rep()
+        super().__init__(self.dimension, fuse_edges)
+
+    def to_ufl(self, name=None):
+        return CellComplexToUFL(self, "quadrilateral")
+
+    def to_fiat(self, name=None):
+        # TODO this should check if it actually is a hypercube
+        fiat = CellComplexToFiatHypercube(self, CellComplexToFiatTensorProduct(self, name))
+        return fiat
+
+    def construct_fuse_rep(self):
+        sub_cells = [self.A, self.B]
+        dims = (self.A.dimension, self.B.dimension)
+
+        points = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells}
+        attachments = {cell: {i: [] for i in range(max(dims) + 1)} for cell in sub_cells}
+
+        for d in range(max(dims) + 1):
+            for cell in sub_cells:
+                if d <= cell.dimension:
+                    sub_ent = cell.d_entities(d, get_class=True)
+                    points[cell][d].extend(sub_ent)
+                    for s in sub_ent:
+                        attachments[cell][d].extend(s.connections)
+
+        prod_points = list(itertools.product(*[points[cell][0] for cell in sub_cells]))
+        # temp = prod_points[1]
+        # prod_points[1] = prod_points[2]
+        # prod_points[2] = temp
+        point_cls = [Point(0) for i in range(len(prod_points))]
+        edges = []
+
+        # generate edges of tensor product result
+        for a in prod_points:
+            for b in prod_points:
+                # of all combinations of point, take those where at least one changes and at least one is the same
+                if any(a[i] == b[i] for i in range(len(a))) and any(a[i] != b[i] for i in range(len(sub_cells))):
+                    # ensure if they change, that edge exists in the existing topology
+                    if all([a[i] == b[i] or (sub_cells[i].local_id(a[i]), sub_cells[i].local_id(b[i])) in list(sub_cells[i].topology[1].values()) for i in range(len(sub_cells))]):
+                        edges.append((a, b))
+        # hasse level 1
+        edge_cls1 = {e: None for e in edges}
+        for i in range(len(sub_cells)):
+            for (a, b) in edges:
+                a_idx = prod_points.index(a)
+                b_idx = prod_points.index(b)
+                if a[i] != b[i]:
+                    a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0]
+                    b_edge = [att for att in attachments[sub_cells[i]][1] if att.point == b[i]][0]
+                    edge_cls1[(a, b)] = Point(1, [Edge(point_cls[a_idx], a_edge.attachment, a_edge.o),
+                                                  Edge(point_cls[b_idx], b_edge.attachment, b_edge.o)])
+        edge_cls2 = []
+        # hasse level 2
+        for i in range(len(sub_cells)):
+            for (a, b) in edges:
+                if a[i] == b[i]:
+                    x = sp.Symbol("x")
+                    a_edge = [att for att in attachments[sub_cells[i]][1] if att.point == a[i]][0]
+                    if i == 0:
+                        attach = (x,) + a_edge.attachment
+                    else:
+                        attach = a_edge.attachment + (x,)
+                    edge_cls2.append(Edge(edge_cls1[(a, b)], attach, a_edge.o))
+        return edge_cls2
+
+    def flatten(self):
+        return self
 
 
 class CellComplexToFiatSimplex(Simplex):
@@ -1011,8 +1111,7 @@ def constructCellComplex(name):
         return polygon(3).to_ufl(name)
         # return firedrake_triangle().to_ufl(name)
     elif name == "quadrilateral":
-        interval = Point(1, [Point(0), Point(0)], vertex_num=2)
-        return TensorProductPoint(interval, interval).flatten().to_ufl(name)
+        return TensorProductPoint(line(), line()).flatten().to_ufl(name)
         # return firedrake_quad().to_ufl(name)
         # return polygon(4).to_ufl(name)
     elif name == "tetrahedron":

fuse/groups.py

@@ -55,13 +55,18 @@ def compute_perm(self, base_val=None):
         return val, val_list
 
     def numeric_rep(self):
+        """ Uses a standard formula to number permutations in the group.
+        For the case where this doesn't automatically number from 0..n (ie the group is not the full symmetry group),
+        a mapping is constructed on group creation"""
         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)
+        if self.group.group_rep_numbering is not None:
+            return self.group.group_rep_numbering[val]
         return val
 
     def __eq__(self, x):
@@ -134,6 +139,11 @@ def __init__(self, perm_list, cell=None):
                 self._members.append(p_rep)
                 counter += 1
 
+            self.group_rep_numbering = None
+            numeric_reps = [m.numeric_rep() for m in self.members()]
+            if sorted(numeric_reps) != list(range(len(numeric_reps))):
+                self.group_rep_numbering = {a: b for a, b in zip(sorted(numeric_reps), list(range(len(numeric_reps))))}
+
     def add_cell(self, cell):
         return PermutationSetRepresentation(self.perm_list, cell=cell)
 
@@ -224,6 +234,11 @@ def __init__(self, base_group, cell=None):
                     self.identity = p_rep
                 counter += 1
 
+            self.group_rep_numbering = None
+            numeric_reps = [m.numeric_rep() for m in self.members()]
+            if sorted(numeric_reps) != list(range(len(numeric_reps))):
+                self.group_rep_numbering = {a: b for a, b in zip(sorted(numeric_reps), list(range(len(numeric_reps))))}
+
             # this order produces simpler generator lists
             # self.generators.reverse()
 

fuse/spaces/polynomial_spaces.py

@@ -1,13 +1,16 @@
 from FIAT.polynomial_set import ONPolynomialSet
+from FIAT.expansions import morton_index2, morton_index3
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.reference_element import cell_to_simplex
 from FIAT import expansions, polynomial_set, reference_element
 from itertools import chain
-from fuse.utils import tabulate_sympy, max_deg_sp_mat
+from fuse.utils import tabulate_sympy, max_deg_sp_expr
 import sympy as sp
 import numpy as np
 from functools import total_ordering
 
+morton_index = {2: morton_index2, 3: morton_index3}
+
 
 @total_ordering
 class PolynomialSpace(object):
@@ -47,7 +50,6 @@ def degree(self):
         return self.maxdegree
 
     def to_ON_polynomial_set(self, ref_el, k=None):
-        # how does super/sub degrees work here
         if not isinstance(ref_el, reference_element.Cell):
             ref_el = ref_el.to_fiat()
         sd = ref_el.get_spatial_dimension()
@@ -56,18 +58,25 @@ def to_ON_polynomial_set(self, ref_el, k=None):
             shape = (sd,)
         else:
             shape = tuple()
+        base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal")
+        indices = None
 
         if self.mindegree > 0:
-            base_ON = ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal")
             dimPmin = expansions.polynomial_dimension(ref_el, self.mindegree)
             dimPmax = expansions.polynomial_dimension(ref_el, self.maxdegree)
             if self.set_shape:
                 indices = list(chain(*(range(i * dimPmin, i * dimPmax) for i in range(sd))))
             else:
                 indices = list(range(dimPmin, dimPmax))
-            restricted_ON = base_ON.take(indices)
-            return restricted_ON
-        return ONPolynomialSet(ref_el, self.maxdegree, shape, scale="orthonormal")
+
+        if self.contains != self.maxdegree and self.contains != -1:
+            indices = [morton_index[sd](p, q) for p in range(self.contains + 1) for q in range(self.contains + 1)]
+
+        if indices is None:
+            return base_ON
+
+        restricted_ON = base_ON.take(indices)
+        return restricted_ON
 
     def __repr__(self):
         res = ""
@@ -161,37 +170,49 @@ def to_ON_polynomial_set(self, ref_el):
         if not isinstance(ref_el, reference_element.Cell):
             ref_el = ref_el.to_fiat()
         k = max([s.maxdegree for s in self.spaces])
-        space_poly_sets = [s.to_ON_polynomial_set(ref_el) for s in self.spaces]
         sd = ref_el.get_spatial_dimension()
         ref_el = cell_to_simplex(ref_el)
 
-        if all([w == 1 for w in self.weights]):
-            weighted_sets = space_poly_sets
-
         # otherwise have to work on this through tabulation
 
         Q = create_quadrature(ref_el, 2 * (k + 1))
         Qpts, Qwts = Q.get_points(), Q.get_weights()
         weighted_sets = []
 
-        for (space, w) in zip(space_poly_sets, self.weights):
+        for (s, w) in zip(self.spaces, self.weights):
+            space = s.to_ON_polynomial_set(ref_el)
+            if s.set_shape:
+                shape = (sd,)
+            else:
+                shape = tuple()
             if not (isinstance(w, sp.Expr) or isinstance(w, sp.Matrix)):
                 weighted_sets.append(space)
             else:
-                w_deg = max_deg_sp_mat(w)
-                Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, scale="orthonormal")
-                vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal")
+                if isinstance(w, sp.Expr):
+                    w = sp.Matrix([[w]])
+                    vec = False
+                else:
+                    vec = True
+                w_deg = max_deg_sp_expr(w)
+                Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, shape, scale="orthonormal")
+                # vec_Pkpw = ONPolynomialSet(ref_el, space.degree + w_deg, (sd,), scale="orthonormal")
 
                 space_at_Qpts = space.tabulate(Qpts)[(0,) * sd]
                 Pkpw_at_Qpts = Pkpw.tabulate(Qpts)[(0,) * sd]
 
                 tabulated_expr = tabulate_sympy(w, Qpts).T
-                scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :]
+                if s.set_shape or vec:
+                    scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :]
+                else:
+                    # breakpoint()
+                    scaled_at_Qpts = space_at_Qpts[:, None, :] * tabulated_expr[None, :, :]
+                    scaled_at_Qpts = scaled_at_Qpts.squeeze()
                 PkHw_coeffs = np.dot(np.multiply(scaled_at_Qpts, Qwts), Pkpw_at_Qpts.T)
+                # breakpoint()
                 weighted_sets.append(polynomial_set.PolynomialSet(ref_el,
                                                                   space.degree + w_deg,
                                                                   space.degree + w_deg,
-                                                                  vec_Pkpw.get_expansion_set(),
+                                                                  Pkpw.get_expansion_set(),
                                                                   PkHw_coeffs))
         combined_sets = weighted_sets[0]
         for i in range(1, len(weighted_sets)):