Fix RaviartThomas on the interval (#71)

Author: pbrubeck

FIAT/expansions.py

@@ -309,6 +309,8 @@ def get_cell_node_map(self, n):
             return self._cell_node_map_cache.setdefault(n, cell_node_map)
 
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
+        """Returns a dict of tabulations such that
+        tabulations[alpha][i, j] = D^alpha phi_i(pts[j])."""
         from FIAT.polynomial_set import mis
         lorder = min(order, self.recurrence_order)
         A, b = self.affine_mappings[cell]
@@ -468,10 +470,11 @@ def __init__(self, ref_el, **kwargs):
             raise ValueError("Must have a point")
         super(PointExpansionSet, self).__init__(ref_el, **kwargs)
 
-    def tabulate(self, n, pts):
-        """Returns a numpy array A[i,j] = phi_i(pts[j]) = 1.0."""
-        assert n == 0
-        return numpy.ones((1, len(pts)))
+    def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
+        """Returns a dict of tabulations such that
+        tabulations[alpha][i, j] = D^alpha phi_i(pts[j])."""
+        assert n == 0 and order == 0
+        return {(): numpy.ones((1, len(pts)))}
 
 
 class LineExpansionSet(ExpansionSet):
@@ -482,8 +485,8 @@ def __init__(self, ref_el, **kwargs):
         super(LineExpansionSet, self).__init__(ref_el, **kwargs)
 
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
-        """Returns a tuple of (vals, derivs) such that
-        vals[i,j] = phi_i(pts[j]), derivs[i,j] = D vals[i,j]."""
+        """Returns a dict of tabulations such that
+        tabulations[alpha][i, j] = D^alpha phi_i(pts[j])."""
         if self.variant is not None:
             return super(LineExpansionSet, self)._tabulate_on_cell(n, pts, order=order, cell=cell, direction=direction)
 

FIAT/quadrature.py

@@ -24,8 +24,12 @@ def map_quadrature(pts_ref, wts_ref, source_cell, target_cell, jacobian=False):
     A, b = reference_element.make_affine_mapping(source_cell.get_vertices(),
                                                  target_cell.get_vertices())
     scale = pseudo_determinant(A)
-    pts = numpy.dot(pts_ref.reshape((-1, A.shape[1])), A.T) + b[None, :]
     wts = scale * wts_ref
+    if pts_ref.size == 0:
+        pts = b[None, :]
+    else:
+        pts = numpy.dot(pts_ref.reshape((-1, A.shape[1])), A.T) + b[None, :]
+
     # return immutable types
     pts = tuple(map(tuple, pts))
     wts = tuple(wts.flat)

FIAT/raviart_thomas.py

@@ -72,10 +72,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
 
         if variant == "integral":
             facet = ref_el.get_facet_element()
-            # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1}
-            # degree is q - 1
-            Q_ref = create_quadrature(facet, interpolant_deg + degree - 1)
-            Pq = polynomial_set.ONPolynomialSet(facet, degree - 1)
+            # Facet nodes are \int_F v\cdot n p ds where p \in P_q
+            q = degree - 1
+            Q_ref = create_quadrature(facet, interpolant_deg + q)
+            Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0)
             Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
             for f in top[sd - 1]:
                 cur = len(nodes)
@@ -87,15 +87,15 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                              for phi in phis)
                 entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
 
-            # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-2}^d
-            if degree > 1:
+            # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d
+            if q > 0:
                 cur = len(nodes)
-                Q = create_quadrature(ref_el, interpolant_deg + degree - 2)
-                Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2)
-                Pkm1_at_qpts = Pkm1.tabulate(Q.get_points())[(0,) * sd]
+                Q = create_quadrature(ref_el, interpolant_deg + q - 1)
+                Pqm1 = polynomial_set.ONPolynomialSet(ref_el, q - 1)
+                Pqm1_at_qpts = Pqm1.tabulate(Q.get_points())[(0,) * sd]
                 nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
                              for d in range(sd)
-                             for phi in Pkm1_at_qpts)
+                             for phi in Pqm1_at_qpts)
                 entity_ids[sd][0] = list(range(cur, len(nodes)))
 
         elif variant == "point":

test/unit/test_fiat.py

@@ -132,24 +132,36 @@ def __init__(self, a, b):
     "CrouzeixRaviart(I, 1)",
     "CrouzeixRaviart(T, 1)",
     "CrouzeixRaviart(S, 1)",
+    "RaviartThomas(I, 1)",
+    "RaviartThomas(I, 2)",
+    "RaviartThomas(I, 3)",
     "RaviartThomas(T, 1)",
     "RaviartThomas(T, 2)",
     "RaviartThomas(T, 3)",
     "RaviartThomas(S, 1)",
     "RaviartThomas(S, 2)",
     "RaviartThomas(S, 3)",
+    'RaviartThomas(I, 1, variant="integral")',
+    'RaviartThomas(I, 2, variant="integral")',
+    'RaviartThomas(I, 3, variant="integral")',
     'RaviartThomas(T, 1, variant="integral")',
     'RaviartThomas(T, 2, variant="integral")',
     'RaviartThomas(T, 3, variant="integral")',
     'RaviartThomas(S, 1, variant="integral")',
     'RaviartThomas(S, 2, variant="integral")',
     'RaviartThomas(S, 3, variant="integral")',
+    'RaviartThomas(I, 1, variant="integral(1)")',
+    'RaviartThomas(I, 2, variant="integral(1)")',
+    'RaviartThomas(I, 3, variant="integral(1)")',
     'RaviartThomas(T, 1, variant="integral(1)")',
     'RaviartThomas(T, 2, variant="integral(1)")',
     'RaviartThomas(T, 3, variant="integral(1)")',
     'RaviartThomas(S, 1, variant="integral(1)")',
     'RaviartThomas(S, 2, variant="integral(1)")',
     'RaviartThomas(S, 3, variant="integral(1)")',
+    'RaviartThomas(I, 1, variant="point")',
+    'RaviartThomas(I, 2, variant="point")',
+    'RaviartThomas(I, 3, variant="point")',
     'RaviartThomas(T, 1, variant="point")',
     'RaviartThomas(T, 2, variant="point")',
     'RaviartThomas(T, 3, variant="point")',