Mimetic spectral elements

FIAT/__init__.py

@@ -62,6 +62,7 @@
 from FIAT.bubble import Bubble, FacetBubble
 from FIAT.hdiv_trace import HDivTrace
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen
+from FIAT.histopolation import Histopolation
 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
 
 # List of supported elements and mapping to element classes
@@ -101,6 +102,7 @@
                       "Gauss-Lobatto-Legendre": GaussLobattoLegendre,
                       "Gauss-Legendre": GaussLegendre,
                       "Gauss-Radau": GaussRadau,
+                      "Histopolation": Histopolation,
                       "Legendre": Legendre,
                       "Integrated Legendre": IntegratedLegendre,
                       "Morley": Morley,

FIAT/histopolation.py

@@ -0,0 +1,70 @@
+# Copyright (C) 2025 Imperial College London and others.
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+#
+# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2025
+
+import numpy
+from FIAT import finite_element, dual_set, functional, quadrature
+from FIAT.reference_element import LINE
+from FIAT.orientation_utils import make_entity_permutations_simplex
+from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points
+from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre
+
+
+class HistopolationDualSet(dual_set.DualSet):
+    r"""The dual basis for 1D histopolation elements.
+
+    We define window functions w_j that satisfy
+
+    \int_K w_j v dx = \ell_j(v)   for all v in P_{k}
+
+    where
+
+    \ell_j(v) = 1/h_j \int_{[x_j, x_{j+1}]} v dx
+
+    is the usual histopolation dual basis.
+
+    The DOFs are defined as integral moments against w_j.
+    """
+    def __init__(self, ref_el, degree):
+        entity_ids = {0: {0: [], 1: []},
+                      1: {0: list(range(0, degree+1))}}
+
+        fe = GaussLobattoLegendre(ref_el, degree+1)
+        points = get_lagrange_points(fe.dual_basis())
+        h = numpy.diff(numpy.reshape(points, (-1,)))
+        B = numpy.diag(1.0 / h[:-1], k=-1)
+        numpy.fill_diagonal(B, -1.0 / h)
+
+        rule = quadrature.GaussLegendreQuadratureLineRule(ref_el, degree+1)
+        self.rule = rule
+
+        phi = fe.tabulate(1, rule.get_points())
+        wts = rule.get_weights()
+        D = phi[(1, )][:-1]
+        A = numpy.dot(numpy.multiply(D, wts), D.T)
+
+        C = numpy.linalg.solve(A, B)
+        F = numpy.dot(C.T, D)
+        nodes = [functional.IntegralMoment(ref_el, rule, f) for f in F]
+
+        entity_permutations = {}
+        entity_permutations[0] = {0: {0: []}, 1: {0: []}}
+        entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree + 1)}
+
+        super().__init__(nodes, ref_el, entity_ids, entity_permutations)
+
+
+class Histopolation(finite_element.CiarletElement):
+    """1D discontinuous element with integral DOFs on GLL subgrid."""
+    def __init__(self, ref_el, degree):
+        if ref_el.shape != LINE:
+            raise ValueError("Histopolation elements are only defined in one dimension.")
+
+        dual = HistopolationDualSet(ref_el, degree)
+        poly_set = LagrangePolynomialSet(ref_el, dual.rule.pts)
+        formdegree = ref_el.get_spatial_dimension()  # n-form
+        super().__init__(poly_set, dual, degree, formdegree)

finat/__init__.py

@@ -3,8 +3,8 @@
                             BrezziDouglasMariniCubeFace, CrouzeixRaviart,                    # noqa: F401
                             DiscontinuousLagrange, DiscontinuousTaylor, DPC,                 # noqa: F401
                             FacetBubble, GopalakrishnanLedererSchoberlFirstKind,             # noqa: F401
-                            GopalakrishnanLedererSchoberlSecondKind, HDivTrace,              # noqa: F401
-                            HellanHerrmannJohnson, KongMulderVeldhuizen,                     # noqa: F401
+                            GopalakrishnanLedererSchoberlSecondKind, Histopolation,          # noqa: F401
+                            HDivTrace, HellanHerrmannJohnson, KongMulderVeldhuizen,          # noqa: F401
                             Lagrange, Real, Serendipity,                                     # noqa: F401
                             TrimmedSerendipityCurl, TrimmedSerendipityDiv,                   # noqa: F401
                             TrimmedSerendipityEdge, TrimmedSerendipityFace,                  # noqa: F401

finat/element_factory.py

@@ -179,7 +179,7 @@ def convert_finiteelement(element, **kwargs):
         kind = 'spectral'  # default variant
 
     if element.family() == "Lagrange":
-        if kind == 'spectral':
+        if kind in ['spectral', 'mimetic']:
             lmbda = finat.GaussLobattoLegendre
         elif element.cell.cellname() == "interval" and kind in cg_interval_variants:
             lmbda = cg_interval_variants[kind]
@@ -200,6 +200,8 @@ def convert_finiteelement(element, **kwargs):
     elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]:
         if kind == 'spectral':
             lmbda = finat.GaussLegendre
+        elif kind == 'mimetic':
+            lmbda = finat.Histopolation
         elif element.cell.cellname() == "interval" and kind in dg_interval_variants:
             lmbda = dg_interval_variants[kind]
         elif any(map(kind.startswith, ['integral', 'demkowicz', 'fdm'])):

finat/fiat_elements.py

@@ -378,6 +378,11 @@ def __init__(self, cell, degree, variant=None):
         super().__init__(FIAT.DiscontinuousLagrange(cell, degree, variant=variant))
 
 
+class Histopolation(ScalarFiatElement):
+    def __init__(self, cell, degree):
+        super().__init__(FIAT.Histopolation(cell, degree))
+
+
 class Real(DiscontinuousLagrange):
     ...
 

test/FIAT/unit/test_fiat.py

@@ -63,6 +63,7 @@
 from FIAT.enriched import EnrichedElement                       # noqa: F401
 from FIAT.nodal_enriched import NodalEnrichedElement
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen    # noqa: F401
+from FIAT.histopolation import Histopolation                    # noqa: F401
 
 P = Point()
 I = UFCInterval()  # noqa: E741
@@ -315,6 +316,9 @@ def __init__(self, a, b):
     "GaussLobattoLegendre(S, 1)",
     "GaussLobattoLegendre(S, 2)",
     "GaussLobattoLegendre(S, 3)",
+    "Histopolation(I, 0)",
+    "Histopolation(I, 1)",
+    "Histopolation(I, 2)",
     "Bubble(I, 2)",
     "Bubble(T, 3)",
     "Bubble(S, 4)",

test/FIAT/unit/test_gauss_legendre.py

@@ -18,6 +18,7 @@
 # Authors:
 #
 # David Ham
+# Pablo Brubeck
 
 import pytest
 import numpy as np

test/FIAT/unit/test_gauss_lobatto_legendre.py

@@ -18,6 +18,7 @@
 # Authors:
 #
 # David Ham
+# Pablo Brubeck
 
 import pytest
 import numpy as np

test/FIAT/unit/test_histopolation.py

@@ -0,0 +1,74 @@
+# Copyright (C) 2025 Imperial College London and others
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+#
+# Authors:
+#
+# Pablo Brubeck
+
+import pytest
+import numpy as np
+
+from FIAT import ufc_simplex, Histopolation, GaussLobattoLegendre
+from FIAT.barycentric_interpolation import get_lagrange_points
+from FIAT.macro import IsoSplit
+from FIAT.quadrature_schemes import create_quadrature
+
+
+@pytest.fixture
+def interval():
+    return ufc_simplex(1)
+
+
+@pytest.mark.parametrize("degree", range(7))
+def test_gll_basis_values(interval, degree):
+    """Ensure that integrating a simple monomial produces the expected results."""
+
+    s = interval
+    q = create_quadrature(s, 2*degree)
+    fe = Histopolation(s, degree)
+    tab = fe.tabulate(0, q.pts)[(0,)]
+
+    for test_degree in range(degree + 1):
+        v = lambda x: sum(x)**test_degree
+        coefs = [n(v) for n in fe.dual.nodes]
+        integral = np.dot(coefs, np.dot(tab, q.wts))
+        reference = q.integrate(v)
+        assert np.allclose(integral, reference, rtol=1e-14)
+
+
+@pytest.mark.parametrize("degree", range(4))
+def test_histopolation_dofs(interval, degree):
+    """Ensure that basis functions integrate to either 1 or 0 on the GLL subgrid."""
+
+    fe = Histopolation(interval, degree)
+    gll = GaussLobattoLegendre(interval, degree+1)
+    pts = get_lagrange_points(gll.dual_basis())
+    x = np.reshape(pts, (-1, ))
+
+    macrocell = IsoSplit(interval, degree+1, variant="gll")
+    assert np.allclose(macrocell.vertices, pts)
+
+    Q = create_quadrature(macrocell, 2*degree)
+    qpts = Q.get_points()
+    qwts = Q.get_weights().reshape((degree+1, -1))
+
+    tab = fe.tabulate(0, qpts)[(0,)]
+    tab = tab.reshape((fe.space_dimension(), degree+1, -1))
+    expected = np.eye(degree+1)
+    for i in range(degree+1):
+        result = np.dot(tab[:, i, :], qwts[i] / (x[i+1] - x[i]))
+        assert np.allclose(result, expected[i])