Extending the types of DOFs that can be converted to FIAT

Makefile

@@ -11,6 +11,9 @@ docs:
 	# put -n back in when things are better documented
 	@(cd docs/ && (make html SPHINXOPTS="-W --keep-going"))
 
+clean:
+	@(cd docs/ && make clean)
+
 lint:
 	@echo "    Linting FUSE codebase"
 	@python3 -m flake8 $(FLAKE8_FORMAT) fuse
@@ -38,8 +41,6 @@ test_cells:
 	@firedrake-clean
 	@python3 -m pytest -rPx --run-cleared test/test_cells.py::test_ref_els[expect1]
 
-clean:
-	@(cd docs/ && make clean)
 
 prepush: lint tests
 	@rm .coverage.*

fuse/__init__.py

@@ -6,7 +6,7 @@
 from fuse.tensor_products import tensor_product
 
 from fuse.spaces.element_sobolev_spaces import CellH1, CellL2, CellHDiv, CellHCurl, CellH2
-from fuse.spaces.polynomial_spaces import P0, P1, P2, P3, Q2, PolynomialSpace
+from fuse.spaces.polynomial_spaces import P0, P1, P2, P3, P5, Q2, Q3, PolynomialSpace
 from fuse.spaces.interpolation_spaces import C0, L2, H1, HDiv
 
 from fuse.element_construction import periodic_table

fuse/cells.py

@@ -767,7 +767,6 @@ def plot3d(self, show=True, plain=False, ax=None, filename=None):
             fig = plt.figure()
             ax = fig.add_subplot(projection='3d')
         xs = np.linspace(-1, 1, 20)
-
         top_level_node = self.d_entities_ids(self.graph_dim())[0]
         nodes = self.d_entities_ids(0)
         min_ids = self.get_starter_ids()

fuse/dof.py

@@ -37,6 +37,10 @@ def __call__(self, kernel, v, cell):
     def tabulate(self):
         return 1
 
+    def get_pts(self, ref_el, total_degree):
+        entity = ref_el.construct_subelement(self.entity.dim())
+        return [(0,) * entity.get_spatial_dimension()], [1], 1
+
     def add_entity(self, entity):
         res = DeltaPairing()
         res.entity = entity
@@ -420,12 +424,13 @@ 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, value_shape=tuple()):
-        # TODO deriv dict needs implementing (currently {})
-        return Functional(ref_el, value_shape, self.to_quadrature(interpolant_degree, value_shape), {}, str(self))
+        pt_dict, deriv_dict = self.to_quadrature(interpolant_degree, value_shape)
+        return Functional(ref_el, value_shape, pt_dict, deriv_dict, str(self))
 
     def to_quadrature(self, arg_degree, value_shape):
         Qpts, Qwts = self.cell_defined_on.quadrature(self.kernel.degree(arg_degree))
         Qwts = Qwts.reshape(Qwts.shape + (1,))
+        deriv_terms = None
         dim = self.cell_defined_on.get_spatial_dimension()
         if dim > 0:
             bvs = np.array(self.cell_defined_on.basis_vectors())
@@ -439,23 +444,13 @@ def immersed(pt):
                 basis = np.array(self.cell_defined_on.basis_vectors()).T
                 basis_coeffs = np.matmul(np.linalg.inv(basis), np.array(pt))
                 J = np.array(self.cell.basis_vectors(entity=self.cell_defined_on)).T
-                # J2 = self.cell.attachment_J(self.cell.id, self.cell_defined_on.id)
-                # if not np.allclose(J2 @ np.array(pt), J @ basis_coeffs):
-                #     breakpoint()
                 return np.matmul(J, basis_coeffs)
         else:
             immersed = self.immersed
 
         if isinstance(self.kernel, BarycentricPolynomialKernel):
-            # if self.pairing.orientation is not None and
-            # self.pairing.orientation.numeric_rep() == 1:
-            #     breakpoint()
-            # print(self)
-            # print(self.cell_defined_on.cartesian_to_barycentric(Qpts))
             pts = [np.matmul(basis_change.T, pt) for pt in Qpts]
             bary_pts = self.cell_defined_on.cartesian_to_barycentric(pts)
-            # print(bary_pts)
-            # print(basis_change)
             pts, wts, comps = self.kernel.evaluate(Qpts, bary_pts, Qwts, basis_change, immersed, self.cell.dimension, value_shape)
         else:
             pts, wts, comps = self.kernel.evaluate(Qpts, Qwts, basis_change, immersed, self.cell.dimension, value_shape)
@@ -472,19 +467,32 @@ def immersed(pt):
                 new_wts = wts
             else:
                 new_wts = np.outer(wts * J_det, immersion)
-                # shape is wrong for 2d face on tet
-            # if isinstance(self.kernel, BarycentricPolynomialKernel) and self.kernel.shape > 1:
-            #     new_wts = np.array([self.cell.attachment(self.cell.id, self.cell_defined_on.id)(*pt) for pt in new_wts])
+
+            # a derivative dof is described as sum_i coeff_i * D^{alpha_i},
+            # already expressed w.r.t. the ambient cell's reference frame
+            # (matching FIAT's own PointDerivative/PointDirectionalDerivative)
+            deriv_terms = self.target_space.tabulate_derivs(pts, self.cell_defined_on)
         else:
             new_wts = wts
         # pt dict is { pt: [(weight, component)]}
-        pt_dict = {tuple(pt): [(w, c) for w, c in zip(wt, cp)] for pt, wt, cp in zip(pts, new_wts, comps)}
-        # if self.cell_defined_on.dimension >= 2:
-        #     print(self)
-        #     np.set_printoptions(linewidth=90, precision=4, suppress=True)
-        #     for key, val in pt_dict.items():
-        #         print(np.array(key), ":", np.array([v[0] for v in val]))
-        return pt_dict
+        # a pure derivative dof has no point-value contribution, so FIAT expects
+        # those points to be absent from pt_dict entirely rather than mapped to [] if (list(zip(wt, cp)) ensures this
+        pt_dict = {tuple(pt): [(w, c) for w, c in zip(wt, cp)] for pt, wt, cp in zip(pts, new_wts, comps) if list(zip(wt, cp))}
+        # deriv dict is {pt: [(weight, alpha, component)]}
+        if deriv_terms is None:
+            deriv_dict = {}
+        else:
+            deriv_dict = {tuple(pt): [(w[0] * coeff * J_det, alpha, cp[0]) for coeff, alpha in deriv_terms]
+                          for pt, w, cp in zip(pts, wts, comps)}
+        # print(self)
+        # np.set_printoptions(linewidth=90, precision=4, suppress=True)
+        # print("pt")
+        # for key, val in pt_dict.items():
+        #     print(np.array(key), ":", np.array([v[0] for v in val]))
+        # print("deriv")
+        # for key, val in deriv_dict.items():
+        #     print(np.array(key), ":", np.array([v[0] for v in val]))
+        return pt_dict, deriv_dict
 
     def __repr__(self, fn="v"):
         return str(self.pairing).format(fn=fn, kernel=self.kernel)

fuse/spaces/polynomial_spaces.py

@@ -230,6 +230,7 @@ def _from_dict(obj_dict):
 P2 = PolynomialSpace(2)
 P3 = PolynomialSpace(3)
 P4 = PolynomialSpace(4)
+P5 = PolynomialSpace(5)
 
 Q1 = PolynomialSpace(1, 2)
 Q2 = PolynomialSpace(2, 3)

fuse/tensor_products.py

@@ -43,8 +43,6 @@ def to_ufl(self):
         if self.flat:
             return FuseElement(self, self.cell.flatten().to_ufl())
         ufl_sub_elements = [e.to_ufl() for e in self.sub_elements()]
-        # self.setup_matrices()
-        # breakpoint()
         return TensorProductElement(*ufl_sub_elements, cell=self.cell.to_ufl())
 
     def flatten(self):

fuse/traces.py

@@ -1,13 +1,63 @@
 import numpy as np
 import sympy as sp
 import matplotlib.pyplot as plt
+from collections import defaultdict
+from functools import reduce
 from fuse.utils import sympy_to_numpy, numpy_to_str_tuple
 
 
+def _resolve_direction(spec, domain, trace_entity):
+    """Resolve a direction spec - a fixed ambient vector, or one of the
+    keywords "tangent"/"normal" - to a concrete vector, given the facet
+    (trace_entity) immersed within domain. "tangent" mirrors
+    TrHCurl.tabulate; "normal" mirrors TrHDiv.tabulate."""
+    if not isinstance(spec, str):
+        return np.asarray(spec, dtype=float)
+    sd = domain.get_spatial_dimension()
+    basis = np.array(domain.basis_vectors(entity=trace_entity))
+    if spec == "tangent":
+        if trace_entity.dimension != 1:
+            raise ValueError('"tangent" direction requires a 1D (edge) entity')
+        return basis[0]
+    if spec == "normal":
+        if trace_entity.dimension != sd - 1:
+            raise ValueError('"normal" direction is only defined on facets (codimension 1 entities)')
+        if sd == 2:
+            return np.matmul(basis, np.array([[0, -1], [1, 0]]))[0]
+        if sd == 3:
+            return np.cross(basis[0], basis[1])
+        raise ValueError("normal direction not implemented in dimension > 3")
+    raise ValueError(f"Unknown direction keyword {spec!r}")
+
+
+def _directional_deriv_terms(directions, domain, trace_entity):
+    """Expand the order-k mixed directional derivative d/dv_1 ... d/dv_k into
+    FIAT-style [(coeff, alpha)] terms, by taking the outer product of the k
+    resolved direction vectors and accumulating entries that land on the same
+    multi-index - generalizing FIAT's own PointSecondDerivative
+    (FIAT/functional.py, which does exactly this for k=2 via numpy.outer and
+    a defaultdict keyed by alpha) to arbitrary k."""
+    sd = domain.get_spatial_dimension()
+    vectors = [_resolve_direction(s, domain, trace_entity) for s in directions]
+    tensor = reduce(np.multiply.outer, vectors)
+    tau = defaultdict(float)
+    for index in np.ndindex(tensor.shape):
+        alpha = [0] * sd
+        for i in index:
+            alpha[i] += 1
+        tau[tuple(alpha)] += tensor[index]
+    return [(coeff, alpha) for alpha, coeff in tau.items()]
+
+
 class Trace():
 
-    def __init__(self, cell):
+    def __init__(self, cell=None, alpha=None, directions=None):
         self.domain = cell
+        self.alpha = alpha
+        self.directions = directions
+
+    def add_cell(self, cell):
+        return type(self)(cell=cell, alpha=self.alpha, directions=self.directions)
 
     def __call__(self, trace_entity):
         raise NotImplementedError("Trace uninstanitated")
@@ -18,6 +68,15 @@ def plot(self, ax, coord, trace_entity, **kwargs):
     def tabulate(self, Qwts, trace_entity):
         raise NotImplementedError("Tabulation uninstantiated")
 
+    def tabulate_derivs(self, Qwts, trace_entity):
+        if self.alpha is not None and self.directions is not None:
+            raise ValueError("Specify either alpha or directions, not both")
+        if self.directions is not None:
+            return _directional_deriv_terms(self.directions, self.domain, trace_entity)
+        if self.alpha is None:
+            return None
+        return [(1.0, self.alpha)]
+
     def _to_dict(self):
         return {"trace": str(self)}
 
@@ -42,9 +101,6 @@ def _from_dict(obj_dict):
 
 class TrH1(Trace):
 
-    def __init__(self, cell):
-        super(TrH1, self).__init__(cell)
-
     def __call__(self, v, trace_entity):
         return v
 
@@ -67,9 +123,6 @@ def __repr__(self):
 
 class TrHDiv(Trace):
 
-    def __init__(self, cell):
-        super(TrHDiv, self).__init__(cell)
-
     def __call__(self, v, trace_entity):
         def apply(*x):
             result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze())
@@ -125,9 +178,6 @@ def __repr__(self):
 
 class TrHCurl(Trace):
 
-    def __init__(self, cell):
-        super(TrHCurl, self).__init__(cell)
-
     def __call__(self, v, trace_entity):
         def apply(*x):
             result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze())
@@ -162,32 +212,38 @@ def __repr__(self):
 
 class TrGrad(Trace):
 
-    def __init__(self, cell):
-        super(TrGrad, self).__init__(cell)
-
     def __call__(self, v, trace_entity):
         # Compute grad v and then dot with tangent rotated according to the group member
-        raise NotImplementedError("Gradient immersions are under development")
-        g = None
-        tangent = np.array(g(np.array(self.domain.basis_vectors())[0]))
-
+        # raise NotImplementedError("Gradient immersions are under development")
         def apply(*x):
-            X = sp.DeferredVector('x')
-            dX = tuple([X[i] for i in range(self.domain.dim())])
-            compute_v = v(*dX, sym=True)
-            grad_v = sp.Matrix([sp.diff(compute_v, dX[i]) for i in range(len(dX))])
-            eval_grad_v = sympy_to_numpy(grad_v, dX, v.attach_func(*x))
-            result = np.dot(tangent, np.array(eval_grad_v))
-
-            if not hasattr(result, "__iter__"):
+            result = np.dot(self.tabulate(None, trace_entity), np.array(v(*x)).squeeze())
+            if isinstance(result, np.float64):
                 return (result,)
             return tuple(result)
         return apply
 
+    def convert_to_fiat(self, qpts, pts, wts):
+        shp = (self.domain.get_spatial_dimension(),)
+        alphas = []
+        for i in range(pts.shape[0]):
+            new = np.zeros(shp, dtype=int)
+            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, **kwargs):
         circle1 = plt.Circle(coord, 0.075, fill=False, **kwargs)
         ax.add_patch(circle1)
 
+    def tabulate(self, Qpts, trace_entity):
+        return np.array([])
+
     def to_tikz(self, coord, trace_entity, scale, color="black"):
         return f"\\draw[{color}] {numpy_to_str_tuple(coord, scale)} circle (4pt) node[anchor = south] {{}};"
 
@@ -197,9 +253,6 @@ def __repr__(self):
 
 class TrHess(Trace):
 
-    def __init__(self, cell):
-        super(TrHess, self).__init__(cell)
-
     def __call__(self, v, trace_entity):
         raise NotImplementedError("Hessian trace needs reviewing")
         g = None
@@ -219,6 +272,9 @@ def apply(*x):
             return tuple(result)
         return apply
 
+    def tabulate(self, Qpts, trace_entity):
+        return np.array([])
+
     def plot(self, ax, coord, trace_entity, **kwargs):
         circle1 = plt.Circle(coord, 0.15, fill=False, **kwargs)
         ax.add_patch(circle1)