convert derivative nodes to fiat def, lint

Author: indiamai

fuse/cells.py

@@ -164,17 +164,15 @@ def polygon(n):
 
     return Point(2, edges, vertex_num=n)
 
+
 def firedrake_triangle():
     vertices = []
     for i in range(3):
         vertices.append(Point(0))
     edges = []
-    edges.append(
-            Point(1, [vertices[1], vertices[2]], vertex_num=2))
-    edges.append(
-            Point(1, [vertices[0], vertices[2]], vertex_num=2))
-    edges.append(
-            Point(1, [vertices[0], vertices[1]], vertex_num=2))
+    edges.append(Point(1, [vertices[1], vertices[2]], vertex_num=2))
+    edges.append(Point(1, [vertices[0], vertices[2]], vertex_num=2))
+    edges.append(Point(1, [vertices[0], vertices[1]], vertex_num=2))
     tri = Point(2, edges, vertex_num=3, edge_orientations={1: [1, 0]})
     # tri = polygon(3)
     s3 = tri.group
@@ -360,7 +358,6 @@ def get_topology(self):
             self.topology[i] = {}
             self.topology_unrelabelled[i] = {}
             for node in dimension:
-                neighbours = list(self.G.neighbors(node))
                 self.topology[i][node - min_ids[i]] = tuple([relabelled_verts[vert] for vert in self.get_node(node).ordered_vertices()])
                 self.topology_unrelabelled[i][node - min_ids[i]] = tuple([vert - min_ids[0] for vert in self.get_node(node).ordered_vertices()])
         return self.topology_unrelabelled
@@ -728,11 +725,9 @@ class CellComplexToFiatSimplex(Simplex):
     def __init__(self, cell, name=None):
         self.fe_cell = cell
         if name is not None:
-            name = "IndiaDefCell"
+            name = "FuseCell"
         self.name = name
 
-        
-
         # verts = [cell.get_node(v, return_coords=True) for v in cell.ordered_vertices()]
         verts = cell.vertices(return_coords=True)
         topology = cell.get_topology()
@@ -838,7 +833,7 @@ def to_fiat(self):
         return self.cell_complex.to_fiat(name=self.cellname())
 
     def __repr__(self):
-        return super(CellComplexToUFL, self).__repr__() 
+        return super(CellComplexToUFL, self).__repr__()
 
     def reconstruct(self, **kwargs):
         """Reconstruct this cell, overwriting properties by those in kwargs."""
@@ -866,4 +861,3 @@ def constructCellComplex(name):
         return make_tetrahedron().to_ufl(name)
     else:
         raise TypeError("Cell complex construction undefined for {}".format(str(name)))
-

fuse/dof.py

@@ -37,7 +37,7 @@ def __call__(self, kernel, v, cell):
     def convert_to_fiat(self, ref_el, dof, interpolant_deg):
         pt = dof.eval(MyTestFunction(lambda *x: x))
         return PointEvaluation(ref_el, pt)
-    
+
     def get_pts(self, ref_el, total_degree):
         entity = ref_el.construct_subelement(self.entity.dim())
         return [(0,) * entity.get_spatial_dimension()], [1], 1
@@ -272,24 +272,20 @@ def add_context(self, dof_gen, cell, space, g, overall_id=None, generator_id=Non
             self.sub_id = generator_id
 
     def convert_to_fiat(self, ref_el, interpolant_degree):
-        return self.pairing.convert_to_fiat(ref_el, self, interpolant_degree)
-    
-    def convert_to_fiat_new(self, ref_el, interpolant_degree):
         total_degree = self.kernel.degree() + interpolant_degree
         pts, wts, jdet = self.pairing.get_pts(ref_el, total_degree)
         f_pts = self.kernel.tabulate(pts).T / jdet
         # TODO need to work out how i can discover the shape in a better way
         if isinstance(self.pairing, DeltaPairing):
             shp = tuple()
-            pt_dict = {tuple(p) : [(w, tuple())] for (p, w) in zip(f_pts.T, wts)}
+            pt_dict = {tuple(p): [(w, tuple())] for (p, w) in zip(f_pts.T, wts)}
         else:
             shp = tuple(f_pts.shape[:-1])
             weights = np.transpose(np.multiply(f_pts, wts), (-1,) + tuple(range(len(shp))))
             alphas = list(np.ndindex(shp))
             pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(pts, weights)}
 
-        return Functional(ref_el, shp, pt_dict, {}, self.__repr__())
-
+        return [Functional(ref_el, shp, pt_dict, {}, self.__repr__())]
 
     def __repr__(self, fn="v"):
         return str(self.pairing).format(fn=fn, kernel=self.kernel)
@@ -329,24 +325,24 @@ def tabulate(self, Qpts):
         immersion = self.target_space.tabulate(Qpts, self.trace_entity, self.g)
         res = self.kernel.tabulate(Qpts)
         return immersion*res
-    
-    def convert_to_fiat_new(self, ref_el, interpolant_degree):
+
+    def convert_to_fiat(self, ref_el, interpolant_degree):
         total_degree = self.kernel.degree() + interpolant_degree
         pts, wts, jdet = self.pairing.get_pts(ref_el, total_degree)
         f_pts = self.kernel.tabulate(pts, self.attachment)
         attached_pts = [self.attachment(*p) for p in pts]
         immersion = self.target_space.tabulate(f_pts, self.trace_entity, self.g)
-        
+
         f_pts = (f_pts * immersion).T / jdet
-        pt_dict, deriv_dict = self.target_space.convert_to_fiat(attached_pts, f_pts, wts)
+        dict_list = self.target_space.convert_to_fiat(attached_pts, f_pts, wts)
 
         # breakpoint()
         # TODO need to work out how i can discover the shape in a better way
         if isinstance(self.pairing, DeltaPairing):
             shp = tuple()
         else:
             shp = tuple(f_pts.shape[:-1])
-        return Functional(ref_el, shp, pt_dict, deriv_dict, self.__repr__())
+        return [Functional(ref_el, shp, pt_dict, deriv_dict, self.__repr__()) for (pt_dict, deriv_dict) in dict_list]
 
     def __call__(self, g):
         permuted = self.cell.permute_entities(g, self.trace_entity.dim())

fuse/traces.py

@@ -55,8 +55,8 @@ def tabulate(self, Qpts, trace_entity, g):
         return np.ones_like(Qpts)
 
     def convert_to_fiat(self, qpts, pts, wts):
-        pt_dict = {tuple(p) : [(w, tuple())] for (p, w) in zip(pts.T, wts)}
-        return pt_dict, {}
+        pt_dict = {tuple(p): [(w, tuple())] for (p, w) in zip(pts.T, wts)}
+        return [(pt_dict, {})]
 
     def __repr__(self):
         return "H1"
@@ -96,8 +96,7 @@ def convert_to_fiat(self, qpts, pts, wts):
         weights = np.transpose(np.multiply(f_at_qpts, wts), (-1,) + tuple(range(len(shp))))
         alphas = list(np.ndindex(shp))
         pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)}
-        return pt_dict, {}
-
+        return [(pt_dict, {})]
 
     def tabulate(self, Qpts, trace_entity, g):
         entityBasis = np.array(trace_entity.basis_vectors())
@@ -135,14 +134,14 @@ def tabulate(self, Qpts, trace_entity, g):
         subEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity))
         result = np.matmul(tangent, subEntityBasis)
         return result
-    
+
     def convert_to_fiat(self, qpts, pts, wts):
         f_at_qpts = pts
         shp = tuple(f_at_qpts.shape[:-1])
         weights = np.transpose(np.multiply(f_at_qpts, wts), (-1,) + tuple(range(len(shp))))
         alphas = list(np.ndindex(shp))
         pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)}
-        return pt_dict, {}
+        return [(pt_dict, {})]
 
     def plot(self, ax, coord, trace_entity, g, **kwargs):
         permuted = self.domain.permute_entities(g, trace_entity.dimension)
@@ -179,10 +178,28 @@ def apply(*x):
             return tuple(result)
         return apply
 
+    def convert_to_fiat(self, qpts, pts, wts):
+        shp = tuple(pts.shape[0:])
+        alphas = []
+        for i in range(pts.shape[0]):
+            new = np.zeros_like(shp)
+            new[i] = 1
+            alphas += [tuple(new)]
+        deriv_dicts = []
+        for alpha in alphas:
+            deriv_dicts += [{tuple(p): [(1.0, tuple(alpha), tuple())] for p in pts.T}]
+
+        # self.alpha = tuple(alpha)
+        # self.order = sum(self.alpha)
+        return [({}, d) for d in deriv_dicts]
+
     def plot(self, ax, coord, trace_entity, g, **kwargs):
         circle1 = plt.Circle(coord, 0.075, fill=False, **kwargs)
         ax.add_patch(circle1)
 
+    def tabulate(self, Qpts, trace_entity, g):
+        return np.ones_like(Qpts)
+
     def __repr__(self):
         return "Grad"
 

fuse/triples.py

@@ -112,21 +112,10 @@ def to_fiat(self):
             dim = entity[0]
             for i in range(len(dofs)):
                 if entity[1] == dofs[i].trace_entity.id - min_ids[dim]:
-                    entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].append(counter)
-                    if hasattr(dofs[i], "convert_to_fiat_new"):    
-                        nodes.append(dofs[i].convert_to_fiat_new(ref_el, degree))
-                        print("old")
-
-                        print(dofs[i].convert_to_fiat(ref_el, degree).pt_dict)
-                        print(dofs[i].convert_to_fiat(ref_el, degree).target_shape)
-                        print("new")
-
-                        print(dofs[i].convert_to_fiat_new(ref_el, degree).pt_dict)
-                        print(dofs[i].convert_to_fiat_new(ref_el, degree).target_shape)
-                    else:
-                        raise ValueError("using old")
-                        nodes.append(dofs[i].convert_to_fiat(ref_el, degree))
-                    counter += 1
+                    fiat_dofs = dofs[i].convert_to_fiat(ref_el, degree)
+                    nodes.extend(fiat_dofs)
+                    entity_ids[dim][dofs[i].trace_entity.id - min_ids[dim]].extend([counter + i for i in range(len(fiat_dofs))])
+                    counter += len(fiat_dofs)
         entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set)
 
         form_degree = 1 if self.spaces[0].set_shape else 0

test/test_cells.py

@@ -111,10 +111,11 @@ def test_oriented_verts():
         print(g, permuted)
         assert g.permute(tetra.ordered_vertices()) == oriented.ordered_vertices()
 
+
 def test_compare_cell_to_firedrake():
     tri1 = polygon(3)
     tri2 = default_simplex(2)
-    
+
     n = 3
     vertices = []
     for i in range(n):
@@ -123,40 +124,38 @@ def test_compare_cell_to_firedrake():
     for i in range(n):
         edges.append(
             Point(1, [vertices[(i) % n], vertices[(i+1) % n]], vertex_num=2))
-    from sympy.combinatorics.named_groups import SymmetricGroup
-    s3 = SymmetricGroup(3)
+
     cellS3 = S3.add_cell(tri1)
     for g in cellS3.members():
         print(g.perm.array_form)
-        try:
-            p = g.perm.array_form
-            # tri3 = Point(2, [edges[p[0]], edges[p[1]], edges[p[2]]], vertex_num=n)
-            print(tri1.orient(g).get_topology())
-        except:
-            print('FAIL')
-
+        p = g.perm.array_form
+        tri3 = Point(2, [edges[p[0]], edges[p[1]], edges[p[2]]], vertex_num=n)
+        print(tri1.orient(g).get_topology())
+        print(tri3.get_topology())
 
     # print(tri1.get_topology())
     print(tri2.get_topology())
     tri3 = firedrake_triangle()
     print(tri3.get_topology())
 
+
 @pytest.fixture
 def mock_cell_complex(mocker, expect):
     mocker.patch('firedrake.mesh.constructCellComplex', return_value=expect.to_ufl("triangle"))
 
+
 @pytest.mark.usefixtures("mock_cell_complex")
-@pytest.mark.parametrize(["expect"],[(firedrake_triangle(),), (polygon(3),)])
+@pytest.mark.parametrize(["expect"], [(firedrake_triangle(),), (polygon(3),)])
 def test_ref_els(expect):
     scale_range = range(3, 6)
 
     diff2 = [0 for i in scale_range]
     for i in scale_range:
-            mesh = UnitSquareMesh(2 ** i, 2 ** i)
+        mesh = UnitSquareMesh(2 ** i, 2 ** i)
 
-            V = FunctionSpace(mesh, "CG", 3)
-            res1 = helmholtz_solve(mesh, V)
-            diff2[i-3] = res1
+        V = FunctionSpace(mesh, "CG", 3)
+        res1 = helmholtz_solve(mesh, V)
+        diff2[i-3] = res1
 
     print("firedrake l2 error norms:", diff2)
     diff2 = np.array(diff2)

test/test_convert_to_fiat.py

@@ -133,6 +133,7 @@ def create_cg2_tri(cell):
                                                 DOFGenerator(edge_xs, C3, S1)])
     return cg
 
+
 def create_hermite(tri):
     vert = tri.vertices()[0]
 
@@ -145,14 +146,11 @@ def create_hermite(tri):
     v_derv_xs = [immerse(tri, dg0, TrGrad)]
     v_derv_dofs = DOFGenerator(v_derv_xs, S3/S2, S1)
 
-    v_derv2_xs = [immerse(tri, dg0, TrHess)]
-    v_derv2_dofs = DOFGenerator(v_derv2_xs, S3/S2, S1)
-
     i_xs = [DOF(DeltaPairing(), PointKernel((0, 0)))]
     i_dofs = DOFGenerator(i_xs, S1, S1)
 
     her = ElementTriple(tri, (P3, CellH2, C0),
-                        [v_dofs, v_derv_dofs, v_derv2_dofs, i_dofs])
+                        [v_dofs, v_derv_dofs, i_dofs])
     return her
 
 
@@ -323,7 +321,6 @@ def test_helmholtz(elem_gen, elem_code, deg, conv_rate):
     diff2 = np.array(diff2)
     conv1 = np.log2(diff2[:-1] / diff2[1:])
     print("firedrake convergence order:", conv1)
-    
 
     print("fuse l2 error norms:", diff)
     diff = np.array(diff)
@@ -351,10 +348,10 @@ def helmholtz_solve(mesh, V):
     a = (inner(grad(u), grad(v)) + inner(u, v)) * dx
     L = inner(f, v) * dx
     u = Function(V)
-    l_a = assemble(L)
-    elem = V.finat_element.fiat_equivalent
-    W = VectorFunctionSpace(mesh, V.ufl_element())
-    X = assemble(interpolate(mesh.coordinates, W))
+    # l_a = assemble(L)
+    # elem = V.finat_element.fiat_equivalent
+    # W = VectorFunctionSpace(mesh, V.ufl_element())
+    # X = assemble(interpolate(mesh.coordinates, W))
     solve(a == L, u)
     f.interpolate(cos(x*pi*2)*cos(y*pi*2))
     return sqrt(assemble(dot(u - f, u - f) * dx))
@@ -564,7 +561,8 @@ def test_project_3d(elem_gen, elem_code, deg):
 
     assert np.allclose(out.dat.data, f.dat.data, rtol=1e-5)
 
-@pytest.mark.xfail(reason='Derivative nodes to fiat')
+
+@pytest.mark.xfail(reason='Handling generation of multiple fiat nodes from one in permutations')
 def test_create_hermite():
     deg = 3
     cell = polygon(3)
@@ -581,13 +579,13 @@ def test_create_hermite():
     Qpts, _ = Q.get_points(), Q.get_weights()
 
     fiat_vals = fiat_elem.tabulate(0, Qpts)
-    # my_vals = my_elem.tabulate(0, Qpts)
+    my_vals = my_elem.tabulate(0, Qpts)
 
     fiat_vals = flatten(fiat_vals[(0,) * sd])
-    # my_vals = flatten(my_vals[(0,) * sd])
+    my_vals = flatten(my_vals[(0,) * sd])
 
-    # (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T)
-    # x1 = np.linalg.inv(x)
-    # assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0)
-    # assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0)
-    # assert np.allclose(res, 0)
+    (x, res, _, _) = np.linalg.lstsq(fiat_vals.T, my_vals.T)
+    x1 = np.linalg.inv(x)
+    assert np.allclose(np.linalg.norm(my_vals.T - fiat_vals.T @ x), 0)
+    assert np.allclose(np.linalg.norm(fiat_vals.T - my_vals.T @ x1), 0)
+    assert np.allclose(res, 0)

test/test_dofs.py

@@ -1,5 +1,5 @@
 from fuse import *
-from test_convert_to_fiat import create_cg1, create_dg1, construct_cg3, construct_rt, construct_nd
+from test_convert_to_fiat import create_cg1, create_dg1, construct_cg3, construct_rt, construct_nd, create_hermite
 import sympy as sp
 import numpy as np
 
@@ -106,3 +106,26 @@ def test_permute_nd():
     #     print(g)
     #     print(nd.cell.permute_entities(g, 0))
     #     print(nd.cell.permute_entities(g, 1))
+
+
+def test_convert_dofs():
+    cell = polygon(3)
+
+    cg3 = create_hermite(cell)
+
+    for dof in cg3.generate():
+        print(dof)
+        # print("old")
+        # # old = dof.convert_to_fiat(cell.to_fiat(), 5).pt_dict
+        # print(old)
+        # print("new")
+        new = dof.convert_to_fiat_new(cell.to_fiat(), 5)[0].pt_dict
+        print(new)
+        new = dof.convert_to_fiat_new(cell.to_fiat(), 5)[0].deriv_dict
+        print(new)
+    from FIAT.hermite import CubicHermite
+    fiat_elem = CubicHermite(cell.to_fiat(), 3)
+
+    print(len([n.pt_dict for n in fiat_elem.dual.nodes]))
+    print([n.pt_dict for n in fiat_elem.dual.nodes])
+    print([n.deriv_dict for n in fiat_elem.dual.nodes])