From 3b02c0592845c59c10840e5f75f7d434f24b0a38 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 6 Feb 2026 18:20:39 +0000 Subject: [PATCH 01/21] test GN trace --- finat/mtw.py | 13 ++++++-- test/FIAT/unit/test_stokes_complex.py | 45 +++++++++++++++++++++++++-- 2 files changed, 53 insertions(+), 5 deletions(-) diff --git a/finat/mtw.py b/finat/mtw.py index 9e4ebaa2..542ae38c 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -16,7 +16,11 @@ def __init__(self, cell, degree=3): reduced_dofs = deepcopy(self._element.entity_dofs()) sd = cell.get_spatial_dimension() - fdofs = sd + 1 + + dim_P1 = sd + dim_Ned1 = (sd*(sd-1))//2 + fdofs = dim_P1 + dim_Ned1 + reduced_dofs[sd][0] = [] for f in reduced_dofs[sd-1]: reduced_dofs[sd-1][f] = reduced_dofs[sd-1][f][:fdofs] @@ -29,13 +33,18 @@ def basis_transformation(self, coordinate_mapping): V = identity(numbf, ndof) sd = self.cell.get_spatial_dimension() + if sd == 2: + facet_transform = normal_tangential_edge_transform + else: + raise NotImplementedError + bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) detJ = coordinate_mapping.detJ_at(bary) entity_dofs = self.entity_dofs() for f in sorted(entity_dofs[sd-1]): cur = entity_dofs[sd-1][f][0] - V[cur+1, cur:cur+sd] = normal_tangential_edge_transform(self.cell, J, detJ, f) + V[cur+1, cur:cur+sd] = facet_transform(self.cell, J, detJ, f) return ListTensor(V.T) diff --git a/test/FIAT/unit/test_stokes_complex.py b/test/FIAT/unit/test_stokes_complex.py index 4e973588..a7725ed7 100644 --- a/test/FIAT/unit/test_stokes_complex.py +++ b/test/FIAT/unit/test_stokes_complex.py @@ -14,6 +14,7 @@ from FIAT.macro import CkPolynomialSet from FIAT.alfeld_sorokina import AlfeldSorokinaSpace from FIAT.arnold_qin import ArnoldQinSpace +from FIAT.bernardi_raugel import BernardiRaugel from FIAT.christiansen_hu import ChristiansenHuSpace from FIAT.guzman_neilan import ( GuzmanNeilanFirstKindH1, @@ -21,12 +22,16 @@ GuzmanNeilanH1div, GuzmanNeilanSpace) from FIAT.restricted import RestrictedElement - +from FIAT.quadrature import FacetQuadratureRule T = ufc_simplex(2) S = ufc_simplex(3) +def inner(u, v, qwts): + return numpy.tensordot(numpy.multiply(u, qwts), v, axes=(tuple(range(1, u.ndim)), )*2) + + def rHCT(cell): return RestrictedElement(HCT(cell, reduced=True), restriction_domain="vertex") @@ -182,12 +187,13 @@ def test_gn_stokes_pairs(cell, kind): check_stokes_complex(spaces, degree) +@pytest.mark.parametrize("element", (GuzmanNeilanFirstKindH1, BernardiRaugel)) @pytest.mark.parametrize("quad_scheme", (None, "KMV(2),powell-sabin")) @pytest.mark.parametrize("sd", (2, 3)) -def test_gn_dofs(sd, quad_scheme): +def test_gn_dofs(element, sd, quad_scheme): cell = symmetric_simplex(sd) - fe = GuzmanNeilanFirstKindH1(cell, 1, quad_scheme=quad_scheme) + fe = element(cell, 1, quad_scheme=quad_scheme) degree = fe.degree() assert degree == sd ref_complex = fe.get_reference_complex() @@ -206,6 +212,39 @@ def test_gn_dofs(sd, quad_scheme): assert numpy.allclose(div_moments, expected) +@pytest.mark.parametrize("sd", (2, 3)) +def test_gn_trace(sd): + degree = 1 + cell = symmetric_simplex(sd) + + gn = GuzmanNeilanFirstKindH1(cell, degree) + br = BernardiRaugel(cell, degree) + + sd = cell.get_spatial_dimension() + ref_face = cell.construct_subelement(sd-1) + Q_face = create_quadrature(ref_face, 2*degree) + + Phis = ONPolynomialSet(ref_face, degree) + phis_at_qpts = Phis.tabulate(Q_face.get_points())[(0,)*(sd-1)] + + for f in cell.topology[sd-1]: + Q = FacetQuadratureRule(cell, sd-1, f, Q_face) + vals = gn.tabulate(0, Q.get_points())[(0,)*sd] + vals -= br.tabulate(0, Q.get_points())[(0,)*sd] + + # Assert that the normal component is the same + n = cell.compute_normal(f) + normal_trace = numpy.tensordot(vals, n, axes=(1, 0)) + normal_moments = inner(normal_trace, phis_at_qpts, Q.get_weights()) + assert numpy.allclose(normal_moments, 0) + + # Assert that the tangential components are the same + for t in cell.compute_tangents(sd-1, f): + tangential_trace = numpy.tensordot(vals, t, axes=(1, 0)) + tangential_moments = inner(tangential_trace, phis_at_qpts, Q.get_weights()) + assert numpy.allclose(tangential_moments, 0) + + @pytest.mark.parametrize("cell", (T, S)) @pytest.mark.parametrize("family", ("AQ", "CH", "GN", "GN2")) def test_minimal_stokes_space(cell, family): From ce787a16d415f4ec6329fdb65afb6618d574fb66 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 11:14:32 +0000 Subject: [PATCH 02/21] WIP --- FIAT/functional.py | 3 +- FIAT/mardal_tai_winther.py | 124 +++++++++++++++++++----------------- FIAT/polynomial_set.py | 18 ++++++ finat/mtw.py | 41 ++++-------- test/FIAT/unit/test_fiat.py | 1 + 5 files changed, 99 insertions(+), 88 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index 82509058..9a6dbdf6 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -459,6 +459,7 @@ def __init__(self, ref_el, Q, f_at_qpts): self.Q = Q sd = ref_el.get_spatial_dimension() + shp = f_at_qpts.shape[1:-1] + (sd,) points = Q.get_points() self.dpts = points @@ -468,7 +469,7 @@ def __init__(self, ref_el, Q, f_at_qpts): dpt_dict = {tuple(pt): [(wt, alphas[i], (i,)) for i in range(sd)] for pt, wt in zip(points, weights)} - super().__init__(ref_el, tuple(), {}, dpt_dict, + super().__init__(ref_el, shp, {}, dpt_dict, "IntegralMomentOfDivergence") diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 2c623c64..4ccdb22d 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -8,26 +8,57 @@ # SPDX-License-Identifier: LGPL-3.0-or-later +import numpy from FIAT import dual_set, expansions, finite_element, polynomial_set -from FIAT.functional import (IntegralMomentOfNormalEvaluation, - IntegralMomentOfTangentialEvaluation, - IntegralLegendreNormalMoment, - IntegralMomentOfDivergence) - +from FIAT.functional import FrobeniusIntegralMoment, IntegralMomentOfDivergence +from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature +from FIAT.nedelec import Nedelec -def DivergenceDubinerMoments(ref_el, start_deg, stop_deg, comp_deg): - sd = ref_el.get_spatial_dimension() - P = polynomial_set.ONPolynomialSet(ref_el, stop_deg) - Q = create_quadrature(ref_el, comp_deg + stop_deg) +def MardalTaiWintherSpace(ref_el): + # Generate constraint nodes on the cell and facets + # * div(v) must be constant on the cell. Since v is a cubic and + # div(v) is quadratic, we need the integral of div(v) against the + # linear and quadratic Dubiner polynomials to vanish. + # There are two linear and three quadratics, so these are five + # constraints + # * v.n must be linear on each facet. Since v.n is cubic, we need + # the integral of v.n against the cubic and quadratic Legendre + # polynomial to vanish on each facet. - dim0 = expansions.polynomial_dimension(ref_el, start_deg-1) - dim1 = expansions.polynomial_dimension(ref_el, stop_deg) - indices = list(range(dim0, dim1)) - phis = P.take(indices).tabulate(Q.get_points())[(0,)*sd] - for phi in phis: - yield IntegralMomentOfDivergence(ref_el, Q, phi) + # So we introduce functionals whose kernel describes this property, + # as described in the FIAT paper. + + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + # Polynomials of degree sd+1 + degree = sd + 1 + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,)) + + constraints = [] + # Normal component in P1 + ref_facet = ref_el.construct_subelement(sd-1) + P = polynomial_set.ONPolynomialSet(ref_facet, degree) + start = expansions.polynomial_dimension(ref_facet, 1) + stop = expansions.polynomial_dimension(ref_facet, P.degree) + Q = create_quadrature(ref_facet, degree + P.degree) + Phis = P.take(range(start, stop)).tabulate(Q.get_points())[(0,)*(sd-1)] + for f in sorted(top[sd-1]): + n = ref_el.compute_normal(f) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + phis = n[None, :, None] * Phis[:, None, :] + constraints.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + + # Divergence in P0 + P = polynomial_set.ONPolynomialSet(ref_el, degree-1) + start = expansions.polynomial_dimension(ref_el, 0) + stop = expansions.polynomial_dimension(ref_el, P.degree) + Q = create_quadrature(ref_el, degree-1 + P.degree) + phis = P.take(range(start, stop)).tabulate(Q.get_points())[(0,)*sd] + constraints.extend(IntegralMomentOfDivergence(ref_el, Q, phi) for phi in phis) + + return polynomial_set.ConstrainedPolynomialSet(constraints, poly_set) class MardalTaiWintherDual(dual_set.DualSet): @@ -36,11 +67,11 @@ def __init__(self, ref_el, degree): sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - if sd != 2: + if sd not in (2, 3): raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2.") - if degree != 3: - raise ValueError("Mardal-Tai-Winther elements are only defined for degree 3.") + if degree != sd+1: + raise ValueError("Mardal-Tai-Winther elements are only defined for degree = dim+1.") entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} nodes = [] @@ -50,44 +81,27 @@ def __init__(self, ref_el, degree): # On each facet, let n be its normal. We need to integrate # u.n and u.t against the first Legendre polynomial (constant) # and u.n against the second (linear). - facet = ref_el.get_facet_element() + ref_facet = ref_el.get_facet_element() # Facet nodes are \int_F v.n p ds where p \in P_{q-1} # degree is q - 1 - Q = create_quadrature(facet, degree+1) - Pq = polynomial_set.ONPolynomialSet(facet, 1) - phis = Pq.tabulate(Q.get_points())[(0,)*(sd - 1)] + Q = create_quadrature(ref_facet, degree+1) + P1 = polynomial_set.ONPolynomialSet(ref_facet, 1) + P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] + if sd == 2: + Phis = P1_at_qpts[:1, None, :] + else: + Ned1 = Nedelec(ref_facet, 1) + Phis = Ned1.tabulate(0, Q.get_points())[(0,)*(sd - 1)] for f in sorted(top[sd-1]): cur = len(nodes) - nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[0], f)) - nodes.append(IntegralMomentOfTangentialEvaluation(ref_el, Q, phis[0], f)) - nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[1], f)) - entity_ids[sd-1][f].extend(range(cur, len(nodes))) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + Jf = Qf.jacobian() + n = ref_el.compute_scaled_normal(f) + phis = numpy.tensordot(Jf, Phis.transpose(1, 0, 2), axes=(-1, 0)).transpose(1, 0, 2) - # Generate constraint nodes on the cell and facets - # * div(v) must be constant on the cell. Since v is a cubic and - # div(v) is quadratic, we need the integral of div(v) against the - # linear and quadratic Dubiner polynomials to vanish. - # There are two linear and three quadratics, so these are five - # constraints - # * v.n must be linear on each facet. Since v.n is cubic, we need - # the integral of v.n against the cubic and quadratic Legendre - # polynomial to vanish on each facet. - - # So we introduce functionals whose kernel describes this property, - # as described in the FIAT paper. - start_order = 2 - stop_order = 3 - qdegree = degree + stop_order - for f in sorted(top[sd-1]): - cur = len(nodes) - nodes.extend(IntegralLegendreNormalMoment(ref_el, f, order, qdegree) - for order in range(start_order, stop_order+1)) + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) entity_ids[sd-1][f].extend(range(cur, len(nodes))) - - cur = len(nodes) - nodes.extend(DivergenceDubinerMoments(ref_el, start_order-1, stop_order-1, degree)) - entity_ids[sd][0].extend(range(cur, len(nodes))) - super().__init__(nodes, ref_el, entity_ids) @@ -95,11 +109,7 @@ class MardalTaiWinther(finite_element.CiarletElement): """The definition of the Mardal-Tai-Winther element. """ def __init__(self, ref_el, degree=3): - sd = ref_el.get_spatial_dimension() - assert degree == 3, "Only defined for degree 3" - assert sd == 2, "Only defined for dimension 2" - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd,)) dual = MardalTaiWintherDual(ref_el, degree) - formdegree = sd-1 - mapping = "contravariant piola" - super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) + poly_set = MardalTaiWintherSpace(ref_el) + formdegree = ref_el.get_spatial_dimension() - 1 + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index b3aee67f..e55fa92b 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -282,6 +282,24 @@ def __init__(self, ref_el, degree, size=None, **kwargs): super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) +class ConstrainedPolynomialSet(PolynomialSet): + + def __init__(self, nodes, poly_set): + from FIAT.dual_set import DualSet + ref_el = poly_set.get_reference_element() + top = ref_el.get_topology() + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + entity_ids[max(top)][0] = list(range(len(nodes))) + + dual = DualSet(nodes, ref_el, entity_ids) + A = dual.to_riesz(poly_set) + B = poly_set.get_coeffs() + dualmat = numpy.tensordot(A, B, axes=(range(1, A.ndim), range(1, B.ndim))) + coeffs = spanning_basis(dualmat, nullspace=True) + degree = poly_set.degree + super().__init__(ref_el, degree, degree, poly_set.get_expansion_set(), coeffs) + + def make_bubbles(ref_el, degree, codim=0, shape=(), scale="L2 piola"): """Construct a polynomial set with codim bubbles up to the given degree. """ diff --git a/finat/mtw.py b/finat/mtw.py index 542ae38c..6fd95939 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -1,12 +1,10 @@ import FIAT - from gem import ListTensor from finat.citations import cite from finat.fiat_elements import FiatElement from finat.physically_mapped import identity, PhysicallyMappedElement from finat.piola_mapped import normal_tangential_edge_transform -from copy import deepcopy class MardalTaiWinther(PhysicallyMappedElement, FiatElement): @@ -14,42 +12,25 @@ def __init__(self, cell, degree=3): cite("Mardal2002") super().__init__(FIAT.MardalTaiWinther(cell, degree)) - reduced_dofs = deepcopy(self._element.entity_dofs()) - sd = cell.get_spatial_dimension() - - dim_P1 = sd - dim_Ned1 = (sd*(sd-1))//2 - fdofs = dim_P1 + dim_Ned1 - - reduced_dofs[sd][0] = [] - for f in reduced_dofs[sd-1]: - reduced_dofs[sd-1][f] = reduced_dofs[sd-1][f][:fdofs] - self._entity_dofs = reduced_dofs - self._space_dimension = fdofs * len(reduced_dofs[sd-1]) - def basis_transformation(self, coordinate_mapping): - numbf = self._element.space_dimension() - ndof = self.space_dimension() - V = identity(numbf, ndof) - sd = self.cell.get_spatial_dimension() + bary, = self.cell.make_points(sd, 0, sd+1) + J = coordinate_mapping.jacobian_at(bary) + detJ = coordinate_mapping.detJ_at(bary) + entity_dofs = self.entity_dofs() if sd == 2: facet_transform = normal_tangential_edge_transform else: raise NotImplementedError - bary, = self.cell.make_points(sd, 0, sd+1) - J = coordinate_mapping.jacobian_at(bary) - detJ = coordinate_mapping.detJ_at(bary) - entity_dofs = self.entity_dofs() + # dim_P1 = sd + # dim_Ned1 = (sd*(sd-1))//2 + # fdofs = dim_P1 + dim_Ned1 + + ndof = self.space_dimension() + V = identity(ndof, ndof) for f in sorted(entity_dofs[sd-1]): cur = entity_dofs[sd-1][f][0] - V[cur+1, cur:cur+sd] = facet_transform(self.cell, J, detJ, f) + V[cur, cur:cur+sd] = facet_transform(self.cell, J, detJ, f)[::-1] return ListTensor(V.T) - - def entity_dofs(self): - return self._entity_dofs - - def space_dimension(self): - return self._space_dimension diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index 361cd510..a4b9f0ea 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -361,6 +361,7 @@ def __init__(self, a, b): "BernardiRaugel(T)", "BernardiRaugel(S)", "MardalTaiWinther(T, 3)", + "MardalTaiWinther(S, 4)", "ArnoldWintherNC(T, 2)", "ArnoldWinther(T, 3)", "HuZhang(T, 3)", From 8b7c738fdc5cc1b476b56e31afbc053bdddf25d1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 11:22:06 +0000 Subject: [PATCH 03/21] WIP --- FIAT/mardal_tai_winther.py | 95 ++++++++++++++++----------------- finat/mtw.py | 44 ++++++++++----- test/finat/conftest.py | 6 +++ test/finat/test_zany_mapping.py | 4 +- 4 files changed, 85 insertions(+), 64 deletions(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 4ccdb22d..2a214029 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -10,55 +10,50 @@ import numpy from FIAT import dual_set, expansions, finite_element, polynomial_set -from FIAT.functional import FrobeniusIntegralMoment, IntegralMomentOfDivergence +from FIAT.functional import FrobeniusIntegralMoment from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -from FIAT.nedelec import Nedelec + + +def curl(table_u): + grad_u = {alpha.index(1): table_u[alpha] for alpha in table_u if sum(alpha) == 1} + d = len(grad_u) + if d == 2: + curl_u = [grad_u[1], -grad_u[0]] + return numpy.concatenate(curl_u, axis=1) + else: + indices = ((i, j) for i in reversed(range(d)) for j in reversed(range(i+1, d))) + curl_u = [((-1)**k) * (grad_u[j][:, i, :] - grad_u[i][:, j, :]) for k, (i, j) in enumerate(indices)] + return numpy.transpose(curl_u, (1, 0, 2)) def MardalTaiWintherSpace(ref_el): - # Generate constraint nodes on the cell and facets - # * div(v) must be constant on the cell. Since v is a cubic and - # div(v) is quadratic, we need the integral of div(v) against the - # linear and quadratic Dubiner polynomials to vanish. - # There are two linear and three quadratics, so these are five - # constraints - # * v.n must be linear on each facet. Since v.n is cubic, we need - # the integral of v.n against the cubic and quadratic Legendre - # polynomial to vanish on each facet. - - # So we introduce functionals whose kernel describes this property, - # as described in the FIAT paper. - - top = ref_el.get_topology() + """Construct the MTW space [P1]^d + curl(B [P1]^d).""" sd = ref_el.get_spatial_dimension() # Polynomials of degree sd+1 degree = sd + 1 - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,)) - - constraints = [] - # Normal component in P1 - ref_facet = ref_el.construct_subelement(sd-1) - P = polynomial_set.ONPolynomialSet(ref_facet, degree) - start = expansions.polynomial_dimension(ref_facet, 1) - stop = expansions.polynomial_dimension(ref_facet, P.degree) - Q = create_quadrature(ref_facet, degree + P.degree) - Phis = P.take(range(start, stop)).tabulate(Q.get_points())[(0,)*(sd-1)] - for f in sorted(top[sd-1]): - n = ref_el.compute_normal(f) - Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) - phis = n[None, :, None] * Phis[:, None, :] - constraints.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) - - # Divergence in P0 - P = polynomial_set.ONPolynomialSet(ref_el, degree-1) - start = expansions.polynomial_dimension(ref_el, 0) - stop = expansions.polynomial_dimension(ref_el, P.degree) - Q = create_quadrature(ref_el, degree-1 + P.degree) - phis = P.take(range(start, stop)).tabulate(Q.get_points())[(0,)*sd] - constraints.extend(IntegralMomentOfDivergence(ref_el, Q, phi) for phi in phis) - - return polynomial_set.ConstrainedPolynomialSet(constraints, poly_set) + Pk = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), scale="orthonormal") + + # Grab [P1]^d from [Pk]^d + dimP1 = expansions.polynomial_dimension(ref_el, 1) + dimPk = expansions.polynomial_dimension(ref_el, degree) + ids = [i+dimPk*j for i in range(dimP1) for j in range(sd)] + P1 = Pk.take(ids) + # Project curl(B [P1]^d) into [Pk]^d + BP1 = polynomial_set.make_bubbles(ref_el, degree+1, shape=((sd*(sd-1))//2,)) + + Q = create_quadrature(ref_el, degree*2) + qpts = Q.get_points() + qwts = Q.get_weights() + Pk_at_qpts = Pk.tabulate(qpts) + BP1_at_qpts = BP1.tabulate(qpts, 1) + + inner = lambda u, v, qwts: numpy.tensordot(u, numpy.multiply(v, qwts), axes=(range(1, u.ndim),)*2) + C = inner(curl(BP1_at_qpts), Pk_at_qpts[(0,)*sd], qwts) + coeffs = numpy.tensordot(C, Pk.get_coeffs(), axes=(1, 0)) + curlBP1 = polynomial_set.PolynomialSet(ref_el, degree, degree, Pk.get_expansion_set(), coeffs) + + return polynomial_set.polynomial_set_union_normalized(P1, curlBP1) class MardalTaiWintherDual(dual_set.DualSet): @@ -68,7 +63,7 @@ def __init__(self, ref_el, degree): top = ref_el.get_topology() if sd not in (2, 3): - raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2.") + raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2 and 3.") if degree != sd+1: raise ValueError("Mardal-Tai-Winther elements are only defined for degree = dim+1.") @@ -90,16 +85,20 @@ def __init__(self, ref_el, degree): if sd == 2: Phis = P1_at_qpts[:1, None, :] else: - Ned1 = Nedelec(ref_facet, 1) - Phis = Ned1.tabulate(0, Q.get_points())[(0,)*(sd - 1)] + Phis = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) + Phis[0, 0, :] = P1_at_qpts[0, None, :] + Phis[1, 1, :] = P1_at_qpts[0, None, :] + Phis[2, 0, :] = P1_at_qpts[1, None, :] + Phis[2, 1, :] = P1_at_qpts[2, None, :] + for f in sorted(top[sd-1]): cur = len(nodes) - Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) - Jf = Qf.jacobian() n = ref_el.compute_scaled_normal(f) - phis = numpy.tensordot(Jf, Phis.transpose(1, 0, 2), axes=(-1, 0)).transpose(1, 0, 2) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + Jf = ref_el.compute_tangents(sd-1, f) + phis = numpy.tensordot(Jf.T, Phis.transpose((1, 0, 2)), (1, 0)).transpose((1, 0, 2)) - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.cross(n, phi, axis=0) if i == 2 else phi) for i, phi in enumerate(phis)) nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) entity_ids[sd-1][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/finat/mtw.py b/finat/mtw.py index 6fd95939..fa391d94 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -1,14 +1,17 @@ import FIAT -from gem import ListTensor +import gem +import numpy from finat.citations import cite from finat.fiat_elements import FiatElement from finat.physically_mapped import identity, PhysicallyMappedElement -from finat.piola_mapped import normal_tangential_edge_transform +from finat.piola_mapped import normal_tangential_edge_transform, normal_tangential_face_transform class MardalTaiWinther(PhysicallyMappedElement, FiatElement): - def __init__(self, cell, degree=3): + def __init__(self, cell, degree=None): + if degree is None: + degree = cell.get_spatial_dimension()+1 cite("Mardal2002") super().__init__(FIAT.MardalTaiWinther(cell, degree)) @@ -18,19 +21,32 @@ def basis_transformation(self, coordinate_mapping): J = coordinate_mapping.jacobian_at(bary) detJ = coordinate_mapping.detJ_at(bary) entity_dofs = self.entity_dofs() + + ndof = self.space_dimension() + V = identity(ndof, ndof) if sd == 2: - facet_transform = normal_tangential_edge_transform + for f in sorted(entity_dofs[sd-1]): + Bnt = normal_tangential_edge_transform(self.cell, J, detJ, f) + cur = entity_dofs[sd-1][f][0] + V[cur, cur:cur+sd] = Bnt[::-1] else: - raise NotImplementedError + adjugate = lambda A: [[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]] + dim_Ned1 = (sd*(sd-1))//2 + for f in sorted(entity_dofs[sd-1]): + Bnt = normal_tangential_face_transform(self.cell, J, detJ, f) + tdofs = entity_dofs[sd-1][f][:dim_Ned1] + ndofs = entity_dofs[sd-1][f][dim_Ned1:] - # dim_P1 = sd - # dim_Ned1 = (sd*(sd-1))//2 - # fdofs = dim_P1 + dim_Ned1 + thats = self.cell.compute_tangents(sd-1, f) + Jt = J @ gem.Literal(thats.T) + Gtt = Jt.T @ Jt + detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] - ndof = self.space_dimension() - V = identity(ndof, ndof) - for f in sorted(entity_dofs[sd-1]): - cur = entity_dofs[sd-1][f][0] - V[cur, cur:cur+sd] = facet_transform(self.cell, J, detJ, f)[::-1] + V[numpy.ix_(tdofs[:2], tdofs[:2])] = adjugate(Gtt) + V[numpy.ix_(tdofs[:2], tdofs[:2])] *= detJ / detG + #V[tdofs[2], tdofs[2]] = detJ # / detG + + V[tdofs[0], ndofs[0]] = -1*Bnt[0][0] + V[tdofs[1], ndofs[0]] = -1*Bnt[1][0] - return ListTensor(V.T) + return gem.ListTensor(V.T) diff --git a/test/finat/conftest.py b/test/finat/conftest.py index c31ef735..b24e0703 100644 --- a/test/finat/conftest.py +++ b/test/finat/conftest.py @@ -116,6 +116,12 @@ def phys_el(): (-0.1, -0.2, 1.38)) return K + K[3].vertices = ((0., 0., 0.), + (1., 0., 0.), + (0., 1., 0.), + (0., 0., 2.)) + return K + @pytest.fixture def ref_to_phys(ref_el, phys_el): diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 030d2a5a..0100ac65 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -96,7 +96,7 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): inds = tuple(map(np.unique, np.nonzero(error))) error = error[np.ix_(*inds)] - error[error != 0] += 1 + #error[error != 0] += 1 pp = pprint.PrettyPrinter(width=140, compact=True) assert np.allclose(residual, 0), pp.pformat((np.round(error, 8).tolist(), *inds)) @@ -150,12 +150,12 @@ def test_argyris_point(ref_to_phys): zany_piola_elements = { 2: [ - finat.MardalTaiWinther, finat.ReducedArnoldQin, finat.ArnoldWinther, finat.ArnoldWintherNC, ], 3: [ + finat.MardalTaiWinther, finat.BernardiRaugel, finat.BernardiRaugelBubble, finat.AlfeldSorokina, From 24ae7c2b30bbb40b4763e19cb90b58bf6284cbb4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 21:26:49 +0000 Subject: [PATCH 04/21] wip --- finat/mtw.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/finat/mtw.py b/finat/mtw.py index fa391d94..23e970f7 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -38,15 +38,24 @@ def basis_transformation(self, coordinate_mapping): ndofs = entity_dofs[sd-1][f][dim_Ned1:] thats = self.cell.compute_tangents(sd-1, f) - Jt = J @ gem.Literal(thats.T) - Gtt = Jt.T @ Jt - detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] + nhat = numpy.cross(*thats) + nhat /= numpy.dot(nhat, nhat) + orths = numpy.array([numpy.cross(thats[1], nhat), + numpy.cross(nhat, thats[0])]) + Jts = J @ gem.Literal(thats.T) + Jorths = J @ gem.Literal(orths.T) + A = Jorths.T @ Jts - V[numpy.ix_(tdofs[:2], tdofs[:2])] = adjugate(Gtt) - V[numpy.ix_(tdofs[:2], tdofs[:2])] *= detJ / detG - #V[tdofs[2], tdofs[2]] = detJ # / detG + + detA = A[0, 0]*A[1, 1] - A[0, 1]*A[1, 0] + + V[numpy.ix_(tdofs[:2], tdofs[:2])] = adjugate(A) + V[numpy.ix_(tdofs[:2], tdofs[:2])] *= detJ / detA + V[tdofs[2], tdofs[2]] = detJ / detA V[tdofs[0], ndofs[0]] = -1*Bnt[0][0] V[tdofs[1], ndofs[0]] = -1*Bnt[1][0] + V[tdofs[2], ndofs[1:]] += Bnt[0][1:] + return gem.ListTensor(V.T) From dd7590a94a1dc1092986a46d8fbe171c4a01a3f9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 22:25:00 +0000 Subject: [PATCH 05/21] Done --- FIAT/mardal_tai_winther.py | 6 ++++-- finat/mtw.py | 19 +++++++++---------- test/finat/test_zany_mapping.py | 2 +- 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 2a214029..b0bff831 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -97,8 +97,10 @@ def __init__(self, ref_el, degree): Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) Jf = ref_el.compute_tangents(sd-1, f) phis = numpy.tensordot(Jf.T, Phis.transpose((1, 0, 2)), (1, 0)).transpose((1, 0, 2)) - - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.cross(n, phi, axis=0) if i == 2 else phi) for i, phi in enumerate(phis)) + if sd == 2: + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + else: + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.cross(n, phi, axis=0)) for phi in phis) nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) entity_ids[sd-1][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/finat/mtw.py b/finat/mtw.py index 23e970f7..f2a6fad3 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -30,10 +30,10 @@ def basis_transformation(self, coordinate_mapping): cur = entity_dofs[sd-1][f][0] V[cur, cur:cur+sd] = Bnt[::-1] else: - adjugate = lambda A: [[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]] dim_Ned1 = (sd*(sd-1))//2 for f in sorted(entity_dofs[sd-1]): Bnt = normal_tangential_face_transform(self.cell, J, detJ, f) + tdofs = entity_dofs[sd-1][f][:dim_Ned1] ndofs = entity_dofs[sd-1][f][dim_Ned1:] @@ -42,20 +42,19 @@ def basis_transformation(self, coordinate_mapping): nhat /= numpy.dot(nhat, nhat) orths = numpy.array([numpy.cross(thats[1], nhat), numpy.cross(nhat, thats[0])]) + Jts = J @ gem.Literal(thats.T) Jorths = J @ gem.Literal(orths.T) A = Jorths.T @ Jts - - detA = A[0, 0]*A[1, 1] - A[0, 1]*A[1, 0] + V[tdofs, tdofs] = detJ / detA - V[numpy.ix_(tdofs[:2], tdofs[:2])] = adjugate(A) - V[numpy.ix_(tdofs[:2], tdofs[:2])] *= detJ / detA - V[tdofs[2], tdofs[2]] = detJ / detA - - V[tdofs[0], ndofs[0]] = -1*Bnt[0][0] - V[tdofs[1], ndofs[0]] = -1*Bnt[1][0] + Q = numpy.dot(thats, thats.T) + a = Bnt[1][0] + b = -1*Bnt[0][0] + a, b = Q @ (a, b) - V[tdofs[2], ndofs[1:]] += Bnt[0][1:] + V[tdofs[:2], ndofs[0]] += (a, b) + V[tdofs[2], ndofs[1:]] += (a, b) return gem.ListTensor(V.T) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 0100ac65..bedc5336 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -96,7 +96,7 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): inds = tuple(map(np.unique, np.nonzero(error))) error = error[np.ix_(*inds)] - #error[error != 0] += 1 + error[error != 0] += 1 pp = pprint.PrettyPrinter(width=140, compact=True) assert np.allclose(residual, 0), pp.pformat((np.round(error, 8).tolist(), *inds)) From e2c346f602a9aceb39bb646b095a6d754c1277c9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 22:29:32 +0000 Subject: [PATCH 06/21] cleanup --- finat/mtw.py | 9 +++------ test/finat/conftest.py | 6 ------ 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/finat/mtw.py b/finat/mtw.py index f2a6fad3..94047a0f 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -50,11 +50,8 @@ def basis_transformation(self, coordinate_mapping): V[tdofs, tdofs] = detJ / detA Q = numpy.dot(thats, thats.T) - a = Bnt[1][0] - b = -1*Bnt[0][0] - a, b = Q @ (a, b) - - V[tdofs[:2], ndofs[0]] += (a, b) - V[tdofs[2], ndofs[1:]] += (a, b) + Bnt = Q @ (Bnt[1][0], -1*Bnt[0][0]) + V[tdofs[:2], ndofs[0]] += Bnt + V[tdofs[2], ndofs[1:]] += Bnt return gem.ListTensor(V.T) diff --git a/test/finat/conftest.py b/test/finat/conftest.py index b24e0703..c31ef735 100644 --- a/test/finat/conftest.py +++ b/test/finat/conftest.py @@ -116,12 +116,6 @@ def phys_el(): (-0.1, -0.2, 1.38)) return K - K[3].vertices = ((0., 0., 0.), - (1., 0., 0.), - (0., 1., 0.), - (0., 0., 2.)) - return K - @pytest.fixture def ref_to_phys(ref_el, phys_el): From d8641c6d3fd67951f99c5797f26d70ce0e36639f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 22:37:17 +0000 Subject: [PATCH 07/21] fix test --- test/FIAT/unit/test_mtw.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/FIAT/unit/test_mtw.py b/test/FIAT/unit/test_mtw.py index c1a71120..67cd7854 100644 --- a/test/FIAT/unit/test_mtw.py +++ b/test/FIAT/unit/test_mtw.py @@ -25,7 +25,7 @@ def test_dofs(): for m in (0, 1): normal_moments[bf, m] += wts[k] * nvals[bf, k] * linevals[m, k] right = np.zeros((9, 2)) - right[3*ed, 0] = 1.0 + right[3*ed+1, 0] = 1.0 right[3*ed+2, 1] = 1.0 assert np.allclose(normal_moments, right) for ed in range(3): @@ -39,5 +39,5 @@ def test_dofs(): for k in range(len(Qline.wts)): tangent_moments[bf] += wts[k] * tvals[bf, k] * linevals[0, k] right = np.zeros(9) - right[3*ed + 1] = 1.0 + right[3*ed] = 1.0 assert np.allclose(tangent_moments, right) From 40e89b5058f6e98e6ec19a7f0eb7a5f9dabbf87f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 22:44:10 +0000 Subject: [PATCH 08/21] more sensible dof ordering --- FIAT/mardal_tai_winther.py | 4 +++- finat/mtw.py | 14 ++++++++------ test/FIAT/unit/test_mtw.py | 6 +++--- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index b0bff831..a7d033bb 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -95,13 +95,15 @@ def __init__(self, ref_el, degree): cur = len(nodes) n = ref_el.compute_scaled_normal(f) Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + # Normal moments against P1 + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) + # Tangential moments against RT0 Jf = ref_el.compute_tangents(sd-1, f) phis = numpy.tensordot(Jf.T, Phis.transpose((1, 0, 2)), (1, 0)).transpose((1, 0, 2)) if sd == 2: nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) else: nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.cross(n, phi, axis=0)) for phi in phis) - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) entity_ids[sd-1][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/finat/mtw.py b/finat/mtw.py index 94047a0f..cae5399f 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -24,18 +24,20 @@ def basis_transformation(self, coordinate_mapping): ndof = self.space_dimension() V = identity(ndof, ndof) + dimP1 = sd if sd == 2: for f in sorted(entity_dofs[sd-1]): Bnt = normal_tangential_edge_transform(self.cell, J, detJ, f) - cur = entity_dofs[sd-1][f][0] - V[cur, cur:cur+sd] = Bnt[::-1] + ndofs = entity_dofs[sd-1][f][:dimP1] + tdofs = entity_dofs[sd-1][f][dimP1:] + + V[tdofs[0], ndofs[0]] = Bnt[0] + V[tdofs[0], tdofs[0]] = Bnt[1] else: - dim_Ned1 = (sd*(sd-1))//2 for f in sorted(entity_dofs[sd-1]): Bnt = normal_tangential_face_transform(self.cell, J, detJ, f) - - tdofs = entity_dofs[sd-1][f][:dim_Ned1] - ndofs = entity_dofs[sd-1][f][dim_Ned1:] + ndofs = entity_dofs[sd-1][f][:dimP1] + tdofs = entity_dofs[sd-1][f][dimP1:] thats = self.cell.compute_tangents(sd-1, f) nhat = numpy.cross(*thats) diff --git a/test/FIAT/unit/test_mtw.py b/test/FIAT/unit/test_mtw.py index 67cd7854..3ca12153 100644 --- a/test/FIAT/unit/test_mtw.py +++ b/test/FIAT/unit/test_mtw.py @@ -25,8 +25,8 @@ def test_dofs(): for m in (0, 1): normal_moments[bf, m] += wts[k] * nvals[bf, k] * linevals[m, k] right = np.zeros((9, 2)) - right[3*ed+1, 0] = 1.0 - right[3*ed+2, 1] = 1.0 + right[3*ed, 0] = 1.0 + right[3*ed+1, 1] = 1.0 assert np.allclose(normal_moments, right) for ed in range(3): t = T.compute_edge_tangent(ed) @@ -39,5 +39,5 @@ def test_dofs(): for k in range(len(Qline.wts)): tangent_moments[bf] += wts[k] * tvals[bf, k] * linevals[0, k] right = np.zeros(9) - right[3*ed] = 1.0 + right[3*ed+2] = 1.0 assert np.allclose(tangent_moments, right) From 5e7765021b7069f25b11b01d0a5ea9a010b83192 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 7 Feb 2026 22:49:51 +0000 Subject: [PATCH 09/21] enable MTW tetrahedron --- finat/ufl/elementlist.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index 4d7ae92d..f3538622 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -111,7 +111,7 @@ def show_elements(): # Zany elements register_element("Bernardi-Raugel", "BR", 1, H1, "contravariant Piola", (1, None), simplices[1:]) register_element("Bernardi-Raugel Bubble", "BRB", 1, H1, "contravariant Piola", (None, None), simplices[1:]) -register_element("Mardal-Tai-Winther", "MTW", 1, H1, "contravariant Piola", (3, 3), ("triangle",)) +register_element("Mardal-Tai-Winther", "MTW", 1, H1, "contravariant Piola", (3, 4), ("triangle", "tetrahedron")) register_element("Hermite", "HER", 0, H1, "custom", (3, 3), simplices) register_element("Argyris", "ARG", 0, H2, "custom", (5, None), ("triangle",)) register_element("Bell", "BELL", 0, H2, "custom", (5, 5), ("triangle",)) From f63ede206b4e33e360bd958ce8ff5cedddca02a4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 8 Feb 2026 08:58:13 +0000 Subject: [PATCH 10/21] cleanup --- FIAT/mardal_tai_winther.py | 38 ++++++++++++++++++-------------------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index a7d033bb..7c2cf58c 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -42,7 +42,7 @@ def MardalTaiWintherSpace(ref_el): # Project curl(B [P1]^d) into [Pk]^d BP1 = polynomial_set.make_bubbles(ref_el, degree+1, shape=((sd*(sd-1))//2,)) - Q = create_quadrature(ref_el, degree*2) + Q = create_quadrature(ref_el, 2*degree) qpts = Q.get_points() qwts = Q.get_weights() Pk_at_qpts = Pk.tabulate(qpts) @@ -71,25 +71,23 @@ def __init__(self, ref_el, degree): entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} nodes = [] - # no vertex dofs - # On each facet, let n be its normal. We need to integrate - # u.n and u.t against the first Legendre polynomial (constant) - # and u.n against the second (linear). + # u.n against a Dubiner basis for P1 + # and u x n against a basis for lowest-order RT. ref_facet = ref_el.get_facet_element() - # Facet nodes are \int_F v.n p ds where p \in P_{q-1} - # degree is q - 1 Q = create_quadrature(ref_facet, degree+1) P1 = polynomial_set.ONPolynomialSet(ref_facet, 1) P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] if sd == 2: - Phis = P1_at_qpts[:1, None, :] + # For 2D just take the constant + RT_at_qpts = P1_at_qpts[:1, None, :] else: - Phis = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) - Phis[0, 0, :] = P1_at_qpts[0, None, :] - Phis[1, 1, :] = P1_at_qpts[0, None, :] - Phis[2, 0, :] = P1_at_qpts[1, None, :] - Phis[2, 1, :] = P1_at_qpts[2, None, :] + # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] + RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) + RT_at_qpts[0, 0, :] = P1_at_qpts[0, None, :] + RT_at_qpts[1, 1, :] = P1_at_qpts[0, None, :] + RT_at_qpts[2, 0, :] = P1_at_qpts[1, None, :] + RT_at_qpts[2, 1, :] = P1_at_qpts[2, None, :] for f in sorted(top[sd-1]): cur = len(nodes) @@ -97,13 +95,13 @@ def __init__(self, ref_el, degree): Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) # Normal moments against P1 nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) - # Tangential moments against RT0 - Jf = ref_el.compute_tangents(sd-1, f) - phis = numpy.tensordot(Jf.T, Phis.transpose((1, 0, 2)), (1, 0)).transpose((1, 0, 2)) - if sd == 2: - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) - else: - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.cross(n, phi, axis=0)) for phi in phis) + # Map the RT basis into the facet + Jf = Qf.jacobian() + phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + # Moments against cross(n, RT) + phis = numpy.cross(n[None, :, None], phis, axis=1) + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) entity_ids[sd-1][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) From eba8ad47129902fa77a633bf0e3a3d0e26384b54 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 8 Feb 2026 14:26:42 +0000 Subject: [PATCH 11/21] Use n x t for tangential dofs --- FIAT/bernardi_raugel.py | 3 +++ FIAT/johnson_mercier.py | 7 ++++++- finat/aw.py | 8 ++++++-- finat/mtw.py | 43 ++++++++++++++--------------------------- finat/piola_mapped.py | 14 ++++++++------ 5 files changed, 37 insertions(+), 38 deletions(-) diff --git a/FIAT/bernardi_raugel.py b/FIAT/bernardi_raugel.py index 852f2b0d..ea2fac3a 100644 --- a/FIAT/bernardi_raugel.py +++ b/FIAT/bernardi_raugel.py @@ -116,6 +116,9 @@ def __init__(self, ref_el, order=1, degree=None, reduced=False, ref_complex=None Q = Qt[f] phi = ft_at_qpts udir = thats[f][i-1] + if sd == 3: + udir = numpy.cross(perp(*thats[f]), udir) + nodes.append(FrobeniusIntegralMoment(ref_el, Q, numpy.outer(udir, phi))) entity_ids[sd-1][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 4521c016..3cd9c9b8 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -28,7 +28,12 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): cur = len(nodes) Q = FacetQuadratureRule(ref_el, dim, f, Qref, avg=True) thats = ref_el.compute_tangents(dim, f) - nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats) + if sd == 2: + nhat = numpy.dot(R, *thats) + else: + nhat = numpy.cross(*thats) + thats = numpy.cross(nhat[None, :], thats, axis=1) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, nhat, comp, Q, phi) for phi in phis for comp in (nhat, *thats)) entity_ids[dim][f].extend(range(cur, len(nodes))) diff --git a/finat/aw.py b/finat/aw.py index 97ad614b..f498a61b 100644 --- a/finat/aw.py +++ b/finat/aw.py @@ -28,10 +28,14 @@ def _facet_transform(fiat_cell, facet_moment_degree, coordinate_mapping): transform = normal_tangential_face_transform for f in range(num_facets): - rows = transform(fiat_cell, J, detJ, f) + *Bnt, Btt = transform(fiat_cell, J, detJ, f) for i in range(dimPk_facet): s = dofs_per_facet*f + i * sd - V[s+1:s+sd, s:s+sd] = rows + ndof = s + tdofs = range(s+1, s+sd) + V[tdofs, ndof] = Bnt + V[tdofs, tdofs] = Btt + return V diff --git a/finat/mtw.py b/finat/mtw.py index cae5399f..588ebf15 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -1,6 +1,5 @@ import FIAT import gem -import numpy from finat.citations import cite from finat.fiat_elements import FiatElement @@ -25,35 +24,21 @@ def basis_transformation(self, coordinate_mapping): ndof = self.space_dimension() V = identity(ndof, ndof) dimP1 = sd - if sd == 2: - for f in sorted(entity_dofs[sd-1]): - Bnt = normal_tangential_edge_transform(self.cell, J, detJ, f) - ndofs = entity_dofs[sd-1][f][:dimP1] - tdofs = entity_dofs[sd-1][f][dimP1:] - V[tdofs[0], ndofs[0]] = Bnt[0] - V[tdofs[0], tdofs[0]] = Bnt[1] + if sd == 2: + transform = normal_tangential_edge_transform else: - for f in sorted(entity_dofs[sd-1]): - Bnt = normal_tangential_face_transform(self.cell, J, detJ, f) - ndofs = entity_dofs[sd-1][f][:dimP1] - tdofs = entity_dofs[sd-1][f][dimP1:] - - thats = self.cell.compute_tangents(sd-1, f) - nhat = numpy.cross(*thats) - nhat /= numpy.dot(nhat, nhat) - orths = numpy.array([numpy.cross(thats[1], nhat), - numpy.cross(nhat, thats[0])]) - - Jts = J @ gem.Literal(thats.T) - Jorths = J @ gem.Literal(orths.T) - A = Jorths.T @ Jts - detA = A[0, 0]*A[1, 1] - A[0, 1]*A[1, 0] - V[tdofs, tdofs] = detJ / detA - - Q = numpy.dot(thats, thats.T) - Bnt = Q @ (Bnt[1][0], -1*Bnt[0][0]) - V[tdofs[:2], ndofs[0]] += Bnt - V[tdofs[2], ndofs[1:]] += Bnt + transform = normal_tangential_face_transform + + for f in sorted(entity_dofs[sd-1]): + *Bnt, Btt = transform(self.cell, J, detJ, f) + ndofs = entity_dofs[sd-1][f][:dimP1] + tdofs = entity_dofs[sd-1][f][dimP1:] + V[tdofs, tdofs] = Btt + if sd == 2: + V[tdofs, ndofs[0]] = Bnt + else: + V[tdofs[-1], ndofs[1:]] = Bnt + V[tdofs[:-1], ndofs[0]] = Bnt return gem.ListTensor(V.T) diff --git a/finat/piola_mapped.py b/finat/piola_mapped.py index d5fbf8a2..62fefea0 100644 --- a/finat/piola_mapped.py +++ b/finat/piola_mapped.py @@ -49,9 +49,11 @@ def normal_tangential_face_transform(fiat_cell, J, detJ, f): det0 = A[1, 0] * A[2, 1] - A[1, 1] * A[2, 0] det1 = A[2, 0] * A[0, 1] - A[2, 1] * A[0, 0] det2 = A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0] - scale = detJ / det0 - rows = ((-1 * det1 / det0, -1 * scale * A[2, 1], scale * A[2, 0]), - (-1 * det2 / det0, scale * A[1, 1], -1 * scale * A[1, 0])) + rows = numpy.array((-1 * det1 / det0, -1 * det2 / det0, detJ / det0)) + + Q = numpy.dot(thats, thats.T) + R = numpy.array([[0, 1], [-1, 0]]) + rows[:2] = Q @ R @ rows[:2] return rows @@ -118,10 +120,10 @@ def basis_transformation(self, coordinate_mapping): transform = normal_tangential_face_transform for f in sorted(dofs[sd-1]): - rows = numpy.asarray(transform(self.cell, J, detJ, f)) - cur_dofs = dofs[sd-1][f] + *Bnt, Btt = transform(self.cell, J, detJ, f) + cur_dof = dofs[sd-1][f][0] cur_bfs = bfs[sd-1][f][1:] - V[numpy.ix_(cur_bfs, cur_dofs)] = rows[..., :len(cur_dofs)] + V[cur_bfs, cur_dof] = Bnt # Fix discrepancy between normal and tangential moments needs_facet_vertex_coupling = len(dofs[0][0]) > 0 and numbf > ndof From 5e206b06590f31cbbb74489288aa587a42546649 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 8 Feb 2026 17:12:54 +0000 Subject: [PATCH 12/21] Reduced Johnson-Mercier --- FIAT/__init__.py | 3 +- FIAT/functional.py | 2 +- FIAT/johnson_mercier.py | 133 +++++++++++++++++++++++++++++++- finat/__init__.py | 2 +- finat/element_factory.py | 1 + finat/johnson_mercier.py | 44 ++++++++--- finat/mtw.py | 9 +-- finat/piola_mapped.py | 2 +- finat/ufl/elementlist.py | 1 + test/finat/test_zany_mapping.py | 1 + 10 files changed, 179 insertions(+), 19 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 4b90bcb2..5ecfba55 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -27,7 +27,7 @@ from FIAT.arnold_qin import ArnoldQin from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div from FIAT.christiansen_hu import ChristiansenHu -from FIAT.johnson_mercier import JohnsonMercier +from FIAT.johnson_mercier import JohnsonMercier, ReducedJohnsonMercier from FIAT.brezzi_douglas_marini import BrezziDouglasMarini from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace from FIAT.SminusDiv import TrimmedSerendipityDiv @@ -104,6 +104,7 @@ "Guzman-Neilan 2nd kind H1": GuzmanNeilanSecondKindH1, "Guzman-Neilan H1(div)": GuzmanNeilanH1div, "Johnson-Mercier": JohnsonMercier, + "Reduced Johnson-Mercier": ReducedJohnsonMercier, "Lagrange": Lagrange, "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen, "Gauss-Lobatto-Legendre": GaussLobattoLegendre, diff --git a/FIAT/functional.py b/FIAT/functional.py index 9a6dbdf6..4041b73e 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -493,7 +493,7 @@ def __init__(self, ref_el, Q, f_at_qpts): dpt_dict = {tuple(pt): [(wt[i], alphas[j], (i, j)) for i, j in numpy.ndindex(shp)] for pt, wt in zip(points, weights)} - super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence") + super().__init__(ref_el, shp, {}, dpt_dict, "IntegralMomentOfDivergence") class PointNormalEvaluation(Functional): diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 3cd9c9b8..24979e4f 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -1,7 +1,12 @@ from FIAT import finite_element, dual_set, macro, polynomial_set from FIAT.check_format_variant import parse_quadrature_scheme -from FIAT.functional import TensorBidirectionalIntegralMoment +from FIAT.functional import (FrobeniusIntegralMoment, IntegralMomentOfTensorDivergence, + TensorBidirectionalIntegralMoment) from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature +from FIAT.nedelec_second_kind import NedelecSecondKind as N2curl +from FIAT.brezzi_douglas_marini import BrezziDouglasMarini as N2div +from FIAT.restricted import RestrictedElement import numpy @@ -62,3 +67,129 @@ def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) + + +def ReducedJohnsonMercierSpace(ref_el, degree): + ref_complex = macro.AlfeldSplit(ref_el) + poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + + # Interior constraints: moments of divergence against the orthogonal complement of RBMs + Q = create_quadrature(ref_complex, 2*(degree-1)) + qpts, qwts = Q.get_points(), Q.get_weights() + + N2 = N2curl(ref_el, degree) + edofs = N2.entity_dofs() + rbm_indices = [edofs[1][entity][0] for entity in edofs[1]] + indices = numpy.setdiff1d(range(N2.space_dimension()), rbm_indices) + + N2_at_qpts = N2.tabulate(1, qpts)[(0,) * sd] + rbms = N2_at_qpts[rbm_indices] + ells = rbms * qwts[None, None, :] + M = numpy.tensordot(ells, N2_at_qpts, axes=((1, 2), (1, 2))) + C = numpy.linalg.solve(M[:, rbm_indices], M[:, indices]) + + phis = N2_at_qpts[indices] + phis -= numpy.tensordot(C, rbms, axes=(0, 0)) + nodes = [IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis] + + # Facet constraints: moments against the orthogonal complement of RBMs + ref_facet = ref_el.get_facet_element() + Q = create_quadrature(ref_facet, 2*(degree-1)) + qpts, qwts = Q.get_points(), Q.get_weights() + + if sd == 2: + P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) + phis = numpy.zeros((sd-1, sd-1, qpts.shape[-1])) + phis[0] = P1.take([1]).tabulate(Q.get_points())[(0,)*(sd-1)] + else: + N2 = N2curl(ref_facet, degree) + edofs = N2.entity_dofs() + rbm_indices = [edofs[1][entity][0] for entity in edofs[1]] + indices = numpy.setdiff1d(range(N2.space_dimension()), rbm_indices) + + N2_at_qpts = N2.tabulate(1, qpts)[(0,) * (sd-1)] + rbms = N2_at_qpts[rbm_indices] + ells = rbms * qwts[None, None, :] + M = numpy.tensordot(ells, N2_at_qpts, axes=((1, 2), (1, 2))) + C = numpy.linalg.solve(M[:, rbm_indices], M[:, indices]) + + phis = N2_at_qpts[indices] + phis -= numpy.tensordot(C, rbms, axes=(0, 0)) + + Phis = phis + for f in top[sd-1]: + n = ref_el.compute_scaled_normal(f) + # Map the BDM basis into the facet + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + Jf = Qf.jacobian() + phis = numpy.tensordot(Jf, Phis.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + # Normal moments against cross(n, RT) + phis = numpy.cross(n[None, :, None], phis, axis=1) + phis = phis[:, None, :] * n[None, :, None] + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + + return polynomial_set.ConstrainedPolynomialSet(nodes, poly_set) + + +class ReducedJohnsonMercierDualSet(dual_set.DualSet): + def __init__(self, ref_el, degree, variant=None, quad_scheme=None): + if degree != 1: + raise ValueError("Reduced Johnson-Mercier only defined for degree=1") + if variant is not None: + raise ValueError(f"Reduced Johnson-Mercier does not have the {variant} variant") + + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + nodes = [] + # On each facet, let n be its normal. We need to integrate + # u.n against a Dubiner basis for P1 + # and u x n against a basis for lowest-order RT. + ref_facet = ref_el.get_facet_element() + Q = create_quadrature(ref_facet, degree+1) + P1 = polynomial_set.ONPolynomialSet(ref_facet, 1) + P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] + if sd == 2: + # For 2D just take the constant + RT_at_qpts = P1_at_qpts[:1, None, :] + else: + # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] + RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) + RT_at_qpts[0, 0, :] = P1_at_qpts[0, None, :] + RT_at_qpts[1, 1, :] = P1_at_qpts[0, None, :] + RT_at_qpts[2, 0, :] = P1_at_qpts[1, None, :] + RT_at_qpts[2, 1, :] = P1_at_qpts[2, None, :] + + for f in sorted(top[sd-1]): + cur = len(nodes) + n = ref_el.compute_scaled_normal(f) + nn = numpy.outer(n, n) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + # Normal-normal moments against P1 + phis = P1_at_qpts[:, None, None, :] * nn[None, :, :, None] + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + # Map the RT basis into the facet + Jf = Qf.jacobian() + phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + # Normal moments against cross(n, RT) + phis = numpy.cross(n[None, :, None], phis, axis=1) + phis = phis[:, None, :] * n[None, :, None] + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) + + +class ReducedJohnsonMercier(finite_element.CiarletElement): + """The Reduced Johnson-Mercier finite element.""" + + def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): + poly_set = ReducedJohnsonMercierSpace(ref_el, degree) + dual = ReducedJohnsonMercierDualSet(ref_el, degree, variant=variant, quad_scheme=quad_scheme) + formdegree = ref_el.get_spatial_dimension() - 1 + mapping = "double contravariant piola" + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/finat/__init__.py b/finat/__init__.py index bde435a1..186816d8 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -28,7 +28,7 @@ from .guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanBubble, GuzmanNeilanH1div # noqa: F401 from .powell_sabin import QuadraticPowellSabin6, QuadraticPowellSabin12 # noqa: F401 from .hermite import Hermite # noqa: F401 -from .johnson_mercier import JohnsonMercier # noqa: F401 +from .johnson_mercier import JohnsonMercier, ReducedJohnsonMercier # noqa: F401 from .mtw import MardalTaiWinther # noqa: F401 from .morley import Morley # noqa: F401 from .walkington import Walkington # noqa: F401 diff --git a/finat/element_factory.py b/finat/element_factory.py index 8563aa60..7129200b 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -66,6 +66,7 @@ "Guzman-Neilan H1(div)": finat.GuzmanNeilanH1div, "Guzman-Neilan Bubble": finat.GuzmanNeilanBubble, "Johnson-Mercier": finat.JohnsonMercier, + "Reduced Johnson-Mercier": finat.ReducedJohnsonMercier, "Lagrange": finat.Lagrange, "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, "Gauss-Lobatto-Legendre": finat.GaussLobattoLegendre, diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index 7b6d485c..4bccc706 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -5,25 +5,51 @@ from finat.citations import cite from finat.fiat_elements import FiatElement from finat.physically_mapped import identity, PhysicallyMappedElement +from finat.piola_mapped import normal_tangential_edge_transform, normal_tangential_face_transform class JohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued def __init__(self, cell, degree=1, variant=None, quad_scheme=None): cite("Gopalakrishnan2024") - self._indices = slice(None, None) super().__init__(FIAT.JohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) def basis_transformation(self, coordinate_mapping): - numbf = self._element.space_dimension() - ndof = self.space_dimension() - - V = identity(numbf, ndof) + V = identity(self.space_dimension()) Vsub = _facet_transform(self.cell, 1, coordinate_mapping) - Vsub = Vsub[:, self._indices] m, n = Vsub.shape V[:m, :n] = Vsub + return ListTensor(V.T) + + +class ReducedJohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued + def __init__(self, cell, degree=1, variant=None, quad_scheme=None): + cite("Gopalakrishnan2024") + super().__init__(FIAT.ReducedJohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) + + def basis_transformation(self, coordinate_mapping): + sd = self.cell.get_spatial_dimension() + bary, = self.cell.make_points(sd, 0, sd+1) + J = coordinate_mapping.jacobian_at(bary) + detJ = coordinate_mapping.detJ_at(bary) + + V = identity(self.space_dimension()) + dimP1 = sd + + if sd == 2: + transform = normal_tangential_edge_transform + else: + transform = normal_tangential_face_transform + + entity_dofs = self.entity_dofs() + for f in sorted(entity_dofs[sd-1]): + *Bnt, Btt = transform(self.cell, J, detJ, f) + ndofs = entity_dofs[sd-1][f][:dimP1] + tdofs = entity_dofs[sd-1][f][dimP1:] + V[tdofs, tdofs] = Btt + if sd == 2: + V[tdofs, ndofs[0]] = Bnt + else: + V[tdofs[-1], ndofs[1:]] = Bnt + V[tdofs[:-1], ndofs[0]] = Bnt - # Note: that the edge DOFs are scaled by edge lengths in FIAT implies - # that they are already have the necessary rescaling to improve - # conditioning. return ListTensor(V.T) diff --git a/finat/mtw.py b/finat/mtw.py index 588ebf15..c4173e31 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -1,5 +1,5 @@ import FIAT -import gem +from gem import ListTensor from finat.citations import cite from finat.fiat_elements import FiatElement @@ -19,10 +19,8 @@ def basis_transformation(self, coordinate_mapping): bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) detJ = coordinate_mapping.detJ_at(bary) - entity_dofs = self.entity_dofs() - ndof = self.space_dimension() - V = identity(ndof, ndof) + V = identity(self.space_dimension()) dimP1 = sd if sd == 2: @@ -30,6 +28,7 @@ def basis_transformation(self, coordinate_mapping): else: transform = normal_tangential_face_transform + entity_dofs = self.entity_dofs() for f in sorted(entity_dofs[sd-1]): *Bnt, Btt = transform(self.cell, J, detJ, f) ndofs = entity_dofs[sd-1][f][:dimP1] @@ -41,4 +40,4 @@ def basis_transformation(self, coordinate_mapping): V[tdofs[-1], ndofs[1:]] = Bnt V[tdofs[:-1], ndofs[0]] = Bnt - return gem.ListTensor(V.T) + return ListTensor(V.T) diff --git a/finat/piola_mapped.py b/finat/piola_mapped.py index 62fefea0..eb3adbf2 100644 --- a/finat/piola_mapped.py +++ b/finat/piola_mapped.py @@ -53,7 +53,7 @@ def normal_tangential_face_transform(fiat_cell, J, detJ, f): Q = numpy.dot(thats, thats.T) R = numpy.array([[0, 1], [-1, 0]]) - rows[:2] = Q @ R @ rows[:2] + rows[:2] = (Q @ R) @ rows[:2] return rows diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index f3538622..b4743094 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -127,6 +127,7 @@ def show_elements(): register_element("Hsieh-Clough-Tocher", "HCT", 0, H2, "custom", (3, None), ("triangle",)) register_element("Reduced-Hsieh-Clough-Tocher", "HCT-red", 0, H2, "custom", (3, 3), ("triangle",)) register_element("Johnson-Mercier", "JM", 2, HDiv, "double contravariant Piola", (1, 1), simplices[1:]) +register_element("Reduced Johnson-Mercier", "JM-red", 2, HDiv, "double contravariant Piola", (1, 1), simplices[1:]) register_element("Walkington", "WALK", 0, H2, "custom", (5, 5), ("tetrahedron",)) register_element("Alfeld C2", "ALF-C2", 0, H3, "custom", (5, None), ("triangle",)) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index bedc5336..b115afde 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -161,6 +161,7 @@ def test_argyris_point(ref_to_phys): finat.AlfeldSorokina, finat.ChristiansenHu, finat.JohnsonMercier, + finat.ReducedJohnsonMercier, finat.GuzmanNeilanFirstKindH1, finat.GuzmanNeilanSecondKindH1, finat.GuzmanNeilanBubble, From cd06febe0a371a7d656821ed257333f4f96414ae Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 8 Feb 2026 18:09:23 +0000 Subject: [PATCH 13/21] cleanup --- finat/piola_mapped.py | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/finat/piola_mapped.py b/finat/piola_mapped.py index 62fefea0..40555cef 100644 --- a/finat/piola_mapped.py +++ b/finat/piola_mapped.py @@ -1,7 +1,7 @@ import numpy from finat.fiat_elements import FiatElement -from finat.physically_mapped import adjugate, identity, PhysicallyMappedElement +from finat.physically_mapped import adjugate, determinant, identity, PhysicallyMappedElement from gem import Literal, ListTensor, Zero from copy import deepcopy from itertools import chain @@ -38,23 +38,21 @@ def normal_tangential_face_transform(fiat_cell, J, detJ, f): thats = fiat_cell.compute_tangents(2, f) nhat = numpy.cross(*thats) nhat /= numpy.dot(nhat, nhat) - orth_vecs = numpy.array([nhat, - numpy.cross(nhat, thats[1]), - numpy.cross(thats[0], nhat)]) - # Compute A = (alpha, beta, gamma) - Jts = J @ Literal(thats.T) - Jorths = J @ Literal(orth_vecs.T) - A = Jorths.T @ Jts - # Compute the last two rows of inv([[1, 0, 0], A.T/detJ]) - det0 = A[1, 0] * A[2, 1] - A[1, 1] * A[2, 0] - det1 = A[2, 0] * A[0, 1] - A[2, 1] * A[0, 0] - det2 = A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0] - rows = numpy.array((-1 * det1 / det0, -1 * det2 / det0, detJ / det0)) + orths = numpy.cross(thats, nhat[None, :], axis=1) + + Jn = J @ Literal(nhat) + Jthats = J @ Literal(thats.T) + Jorths = J @ Literal(orths.T) + A = Jthats.T @ Jorths + B = Jn @ Jthats + A = numpy.array([[A[i, j] for j in range(A.shape[1])] for i in range(A.shape[0])]) + B = numpy.array([B[i] for i in range(B.shape[0])]) Q = numpy.dot(thats, thats.T) - R = numpy.array([[0, 1], [-1, 0]]) - rows[:2] = Q @ R @ rows[:2] - return rows + beta = determinant(A) + alpha = Q @ (adjugate(A) @ B) + row = (alpha[0] / beta, alpha[1] / beta, detJ / beta) + return row class PiolaBubbleElement(PhysicallyMappedElement, FiatElement): From 0dd59b72aa082b5f0fc9f97f19487f5f65fe0966 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 9 Feb 2026 14:03:54 +0000 Subject: [PATCH 14/21] WIP --- FIAT/johnson_mercier.py | 109 +++++++++++++++++++-------------------- FIAT/polynomial_set.py | 5 +- finat/johnson_mercier.py | 44 +++++++++++++++- 3 files changed, 97 insertions(+), 61 deletions(-) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 24979e4f..f45e61af 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -6,7 +6,6 @@ from FIAT.quadrature_schemes import create_quadrature from FIAT.nedelec_second_kind import NedelecSecondKind as N2curl from FIAT.brezzi_douglas_marini import BrezziDouglasMarini as N2div -from FIAT.restricted import RestrictedElement import numpy @@ -69,69 +68,61 @@ def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) -def ReducedJohnsonMercierSpace(ref_el, degree): - ref_complex = macro.AlfeldSplit(ref_el) - poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) - sd = ref_el.get_spatial_dimension() - top = ref_el.get_topology() - - # Interior constraints: moments of divergence against the orthogonal complement of RBMs - Q = create_quadrature(ref_complex, 2*(degree-1)) - qpts, qwts = Q.get_points(), Q.get_weights() - - N2 = N2curl(ref_el, degree) - edofs = N2.entity_dofs() - rbm_indices = [edofs[1][entity][0] for entity in edofs[1]] - indices = numpy.setdiff1d(range(N2.space_dimension()), rbm_indices) +def rbm_complement(ref_el, fe=None): + if ref_el.get_spatial_dimension() == 1: + P1 = polynomial_set.ONPolynomialSet(ref_el, 1, shape=(1,)) + return P1.take(range(1, len(P1))) + else: + if fe is None: + fe = N2curl + P1 = fe(ref_el, 1) + entity_ids = P1.entity_dofs() + ids = [] + for entity in entity_ids[1]: + ids.extend(entity_ids[1][entity][1:]) + return P1.get_nodal_basis().take(ids) - N2_at_qpts = N2.tabulate(1, qpts)[(0,) * sd] - rbms = N2_at_qpts[rbm_indices] - ells = rbms * qwts[None, None, :] - M = numpy.tensordot(ells, N2_at_qpts, axes=((1, 2), (1, 2))) - C = numpy.linalg.solve(M[:, rbm_indices], M[:, indices]) - phis = N2_at_qpts[indices] - phis -= numpy.tensordot(C, rbms, axes=(0, 0)) - nodes = [IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis] +def constraints(ref_el, degree): + sd = ref_el.get_spatial_dimension() + nodes = [] # Facet constraints: moments against the orthogonal complement of RBMs - ref_facet = ref_el.get_facet_element() - Q = create_quadrature(ref_facet, 2*(degree-1)) - qpts, qwts = Q.get_points(), Q.get_weights() + ref_facet = ref_el.construct_subelement(sd-1) + Q = create_quadrature(ref_facet, 2*degree) + P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) + P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] if sd == 2: - P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) - phis = numpy.zeros((sd-1, sd-1, qpts.shape[-1])) - phis[0] = P1.take([1]).tabulate(Q.get_points())[(0,)*(sd-1)] + # For 2D just take the constant + RT_at_qpts = P1_at_qpts[1:, None, :] else: - N2 = N2curl(ref_facet, degree) - edofs = N2.entity_dofs() - rbm_indices = [edofs[1][entity][0] for entity in edofs[1]] - indices = numpy.setdiff1d(range(N2.space_dimension()), rbm_indices) - - N2_at_qpts = N2.tabulate(1, qpts)[(0,) * (sd-1)] - rbms = N2_at_qpts[rbm_indices] - ells = rbms * qwts[None, None, :] - M = numpy.tensordot(ells, N2_at_qpts, axes=((1, 2), (1, 2))) - C = numpy.linalg.solve(M[:, rbm_indices], M[:, indices]) - - phis = N2_at_qpts[indices] - phis -= numpy.tensordot(C, rbms, axes=(0, 0)) - - Phis = phis - for f in top[sd-1]: + # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] + RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) + RT_at_qpts[0, 0, :] = P1_at_qpts[1, None, :] - P1_at_qpts[2, None, :] + RT_at_qpts[1, 1, :] = P1_at_qpts[2, None, :] - P1_at_qpts[1, None, :] + RT_at_qpts[2, 0, :] = P1_at_qpts[2, None, :] + RT_at_qpts[2, 1, :] = -P1_at_qpts[1, None, :] + + Phis = RT_at_qpts + for f in ref_el.topology[sd-1]: n = ref_el.compute_scaled_normal(f) - # Map the BDM basis into the facet Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) Jf = Qf.jacobian() phis = numpy.tensordot(Jf, Phis.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) if sd == 3: - # Normal moments against cross(n, RT) phis = numpy.cross(n[None, :, None], phis, axis=1) - phis = phis[:, None, :] * n[None, :, None] + phis = phis[:, :, None, :] * n[None, None, :, None] nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) - return polynomial_set.ConstrainedPolynomialSet(nodes, poly_set) + # Interior constraints: moments of divergence against the orthogonal complement of RBMs + ref_complex = macro.AlfeldSplit(ref_el) + Q = create_quadrature(ref_complex, 2*degree-1) + + comp_space = rbm_complement(ref_el) + phis = comp_space.tabulate(Q.get_points())[(0,)*sd] + nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) + return nodes class ReducedJohnsonMercierDualSet(dual_set.DualSet): @@ -146,11 +137,11 @@ def __init__(self, ref_el, degree, variant=None, quad_scheme=None): entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} nodes = [] # On each facet, let n be its normal. We need to integrate - # u.n against a Dubiner basis for P1 - # and u x n against a basis for lowest-order RT. + # sigma.nn against a Dubiner basis for P1 + # and (sigma.n) x n against a basis for lowest-order RT. ref_facet = ref_el.get_facet_element() Q = create_quadrature(ref_facet, degree+1) - P1 = polynomial_set.ONPolynomialSet(ref_facet, 1) + P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] if sd == 2: # For 2D just take the constant @@ -166,21 +157,23 @@ def __init__(self, ref_el, degree, variant=None, quad_scheme=None): for f in sorted(top[sd-1]): cur = len(nodes) n = ref_el.compute_scaled_normal(f) - nn = numpy.outer(n, n) Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) # Normal-normal moments against P1 - phis = P1_at_qpts[:, None, None, :] * nn[None, :, :, None] - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n, n, Qf, phi) for phi in P1_at_qpts) # Map the RT basis into the facet Jf = Qf.jacobian() phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) if sd == 3: # Normal moments against cross(n, RT) phis = numpy.cross(n[None, :, None], phis, axis=1) - phis = phis[:, None, :] * n[None, :, None] + phis = phis[:, :, None, :] * n[None, None, :, None] nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) entity_ids[sd-1][f].extend(range(cur, len(nodes))) + cur = len(nodes) + nodes.extend(constraints(ref_el, degree)) + entity_ids[sd][0].extend(range(cur, len(nodes))) + super().__init__(nodes, ref_el, entity_ids) @@ -188,7 +181,9 @@ class ReducedJohnsonMercier(finite_element.CiarletElement): """The Reduced Johnson-Mercier finite element.""" def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): - poly_set = ReducedJohnsonMercierSpace(ref_el, degree) + + ref_complex = macro.AlfeldSplit(ref_el) + poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) dual = ReducedJohnsonMercierDualSet(ref_el, degree, variant=variant, quad_scheme=quad_scheme) formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index e55fa92b..91dbf46c 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -284,7 +284,7 @@ def __init__(self, ref_el, degree, size=None, **kwargs): class ConstrainedPolynomialSet(PolynomialSet): - def __init__(self, nodes, poly_set): + def __init__(self, nodes, poly_set, nullspace=True): from FIAT.dual_set import DualSet ref_el = poly_set.get_reference_element() top = ref_el.get_topology() @@ -295,7 +295,8 @@ def __init__(self, nodes, poly_set): A = dual.to_riesz(poly_set) B = poly_set.get_coeffs() dualmat = numpy.tensordot(A, B, axes=(range(1, A.ndim), range(1, B.ndim))) - coeffs = spanning_basis(dualmat, nullspace=True) + nsp = spanning_basis(dualmat, nullspace=nullspace) + coeffs = numpy.tensordot(nsp, B, (1, 0)) degree = poly_set.degree super().__init__(ref_el, degree, degree, poly_set.get_expansion_set(), coeffs) diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index 4bccc706..7fe27cf1 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -1,4 +1,5 @@ import FIAT +import numpy from gem import ListTensor from finat.aw import _facet_transform @@ -6,6 +7,7 @@ from finat.fiat_elements import FiatElement from finat.physically_mapped import identity, PhysicallyMappedElement from finat.piola_mapped import normal_tangential_edge_transform, normal_tangential_face_transform +from copy import deepcopy class JohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued @@ -22,9 +24,35 @@ def basis_transformation(self, coordinate_mapping): class ReducedJohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued + def __init__(self, cell, degree=1, variant=None, quad_scheme=None): cite("Gopalakrishnan2024") super().__init__(FIAT.ReducedJohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) + # On each facet we expect the normal dof followed by the tangential ones + # The tangential dofs should be numbered last, and are constrained to be zero + sd = self.cell.get_spatial_dimension() + reduced_dofs = deepcopy(self._element.entity_dofs()) + reduced_dim = 0 + + dimP1 = sd + dimNed1 = (sd*(sd-1))//2 + r = dimP1 + dimNed1 + for entity in sorted(reduced_dofs[sd-1]): + reduced_dim += len(reduced_dofs[sd-1][entity][r:]) + reduced_dofs[sd-1][entity] = reduced_dofs[sd-1][entity][:r] + + for entity in sorted(reduced_dofs[sd]): + reduced_dim += len(reduced_dofs[sd][entity]) + reduced_dofs[sd][entity] = [] + + self._entity_dofs = reduced_dofs + self._space_dimension = self._element.space_dimension() - reduced_dim + + def entity_dofs(self): + return self._entity_dofs + + def space_dimension(self): + return self._space_dimension def basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() @@ -32,8 +60,11 @@ def basis_transformation(self, coordinate_mapping): J = coordinate_mapping.jacobian_at(bary) detJ = coordinate_mapping.detJ_at(bary) - V = identity(self.space_dimension()) + num_dofs = self.space_dimension() + num_bfs = self._element.space_dimension() + V = identity(num_bfs, num_dofs) dimP1 = sd + dimNed1 = (sd*(sd-1))//2 if sd == 2: transform = normal_tangential_edge_transform @@ -41,15 +72,24 @@ def basis_transformation(self, coordinate_mapping): transform = normal_tangential_face_transform entity_dofs = self.entity_dofs() + constraints = {f: list(range(num_dofs+f*dimNed1, num_dofs+(f+1)*dimNed1)) + for f in entity_dofs[sd-1]} for f in sorted(entity_dofs[sd-1]): *Bnt, Btt = transform(self.cell, J, detJ, f) ndofs = entity_dofs[sd-1][f][:dimP1] tdofs = entity_dofs[sd-1][f][dimP1:] + cdofs = constraints[f] + V[tdofs, tdofs] = Btt if sd == 2: V[tdofs, ndofs[0]] = Bnt + V[cdofs, ndofs[1]] = Bnt else: - V[tdofs[-1], ndofs[1:]] = Bnt V[tdofs[:-1], ndofs[0]] = Bnt + V[tdofs[-1], ndofs[1:]] = Bnt + # FIXME + V[cdofs[0], ndofs[1:]] = Bnt + V[cdofs[1], ndofs[1:]] = Bnt + V[cdofs[2], ndofs[1:]] = Bnt return ListTensor(V.T) From 5e831be43c26cbef1e14e40dbdbf12ba90641824 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 9 Feb 2026 21:19:19 +0000 Subject: [PATCH 15/21] Update FIAT/polynomial_set.py --- FIAT/polynomial_set.py | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index e55fa92b..b3aee67f 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -282,24 +282,6 @@ def __init__(self, ref_el, degree, size=None, **kwargs): super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) -class ConstrainedPolynomialSet(PolynomialSet): - - def __init__(self, nodes, poly_set): - from FIAT.dual_set import DualSet - ref_el = poly_set.get_reference_element() - top = ref_el.get_topology() - entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} - entity_ids[max(top)][0] = list(range(len(nodes))) - - dual = DualSet(nodes, ref_el, entity_ids) - A = dual.to_riesz(poly_set) - B = poly_set.get_coeffs() - dualmat = numpy.tensordot(A, B, axes=(range(1, A.ndim), range(1, B.ndim))) - coeffs = spanning_basis(dualmat, nullspace=True) - degree = poly_set.degree - super().__init__(ref_el, degree, degree, poly_set.get_expansion_set(), coeffs) - - def make_bubbles(ref_el, degree, codim=0, shape=(), scale="L2 piola"): """Construct a polynomial set with codim bubbles up to the given degree. """ From c976012763cfd5960e7869c9a476208bcb408104 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Feb 2026 08:58:01 +0000 Subject: [PATCH 16/21] fixup --- FIAT/johnson_mercier.py | 114 ++++++++++++++++----------------------- finat/johnson_mercier.py | 11 ++-- 2 files changed, 52 insertions(+), 73 deletions(-) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index f45e61af..19ef370b 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -4,8 +4,7 @@ TensorBidirectionalIntegralMoment) from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -from FIAT.nedelec_second_kind import NedelecSecondKind as N2curl -from FIAT.brezzi_douglas_marini import BrezziDouglasMarini as N2div +from FIAT.nedelec_second_kind import NedelecSecondKind import numpy @@ -68,14 +67,13 @@ def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) -def rbm_complement(ref_el, fe=None): +def rbm_complement(ref_el): + """Constructs a basis for the complement of the rigid body motions over P1.""" if ref_el.get_spatial_dimension() == 1: P1 = polynomial_set.ONPolynomialSet(ref_el, 1, shape=(1,)) return P1.take(range(1, len(P1))) else: - if fe is None: - fe = N2curl - P1 = fe(ref_el, 1) + P1 = NedelecSecondKind(ref_el, 1) entity_ids = P1.entity_dofs() ids = [] for entity in entity_ids[1]: @@ -83,97 +81,78 @@ def rbm_complement(ref_el, fe=None): return P1.get_nodal_basis().take(ids) -def constraints(ref_el, degree): - sd = ref_el.get_spatial_dimension() - nodes = [] - - # Facet constraints: moments against the orthogonal complement of RBMs - ref_facet = ref_el.construct_subelement(sd-1) - Q = create_quadrature(ref_facet, 2*degree) - - P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) - P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] - if sd == 2: - # For 2D just take the constant - RT_at_qpts = P1_at_qpts[1:, None, :] - else: - # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] - RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) - RT_at_qpts[0, 0, :] = P1_at_qpts[1, None, :] - P1_at_qpts[2, None, :] - RT_at_qpts[1, 1, :] = P1_at_qpts[2, None, :] - P1_at_qpts[1, None, :] - RT_at_qpts[2, 0, :] = P1_at_qpts[2, None, :] - RT_at_qpts[2, 1, :] = -P1_at_qpts[1, None, :] - - Phis = RT_at_qpts - for f in ref_el.topology[sd-1]: - n = ref_el.compute_scaled_normal(f) - Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) - Jf = Qf.jacobian() - phis = numpy.tensordot(Jf, Phis.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) - if sd == 3: - phis = numpy.cross(n[None, :, None], phis, axis=1) - phis = phis[:, :, None, :] * n[None, None, :, None] - nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) - - # Interior constraints: moments of divergence against the orthogonal complement of RBMs - ref_complex = macro.AlfeldSplit(ref_el) - Q = create_quadrature(ref_complex, 2*degree-1) - - comp_space = rbm_complement(ref_el) - phis = comp_space.tabulate(Q.get_points())[(0,)*sd] - nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) - return nodes - - class ReducedJohnsonMercierDualSet(dual_set.DualSet): - def __init__(self, ref_el, degree, variant=None, quad_scheme=None): + def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): if degree != 1: raise ValueError("Reduced Johnson-Mercier only defined for degree=1") if variant is not None: raise ValueError(f"Reduced Johnson-Mercier does not have the {variant} variant") + ref_el = ref_complex.get_parent() sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} nodes = [] - # On each facet, let n be its normal. We need to integrate - # sigma.nn against a Dubiner basis for P1 - # and (sigma.n) x n against a basis for lowest-order RT. + ref_facet = ref_el.get_facet_element() Q = create_quadrature(ref_facet, degree+1) P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] - if sd == 2: - # For 2D just take the constant - RT_at_qpts = P1_at_qpts[:1, None, :] - else: + dimP1 = len(P1)*(sd-1) + dimNed1 = dimP1 // 2 + if sd == 3: # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] - RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) - RT_at_qpts[0, 0, :] = P1_at_qpts[0, None, :] - RT_at_qpts[1, 1, :] = P1_at_qpts[0, None, :] - RT_at_qpts[2, 0, :] = P1_at_qpts[1, None, :] - RT_at_qpts[2, 1, :] = P1_at_qpts[2, None, :] + RT_at_qpts = numpy.zeros((dimP1, sd-1, P1_at_qpts.shape[-1])) + RT_at_qpts[0, 0, :] = P1_at_qpts[0] + RT_at_qpts[1, 1, :] = P1_at_qpts[0] + RT_at_qpts[2, 0, :] = P1_at_qpts[1] + RT_at_qpts[2, 1, :] = P1_at_qpts[2] + # Basis for complement of RT [(y-x, 0), (0, y-x), (y, -x)] + RT_at_qpts[3, 0, :] = P1_at_qpts[2] - P1_at_qpts[1] + RT_at_qpts[4, 1, :] = P1_at_qpts[2] - P1_at_qpts[1] + RT_at_qpts[5, 0, :] = P1_at_qpts[2] + RT_at_qpts[5, 1, :] = -P1_at_qpts[1] + else: + RT_at_qpts = P1_at_qpts[:, None, :] + # Facet degrees of freedom for f in sorted(top[sd-1]): cur = len(nodes) n = ref_el.compute_scaled_normal(f) Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + Jf = Qf.jacobian() # Normal-normal moments against P1 nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n, n, Qf, phi) for phi in P1_at_qpts) - # Map the RT basis into the facet + # Normal-tangential moments against n x RT + phis = numpy.tensordot(Jf, RT_at_qpts[:dimNed1].transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + phis = numpy.cross(n[None, :, None], phis, axis=1) + phis = phis[:, :, None, :] * n[None, None, :, None] + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + # Facet constraints + for f in sorted(top[sd-1]): + cur = len(nodes) + n = ref_el.compute_scaled_normal(f) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) Jf = Qf.jacobian() - phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + # Normal-tangential moments against n x (P1 \ RT) + phis = numpy.tensordot(Jf, RT_at_qpts[dimNed1:].transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) if sd == 3: - # Normal moments against cross(n, RT) phis = numpy.cross(n[None, :, None], phis, axis=1) phis = phis[:, :, None, :] * n[None, None, :, None] nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) entity_ids[sd-1][f].extend(range(cur, len(nodes))) + # Interior constraints: moments of divergence against the orthogonal complement of RBMs + ref_complex = macro.AlfeldSplit(ref_el) + Q = create_quadrature(ref_complex, 2*degree-1) + comp_space = rbm_complement(ref_el) + phis = comp_space.tabulate(Q.get_points())[(0,)*sd] cur = len(nodes) - nodes.extend(constraints(ref_el, degree)) + nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) entity_ids[sd][0].extend(range(cur, len(nodes))) - super().__init__(nodes, ref_el, entity_ids) @@ -181,10 +160,9 @@ class ReducedJohnsonMercier(finite_element.CiarletElement): """The Reduced Johnson-Mercier finite element.""" def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): - ref_complex = macro.AlfeldSplit(ref_el) poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) - dual = ReducedJohnsonMercierDualSet(ref_el, degree, variant=variant, quad_scheme=quad_scheme) + dual = ReducedJohnsonMercierDualSet(ref_complex, degree, variant=variant, quad_scheme=quad_scheme) formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index 7fe27cf1..102c60c7 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -1,5 +1,4 @@ import FIAT -import numpy from gem import ListTensor from finat.aw import _facet_transform @@ -87,9 +86,11 @@ def basis_transformation(self, coordinate_mapping): else: V[tdofs[:-1], ndofs[0]] = Bnt V[tdofs[-1], ndofs[1:]] = Bnt - # FIXME - V[cdofs[0], ndofs[1:]] = Bnt - V[cdofs[1], ndofs[1:]] = Bnt - V[cdofs[2], ndofs[1:]] = Bnt + # (n x t0) * (y - x) + V[cdofs[0], ndofs[1:]] = (-1*Bnt[0], Bnt[0]) + # (n x t1) * (y - x) + V[cdofs[1], ndofs[1:]] = (-1*Bnt[1], Bnt[1]) + # n x (t0 * y - t1 * x) + V[cdofs[2], ndofs[1:]] = (-1*Bnt[1], Bnt[0]) return ListTensor(V.T) From 4e9049144ea0f1e801164950e9d5860d4ebce570 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Feb 2026 09:14:13 +0000 Subject: [PATCH 17/21] cleanup --- FIAT/johnson_mercier.py | 8 +++----- finat/johnson_mercier.py | 13 +++++-------- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 19ef370b..afa200a6 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -3,7 +3,6 @@ from FIAT.functional import (FrobeniusIntegralMoment, IntegralMomentOfTensorDivergence, TensorBidirectionalIntegralMoment) from FIAT.quadrature import FacetQuadratureRule -from FIAT.quadrature_schemes import create_quadrature from FIAT.nedelec_second_kind import NedelecSecondKind import numpy @@ -87,7 +86,6 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): raise ValueError("Reduced Johnson-Mercier only defined for degree=1") if variant is not None: raise ValueError(f"Reduced Johnson-Mercier does not have the {variant} variant") - ref_el = ref_complex.get_parent() sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() @@ -95,7 +93,7 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): nodes = [] ref_facet = ref_el.get_facet_element() - Q = create_quadrature(ref_facet, degree+1) + Q = parse_quadrature_scheme(ref_facet, degree+1, quad_scheme) P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] dimP1 = len(P1)*(sd-1) @@ -145,9 +143,9 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) entity_ids[sd-1][f].extend(range(cur, len(nodes))) - # Interior constraints: moments of divergence against the orthogonal complement of RBMs + # Interior constraints: moments of divergence against (P1 \ Nedelec) ref_complex = macro.AlfeldSplit(ref_el) - Q = create_quadrature(ref_complex, 2*degree-1) + Q = parse_quadrature_scheme(ref_complex, 2*degree-1) comp_space = rbm_complement(ref_el) phis = comp_space.tabulate(Q.get_points())[(0,)*sd] cur = len(nodes) diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index 102c60c7..e1402633 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -23,12 +23,11 @@ def basis_transformation(self, coordinate_mapping): class ReducedJohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued - def __init__(self, cell, degree=1, variant=None, quad_scheme=None): cite("Gopalakrishnan2024") super().__init__(FIAT.ReducedJohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) - # On each facet we expect the normal dof followed by the tangential ones - # The tangential dofs should be numbered last, and are constrained to be zero + # On each facet we expect the normal dofs followed by the tangential ones + # The tangential and interior constraints should be numbered last. sd = self.cell.get_spatial_dimension() reduced_dofs = deepcopy(self._element.entity_dofs()) reduced_dim = 0 @@ -70,14 +69,12 @@ def basis_transformation(self, coordinate_mapping): else: transform = normal_tangential_face_transform - entity_dofs = self.entity_dofs() - constraints = {f: list(range(num_dofs+f*dimNed1, num_dofs+(f+1)*dimNed1)) - for f in entity_dofs[sd-1]} + entity_dofs = self._element.entity_dofs() for f in sorted(entity_dofs[sd-1]): *Bnt, Btt = transform(self.cell, J, detJ, f) ndofs = entity_dofs[sd-1][f][:dimP1] - tdofs = entity_dofs[sd-1][f][dimP1:] - cdofs = constraints[f] + tdofs = entity_dofs[sd-1][f][dimP1:(dimP1+dimNed1)] + cdofs = entity_dofs[sd-1][f][dimP1+dimNed1:] V[tdofs, tdofs] = Btt if sd == 2: From 64020457aaf955ff5c7ab598c784315e4e0b03d7 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Feb 2026 09:16:40 +0000 Subject: [PATCH 18/21] cleanup --- finat/mtw.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/finat/mtw.py b/finat/mtw.py index 588ebf15..c4173e31 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -1,5 +1,5 @@ import FIAT -import gem +from gem import ListTensor from finat.citations import cite from finat.fiat_elements import FiatElement @@ -19,10 +19,8 @@ def basis_transformation(self, coordinate_mapping): bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) detJ = coordinate_mapping.detJ_at(bary) - entity_dofs = self.entity_dofs() - ndof = self.space_dimension() - V = identity(ndof, ndof) + V = identity(self.space_dimension()) dimP1 = sd if sd == 2: @@ -30,6 +28,7 @@ def basis_transformation(self, coordinate_mapping): else: transform = normal_tangential_face_transform + entity_dofs = self.entity_dofs() for f in sorted(entity_dofs[sd-1]): *Bnt, Btt = transform(self.cell, J, detJ, f) ndofs = entity_dofs[sd-1][f][:dimP1] @@ -41,4 +40,4 @@ def basis_transformation(self, coordinate_mapping): V[tdofs[-1], ndofs[1:]] = Bnt V[tdofs[:-1], ndofs[0]] = Bnt - return gem.ListTensor(V.T) + return ListTensor(V.T) From cf759a96ec1e403484122d07693f717eec2d31e4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Feb 2026 09:37:40 +0000 Subject: [PATCH 19/21] tests --- test/FIAT/unit/test_fiat.py | 4 +++- test/FIAT/unit/test_johnson_mercier.py | 21 +++++++-------------- 2 files changed, 10 insertions(+), 15 deletions(-) diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index a4b9f0ea..11f44313 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -59,7 +59,7 @@ from FIAT.christiansen_hu import ChristiansenHu # noqa: F401 from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1 # noqa: F401 from FIAT.guzman_neilan import GuzmanNeilanSecondKindH1 # noqa: F401 -from FIAT.johnson_mercier import JohnsonMercier # noqa: F401 +from FIAT.johnson_mercier import JohnsonMercier, ReducedJohnsonMercier # noqa: F401 from FIAT.bubble import Bubble, FacetBubble # noqa: F401 from FIAT.enriched import EnrichedElement # noqa: F401 from FIAT.nodal_enriched import NodalEnrichedElement @@ -387,6 +387,8 @@ def __init__(self, a, b): "AlfeldC2(T, 6)", "JohnsonMercier(T)", "JohnsonMercier(S)", + "ReducedJohnsonMercier(T)", + "ReducedJohnsonMercier(S)", "AlfeldSorokina(T)", "AlfeldSorokina(S)", "ArnoldQin(T, reduced=False)", diff --git a/test/FIAT/unit/test_johnson_mercier.py b/test/FIAT/unit/test_johnson_mercier.py index dcf60e6c..5a9fc31f 100644 --- a/test/FIAT/unit/test_johnson_mercier.py +++ b/test/FIAT/unit/test_johnson_mercier.py @@ -1,7 +1,7 @@ import pytest import numpy -from FIAT import JohnsonMercier, Nedelec +from FIAT import JohnsonMercier, ReducedJohnsonMercier, Nedelec from FIAT.reference_element import ufc_simplex from FIAT.quadrature_schemes import create_quadrature @@ -20,13 +20,16 @@ def cell(request): return K -def test_johnson_mercier_divergence_rigid_body_motions(cell): +@pytest.mark.parametrize("reduced", (False, True)) +def test_johnson_mercier_divergence_rigid_body_motions(cell, reduced): # test that the divergence of interior JM basis functions is orthogonal to # the rigid-body motions degree = 1 - variant = None sd = cell.get_spatial_dimension() - JM = JohnsonMercier(cell, degree, variant=variant) + if reduced: + JM = ReducedJohnsonMercier(cell, degree) + else: + JM = JohnsonMercier(cell, degree) ref_complex = JM.get_reference_complex() Q = create_quadrature(ref_complex, 2*(degree)-1) @@ -45,13 +48,3 @@ def test_johnson_mercier_divergence_rigid_body_motions(cell): idofs = edofs[sd][0] L = numpy.tensordot(div, ells, axes=((1, 2), (1, 2))) assert numpy.allclose(L[idofs], 0) - - if variant == "divergence": - edofs = JM.entity_dofs() - cdofs = [] - for entity in edofs[sd-1]: - cdofs.extend(edofs[sd-1][entity][:sd]) - D = L[cdofs] - M = numpy.tensordot(rbms, ells, axes=((1, 2), (1, 2))) - X = numpy.linalg.solve(M, D.T) - assert numpy.allclose(numpy.tensordot(X, rbms, axes=(0, 0)), div[cdofs]) From 0cc0d305caa4acb2594fc12aa8b8be9f2baf2088 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Feb 2026 13:56:06 +0000 Subject: [PATCH 20/21] fix --- FIAT/johnson_mercier.py | 22 ++++++++++++---------- finat/johnson_mercier.py | 9 +++++---- finat/physically_mapped.py | 4 ++-- test/FIAT/unit/test_johnson_mercier.py | 13 ++++++++++--- 4 files changed, 29 insertions(+), 19 deletions(-) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index afa200a6..f7eec61b 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -94,22 +94,24 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): ref_facet = ref_el.get_facet_element() Q = parse_quadrature_scheme(ref_facet, degree+1, quad_scheme) - P1 = polynomial_set.ONPolynomialSet(ref_facet, degree) + P1 = polynomial_set.ONPolynomialSet(ref_facet, degree, scale="orthonormal") P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] dimP1 = len(P1)*(sd-1) dimNed1 = dimP1 // 2 if sd == 3: # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] RT_at_qpts = numpy.zeros((dimP1, sd-1, P1_at_qpts.shape[-1])) - RT_at_qpts[0, 0, :] = P1_at_qpts[0] - RT_at_qpts[1, 1, :] = P1_at_qpts[0] - RT_at_qpts[2, 0, :] = P1_at_qpts[1] - RT_at_qpts[2, 1, :] = P1_at_qpts[2] - # Basis for complement of RT [(y-x, 0), (0, y-x), (y, -x)] - RT_at_qpts[3, 0, :] = P1_at_qpts[2] - P1_at_qpts[1] - RT_at_qpts[4, 1, :] = P1_at_qpts[2] - P1_at_qpts[1] - RT_at_qpts[5, 0, :] = P1_at_qpts[2] - RT_at_qpts[5, 1, :] = -P1_at_qpts[1] + RT_at_qpts[0, 0] = P1_at_qpts[0] + RT_at_qpts[1, 1] = P1_at_qpts[0] + RT_at_qpts[2, 0] = P1_at_qpts[1] + RT_at_qpts[2, 1] = P1_at_qpts[2] + # Basis for the complement of RT [(y, x), (x, -y), (y, -x)] + RT_at_qpts[3, 0] = P1_at_qpts[2] + RT_at_qpts[3, 1] = P1_at_qpts[1] + RT_at_qpts[4, 0] = P1_at_qpts[1] + RT_at_qpts[4, 1] = -P1_at_qpts[2] + RT_at_qpts[5, 0] = P1_at_qpts[2] + RT_at_qpts[5, 1] = -P1_at_qpts[1] else: RT_at_qpts = P1_at_qpts[:, None, :] diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index e1402633..de6c3a0d 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -82,11 +82,12 @@ def basis_transformation(self, coordinate_mapping): V[cdofs, ndofs[1]] = Bnt else: V[tdofs[:-1], ndofs[0]] = Bnt + # n x (t0 * x + t1 * y) V[tdofs[-1], ndofs[1:]] = Bnt - # (n x t0) * (y - x) - V[cdofs[0], ndofs[1:]] = (-1*Bnt[0], Bnt[0]) - # (n x t1) * (y - x) - V[cdofs[1], ndofs[1:]] = (-1*Bnt[1], Bnt[1]) + # n x (t0 * y + t1 * x) + V[cdofs[0], ndofs[1:]] = (Bnt[1], Bnt[0]) + # n x (t0 * x - t1 * y) + V[cdofs[1], ndofs[1:]] = (Bnt[0], -1*Bnt[1]) # n x (t0 * y - t1 * x) V[cdofs[2], ndofs[1:]] = (-1*Bnt[1], Bnt[0]) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index 978bbaef..fe004440 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -38,8 +38,8 @@ def matvec(self, table): for i, js in enumerate(self.csr)] val = gem.ListTensor(exprs) - # val = self.M @ table - return gem.optimise.aggressive_unroll(val) + # val = gem.optimise.aggressive_unroll(self.M @ table) + return val def __getitem__(self, alpha): try: diff --git a/test/FIAT/unit/test_johnson_mercier.py b/test/FIAT/unit/test_johnson_mercier.py index 5a9fc31f..0328b078 100644 --- a/test/FIAT/unit/test_johnson_mercier.py +++ b/test/FIAT/unit/test_johnson_mercier.py @@ -43,8 +43,15 @@ def test_johnson_mercier_divergence_rigid_body_motions(cell, reduced): N1_at_qpts = N1.tabulate(0, qpts) rbms = N1_at_qpts[(0,)*sd] ells = rbms * qwts[None, None, :] + L = numpy.tensordot(div, ells, axes=((1, 2), (1, 2))) edofs = JM.entity_dofs() - idofs = edofs[sd][0] - L = numpy.tensordot(div, ells, axes=((1, 2), (1, 2))) - assert numpy.allclose(L[idofs], 0) + idofs = list(edofs[sd][0]) + if reduced: + dimP1 = sd + dimNed1 = (sd*(sd-1))//2 + for f in edofs[sd-1]: + idofs.extend(edofs[sd-1][f][dimP1+dimNed1:]) + + L = L[idofs].reshape(len(idofs), -1) + assert numpy.allclose(L, 0) From c12967298e4e775a1c7f61fa22f045661e2ab106 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 10 Feb 2026 18:39:46 +0000 Subject: [PATCH 21/21] WIP --- FIAT/johnson_mercier.py | 47 ++++++++++++++++++-------- test/FIAT/unit/test_johnson_mercier.py | 6 ---- 2 files changed, 33 insertions(+), 20 deletions(-) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index f7eec61b..2b20fcf4 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -66,18 +66,37 @@ def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) -def rbm_complement(ref_el): - """Constructs a basis for the complement of the rigid body motions over P1.""" - if ref_el.get_spatial_dimension() == 1: - P1 = polynomial_set.ONPolynomialSet(ref_el, 1, shape=(1,)) - return P1.take(range(1, len(P1))) - else: - P1 = NedelecSecondKind(ref_el, 1) - entity_ids = P1.entity_dofs() - ids = [] - for entity in entity_ids[1]: - ids.extend(entity_ids[1][entity][1:]) - return P1.get_nodal_basis().take(ids) +class DG0Alfeld(finite_element.CiarletElement): + + def __init__(self, ref_complex, P1=None, quad_scheme=None): + ref_el = ref_complex.get_parent() + sd = ref_el.get_spatial_dimension() + if P1 is None: + if sd == 1: + P1 = polynomial_set.ONPolynomialSet(ref_el, 1, shape=(sd,)) + else: + P1 = NedelecSecondKind(ref_el, 1) + + Q = parse_quadrature_scheme(ref_complex, P1.degree(), quad_scheme) + phis = P1.tabulate(0, Q.get_points())[(0,)*sd] + nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis] + dual = dual_set.DualSet(nodes, ref_el, P1.entity_dofs()) + + poly_set = polynomial_set.ONPolynomialSet(ref_complex, 0, shape=(sd,)) + formdegree = sd + super().__init__(poly_set, dual, poly_set.degree, formdegree) + + +def rbm_complement(ref_complex): + """Constructs a basis for the complement of the rigid body motions over P0(Alfeld).""" + P0 = DG0Alfeld(ref_complex) + P0 = DG0Alfeld(ref_complex, P1=P0) + + entity_ids = P0.entity_dofs() + ids = [] + for entity in entity_ids[1]: + ids.extend(entity_ids[1][entity][1:]) + return P0.get_nodal_basis().take(ids) class ReducedJohnsonMercierDualSet(dual_set.DualSet): @@ -147,8 +166,8 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): # Interior constraints: moments of divergence against (P1 \ Nedelec) ref_complex = macro.AlfeldSplit(ref_el) - Q = parse_quadrature_scheme(ref_complex, 2*degree-1) - comp_space = rbm_complement(ref_el) + Q = parse_quadrature_scheme(ref_complex, 2*(degree-1)) + comp_space = rbm_complement(ref_complex) phis = comp_space.tabulate(Q.get_points())[(0,)*sd] cur = len(nodes) nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) diff --git a/test/FIAT/unit/test_johnson_mercier.py b/test/FIAT/unit/test_johnson_mercier.py index 0328b078..73d8f536 100644 --- a/test/FIAT/unit/test_johnson_mercier.py +++ b/test/FIAT/unit/test_johnson_mercier.py @@ -47,11 +47,5 @@ def test_johnson_mercier_divergence_rigid_body_motions(cell, reduced): edofs = JM.entity_dofs() idofs = list(edofs[sd][0]) - if reduced: - dimP1 = sd - dimNed1 = (sd*(sd-1))//2 - for f in edofs[sd-1]: - idofs.extend(edofs[sd-1][f][dimP1+dimNed1:]) - L = L[idofs].reshape(len(idofs), -1) assert numpy.allclose(L, 0)