From 9435f1a671d8ae72e8344e93389da70fed3ea9ec Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Feb 2026 17:28:31 +0000 Subject: [PATCH 1/5] NodalEnrichment with mismatching ExpansionSets --- FIAT/nodal_enriched.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index 1b258e6d6..75ff407cf 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -12,6 +12,7 @@ from FIAT.dual_set import DualSet from FIAT.finite_element import CiarletElement from FIAT.barycentric_interpolation import LagrangeLineExpansionSet +from FIAT.quadrature_schemes import create_quadrature __all__ = ['NodalEnrichedElement'] @@ -50,6 +51,7 @@ def __init__(self, *elements): expansion_set = elem.get_nodal_basis().get_expansion_set() mapping = elem.mapping()[0] value_shape = elem.value_shape() + sd = ref_el.get_spatial_dimension() # Sanity check assert all(e.get_reference_complex() == ref_el for e in elements) @@ -61,11 +63,21 @@ def __init__(self, *elements): # Obtain coefficients via interpolation points = expansion_set.get_points() coeffs = np.vstack([e.tabulate(0, points)[(0,)] for e in elements]) - else: - assert all(e.get_nodal_basis().get_expansion_set() == expansion_set - for e in elements) + elif all(e.get_nodal_basis().get_expansion_set() == expansion_set for e in elements): + # All elements have the same ExpansionSet, just merge coefficients coeffs = [e.get_coeffs() for e in elements] coeffs = _merge_coeffs(coeffs, ref_el, embedded_degrees, expansion_set.continuity) + else: + # Obtain coefficients via projection + Q = create_quadrature(ref_el, 2*embedded_degree) + qpts, qwts = Q.get_points(), Q.get_weights() + phis = expansion_set._tabulate(embedded_degree, qpts, 0)[(0,)*sd] + PhiW = phis * qwts[None, :] + M = np.tensordot(PhiW, phis, (-1, -1)) + MinvPhiW = np.linalg.solve(M, PhiW) + coeffs = [np.tensordot(e.tabulate(0, qpts)[(0,)*sd], MinvPhiW, (-1, -1)) for e in elements] + coeffs = np.concatenate(coeffs, axis=0) + poly_set = PolynomialSet(ref_el, degree, embedded_degree, From 2c46a2682918688403fe1e30b9d5111385adf6ed Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Feb 2026 19:12:20 +0000 Subject: [PATCH 2/5] test --- FIAT/nodal_enriched.py | 18 +++++++----- test/FIAT/unit/test_nodal_enriched.py | 40 +++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 7 deletions(-) create mode 100644 test/FIAT/unit/test_nodal_enriched.py diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index 75ff407cf..afd312dce 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -45,8 +45,10 @@ def __init__(self, *elements): order = max(e.get_order() for e in elements) formdegree = None if any(e.get_formdegree() is None for e in elements) \ else max(e.get_formdegree() for e in elements) + + # Grab the ExpansionSet defined on the maximal complex # LagrangeExpansionSet has fixed degree, ensure we grab the highest one - elem = elements[embedded_degrees.index(embedded_degree)] + elem = max(elements, key=lambda e: (e.get_reference_complex(), e.degree())) ref_el = elem.get_reference_complex() expansion_set = elem.get_nodal_basis().get_expansion_set() mapping = elem.mapping()[0] @@ -54,7 +56,7 @@ def __init__(self, *elements): sd = ref_el.get_spatial_dimension() # Sanity check - assert all(e.get_reference_complex() == ref_el for e in elements) + assert all(e.get_reference_complex() <= ref_el for e in elements) assert all(set(e.mapping()) == {mapping, } for e in elements) assert all(e.value_shape() == value_shape for e in elements) @@ -70,13 +72,15 @@ def __init__(self, *elements): else: # Obtain coefficients via projection Q = create_quadrature(ref_el, 2*embedded_degree) - qpts, qwts = Q.get_points(), Q.get_weights() + qpts = Q.get_points() phis = expansion_set._tabulate(embedded_degree, qpts, 0)[(0,)*sd] - PhiW = phis * qwts[None, :] - M = np.tensordot(PhiW, phis, (-1, -1)) + PhiW = phis * Q.get_weights() + M = np.tensordot(phis, PhiW, (-1, -1)) MinvPhiW = np.linalg.solve(M, PhiW) - coeffs = [np.tensordot(e.tabulate(0, qpts)[(0,)*sd], MinvPhiW, (-1, -1)) for e in elements] - coeffs = np.concatenate(coeffs, axis=0) + + tabulations = np.concatenate([e.tabulate(0, qpts)[(0,)*sd] for e in elements], axis=0) + coeffs = np.tensordot(tabulations, MinvPhiW, (-1, -1)) + assert coeffs.shape[1:-1] == value_shape poly_set = PolynomialSet(ref_el, degree, diff --git a/test/FIAT/unit/test_nodal_enriched.py b/test/FIAT/unit/test_nodal_enriched.py new file mode 100644 index 000000000..f2c7f94ac --- /dev/null +++ b/test/FIAT/unit/test_nodal_enriched.py @@ -0,0 +1,40 @@ +import numpy +from FIAT import (NodalEnrichedElement, RestrictedElement, + BernardiRaugel, GuzmanNeilanFirstKindH1, + ufc_simplex) + + +def test_nodal_enriched_mismatching_expansion_set(): + # Add two elements with mismatching expansion sets + sd = 2 + ref_el = ufc_simplex(sd) + + BR = BernardiRaugel(ref_el, 1) + P1 = RestrictedElement(BR, restriction_domain="vertex", take_closure=False) + + GN = GuzmanNeilanFirstKindH1(ref_el, 1) + GNB = RestrictedElement(GN, restriction_domain="facet", take_closure=False) + + fe = NodalEnrichedElement(P1, GNB) + coeffs = fe.poly_set.get_coeffs() + e = numpy.tensordot(GN.dual.to_riesz(fe.poly_set), + coeffs, axes=(range(1, coeffs.ndim),)*2) + e[abs(e) < 1E-12] = 0 + print(e) + + degree = GN.degree() + ref_complex = GN.get_reference_complex() + top = ref_complex.topology + pts = [] + for dim in top: + for entity in top[dim]: + pts.extend(ref_complex.make_points(dim, entity, degree)) + + result = fe.tabulate(0, pts)[(0,)*sd] + expected = GN.tabulate(0, pts)[(0,)*sd] + + result -= expected + expected = 0 + result[abs(result) < 1E-12] = 0 + result = result.reshape((result.shape[0], -1)) + assert numpy.allclose(result, expected) From 953bed5f368ec5a1b9b4fe9526a34f065ee9d603 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Feb 2026 12:29:13 +0000 Subject: [PATCH 3/5] Update test/FIAT/unit/test_nodal_enriched.py --- test/FIAT/unit/test_nodal_enriched.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/FIAT/unit/test_nodal_enriched.py b/test/FIAT/unit/test_nodal_enriched.py index f2c7f94ac..a6c289141 100644 --- a/test/FIAT/unit/test_nodal_enriched.py +++ b/test/FIAT/unit/test_nodal_enriched.py @@ -9,7 +9,7 @@ def test_nodal_enriched_mismatching_expansion_set(): sd = 2 ref_el = ufc_simplex(sd) - BR = BernardiRaugel(ref_el, 1) + BR = BernardiRaugel(ref_el, 1, hierarchical=True) P1 = RestrictedElement(BR, restriction_domain="vertex", take_closure=False) GN = GuzmanNeilanFirstKindH1(ref_el, 1) From 4b16836aca1c84d1752ce965613002f00604b169 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Feb 2026 13:49:29 +0000 Subject: [PATCH 4/5] cleanup --- FIAT/nodal_enriched.py | 16 +++++++--------- test/FIAT/unit/test_nodal_enriched.py | 20 ++++++++++---------- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index afd312dce..51973ed4e 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -39,21 +39,18 @@ def __init__(self, *elements): "of NodalEnrichedElement are nodal") # Extract common data - embedded_degrees = [e.get_nodal_basis().get_embedded_degree() for e in elements] + embedded_degrees = [e.degree() for e in elements] embedded_degree = max(embedded_degrees) - degree = max(e.degree() for e in elements) order = max(e.get_order() for e in elements) formdegree = None if any(e.get_formdegree() is None for e in elements) \ else max(e.get_formdegree() for e in elements) - # Grab the ExpansionSet defined on the maximal complex - # LagrangeExpansionSet has fixed degree, ensure we grab the highest one + # Grab the ExpansionSet defined on the maximal complex with the highest degree elem = max(elements, key=lambda e: (e.get_reference_complex(), e.degree())) ref_el = elem.get_reference_complex() expansion_set = elem.get_nodal_basis().get_expansion_set() mapping = elem.mapping()[0] value_shape = elem.value_shape() - sd = ref_el.get_spatial_dimension() # Sanity check assert all(e.get_reference_complex() <= ref_el for e in elements) @@ -61,7 +58,7 @@ def __init__(self, *elements): assert all(e.value_shape() == value_shape for e in elements) # Merge polynomial sets - if isinstance(expansion_set, LagrangeLineExpansionSet): + if isinstance(expansion_set, LagrangeLineExpansionSet) and expansion_set.degree == embedded_degree: # Obtain coefficients via interpolation points = expansion_set.get_points() coeffs = np.vstack([e.tabulate(0, points)[(0,)] for e in elements]) @@ -71,10 +68,11 @@ def __init__(self, *elements): coeffs = _merge_coeffs(coeffs, ref_el, embedded_degrees, expansion_set.continuity) else: # Obtain coefficients via projection + sd = ref_el.get_spatial_dimension() Q = create_quadrature(ref_el, 2*embedded_degree) qpts = Q.get_points() phis = expansion_set._tabulate(embedded_degree, qpts, 0)[(0,)*sd] - PhiW = phis * Q.get_weights() + PhiW = np.multiply(phis, Q.get_weights()) M = np.tensordot(phis, PhiW, (-1, -1)) MinvPhiW = np.linalg.solve(M, PhiW) @@ -83,7 +81,7 @@ def __init__(self, *elements): assert coeffs.shape[1:-1] == value_shape poly_set = PolynomialSet(ref_el, - degree, + embedded_degree, embedded_degree, expansion_set, coeffs) @@ -151,5 +149,5 @@ def _merge_entity_ids(entity_ids, offsets): for entity in ids[dim]: if entity not in ret[dim]: ret[dim][entity] = [] - ret[dim][entity].extend(np.array(ids[dim][entity]) + offsets[i]) + ret[dim][entity].extend(offsets[i] + dof for dof in ids[dim][entity]) return ret diff --git a/test/FIAT/unit/test_nodal_enriched.py b/test/FIAT/unit/test_nodal_enriched.py index a6c289141..40af1f5a3 100644 --- a/test/FIAT/unit/test_nodal_enriched.py +++ b/test/FIAT/unit/test_nodal_enriched.py @@ -5,23 +5,28 @@ def test_nodal_enriched_mismatching_expansion_set(): - # Add two elements with mismatching expansion sets sd = 2 ref_el = ufc_simplex(sd) + # Extract (non-macro) vector P1 from BR BR = BernardiRaugel(ref_el, 1, hierarchical=True) P1 = RestrictedElement(BR, restriction_domain="vertex", take_closure=False) + # Extract the macro face bubbles from GN GN = GuzmanNeilanFirstKindH1(ref_el, 1) - GNB = RestrictedElement(GN, restriction_domain="facet", take_closure=False) + MFB = RestrictedElement(GN, restriction_domain="facet", take_closure=False) + + # Add two elements with mismatching expansion sets + # The resulting element should be identical to GN + fe = NodalEnrichedElement(P1, MFB) - fe = NodalEnrichedElement(P1, GNB) + # Test nodality coeffs = fe.poly_set.get_coeffs() e = numpy.tensordot(GN.dual.to_riesz(fe.poly_set), coeffs, axes=(range(1, coeffs.ndim),)*2) - e[abs(e) < 1E-12] = 0 - print(e) + assert numpy.allclose(e, numpy.eye(*e.shape)) + # Test that the spaces are equal degree = GN.degree() ref_complex = GN.get_reference_complex() top = ref_complex.topology @@ -32,9 +37,4 @@ def test_nodal_enriched_mismatching_expansion_set(): result = fe.tabulate(0, pts)[(0,)*sd] expected = GN.tabulate(0, pts)[(0,)*sd] - - result -= expected - expected = 0 - result[abs(result) < 1E-12] = 0 - result = result.reshape((result.shape[0], -1)) assert numpy.allclose(result, expected) From 4d1379b5a07b1cae86f28ced3bc518fc6b77e046 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Feb 2026 21:59:23 +0000 Subject: [PATCH 5/5] Update test/FIAT/unit/test_nodal_enriched.py --- test/FIAT/unit/test_nodal_enriched.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/FIAT/unit/test_nodal_enriched.py b/test/FIAT/unit/test_nodal_enriched.py index 40af1f5a3..bf4f449b8 100644 --- a/test/FIAT/unit/test_nodal_enriched.py +++ b/test/FIAT/unit/test_nodal_enriched.py @@ -22,9 +22,9 @@ def test_nodal_enriched_mismatching_expansion_set(): # Test nodality coeffs = fe.poly_set.get_coeffs() - e = numpy.tensordot(GN.dual.to_riesz(fe.poly_set), + V = numpy.tensordot(GN.dual.to_riesz(fe.poly_set), coeffs, axes=(range(1, coeffs.ndim),)*2) - assert numpy.allclose(e, numpy.eye(*e.shape)) + assert numpy.allclose(V, numpy.eye(*V.shape)) # Test that the spaces are equal degree = GN.degree()