diff --git a/FIAT/__init__.py b/FIAT/__init__.py index a68db2bb8..364cf4297 100644 --- a/FIAT/__init__.py +++ b/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, diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py new file mode 100644 index 000000000..ff5f06965 --- /dev/null +++ b/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) diff --git a/finat/__init__.py b/finat/__init__.py index 2bf716c4b..5390fa779 100644 --- a/finat/__init__.py +++ b/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 diff --git a/finat/element_factory.py b/finat/element_factory.py index 64053bbc7..0a7cb2ffc 100644 --- a/finat/element_factory.py +++ b/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'])): diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 0efea5e00..1b2161de2 100644 --- a/finat/fiat_elements.py +++ b/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): ... diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index ba81befa8..ffb192afc 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/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)", diff --git a/test/FIAT/unit/test_gauss_legendre.py b/test/FIAT/unit/test_gauss_legendre.py index 36fd53bc9..d735e3452 100644 --- a/test/FIAT/unit/test_gauss_legendre.py +++ b/test/FIAT/unit/test_gauss_legendre.py @@ -18,6 +18,7 @@ # Authors: # # David Ham +# Pablo Brubeck import pytest import numpy as np diff --git a/test/FIAT/unit/test_gauss_lobatto_legendre.py b/test/FIAT/unit/test_gauss_lobatto_legendre.py index 5a7d0e07d..9edf4a6d7 100644 --- a/test/FIAT/unit/test_gauss_lobatto_legendre.py +++ b/test/FIAT/unit/test_gauss_lobatto_legendre.py @@ -18,6 +18,7 @@ # Authors: # # David Ham +# Pablo Brubeck import pytest import numpy as np diff --git a/test/FIAT/unit/test_histopolation.py b/test/FIAT/unit/test_histopolation.py new file mode 100644 index 000000000..0693eaf74 --- /dev/null +++ b/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 . +# +# 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])