HDivTrace variants

FIAT/hdiv_trace.py

@@ -5,6 +5,7 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 import numpy as np
+from FIAT.barycentric_interpolation import get_lagrange_points
 from FIAT.discontinuous_lagrange import DiscontinuousLagrange
 from FIAT.dual_set import DualSet
 from FIAT.finite_element import FiniteElement
@@ -39,14 +40,15 @@ class HDivTrace(FiniteElement):
     arises in several DG formulations.
     """
 
-    def __init__(self, ref_el, degree):
+    def __init__(self, ref_el, degree, variant=None):
         """Constructor for the HDivTrace element.
 
         :arg ref_el: A reference element, which may be a tensor product
                      cell.
         :arg degree: The degree of approximation. If on a tensor product
                      cell, then provide a tuple of degrees if you want
                      varying degrees.
+        :arg variant: The point distribution variant passed on to recursivenodes.
         """
         sd = ref_el.get_spatial_dimension()
         if sd in (0, 1):
@@ -83,7 +85,7 @@ def __init__(self, ref_el, degree):
 
             # We have a facet entity!
             if cell.get_spatial_dimension() == facet_sd:
-                dg_elements[top_dim] = construct_dg_element(cell, degree)
+                dg_elements[top_dim] = construct_dg_element(cell, degree, variant)
             # Initialize
             for entity in entities:
                 entity_dofs[top_dim][entity] = []
@@ -97,15 +99,14 @@ def __init__(self, ref_el, degree):
             nf = element.space_dimension()
             num_facets = len(topology[facet_dim])
 
+            facet_pts = get_lagrange_points(element.dual_basis())
             for i in range(num_facets):
                 entity_dofs[facet_dim][i] = list(range(offset, offset + nf))
                 offset += nf
 
                 # Run over nodes and collect the points for point evaluations
-                for dof in element.dual_basis():
-                    facet_pt, = dof.get_point_dict()
-                    transform = ref_el.get_entity_transform(facet_dim, i)
-                    pts.append(tuple(transform(facet_pt)))
+                transform = ref_el.get_entity_transform(facet_dim, i)
+                pts.extend(transform(facet_pts))
 
         # Setting up dual basis - only point evaluations
         nodes = [PointEvaluation(ref_el, pt) for pt in pts]
@@ -265,27 +266,26 @@ def is_nodal():
         return True
 
 
-def construct_dg_element(ref_el, degree):
+def construct_dg_element(ref_el, degree, variant):
     """Constructs a discontinuous galerkin element of a given degree
     on a particular reference cell.
     """
     if ref_el.get_shape() in [LINE, TRIANGLE]:
-        dg_element = DiscontinuousLagrange(ref_el, degree)
+        dg_element = DiscontinuousLagrange(ref_el, degree, variant)
 
     # Quadrilateral facets could be on a FiredrakeQuadrilateral.
     # In this case, we treat this as an interval x interval cell:
     elif ref_el.get_shape() == QUADRILATERAL:
-        dg_a = DiscontinuousLagrange(ufc_simplex(1), degree)
-        dg_b = DiscontinuousLagrange(ufc_simplex(1), degree)
-        dg_element = TensorProductElement(dg_a, dg_b)
+        dg_line = DiscontinuousLagrange(ufc_simplex(1), degree, variant)
+        dg_element = TensorProductElement(dg_line, dg_line)
 
     # This handles the more general case for facets:
     elif ref_el.get_shape() == TENSORPRODUCT:
         assert len(degree) == len(ref_el.cells), (
             "Must provide the same number of degrees as the number "
             "of cells that make up the tensor product cell."
         )
-        sub_elements = [construct_dg_element(c, d)
+        sub_elements = [construct_dg_element(c, d, variant)
                         for c, d in zip(ref_el.cells, degree)
                         if c.get_shape() != POINT]
 

finat/__init__.py

@@ -3,8 +3,9 @@
                             BrezziDouglasMariniCubeFace, CrouzeixRaviart,                    # noqa: F401
                             DiscontinuousLagrange, DiscontinuousTaylor, DPC,                 # noqa: F401
                             FacetBubble, GopalakrishnanLedererSchoberlFirstKind,             # noqa: F401
-                            GopalakrishnanLedererSchoberlSecondKind, HellanHerrmannJohnson,  # noqa: F401
-                            KongMulderVeldhuizen, Lagrange, Real, Serendipity,               # noqa: F401
+                            GopalakrishnanLedererSchoberlSecondKind, HDivTrace,              # noqa: F401
+                            HellanHerrmannJohnson, KongMulderVeldhuizen,                     # noqa: F401
+                            Lagrange, Real, Serendipity,                                     # noqa: F401
                             TrimmedSerendipityCurl, TrimmedSerendipityDiv,                   # noqa: F401
                             TrimmedSerendipityEdge, TrimmedSerendipityFace,                  # noqa: F401
                             Nedelec, NedelecSecondKind, RaviartThomas, Regge)                # noqa: F401
@@ -24,9 +25,7 @@
 from .johnson_mercier import JohnsonMercier                        # noqa: F401
 from .mtw import MardalTaiWinther                                  # noqa: F401
 from .morley import Morley                                         # noqa: F401
-from .trace import HDivTrace                                       # noqa: F401
 from .direct_serendipity import DirectSerendipity                  # noqa: F401
-
 from .spectral import GaussLobattoLegendre, GaussLegendre, Legendre, IntegratedLegendre, FDMLagrange, FDMQuadrature, FDMDiscontinuousLagrange, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
 from .tensorfiniteelement import TensorFiniteElement               # noqa: F401
 from .tensor_product import TensorProductElement                   # noqa: F401

finat/element_factory.py

@@ -216,6 +216,9 @@ def convert_finiteelement(element, **kwargs):
             lmbda = finat.DiscontinuousLagrange
             finat_kwargs["variant"] = kind
 
+    elif element.family() == "HDiv Trace":
+        finat_kwargs["variant"] = kind
+
     elif element.variant() is not None:
         finat_kwargs["variant"] = element.variant()
 

finat/fiat_elements.py

@@ -397,6 +397,11 @@ def __init__(self, cell, degree):
         super().__init__(FIAT.DiscontinuousTaylor(cell, degree))
 
 
+class HDivTrace(ScalarFiatElement):
+    def __init__(self, cell, degree, variant=None):
+        super().__init__(FIAT.HDivTrace(cell, degree, variant=variant))
+
+
 class VectorFiatElement(FiatElement):
     @property
     def value_shape(self):

test/FIAT/unit/test_fiat.py

@@ -552,11 +552,11 @@ def test_nodality_tabulate(element):
     "HDivTrace(T, 0)",
     "HDivTrace(T, 1)",
     "HDivTrace(T, 2)",
-    "HDivTrace(T, 3)",
+    "HDivTrace(T, 3, 'spectral')",
     "HDivTrace(S, 0)",
     "HDivTrace(S, 1)",
     "HDivTrace(S, 2)",
-    "HDivTrace(S, 3)",
+    "HDivTrace(S, 3, 'spectral')",
 ])
 def test_facet_nodality_tabulate(element):
     """Check that certain elements (which do no implement get_nodal_basis)