Refactor to allow derivative dof conversion

Author: indiamai

fuse/dof.py

@@ -1,6 +1,6 @@
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.quadrature import FacetQuadratureRule
-from FIAT.functional import PointEvaluation, FrobeniusIntegralMoment
+from FIAT.functional import PointEvaluation, FrobeniusIntegralMoment, Functional
 from fuse.utils import sympy_to_numpy
 import numpy as np
 import sympy as sp
@@ -40,12 +40,7 @@ def convert_to_fiat(self, ref_el, dof, interpolant_deg):
 
     def get_pts(self, ref_el, total_degree):
         entity = ref_el.construct_subelement(self.entity.dim())
-        Q_ref = create_quadrature(entity, total_deg)
-
-        ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()]
-        Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref)
-        qpts, qwts = Q.get_points(), Q.get_weights()
-        return [(1,)], [(1,)]
+        return [(0,) * entity.get_spatial_dimension()], [1], 1
 
     def add_entity(self, entity):
         res = DeltaPairing()
@@ -101,12 +96,15 @@ def add_entity(self, entity):
 
     def get_pts(self, ref_el, total_degree):
         entity = ref_el.construct_subelement(self.entity.dim())
-        Q_ref = create_quadrature(entity, total_deg)
+        Q_ref = create_quadrature(entity, total_degree)
 
         ent_id = self.entity.id - ref_el.fe_cell.get_starter_ids()[self.entity.dim()]
         Q = FacetQuadratureRule(ref_el, self.entity.dim(), ent_id, Q_ref)
-        qpts, qwts = Q.get_points(), Q.get_weights()
-        return qpts, qwts
+        Jdet = Q.jacobian_determinant()
+        # TODO work out how to get J from attachment
+
+        qpts, qwts = Q_ref.get_points(), Q_ref.get_weights()
+        return qpts, qwts, Jdet
 
     def convert_to_fiat(self, ref_el, dof, interpolant_degree):
         total_deg = interpolant_degree + dof.kernel.degree()
@@ -168,8 +166,11 @@ def permute(self, g):
     def __call__(self, *args):
         return self.pt
 
-    def tabulate(self, Qpts):
-        return np.array([self.pt for _ in Qpts]).astype(np.float64)
+    def tabulate(self, Qpts, attachment=None):
+        if attachment is None:
+            return np.array([self.pt for _ in Qpts]).astype(np.float64)
+        else:
+            return np.array([attachment(*self.pt) for _ in Qpts]).astype(np.float64)
 
     def _to_dict(self):
         o_dict = {"pt": self.pt}
@@ -209,8 +210,10 @@ def __call__(self, *args):
             return [res]
         return res
 
-    def tabulate(self, Qpts):
-        return np.array([self(*pt) for pt in Qpts]).astype(np.float64)
+    def tabulate(self, Qpts, attachment=None):
+        if attachment is None:
+            return np.array([self(*pt) for pt in Qpts]).astype(np.float64)
+        return np.array([self(*attachment(*pt)) for pt in Qpts]).astype(np.float64)
 
     def _to_dict(self):
         o_dict = {"fn": self.fn}
@@ -270,7 +273,23 @@ 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 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)}
+        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__())
+
 
     def __repr__(self, fn="v"):
         return str(self.pairing).format(fn=fn, kernel=self.kernel)
@@ -309,13 +328,25 @@ 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)
-        # [self.attachment(*tuple(r)) for r in res]
         return immersion*res
 
-    def convert_to_fiat(self, ref_el):
-        pts, wts = self.pairing.get_pts()
-        self.target_space.convert_to_fiat(tabulated, wts)
+    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, 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)
+
+        # 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__())
 
     def __call__(self, g):
         permuted = self.cell.permute_entities(g, self.trace_entity.dim())

fuse/traces.py

@@ -54,6 +54,10 @@ def plot(self, ax, coord, trace_entity, g, **kwargs):
     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, {}
+
     def __repr__(self):
         return "H1"
 
@@ -86,6 +90,15 @@ def plot(self, ax, coord, trace_entity, g, **kwargs):
             vec = np.cross(basis[0], basis[1])
         ax.quiver(*coord, *vec, **kwargs)
 
+    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, {}
+
+
     def tabulate(self, Qpts, trace_entity, g):
         entityBasis = np.array(trace_entity.basis_vectors())
         cellEntityBasis = np.array(self.domain.basis_vectors(entity=trace_entity))
@@ -122,6 +135,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, {}
 
     def plot(self, ax, coord, trace_entity, g, **kwargs):
         permuted = self.domain.permute_entities(g, trace_entity.dimension)

fuse/triples.py

@@ -113,7 +113,19 @@ def to_fiat(self):
             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)
-                    nodes.append(dofs[i].convert_to_fiat(ref_el, degree))
+                    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
         entity_perms, pure_perm = self.make_dof_perms(ref_el, entity_ids, nodes, poly_set)
 
@@ -189,7 +201,7 @@ def compute_dense_matrix(self, ref_el, entity_ids, nodes, poly_set):
         A = dualmat.reshape((shp[0], -1))
         B = old_coeffs.reshape((shp[0], -1))
         V = np.dot(A, np.transpose(B))
-
+        print(V)
         with warnings.catch_warnings():
             warnings.filterwarnings("error")
             try:

test/test_convert_to_fiat.py

@@ -564,7 +564,7 @@ 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')
 def test_create_hermite():
     deg = 3
     cell = polygon(3)