Merge pull request #47 from firedrakeproject/pbrubeck/dmat-without-sympy

Construct differentiation matrix without sympy

FIAT/barycentric_interpolation.py

@@ -28,7 +28,6 @@ class LagrangeLineExpansionSet(expansions.LineExpansionSet):
     via the second barycentric interpolation formula. See Berrut and Trefethen (2004)
     https://doi.org/10.1137/S0036144502417715 Eq. (4.2) & (9.4)
     """
-
     def __init__(self, ref_el, pts):
         self.nodes = numpy.array(pts).flatten()
         self.dmat, self.weights = make_dmat(self.nodes)
@@ -37,6 +36,9 @@ def __init__(self, ref_el, pts):
     def get_num_members(self, n):
         return len(self.nodes)
 
+    def get_dmats(self, degree):
+        return [self.dmat.T]
+
     def tabulate(self, n, pts):
         assert n == len(self.nodes)-1
         results = numpy.add.outer(-self.nodes, numpy.array(pts).flatten())
@@ -51,8 +53,15 @@ def tabulate(self, n, pts):
             results = numpy.array(list(map(simplify, results)))
         return results
 
-    def tabulate_derivatives(self, n, pts):
-        return numpy.dot(self.dmat, self.tabulate(n, pts))
+    def _tabulate(self, n, pts, order=0):
+        results = [self.tabulate(n, pts)]
+        for r in range(order):
+            results.append(numpy.dot(self.dmat, results[-1]))
+        for r in range(order+1):
+            shape = results[r].shape
+            shape = shape[:1] + (1,)*r + shape[1:]
+            results[r] = numpy.reshape(results[r], shape)
+        return results
 
 
 class LagrangePolynomialSet(polynomial_set.PolynomialSet):
@@ -67,7 +76,10 @@ def __init__(self, ref_el, pts, shape=tuple()):
         num_exp_functions = expansions.polynomial_dimension(ref_el, degree)
         num_members = num_components * num_exp_functions
         embedded_degree = degree
-        expansion_set = get_expansion_set(ref_el, pts)
+        if ref_el.get_shape() == reference_element.LINE:
+            expansion_set = LagrangeLineExpansionSet(ref_el, pts)
+        else:
+            raise ValueError("Invalid reference element type.")
 
         # set up coefficients
         if shape == tuple():
@@ -84,15 +96,5 @@ def __init__(self, ref_el, pts, shape=tuple()):
                     coeffs[cur_idx] = 1.0
                     cur_bf += 1
 
-        dmats = [numpy.transpose(expansion_set.dmat)]
         super(LagrangePolynomialSet, self).__init__(ref_el, degree, embedded_degree,
-                                                    expansion_set, coeffs, dmats)
-
-
-def get_expansion_set(ref_el, pts):
-    """Returns an ExpansionSet instance appopriate for the given
-    reference element."""
-    if ref_el.get_shape() == reference_element.LINE:
-        return LagrangeLineExpansionSet(ref_el, pts)
-    else:
-        raise ValueError("Invalid reference element type.")
+                                                    expansion_set, coeffs)

FIAT/brezzi_douglas_fortin_marini.py

@@ -95,8 +95,7 @@ def BDFMSpace(ref_el, order):
                                               order,
                                               order,
                                               vec_poly_set.get_expansion_set(),
-                                              new_coeffs,
-                                              vec_poly_set.get_dmats())
+                                              new_coeffs)
 
     element_set = polynomial_set.polynomial_set_union_normalized(bubble_set, vec_poly_set)
     return element_set