Merge branch 'meg/fix-dolfin-issue-871'

FIAT/hdiv_trace.py

@@ -275,6 +275,10 @@ def get_num_members(self, arg):
         """Return number of members of the expansion set."""
         raise NotImplementedError("get_num_members not implemented for the trace element.")
 
+    @staticmethod
+    def is_nodal():
+        return True
+
 
 def construct_dg_element(ref_el, degree):
     """Constructs a discontinuous galerkin element of a given degree

test/unit/test_fiat.py

@@ -35,7 +35,7 @@
 from FIAT.nedelec import Nedelec                                # noqa: F401
 from FIAT.nedelec_second_kind import NedelecSecondKind          # noqa: F401
 from FIAT.regge import Regge                                    # noqa: F401
-from FIAT.hdiv_trace import HDivTrace                           # noqa: F401
+from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet   # noqa: F401
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson  # noqa: F401
 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini  # noqa: F401
 from FIAT.gauss_legendre import GaussLegendre                   # noqa: F401
@@ -226,6 +226,14 @@ def __init__(self, a, b):
                "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
                ")"),
     # Following elements are checked using tabulate
+    xfail_impl("HDivTrace(T, 0)"),
+    xfail_impl("HDivTrace(T, 1)"),
+    xfail_impl("HDivTrace(T, 2)"),
+    xfail_impl("HDivTrace(T, 3)"),
+    xfail_impl("HDivTrace(S, 0)"),
+    xfail_impl("HDivTrace(S, 1)"),
+    xfail_impl("HDivTrace(S, 2)"),
+    xfail_impl("HDivTrace(S, 3)"),
     xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"),
     xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"),
     xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"),
@@ -347,24 +355,80 @@ def test_mixed_is_not_nodal():
     "FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))",
     "FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))",
 ])
-def test_nodality_tpe(element):
-    """Check that generated TensorProductElements are nodal, i.e. nodes evaluated
-    on basis functions give Kronecker delta
+def test_nodality_tabulate(element):
+    """Check that certain elements (which do no implement
+    get_nodal_basis) are nodal too, by tabulating at nodes
+    (assuming nodes are point evaluation)
     """
     # Instantiate element
     element = eval(element)
 
     # Get nodes coordinates
-    nodes_coord = list(next(iter(f.get_point_dict().keys())) for f in element.dual_basis())
+    nodes_coords = []
+    for node in element.dual_basis():
+        # Assume point evaluation
+        (coords, weights), = node.get_point_dict().items()
+        assert weights == [(1.0, ())]
+
+        nodes_coords.append(coords)
 
     # Check nodality
-    for j, x in enumerate(nodes_coord):
+    for j, x in enumerate(nodes_coords):
         table = element.tabulate(0, (x,))
         basis = table[list(table.keys())[0]]
         for i in range(len(basis)):
             assert np.isclose(basis[i], 1.0 if i == j else 0.0)
 
 
+@pytest.mark.parametrize('element', [
+    "HDivTrace(T, 0)",
+    "HDivTrace(T, 1)",
+    "HDivTrace(T, 2)",
+    "HDivTrace(T, 3)",
+    "HDivTrace(S, 0)",
+    "HDivTrace(S, 1)",
+    "HDivTrace(S, 2)",
+    "HDivTrace(S, 3)",
+])
+def test_facet_nodality_tabulate(element):
+    """Check that certain elements (which do no implement get_nodal_basis)
+    are nodal too, by tabulating facet-wise at nodes (assuming nodes
+    are point evaluation)
+    """
+    # Instantiate element
+    element = eval(element)
+
+    # Dof/Node coordinates and respective facet
+    nodes_coords = []
+
+    # Iterate over facet degrees of freedom
+    entity_dofs = element.dual.entity_ids
+    facet_dim = sorted(entity_dofs.keys())[-2]
+    facet_dofs = entity_dofs[facet_dim]
+    dofs = element.dual_basis()
+    vertices = element.ref_el.vertices
+
+    for (facet, indices) in facet_dofs.items():
+        for i in indices:
+            node = dofs[i]
+            # Assume point evaluation
+            (coords, weights), = node.get_point_dict().items()
+            assert weights == [(1.0, ())]
+
+            # Map dof coordinates to reference element due to
+            # HdivTrace interface peculiarity
+            ref_coords = map_to_reference_facet((coords,), vertices, facet)
+            assert len(ref_coords) == 1
+            nodes_coords.append((facet, ref_coords[0]))
+
+    # Check nodality
+    for j, (facet, x) in enumerate(nodes_coords):
+        table = element.tabulate(0, (x,), entity=(facet_dim, facet))
+        basis = table[list(table.keys())[0]]
+        for i in range(len(basis)):
+            assert np.isclose(basis[i], 1.0 if i == j else 0.0)
+
+
 if __name__ == '__main__':
     import os
     pytest.main(os.path.abspath(__file__))