Orientations that depends on basis components

Makefile

@@ -26,6 +26,12 @@ tests:
 	@echo "    Running all tests"
 	@FIREDRAKE_USE_FUSE=1 python3 -m coverage run -p -m pytest -rx test
 
+mini_tests:
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_2d_examples_docs.py
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_1d
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_orientations.py::test_surface_vec_rt
+	@FIREDRAKE_USE_FUSE=1 python3 -m pytest test/test_convert_to_fiat.py::test_projection_convergence_3d\[construct_tet_ned-N1curl-1-0.8\]
+
 coverage:
 	@python3 -m coverage combine
 	@python3 -m coverage report -m 

fuse/__init__.py

@@ -1,6 +1,7 @@
-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, diff_C3, tet_edges, tet_faces, sq_edges, GroupRepresentation, PermutationSetRepresentation, get_cyc_group, get_sym_group
-from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, PolynomialKernel, ComponentKernel
+from fuse.dof import DeltaPairing, DOF, L2Pairing, FuseFunction, PointKernel, VectorKernel, BarycentricPolynomialKernel, PolynomialKernel, ComponentKernel
 from fuse.triples import ElementTriple, DOFGenerator, immerse
 from fuse.traces import TrH1, TrGrad, TrHess, TrHCurl, TrHDiv
 from fuse.tensor_products import tensor_product

fuse/cells.py

@@ -151,6 +151,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
@@ -315,8 +322,8 @@ def __init__(self, d, edges=[], vertex_num=None, oriented=False, group=None, edg
         self.group = group
         if not group:
             self.group = self.compute_cell_group()
-
         self.group = self.group.add_cell(self)
+        self.basis_group = self.compute_basis_group()
 
     def compute_attachments(self, n, points, orientations=None):
         """
@@ -373,6 +380,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()]
@@ -388,6 +396,19 @@ def compute_cell_group(self):
                         break
         return fuse_groups.PermutationSetRepresentation(list(accepted_perms))
 
+    def compute_basis_group(self):
+        if self.dimension == 0:
+            return fuse_groups.S1.add_cell(self)
+        bvs = self.basis_vectors()
+        bv_0 = bvs[0]
+        basis_group_members = []
+        for bv in bvs:
+            for g in self.group.members():
+                if np.allclose(np.array(g(bv_0)), np.array(bv)):
+                    basis_group_members += [g.perm]
+                    break
+        return fuse_groups.PermutationSetRepresentation(list(basis_group_members)).add_cell(self)
+
     def get_spatial_dimension(self):
         return self.dimension
 
@@ -475,12 +496,8 @@ def _subentity_traversal(self, sub_ents, min_ids):
                     sub_ents = self.get_node(p)._subentity_traversal(sub_ents, min_ids)
         if dim > 1:
             connections = [(c.point.id, c.point.group.identity) for c in self.connections]
-            # if self.oriented:
-            #     connections = self.permute_entities(self.oriented, dim - 1)
             # if self.dimension == 2:
             #     connections = [connections[-1]] + connections[:-1]
-            #     print([self.get_node(c[0]).id - min_ids[1] for c in connections])
-            #     print([c.point.id - min_ids[1] for c in self.connections])
             for e, o in connections:
                 p = self.get_node(e).orient(o)
                 p_dim = p.get_spatial_dimension()
@@ -502,6 +519,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
@@ -544,6 +570,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)
 
@@ -639,13 +668,12 @@ def permute_entities(self, g, d):
 
         return reordered_entities
 
-    def basis_vectors(self, return_coords=True, entity=None, order=False):
+    def basis_vectors(self, return_coords=True, entity=None, order=False, norm=True):
         if not entity:
             entity = self
         self_levels = [generation for generation in nx.topological_generations(self.G)]
         vertices = entity.ordered_vertices()
         if self.dimension == 0:
-            # return [[]
             raise ValueError("Dimension 0 entities cannot have Basis Vectors")
         if self.oriented:
             # ordered_vertices() handles the orientation so we want to drop the orientation node
@@ -659,7 +687,10 @@ def basis_vectors(self, return_coords=True, entity=None, order=False):
         for v in vertices[1:]:
             if return_coords:
                 v_coords = self.attachment(top_level_node, v)()
-                sub = normalise(np.subtract(v_coords, v_0_coords))
+                if norm:
+                    sub = normalise(np.subtract(v_coords, v_0_coords))
+                else:
+                    sub = np.subtract(v_coords, v_0_coords)
                 if not hasattr(sub, "__iter__"):
                     basis_vecs.append((sub,))
                 else:
@@ -816,12 +847,35 @@ def attachment(self, source, dst):
 
         return lambda *x: fold_reduce(attachments[0], *x)
 
+    def attachment_J_det(self, source, dst):
+        attachment = self.attachment(source, dst)
+        symbol_names = ["x", "y", "z"]
+        symbols = []
+        if self.dim_of_node(dst) == 0:
+            return 1
+        for i in range(self.dim_of_node(dst)):
+            symbols += [sp.Symbol(symbol_names[i])]
+        J = sp.Matrix(attachment(*symbols)).jacobian(sp.Matrix(symbols))
+        return np.sqrt(abs(float(sp.det(J.T * J))))
+
     def quadrature(self, degree):
         fiat_el = self.to_fiat()
         Q = create_quadrature(fiat_el, degree)
         pts, wts = Q.get_points(), Q.get_weights()
         return pts, wts
 
+    def cartesian_to_barycentric(self, pts):
+        verts = np.array(self.ordered_vertex_coords())
+        v_0 = self.ordered_vertex_coords()[0]
+        bvs = np.array(self.basis_vectors(norm=False))
+        bary_coords = []
+        for pt in pts:
+            res = np.matmul(np.linalg.inv(bvs.T), np.array(pt - v_0))
+            assert np.allclose(sum(bvs[i]*res[i] for i in range(len(bvs))), np.array(pt - v_0))
+            bary_coords += [(1 - sum(res),) + tuple(res[i] for i in range(len(res)))]
+            assert np.allclose(np.array(sum(bary_coords[-1][i]*verts[i] for i in range(len(verts)))), pt)
+        return bary_coords
+
     def cell_attachment(self, dst):
         if not isinstance(dst, int):
             raise ValueError
@@ -874,6 +928,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():
     """
@@ -897,7 +958,12 @@ def __call__(self, *x):
             if hasattr(self.attachment, '__iter__'):
                 res = []
                 for attach_comp in self.attachment:
-                    res.append(sympy_to_numpy(attach_comp, syms, x))
+                    if len(attach_comp.atoms(sp.Symbol)) <= len(x):
+                        res.append(sympy_to_numpy(attach_comp, syms, x))
+                    else:
+                        res_val = attach_comp.subs({syms[i]: x[i] for i in range(len(x))})
+                        res.append(res_val)
+
                 return tuple(res)
             return sympy_to_numpy(self.attachment, syms, x)
         return x
@@ -929,11 +995,11 @@ 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 ordered_vertices(self):
         return self.A.ordered_vertices() + self.B.ordered_vertices()
@@ -945,8 +1011,8 @@ def get_sub_entities(self):
         self.A.get_sub_entities()
         self.B.get_sub_entities()
 
-    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)
@@ -961,17 +1027,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):
@@ -1161,9 +1301,8 @@ def constructCellComplex(name):
         return polygon(3).to_ufl(name)
         # return ufc_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 ufc_quad().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":
         # return ufc_tetrahedron().to_ufl(name)