wence/fix/serendipity tensor product (#62)

* pointwise_dual: Get correct spatial dimension

* serendipity: Correct default entity

We need to use get_dimension so that this does the right thing on
tensor product cells.

* serendipity: Use flattened element to construct unisolvent points

* Add test of duals for serendipity in tensor and ufc element cases

Co-authored-by: Rob Kirby <[email protected]>

FIAT/pointwise_dual.py

@@ -32,7 +32,7 @@ def compute_pointwise_dual(el, pts):
     nbf = el.space_dimension()
 
     T = el.ref_el
-    sd = T.get_dimension()
+    sd = T.get_spatial_dimension()
 
     assert np.asarray(pts).shape == (int(nbf / np.prod(el.value_shape())), sd)
 

FIAT/serendipity.py

@@ -143,7 +143,7 @@ def get_coeffs(self):
     def tabulate(self, order, points, entity=None):
 
         if entity is None:
-            entity = (self.ref_el.get_spatial_dimension(), 0)
+            entity = (self.ref_el.get_dimension(), 0)
 
         entity_dim, entity_id = entity
         transform = self.ref_el.get_entity_transform(entity_dim, entity_id)
@@ -248,9 +248,9 @@ def unisolvent_pts(K, deg):
     flat_el = flatten_reference_cube(K)
     dim = flat_el.get_spatial_dimension()
     if dim == 2:
-        return unisolvent_pts_quad(K, deg)
+        return unisolvent_pts_quad(flat_el, deg)
     elif dim == 3:
-        return unisolvent_pts_hex(K, deg)
+        return unisolvent_pts_hex(flat_el, deg)
     else:
         raise ValueError("Serendipity only defined for quads and hexes")
 

test/unit/test_serendipity.py

@@ -1,4 +1,5 @@
-from FIAT.reference_element import UFCQuadrilateral
+from FIAT.reference_element import (
+    UFCQuadrilateral, UFCInterval, TensorProductCell)
 from FIAT import Serendipity
 import numpy as np
 import sympy
@@ -30,3 +31,16 @@ def test_serendipity_derivatives():
         assert actual.shape == (8, 2)
         assert np.allclose(np.asarray(expect, dtype=float),
                            actual.reshape(8, 2))
+
+
+def test_dual_tensor_versus_ufc():
+    K0 = UFCQuadrilateral()
+    ell = UFCInterval()
+    K1 = TensorProductCell(ell, ell)
+    S0 = Serendipity(K0, 2)
+    S1 = Serendipity(K1, 2)
+    # since both elements go through the flattened cell to produce the
+    # dual basis, they ought to do *exactly* the same calculations and
+    # hence form exactly the same nodes.
+    for i in range(S0.space_dimension()):
+        assert S0.dual.nodes[i].pt_dict == S1.dual.nodes[i].pt_dict