From 3634bc59d93a3b038ef4dd3215778526a66ea3d8 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 4 Sep 2025 17:15:35 +0100 Subject: [PATCH 01/11] Define P1 with trivial Vandermonde matrix --- FIAT/expansions.py | 42 +++++++++++++++++++++++------------------- FIAT/lagrange.py | 3 +-- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index b7d846f26..e25f1c2f1 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -81,7 +81,6 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): raise ValueError("Higher order derivatives not supported") if variant not in [None, "bubble", "dual"]: raise ValueError(f"Invalid variant {variant}") - if variant == "bubble": scale = -scale @@ -94,9 +93,10 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): pad_dim = dim + 2 dX = pad_jacobian(Jinv, pad_dim) - phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), scale) + + phi[0] = sum((ref_pts[i] * 0 for i in range(dim)), scale) if dphi is not None: - dphi[0] = (phi[0] - phi[0]) * dX[0] + dphi[0] = (phi[0] * 0) * dX[0] if ddphi is not None: ddphi[0] = outer(dphi[0], dX[0]) if dim == 0 or n == 0: @@ -127,29 +127,32 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): a = 0.5 * (alpha + beta) + 1.0 b = 0.5 * (alpha - beta) - factor = a * fa - b * fb - phi[inext] = factor * phi[icur] + fcur = a * fa - b * fb + phi[inext] = fcur * phi[icur] if dphi is not None: - dfactor = a * dfa - b * dfb - dphi[inext] = factor * dphi[icur] + phi[icur] * dfactor + dfcur = a * dfa - b * dfb + dphi[inext] = fcur * dphi[icur] + phi[icur] * dfcur if ddphi is not None: - ddphi[inext] = factor * ddphi[icur] + sym_outer(dphi[icur], dfactor) + ddphi[inext] = fcur * ddphi[icur] + sym_outer(dphi[icur], dfcur) # general i by recurrence for i in range(1, n - sum(sub_index)): iprev, icur, inext = icur, inext, idx(*sub_index, i + 1) a, b, c = coefficients(alpha, beta, i) - factor = a * fa - b * fb - phi[inext] = factor * phi[icur] - c * (fc * phi[iprev]) + fcur = a * fa - b * fb + fprev = -c * fc + phi[inext] = fcur * phi[icur] + fprev * phi[iprev] if dphi is None: continue - dfactor = a * dfa - b * dfb - dphi[inext] = (factor * dphi[icur] + phi[icur] * dfactor - - c * (fc * dphi[iprev] + phi[iprev] * dfc)) + dfcur = a * dfa - b * dfb + dfprev = -c * dfc + dphi[inext] = (fcur * dphi[icur] + phi[icur] * dfcur + + fprev * dphi[iprev] + phi[iprev] * dfprev) if ddphi is None: continue - ddphi[inext] = (factor * ddphi[icur] + sym_outer(dphi[icur], dfactor) - - c * (fc * ddphi[iprev] + sym_outer(dphi[iprev], dfc) + phi[iprev] * ddfc)) + ddfcur = -c * ddfc + ddphi[inext] = (fcur * ddphi[icur] + sym_outer(dphi[icur], dfcur) + + fprev * ddphi[iprev] + sym_outer(dphi[iprev], dfcur) + phi[iprev] * ddfcur) # normalize d = codim + 1 @@ -159,9 +162,10 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): if variant is not None: p = index[-1] + shift alpha = 2 * (sum(index[:-1]) + d * shift) - 1 - norm2 = (0.5 + d) / d if p > 0 and p + alpha > 0: - norm2 *= (p + alpha) * (2*p + alpha) / p + norm2 = (2*d+1) * (p + alpha) * (2*p + alpha) / (2*d*p) + else: + norm2 = 1.0 else: norm2 = (2*sum(index) + d) / d scale = math.sqrt(norm2) @@ -183,9 +187,9 @@ def C0_basis(dim, n, tabulations): # Recover facet bubbles for phi in tabulations: icur = 0 - phi[icur] *= -1 for inext in range(1, dim+1): - phi[icur] -= phi[inext] + phi[icur] += phi[inext] + phi[icur] *= -1 if dim == 2: for i in range(2, n+1): phi[idx(0, i)] -= phi[idx(1, i-1)] diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index eead0cfa6..cb605ec89 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -83,7 +83,6 @@ def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False): points = get_lagrange_points(dual) poly_set = LagrangePolynomialSet(ref_el, points) else: - poly_variant = "bubble" if ref_el.is_macrocell() else None - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant=poly_variant) + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble", scale=1) formdegree = 0 # 0-form super().__init__(poly_set, dual, degree, formdegree) From 5db64307bc7c5e619353fd485b080dd9bc88f590 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 21 Oct 2025 12:30:53 +0100 Subject: [PATCH 02/11] fixup --- FIAT/expansions.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index e25f1c2f1..142b64854 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -93,8 +93,8 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): pad_dim = dim + 2 dX = pad_jacobian(Jinv, pad_dim) - phi[0] = sum((ref_pts[i] * 0 for i in range(dim)), scale) + if dphi is not None: dphi[0] = (phi[0] * 0) * dX[0] if ddphi is not None: @@ -150,9 +150,9 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): fprev * dphi[iprev] + phi[iprev] * dfprev) if ddphi is None: continue - ddfcur = -c * ddfc + ddfprev = -c * ddfc ddphi[inext] = (fcur * ddphi[icur] + sym_outer(dphi[icur], dfcur) + - fprev * ddphi[iprev] + sym_outer(dphi[iprev], dfcur) + phi[iprev] * ddfcur) + fprev * ddphi[iprev] + sym_outer(dphi[iprev], dfprev) + phi[iprev] * ddfprev) # normalize d = codim + 1 @@ -162,10 +162,9 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): if variant is not None: p = index[-1] + shift alpha = 2 * (sum(index[:-1]) + d * shift) - 1 + norm2 = (0.5 + d) / d if p > 0 and p + alpha > 0: - norm2 = (2*d+1) * (p + alpha) * (2*p + alpha) / (2*d*p) - else: - norm2 = 1.0 + norm2 *= (p + alpha) * (2*p + alpha) / p else: norm2 = (2*sum(index) + d) / d scale = math.sqrt(norm2) From be2f46be372fe3945f69d566e07f8a4fa30d5b44 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 21 Oct 2025 13:16:28 +0100 Subject: [PATCH 03/11] in-place summation --- FIAT/expansions.py | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 142b64854..df25ca809 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -89,7 +89,6 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): phi, dphi, ddphi = results + (None,) * (2-order) outer = lambda x, y: x[:, None, ...] * y[None, ...] - sym_outer = lambda x, y: outer(x, y) + outer(y, x) pad_dim = dim + 2 dX = pad_jacobian(Jinv, pad_dim) @@ -131,28 +130,42 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): phi[inext] = fcur * phi[icur] if dphi is not None: dfcur = a * dfa - b * dfb - dphi[inext] = fcur * dphi[icur] + phi[icur] * dfcur + dphi[inext] = fcur * dphi[icur] + dphi[inext] += phi[icur] * dfcur if ddphi is not None: - ddphi[inext] = fcur * ddphi[icur] + sym_outer(dphi[icur], dfcur) + ddphi[inext] = fcur * ddphi[icur] + ddphi[inext] += outer(dfcur, dphi[icur]) + ddphi[inext] += outer(dphi[icur], dfcur) # general i by recurrence for i in range(1, n - sum(sub_index)): iprev, icur, inext = icur, inext, idx(*sub_index, i + 1) a, b, c = coefficients(alpha, beta, i) + fcur = a * fa - b * fb fprev = -c * fc - phi[inext] = fcur * phi[icur] + fprev * phi[iprev] + phi[inext] = fcur * phi[icur] + phi[inext] += fprev * phi[iprev] if dphi is None: continue + dfcur = a * dfa - b * dfb dfprev = -c * dfc - dphi[inext] = (fcur * dphi[icur] + phi[icur] * dfcur + - fprev * dphi[iprev] + phi[iprev] * dfprev) + dphi[inext] = fcur * dphi[icur] + dphi[inext] += fprev * dphi[iprev] + dphi[inext] += phi[icur] * dfcur + dphi[inext] += phi[iprev] * dfprev if ddphi is None: continue + ddfprev = -c * ddfc - ddphi[inext] = (fcur * ddphi[icur] + sym_outer(dphi[icur], dfcur) + - fprev * ddphi[iprev] + sym_outer(dphi[iprev], dfprev) + phi[iprev] * ddfprev) + ddphi[inext] = fcur * ddphi[icur] + ddphi[inext] += fprev * ddphi[iprev] + ddphi[inext] += outer(dfcur, dphi[icur]) + ddphi[inext] += outer(dphi[icur], dfcur) + ddphi[inext] += outer(dfprev, dphi[iprev]) + ddphi[inext] += outer(dphi[iprev], dfprev) + ddphi[inext] += phi[iprev] * ddfprev # normalize d = codim + 1 From 2412e7547718b10bdc3c0f8ae2ceb32f23f37004 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 21 Oct 2025 13:25:46 +0100 Subject: [PATCH 04/11] Fix BDFM --- FIAT/brezzi_douglas_fortin_marini.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/brezzi_douglas_fortin_marini.py b/FIAT/brezzi_douglas_fortin_marini.py index 4361e70f6..aeffc679b 100644 --- a/FIAT/brezzi_douglas_fortin_marini.py +++ b/FIAT/brezzi_douglas_fortin_marini.py @@ -70,7 +70,7 @@ def BDFMSpace(ref_el, order): # Linear vector valued space. Since the embedding degree of this element # is 2, this is implemented by taking the quadratic space and selecting # the linear polynomials. - vec_poly_set = polynomial_set.ONPolynomialSet(ref_el, order, (sd,)) + vec_poly_set = polynomial_set.ONPolynomialSet(ref_el, order, (sd,), variant="bubble") # Linears are the first three polynomials in each dimension. vec_poly_set = vec_poly_set.take([0, 1, 2, 6, 7, 8]) From 72f5c865659047ef64706cbef22703d73864b2fb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 15:10:38 +0000 Subject: [PATCH 05/11] Apply suggestion from @pbrubeck --- FIAT/expansions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 52df45cda..7df9c4de3 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -200,7 +200,7 @@ def C0_basis(dim, n, tabulations): icur = 0 phi[icur] *= -1.0 for inext in range(1, dim+1): - phi[icur] += phi[inext] + phi[icur] -= phi[inext] if dim == 2: for i in range(2, n+1): phi[idx(0, i)] -= phi[idx(1, i-1)] From b37c52d371203c004432a4e701018a64476848cf Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 15:27:09 +0000 Subject: [PATCH 06/11] fix WuXu --- FIAT/wuxu.py | 70 +++++++++++++++++++--------------------------------- 1 file changed, 25 insertions(+), 45 deletions(-) diff --git a/FIAT/wuxu.py b/FIAT/wuxu.py index 37f187489..b703cd359 100644 --- a/FIAT/wuxu.py +++ b/FIAT/wuxu.py @@ -28,27 +28,30 @@ polydim = expansions.polynomial_dimension -def WuXuRobustH3NCSpace(ref_el): - """Constructs a basis for the the Wu Xu H^3 nonconforming space - P^{(3,2)}(T) = P_3(T) + b_T P_1(T) + b_T^2 P_1(T), +def WuXuH3NCSpace(ref_el, robust=False): + """Constructs a basis for the the Wu Xu H^3 nonconforming spaces + + P^{(3,1)}(T) = P_3(T) + b_T P_1(T), if robust = False + + P^{(3,2)}(T) = P_3(T) + b_T P_1(T) + b_T^2 P_1(T), if robust = True + where b_T is the standard cubic bubble.""" sd = ref_el.get_spatial_dimension() assert sd == 2 - em_deg = 7 - # Unfortunately, b_T^2 P_1 has degree 7 (cubic squared times a linear) # so we need a high embedded degree! - p7 = polynomial_set.ONPolynomialSet(ref_el, 7) + embedded_degree = 7 if robust else 4 + pk = polynomial_set.ONPolynomialSet(ref_el, embedded_degree) dimp1 = polydim(ref_el, 1) dimp3 = polydim(ref_el, 3) - dimp7 = polydim(ref_el, 7) + dimpk = polydim(ref_el, embedded_degree) # Here's the first bit we'll work with. It's already expressed in terms # of the ON basis for P7, so we're golden. - p3fromp7 = p7.take(list(range(dimp3))) + p3frompk = pk.take(list(range(dimp3))) # Rather than creating the barycentric coordinates ourself, let's # reuse the existing bubble functionality @@ -56,56 +59,33 @@ def WuXuRobustH3NCSpace(ref_el): p1 = Lagrange(ref_el, 1) # next, we'll have to project b_T P1 and b_T^2 P1 onto P^7 - Q = create_quadrature(ref_el, 14) + Q = create_quadrature(ref_el, 2*embedded_degree) Qpts = numpy.array(Q.get_points()) Qwts = numpy.array(Q.get_weights()) - zero_index = tuple([0 for i in range(sd)]) - # it's just one bubble function: let's get a 1d array! - bT_at_qpts = bT.tabulate(0, Qpts)[zero_index][0, :] - p1_at_qpts = p1.tabulate(0, Qpts)[zero_index] + bT_at_qpts = bT.tabulate(0, Qpts)[(0,)*sd][0, :] + p1_at_qpts = p1.tabulate(0, Qpts)[(0,)*sd] - # Note: difference in signature because bT, p1 are FE and p7 is a + # Note: difference in signature because bT, p1 are FE and pk is a # polynomial set - p7_at_qpts = p7.tabulate(Qpts)[zero_index] + pk_at_qpts = pk.tabulate(Qpts)[(0,)*sd] - bubble_coeffs = numpy.zeros((6, dimp7), "d") + bubble_coeffs = numpy.zeros((6, dimpk), "d") # first three: bT P1, last three will be bT^2 P1 foo = bT_at_qpts * p1_at_qpts * Qwts - bubble_coeffs[:dimp1, :] = numpy.dot(foo, p7_at_qpts.T) + bubble_coeffs[:dimp1, :] = numpy.dot(foo, pk_at_qpts.T) - foo = bT_at_qpts * foo - bubble_coeffs[dimp1:2*dimp1, :] = numpy.dot(foo, p7_at_qpts.T) + if robust: + foo = bT_at_qpts * foo + bubble_coeffs[dimp1:2*dimp1, :] = numpy.dot(foo, pk_at_qpts.T) - bubbles = polynomial_set.PolynomialSet(ref_el, 3, em_deg, - p7.get_expansion_set(), + bubbles = polynomial_set.PolynomialSet(ref_el, 3, embedded_degree, + pk.get_expansion_set(), bubble_coeffs) - return polynomial_set.polynomial_set_union_normalized(p3fromp7, bubbles) - - -def WuXuH3NCSpace(ref_el): - """Constructs a basis for the the Wu Xu H^3 nonconforming space - P(T) = P_3(T) + b_T P_1(T), - where b_T is the standard cubic bubble.""" - - assert ref_el.get_spatial_dimension() == 2 - - em_deg = 4 - p4 = polynomial_set.ONPolynomialSet(ref_el, em_deg) - - # Here's the first bit we'll work with. It's already expressed in terms - # of the ON basis for P4, so we're golden. - dimp3 = polydim(ref_el, 3) - p3fromp4 = p4.take(list(range(dimp3))) - - # Rather than creating the barycentric coordinates ourself, let's - # reuse the existing bubble functionality - bT = Bubble(ref_el, 4) - - return polynomial_set.polynomial_set_union_normalized(p3fromp4, bT.get_nodal_basis()) + return polynomial_set.polynomial_set_union_normalized(p3frompk, bubbles) class WuXuRobustH3NCDualSet(dual_set.DualSet): @@ -176,7 +156,7 @@ def __init__(self, ref_el, degree): class WuXuRobustH3NC(finite_element.CiarletElement): """The Wu-Xu robust H3 nonconforming finite element""" def __init__(self, ref_el, degree=7): - poly_set = WuXuRobustH3NCSpace(ref_el) + poly_set = WuXuH3NCSpace(ref_el, robust=True) assert degree == poly_set.degree dual = WuXuRobustH3NCDualSet(ref_el, degree) super().__init__(poly_set, dual, degree) From 0ac7d75ac401e6cf9398b289e9f2759ddfa73bb8 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 15:36:07 +0000 Subject: [PATCH 07/11] Also CR --- FIAT/crouzeix_raviart.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index 1fd990cc9..a786c9a9b 100644 --- a/FIAT/crouzeix_raviart.py +++ b/FIAT/crouzeix_raviart.py @@ -89,7 +89,7 @@ def __init__(self, ref_el, degree, variant=None, quad_scheme=None): base_element = type(self)(ref_el.get_parent(), degree) poly_set = macro.MacroPolynomialSet(ref_el, base_element) else: - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble", scale=1) dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg, quad_scheme) super().__init__(poly_set, dual, degree) From 82ef1d322f63cb36e442690de5afdc898eaf754e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 16:00:14 +0000 Subject: [PATCH 08/11] Fix NodalEnrichment across different ExpansionSets --- FIAT/crouzeix_raviart.py | 2 +- FIAT/nodal_enriched.py | 18 +++++++++++++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index a786c9a9b..1fd990cc9 100644 --- a/FIAT/crouzeix_raviart.py +++ b/FIAT/crouzeix_raviart.py @@ -89,7 +89,7 @@ def __init__(self, ref_el, degree, variant=None, quad_scheme=None): base_element = type(self)(ref_el.get_parent(), degree) poly_set = macro.MacroPolynomialSet(ref_el, base_element) else: - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble", scale=1) + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg, quad_scheme) super().__init__(poly_set, dual, degree) 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 2c295a57e655727f751173f6b1625142ae9f179e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 14 Feb 2026 08:47:33 +0000 Subject: [PATCH 09/11] fix --- FIAT/expansions.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 7df9c4de3..1067bf446 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -151,17 +151,17 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): dfcur = a * dfa - b * dfb dfprev = -c * dfc dphi[inext] = fcur * dphi[icur] - dphi[inext] += fprev * dphi[iprev] dphi[inext] += phi[icur] * dfcur + dphi[inext] += fprev * dphi[iprev] dphi[inext] += phi[iprev] * dfprev if ddphi is None: continue ddfprev = -c * ddfc ddphi[inext] = fcur * ddphi[icur] - ddphi[inext] += fprev * ddphi[iprev] ddphi[inext] += outer(dfcur, dphi[icur]) ddphi[inext] += outer(dphi[icur], dfcur) + ddphi[inext] += fprev * ddphi[iprev] ddphi[inext] += outer(dfprev, dphi[iprev]) ddphi[inext] += outer(dphi[iprev], dfprev) ddphi[inext] += phi[iprev] * ddfprev From 3686c6e7089fd305761303258fbb902caeb8b2dd Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 14 Feb 2026 12:44:46 +0000 Subject: [PATCH 10/11] optimise subtraction --- FIAT/expansions.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 1067bf446..132a9dc54 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -243,7 +243,7 @@ def C0_basis(dim, n, tabulations): def xi_triangle(eta): """Maps from [-1,1]^2 to the (-1,1) reference triangle.""" eta1, eta2 = eta - xi1 = 0.5 * (1.0 + eta1) * (1.0 - eta2) - 1.0 + xi1 = -0.5 * (eta1 + 1.0) * (eta2 - 1.0) - 1.0 xi2 = eta2 return (xi1, xi2) @@ -251,9 +251,9 @@ def xi_triangle(eta): def xi_tetrahedron(eta): """Maps from [-1,1]^3 to the -1/1 reference tetrahedron.""" eta1, eta2, eta3 = eta - xi1 = 0.25 * (1. + eta1) * (1. - eta2) * (1. - eta3) - 1. - xi2 = 0.5 * (1. + eta2) * (1. - eta3) - 1. xi3 = eta3 + xi2 = -0.5 * (eta2 + 1.0) * (eta3 - 1.0) - 1.0 + xi1 = 0.25 * (eta1 + 1.0) * (eta2 - 1.0) * (eta3 - 1.0) - 1.0 return xi1, xi2, xi3 From 1ebf6b24a31af3894ee0fd520e82c99d8f2ded2f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Feb 2026 09:57:06 +0000 Subject: [PATCH 11/11] Stable order of operations --- FIAT/expansions.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 4ab246430..3049d44f0 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -96,7 +96,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): for k in range(order+1)] phi, dphi, ddphi = results + [None] * (2-order) - phi[0] += scale + phi[0] = scale if dim == 0 or n == 0: return results if dim > 3 or dim < 0: @@ -129,12 +129,12 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): phi[inext] = fcur * phi[icur] if dphi is not None: dfcur = a * dfa - b * dfb - dphi[inext] = fcur * dphi[icur] - dphi[inext] += phi[icur] * dfcur + dphi[inext] = phi[icur] * dfcur + dphi[inext] += fcur * dphi[icur] if ddphi is not None: - ddphi[inext] = fcur * ddphi[icur] + ddphi[inext] = outer(dphi[icur], dfcur) ddphi[inext] += outer(dfcur, dphi[icur]) - ddphi[inext] += outer(dphi[icur], dfcur) + ddphi[inext] += fcur * ddphi[icur] # general i by recurrence for i in range(1, n - sum(sub_index)): @@ -150,21 +150,21 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None): dfcur = a * dfa - b * dfb dfprev = -c * dfc - dphi[inext] = fcur * dphi[icur] - dphi[inext] += phi[icur] * dfcur - dphi[inext] += fprev * dphi[iprev] + dphi[inext] = phi[icur] * dfcur dphi[inext] += phi[iprev] * dfprev + dphi[inext] += fcur * dphi[icur] + dphi[inext] += fprev * dphi[iprev] if ddphi is None: continue ddfprev = -c * ddfc - ddphi[inext] = fcur * ddphi[icur] - ddphi[inext] += outer(dfcur, dphi[icur]) + ddphi[inext] = phi[iprev] * ddfprev ddphi[inext] += outer(dphi[icur], dfcur) - ddphi[inext] += fprev * ddphi[iprev] - ddphi[inext] += outer(dfprev, dphi[iprev]) + ddphi[inext] += outer(dfcur, dphi[icur]) ddphi[inext] += outer(dphi[iprev], dfprev) - ddphi[inext] += phi[iprev] * ddfprev + ddphi[inext] += outer(dfprev, dphi[iprev]) + ddphi[inext] += fcur * ddphi[icur] + ddphi[inext] += fprev * ddphi[iprev] # normalize d = codim + 1 @@ -243,7 +243,7 @@ def C0_basis(dim, n, tabulations): def xi_triangle(eta): """Maps from [-1,1]^2 to the (-1,1) reference triangle.""" eta1, eta2 = eta - xi1 = -0.5 * (eta1 + 1.0) * (eta2 - 1.0) - 1.0 + xi1 = 0.5 * (1.0 + eta1) * (1.0 - eta2) - 1.0 xi2 = eta2 return (xi1, xi2) @@ -251,9 +251,9 @@ def xi_triangle(eta): def xi_tetrahedron(eta): """Maps from [-1,1]^3 to the -1/1 reference tetrahedron.""" eta1, eta2, eta3 = eta + xi1 = 0.25 * (1. + eta1) * (1. - eta2) * (1. - eta3) - 1. + xi2 = 0.5 * (1. + eta2) * (1. - eta3) - 1. xi3 = eta3 - xi2 = -0.5 * (eta2 + 1.0) * (eta3 - 1.0) - 1.0 - xi1 = 0.25 * (eta1 + 1.0) * (eta2 - 1.0) * (eta3 - 1.0) - 1.0 return xi1, xi2, xi3