From 7dfa2d2ed420830e144c1de56060cbc1a73478a5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 29 Jul 2023 13:13:05 +0100 Subject: [PATCH 1/8] Histopolation DG element --- FIAT/__init__.py | 1 + FIAT/histopolation.py | 63 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+) create mode 100644 FIAT/histopolation.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 534a3216f..789bb2388 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -48,6 +48,7 @@ from FIAT.quadrature_element import QuadratureElement # noqa: F401 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen # noqa: F401 from FIAT.hierarchical import Legendre, IntegratedLegendre # noqa: F401 +from FIAT.histopolation import Histopolation # noqa: F401 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite # noqa: F401 # Important functionality diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py new file mode 100644 index 000000000..67b5e2ad3 --- /dev/null +++ b/FIAT/histopolation.py @@ -0,0 +1,63 @@ +# Copyright (C) 2015 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 David A. Ham (david.ham@imperial.ac.uk), 2015 +# +# Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021 + +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 +from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre + + +class HistopolationDualSet(dual_set.DualSet): + """The dual basis for 1D histopolation elements""" + def __init__(self, ref_el, degree): + entity_ids = {0: {0: [], 1: []}, + 1: {0: list(range(0, degree+1))}} + + embedded = GaussLobattoLegendre(ref_el, degree+1) + points = [] + for node in embedded.dual_basis(): + # Assert singleton point for each node. + pt, = node.get_point_dict().keys() + points.append(pt[0]) + h = numpy.diff(numpy.array(points)) + 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 = embedded.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(HistopolationDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + + +class Histopolation(finite_element.CiarletElement): + """1D discontinuous element with average between consecutive Gauss-Lobatto-Legendre points.""" + 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(Histopolation, self).__init__(poly_set, dual, degree, formdegree) From 781ebab8bff18ad323ac6805af444382c858618b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 29 Jul 2023 13:24:54 +0100 Subject: [PATCH 2/8] documentation --- FIAT/histopolation.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py index 67b5e2ad3..5923bda1c 100644 --- a/FIAT/histopolation.py +++ b/FIAT/histopolation.py @@ -17,7 +17,20 @@ class HistopolationDualSet(dual_set.DualSet): - """The dual basis for 1D histopolation elements""" + """The dual basis for 1D histopolation elements. + + We define window functions w_j that satisfy + + \int_{K} w_j v dx = \phi_j(v) for all v \in P_{k} + + where + + \phi_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 the window functions. + """ def __init__(self, ref_el, degree): entity_ids = {0: {0: [], 1: []}, 1: {0: list(range(0, degree+1))}} @@ -52,7 +65,7 @@ def __init__(self, ref_el, degree): class Histopolation(finite_element.CiarletElement): - """1D discontinuous element with average between consecutive Gauss-Lobatto-Legendre points.""" + """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.") From b58106769025f28b1dcb5352e199b5728e1fb7a6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 29 Jul 2023 13:31:39 +0100 Subject: [PATCH 3/8] flake8 --- FIAT/histopolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py index 5923bda1c..706224591 100644 --- a/FIAT/histopolation.py +++ b/FIAT/histopolation.py @@ -21,11 +21,11 @@ class HistopolationDualSet(dual_set.DualSet): We define window functions w_j that satisfy - \int_{K} w_j v dx = \phi_j(v) for all v \in P_{k} + int_{K} w_j v dx = phi_j(v) for all v in P_{k} where - \phi_j(v) = 1/h_j \int_{x_j}^{x_{j+1}} v dx + phi_j(v) = 1/h_j int_{x_j}^{x_{j+1}} v dx is the usual histopolation dual basis. From 154b0e61b7c22086a6bb3a00d82ee6d44829a9ad Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 21 Mar 2025 14:43:04 +0000 Subject: [PATCH 4/8] Mimetic tensor product elements --- FIAT/histopolation.py | 24 ++++++++++-------------- finat/__init__.py | 4 ++-- finat/element_factory.py | 4 +++- finat/fiat_elements.py | 5 +++++ test/FIAT/unit/test_fiat.py | 4 ++++ 5 files changed, 24 insertions(+), 17 deletions(-) diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py index 706224591..d52c2ad86 100644 --- a/FIAT/histopolation.py +++ b/FIAT/histopolation.py @@ -12,7 +12,7 @@ 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 +from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre @@ -21,34 +21,30 @@ class HistopolationDualSet(dual_set.DualSet): We define window functions w_j that satisfy - int_{K} w_j v dx = phi_j(v) for all v in P_{k} + \int_{K} w_j v dx = \ell_j(v) for all v in P_{k} where - phi_j(v) = 1/h_j int_{x_j}^{x_{j+1}} v dx + \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 the window functions. + 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))}} - embedded = GaussLobattoLegendre(ref_el, degree+1) - points = [] - for node in embedded.dual_basis(): - # Assert singleton point for each node. - pt, = node.get_point_dict().keys() - points.append(pt[0]) - h = numpy.diff(numpy.array(points)) + 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 = embedded.tabulate(1, rule.get_points()) + phi = fe.tabulate(1, rule.get_points()) wts = rule.get_weights() D = phi[(1, )][:-1] A = numpy.dot(numpy.multiply(D, wts), D.T) @@ -61,7 +57,7 @@ def __init__(self, ref_el, degree): entity_permutations[0] = {0: {0: []}, 1: {0: []}} entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree + 1)} - super(HistopolationDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class Histopolation(finite_element.CiarletElement): @@ -73,4 +69,4 @@ def __init__(self, ref_el, degree): dual = HistopolationDualSet(ref_el, degree) poly_set = LagrangePolynomialSet(ref_el, dual.rule.pts) formdegree = ref_el.get_spatial_dimension() # n-form - super(Histopolation, self).__init__(poly_set, dual, degree, formdegree) + 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)", From 15d81498b7dc5a3c9b106f7a85cf132ad3a4c885 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 21 Mar 2025 15:17:21 +0000 Subject: [PATCH 5/8] flake --- FIAT/__init__.py | 1 + FIAT/histopolation.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index e55b62fb1..364cf4297 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -102,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 index d52c2ad86..e086008fb 100644 --- a/FIAT/histopolation.py +++ b/FIAT/histopolation.py @@ -17,11 +17,11 @@ class HistopolationDualSet(dual_set.DualSet): - """The dual basis for 1D histopolation elements. + 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} + \int_K w_j v dx = \ell_j(v) for all v in P_{k} where From cf9c30e237c95ced30d0c99bac2c00859db02879 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 25 Mar 2025 12:04:48 +0000 Subject: [PATCH 6/8] Update Copyright --- FIAT/histopolation.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py index e086008fb..ff5f06965 100644 --- a/FIAT/histopolation.py +++ b/FIAT/histopolation.py @@ -1,12 +1,10 @@ -# Copyright (C) 2015 Imperial College London and others. +# 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 David A. Ham (david.ham@imperial.ac.uk), 2015 -# -# Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021 +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2025 import numpy from FIAT import finite_element, dual_set, functional, quadrature From 7321e846922db81c9b83ae8759eb20ed24afd383 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 25 Mar 2025 12:35:04 +0000 Subject: [PATCH 7/8] Add tests for Histopolation --- FIAT/dual_set.py | 3 +- test/FIAT/unit/test_gauss_legendre.py | 1 + test/FIAT/unit/test_gauss_lobatto_legendre.py | 1 + test/FIAT/unit/test_histopolation.py | 73 +++++++++++++++++++ 4 files changed, 77 insertions(+), 1 deletion(-) create mode 100644 test/FIAT/unit/test_histopolation.py diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index c09db974f..ef95750bf 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -281,7 +281,8 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations): parent_permutations = None parent_to_children = ref_el.get_parent_to_children() - if all(isinstance(node, functional.PointEvaluation) for node in nodes): + sd = parent_cell.get_spatial_dimension() + if sd > 1 and all(isinstance(node, functional.PointEvaluation) for node in nodes): # Merge Lagrange dual with lexicographical reordering parent_nodes = [] for dim in sorted(parent_to_children): 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..dad4cfe4d --- /dev/null +++ b/test/FIAT/unit/test_histopolation.py @@ -0,0 +1,73 @@ +# 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]) From 0cd53b00054d7e41389f77588255caa9aac25445 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 26 Mar 2025 11:21:03 +0000 Subject: [PATCH 8/8] fix tests --- FIAT/dual_set.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index ef95750bf..c09db974f 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -281,8 +281,7 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations): parent_permutations = None parent_to_children = ref_el.get_parent_to_children() - sd = parent_cell.get_spatial_dimension() - if sd > 1 and all(isinstance(node, functional.PointEvaluation) for node in nodes): + if all(isinstance(node, functional.PointEvaluation) for node in nodes): # Merge Lagrange dual with lexicographical reordering parent_nodes = [] for dim in sorted(parent_to_children):