From b68de43a3a9b5bbc403ec8336a0ac48c0de9b3fa Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 19 Jul 2025 15:35:16 -0400 Subject: [PATCH 01/25] Morley Tetrahedron --- FIAT/morley.py | 81 +++++++++++++++++-------------------- test/FIAT/unit/test_fiat.py | 1 + 2 files changed, 39 insertions(+), 43 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index 5100a4663..b5f893cee 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -5,52 +5,43 @@ # SPDX-License-Identifier: LGPL-3.0-or-later from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import TRIANGLE, ufc_simplex +from FIAT.reference_element import TETRAHEDRON, TRIANGLE +from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature import numpy class MorleyDualSet(dual_set.DualSet): - """The dual basis for Lagrange elements. This class works for - simplices of any dimension. Nodes are point evaluation at - equispaced points.""" + """The dual basis for Morley elements. This class works for + simplices of any dimension. Nodes are constant moments on codim=2 entities + and constant moments of normal derivatives on codim=1 entities.""" - def __init__(self, ref_el): - entity_ids = {} - nodes = [] - cur = 0 - - # make nodes by getting points - # need to do this dimension-by-dimension, facet-by-facet + def __init__(self, ref_el, degree): top = ref_el.get_topology() - verts = ref_el.get_vertices() - if ref_el.get_shape() != TRIANGLE: - raise ValueError("Morley only defined on triangles") - - # vertex point evaluations - - entity_ids[0] = {} - for v in sorted(top[0]): - nodes.append(functional.PointEvaluation(ref_el, verts[v])) - - entity_ids[0][v] = [cur] - cur += 1 - - # edge dof -- average of normal derivative at each edge - rline = ufc_simplex(1) - degree = 2 - Q = create_quadrature(rline, degree-1) - qwts = Q.get_weights() - scale = numpy.ones(qwts.shape) - - entity_ids[1] = {} - for e in sorted(top[1]): - n = functional.IntegralMomentOfNormalDerivative(ref_el, e, Q, scale) - nodes.append(n) - entity_ids[1][e] = [cur] - cur += 1 - - entity_ids[2] = {0: []} + sd = ref_el.get_spatial_dimension() + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] + for dim in (sd-2, sd-1): + facet = ref_el.construct_subelement(dim) + Q_ref = create_quadrature(facet, degree-1) + scale = numpy.ones(Q_ref.get_weights().shape) + + for entity in sorted(top[dim]): + cur = len(nodes) + if dim == sd-1: + # codim=1 dof -- average of normal derivative at each facet + ell = functional.IntegralMomentOfNormalDerivative(ref_el, entity, Q_ref, scale) + elif dim == 0: + # codim=2 vertex dof -- point evaluation + pt, = ref_el.make_points(dim, entity, degree) + ell = functional.PointEvaluation(ref_el, pt) + else: + # codim=2 edge dof -- integral average + Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) + ell = functional.IntegralMoment(ref_el, Q, scale) + + nodes.append(ell) + entity_ids[dim][entity].extend(list(range(cur, len(nodes)))) super().__init__(nodes, ref_el, entity_ids) @@ -58,7 +49,11 @@ def __init__(self, ref_el): class Morley(finite_element.CiarletElement): """The Morley finite element.""" - def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(ref_el, 2) - dual = MorleyDualSet(ref_el) - super().__init__(poly_set, dual, 2) + def __init__(self, ref_el, degree=2): + if ref_el.get_shape() not in {TRIANGLE, TETRAHEDRON}: + raise ValueError("Morley only defined on simplices") + if degree != 2: + raise ValueError("{type(self).__name__} only defined for degree=2") + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) + dual = MorleyDualSet(ref_el, degree) + super().__init__(poly_set, dual, degree) diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index 745518a75..bd4a2c1c2 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -350,6 +350,7 @@ def __init__(self, a, b): "CubicHermite(T)", "CubicHermite(S)", "Morley(T)", + "Morley(S)", "BernardiRaugel(T)", "BernardiRaugel(S)", "MardalTaiWinther(T, 3)", From 819f6367175b56f3de86df698ecf77b344773f27 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 19 Jul 2025 16:41:23 -0400 Subject: [PATCH 02/25] Dimension independent zany map --- finat/morley.py | 42 +++++++++++++++++++++++------------------- 1 file changed, 23 insertions(+), 19 deletions(-) diff --git a/finat/morley.py b/finat/morley.py index 986a4471b..a69f538cb 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -1,6 +1,6 @@ import FIAT -from gem import ListTensor +from gem import ListTensor, partial_indexed from finat.fiat_elements import ScalarFiatElement from finat.physically_mapped import Citations, identity, PhysicallyMappedElement @@ -10,37 +10,41 @@ class Morley(PhysicallyMappedElement, ScalarFiatElement): def __init__(self, cell, degree=2): if Citations is not None: Citations().register("Morley1971") - super().__init__(FIAT.Morley(cell)) + super().__init__(FIAT.Morley(cell, degree=degree)) def basis_transformation(self, coordinate_mapping): + sd = self.cell.get_spatial_dimension() + top = self.cell.get_topology() # Jacobians at edge midpoints - J = coordinate_mapping.jacobian_at([1/3, 1/3]) + bary, = self.cell.make_points(sd, 0, sd+1) + J = coordinate_mapping.jacobian_at(bary) rns = coordinate_mapping.reference_normals() pns = coordinate_mapping.physical_normals() - pts = coordinate_mapping.physical_tangents() - pel = coordinate_mapping.physical_edge_lengths() V = identity(self.space_dimension()) - for i in range(3): - V[i+3, i+3] = (rns[i, 0]*(pns[i, 0]*J[0, 0] + pns[i, 1]*J[1, 0]) - + rns[i, 1]*(pns[i, 0]*J[0, 1] + pns[i, 1]*J[1, 1])) - - for i, c in enumerate([(1, 2), (0, 2), (0, 1)]): - B12 = (rns[i, 0]*(pts[i, 0]*J[0, 0] + pts[i, 1]*J[1, 0]) - + rns[i, 1]*(pts[i, 0]*J[0, 1] + pts[i, 1]*J[1, 1])) - V[3+i, c[0]] = -1*B12 / pel[i] - V[3+i, c[1]] = B12 / pel[i] + offset = len(top[sd-2]) + for i in top[sd-1]: + nhat = partial_indexed(rns, (i,)) + n = partial_indexed(pns, (i,)) + t = partial_indexed(pts, (i,)) + Bn = J @ nhat + Bnn = n @ Bn + Btn = t @ Bn + + c = list(top[sd-1][i]) + s = i + offset + V[s, s] = Bnn + V[s, c] = Btn / pel[i] + V[s, c[0]] *= -1 # diagonal post-scaling to patch up conditioning h = coordinate_mapping.cell_size() - - for e in range(3): - v0id, v1id = [i for i in range(3) if i != e] - for i in range(6): - V[i, 3+e] = 2*V[i, 3+e] / (h[v0id] + h[v1id]) + for i in top[sd-1]: + s = i + offset + V[:, s] *= 2 / sum(h[v] for v in top[sd-1][i]) return ListTensor(V.T) From 62d933d1b906f276b7311b5a82641b89d6dd6812 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 19 Jul 2025 16:55:46 -0400 Subject: [PATCH 03/25] Failing zany test --- FIAT/morley.py | 6 +++--- finat/morley.py | 3 ++- test/finat/test_zany_mapping.py | 11 +++++++++-- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index b5f893cee..c9f30ca0c 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -13,8 +13,8 @@ class MorleyDualSet(dual_set.DualSet): """The dual basis for Morley elements. This class works for - simplices of any dimension. Nodes are constant moments on codim=2 entities - and constant moments of normal derivatives on codim=1 entities.""" + simplices of any dimension. Nodes are average on codim=2 entities + and average of normal derivative on codim=1 entities.""" def __init__(self, ref_el, degree): top = ref_el.get_topology() @@ -38,7 +38,7 @@ def __init__(self, ref_el, degree): else: # codim=2 edge dof -- integral average Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) - ell = functional.IntegralMoment(ref_el, Q, scale) + ell = functional.IntegralMoment(ref_el, Q, scale / Q.jacobian_determinant()) nodes.append(ell) entity_ids[dim][entity].extend(list(range(cur, len(nodes)))) diff --git a/finat/morley.py b/finat/morley.py index a69f538cb..9e90032f3 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -15,6 +15,7 @@ def __init__(self, cell, degree=2): def basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() + conn = self.cell.get_connectivity() # Jacobians at edge midpoints bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) @@ -35,8 +36,8 @@ def basis_transformation(self, coordinate_mapping): Bnn = n @ Bn Btn = t @ Bn - c = list(top[sd-1][i]) s = i + offset + c = list(conn[(sd-1, sd-2)][i]) V[s, s] = Bnn V[s, c] = Btn / pel[i] V[s, c[0]] *= -1 diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index d5136e651..7841cb807 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -82,7 +82,7 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): Mh[abs(Mh) < 1E-10] = 0.0 M[abs(M) < 1E-10] = 0.0 assert np.allclose(residual, 0), str(M.T - Mh.T) - assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]) + assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(M.T - Mh.T) @pytest.mark.parametrize("element", [ @@ -90,10 +90,17 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): finat.Hermite, finat.Bell, ]) -def test_C1_elements(ref_to_phys, element): +def test_C1_triangle(ref_to_phys, element): check_zany_mapping(element, ref_to_phys[2]) +@pytest.mark.parametrize("element", [ + finat.Morley, + ]) +def test_C1_tetrahedron(ref_to_phys, element): + check_zany_mapping(element, ref_to_phys[3]) + + @pytest.mark.parametrize("element", [ finat.QuadraticPowellSabin6, finat.QuadraticPowellSabin12, From 51e3930d34f6826a5b6bfae094394c43b305c5af Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 19 Jul 2025 17:09:42 -0400 Subject: [PATCH 04/25] Fix quadrature degree --- FIAT/morley.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index c9f30ca0c..c46accdef 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -21,9 +21,10 @@ def __init__(self, ref_el, degree): sd = ref_el.get_spatial_dimension() entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} nodes = [] - for dim in (sd-2, sd-1): + for codim in (2, 1): + dim = sd - codim facet = ref_el.construct_subelement(dim) - Q_ref = create_quadrature(facet, degree-1) + Q_ref = create_quadrature(facet, degree+codim-2) scale = numpy.ones(Q_ref.get_weights().shape) for entity in sorted(top[dim]): From dc9edc148bd9c535fb3bf712e81a2253e0b7b364 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 21 Jul 2025 11:31:29 +0100 Subject: [PATCH 05/25] Clean up Bell --- FIAT/bell.py | 71 ++++++++++++++++++++-------------------------------- 1 file changed, 27 insertions(+), 44 deletions(-) diff --git a/FIAT/bell.py b/FIAT/bell.py index 1d3fe11b3..e746078e4 100644 --- a/FIAT/bell.py +++ b/FIAT/bell.py @@ -10,59 +10,38 @@ from FIAT import finite_element, polynomial_set, dual_set, functional from FIAT.reference_element import TRIANGLE, ufc_simplex +from FIAT.quadrature_schemes import create_quadrature +from FIAT.jacobi import eval_jacobi class BellDualSet(dual_set.DualSet): - def __init__(self, ref_el): - entity_ids = {} - nodes = [] - cur = 0 - - # make nodes by getting points - # need to do this dimension-by-dimension, facet-by-facet + def __init__(self, ref_el, degree): top = ref_el.get_topology() - verts = ref_el.get_vertices() sd = ref_el.get_spatial_dimension() - if ref_el.get_shape() != TRIANGLE: - raise ValueError("Bell only defined on triangles") - - pd = functional.PointDerivative + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] # get jet at each vertex - - entity_ids[0] = {} for v in sorted(top[0]): - nodes.append(functional.PointEvaluation(ref_el, verts[v])) - - # first derivatives - for i in range(sd): - alpha = [0] * sd - alpha[i] = 1 - nodes.append(pd(ref_el, verts[v], alpha)) - - # second derivatives - alphas = [[2, 0], [1, 1], [0, 2]] - for alpha in alphas: - nodes.append(pd(ref_el, verts[v], alpha)) + cur = len(nodes) + x, = ref_el.make_points(0, v, degree) + nodes.append(functional.PointEvaluation(ref_el, x)) - entity_ids[0][v] = list(range(cur, cur + 6)) - cur += 6 + # first and second derivatives + nodes.extend(functional.PointDerivative(ref_el, x, alpha) + for i in (1, 2) for alpha in polynomial_set.mis(sd, i)) + entity_ids[0][v].extend(range(cur, len(nodes))) # we need an edge quadrature rule for the moment - from FIAT.quadrature_schemes import create_quadrature - from FIAT.jacobi import eval_jacobi - rline = ufc_simplex(1) - q1d = create_quadrature(rline, 8) - q1dpts = q1d.get_points()[:, 0] - leg4_at_qpts = eval_jacobi(0, 0, 4, 2.0*q1dpts - 1) + facet = ufc_simplex(1) + Q_ref = create_quadrature(facet, 8) + qpts = Q_ref.get_points()[:, 0] + leg4_at_qpts = eval_jacobi(0, 0, 4, 2.0*qpts - 1) - imond = functional.IntegralMomentOfNormalDerivative - entity_ids[1] = {} for e in sorted(top[1]): - entity_ids[1][e] = [18+e] - nodes.append(imond(ref_el, e, q1d, leg4_at_qpts)) - - entity_ids[2] = {0: []} + cur = len(nodes) + nodes.append(functional.IntegralMomentOfNormalDerivative(ref_el, e, Q_ref, leg4_at_qpts)) + entity_ids[1][e].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -70,7 +49,11 @@ def __init__(self, ref_el): class Bell(finite_element.CiarletElement): """The Bell finite element.""" - def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(ref_el, 5) - dual = BellDualSet(ref_el) - super().__init__(poly_set, dual, 5) + def __init__(self, ref_el, degree=5): + if ref_el.get_shape() != TRIANGLE: + raise ValueError(f"{type(self).__name__} only defined on triangles") + if degree != 5: + raise ValueError(f"{type(self).__name__} only defined for degree=5.") + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) + dual = BellDualSet(ref_el, degree) + super().__init__(poly_set, dual, degree) From cbae4f180fe8bc1433a8cb085bb47ac38edd87ea Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 21 Jul 2025 18:04:10 +0100 Subject: [PATCH 06/25] Walkington C1 macroelement --- FIAT/__init__.py | 4 +- FIAT/bell.py | 12 ++-- FIAT/functional.py | 5 +- FIAT/walkington.py | 106 ++++++++++++++++++++++++++++++++ finat/__init__.py | 1 + finat/argyris.py | 14 +++-- finat/bell.py | 2 +- finat/walkington.py | 90 +++++++++++++++++++++++++++ test/finat/test_zany_mapping.py | 12 ++-- 9 files changed, 228 insertions(+), 18 deletions(-) create mode 100644 FIAT/walkington.py create mode 100644 finat/walkington.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 364cf4297..fcfdea88a 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -62,6 +62,7 @@ from FIAT.bubble import Bubble, FacetBubble from FIAT.hdiv_trace import HDivTrace from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen +from FIAT.walkington import Walkington from FIAT.histopolation import Histopolation from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite # noqa: F401 @@ -117,7 +118,8 @@ "Conforming Arnold-Winther": ArnoldWinther, "Nonconforming Arnold-Winther": ArnoldWintherNC, "Hu-Zhang": HuZhang, - "Mardal-Tai-Winther": MardalTaiWinther} + "Mardal-Tai-Winther": MardalTaiWinther, + "Walkington": Walkington} # List of extra elements extra_elements = {"P0": P0} diff --git a/FIAT/bell.py b/FIAT/bell.py index e746078e4..8fccef82e 100644 --- a/FIAT/bell.py +++ b/FIAT/bell.py @@ -9,7 +9,7 @@ # bfs, but the extra three are used in the transformation theory. from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import TRIANGLE, ufc_simplex +from FIAT.reference_element import TRIANGLE from FIAT.quadrature_schemes import create_quadrature from FIAT.jacobi import eval_jacobi @@ -33,10 +33,10 @@ def __init__(self, ref_el, degree): entity_ids[0][v].extend(range(cur, len(nodes))) # we need an edge quadrature rule for the moment - facet = ufc_simplex(1) - Q_ref = create_quadrature(facet, 8) - qpts = Q_ref.get_points()[:, 0] - leg4_at_qpts = eval_jacobi(0, 0, 4, 2.0*qpts - 1) + facet = ref_el.construct_subelement(1) + Q_ref = create_quadrature(facet, 2*(degree-1)) + x = facet.compute_barycentric_coordinates(Q_ref.get_points()) + leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0]) for e in sorted(top[1]): cur = len(nodes) @@ -53,7 +53,7 @@ def __init__(self, ref_el, degree=5): if ref_el.get_shape() != TRIANGLE: raise ValueError(f"{type(self).__name__} only defined on triangles") if degree != 5: - raise ValueError(f"{type(self).__name__} only defined for degree=5.") + raise ValueError(f"{type(self).__name__} only defined for degree = 5.") poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = BellDualSet(ref_el, degree) super().__init__(poly_set, dual, degree) diff --git a/FIAT/functional.py b/FIAT/functional.py index 18de29942..18ab6e28a 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -349,8 +349,9 @@ def __init__(self, ref_el, s, Q, f_at_qpts, comp=(), shp=()): class IntegralMomentOfNormalDerivative(Functional): """Functional giving normal derivative integrated against some function on a facet.""" - def __init__(self, ref_el, facet_no, Q, f_at_qpts): - n = ref_el.compute_normal(facet_no) + def __init__(self, ref_el, facet_no, Q, f_at_qpts, n=None): + if n is None: + n = ref_el.compute_normal(facet_no) self.n = n self.f_at_qpts = f_at_qpts self.Q = Q diff --git a/FIAT/walkington.py b/FIAT/walkington.py new file mode 100644 index 000000000..69b02b626 --- /dev/null +++ b/FIAT/walkington.py @@ -0,0 +1,106 @@ +# Copyright (C) 2025 Pablo D. Brubeck +# +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +# This is not quite Bell, but is 65-dofs and includes 20 extra constraint +# functionals. The first 45 basis functions are the reference element +# bfs, but the extra three are used in the transformation theory. + +from FIAT import finite_element, polynomial_set, dual_set, functional, macro +from FIAT.reference_element import TETRAHEDRON +from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature +from FIAT.jacobi import eval_jacobi +from FIAT.hierarchical import make_dual_bubbles +import numpy + + +class WalkingtonDualSet(dual_set.DualSet): + def __init__(self, ref_el, degree): + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] + + # Vertex dofs: second order jet + for v in sorted(top[0]): + cur = len(nodes) + x, = ref_el.make_points(0, v, degree) + nodes.append(functional.PointEvaluation(ref_el, x)) + + # first and second derivatives + nodes.extend(functional.PointDerivative(ref_el, x, alpha) + for i in (1, 2) for alpha in polynomial_set.mis(sd, i)) + entity_ids[0][v].extend(range(cur, len(nodes))) + + # Face dofs: moments against cubic bubble + ref_face = ref_el.construct_subelement(2) + Q_face, phis = make_dual_bubbles(ref_face, degree-1) + + for face in sorted(top[sd-1]): + cur = len(nodes) + nface = ref_el.compute_scaled_normal(face) + nface /= numpy.dot(nface, nface) + + nodes.append(functional.IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phis[0], n=nface)) + entity_ids[sd-1][face].extend(range(cur, len(nodes))) + + # Interior dof: point evaluation at barycenter + for entity in top[sd]: + cur = len(nodes) + x, = ref_el.make_points(sd, entity, sd+1) + nodes.append(functional.PointEvaluation(ref_el, x)) + entity_ids[sd][entity].extend(range(cur, len(nodes))) + + # Constraint dofs + # Face constraint: normal derivative is cubic + # Face-edge constraint: normal derivative along edge is cubic + edges = ref_el.get_connectivity()[(2, 1)] + ref_edge = ref_el.construct_subelement(1) + Q_edge = create_quadrature(ref_edge, 2*(degree-1)) + x = ref_edge.compute_barycentric_coordinates(Q_edge.get_points()) + leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0]) + + for face in sorted(top[2]): + cur = len(nodes) + nface = ref_el.compute_scaled_normal(face) + nface /= numpy.dot(nface, nface) + + for i, e in enumerate(edges[face]): + Q = FacetQuadratureRule(ref_face, 1, i, Q_edge) + + te = ref_el.compute_edge_tangent(e) + nfe = numpy.cross(te, nface) + nfe /= numpy.linalg.norm(nfe) + nfe /= Q.jacobian_determinant() + + nodes.append(functional.IntegralMomentOfNormalDerivative(ref_el, face, Q, leg4_at_qpts, n=nfe)) + + nodes.extend(functional.IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phi, n=nface) for phi in phis[1:]) + entity_ids[sd-1][face].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) + + +class Walkington(finite_element.CiarletElement): + """The Walkington C1 macroelement.""" + + def __init__(self, ref_el, degree=5): + if ref_el.get_shape() != TETRAHEDRON: + raise ValueError(f"{type(self).__name__} only defined on tetrahedron") + if degree != 5: + raise ValueError(f"{type(self).__name__} only defined for degree = 5.") + + dual = WalkingtonDualSet(ref_el, degree) + ref_complex = macro.AlfeldSplit(ref_el) + poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=4, variant="bubble") + super().__init__(poly_set, dual, degree) + + +if __name__ == "__main__": + from FIAT import ufc_simplex + cell = ufc_simplex(3) + fe = Walkington(cell) + print(fe.entity_dofs()) diff --git a/finat/__init__.py b/finat/__init__.py index 5390fa779..30caef724 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -25,6 +25,7 @@ from .johnson_mercier import JohnsonMercier # noqa: F401 from .mtw import MardalTaiWinther # noqa: F401 from .morley import Morley # noqa: F401 +from .walkington import Walkington # noqa: F401 from .direct_serendipity import DirectSerendipity # noqa: F401 from .spectral import GaussLobattoLegendre, GaussLegendre, Legendre, IntegratedLegendre, FDMLagrange, FDMQuadrature, FDMDiscontinuousLagrange, FDMBrokenH1, FDMBrokenL2, FDMHermite # noqa: F401 from .tensorfiniteelement import TensorFiniteElement # noqa: F401 diff --git a/finat/argyris.py b/finat/argyris.py index 562ae41b2..2a7edc6e5 100644 --- a/finat/argyris.py +++ b/finat/argyris.py @@ -41,10 +41,16 @@ def _vertex_transform(V, vorder, fiat_cell, coordinate_mapping): return V -def _normal_tangential_transform(fiat_cell, J, detJ, f): - R = numpy.array([[0, 1], [-1, 0]]) - that = fiat_cell.compute_edge_tangent(f) - nhat = R @ that +def _normal_tangential_transform(fiat_cell, J, detJ, edge, face=None): + that = fiat_cell.compute_edge_tangent(edge) + if fiat_cell.get_spatial_dimension() == 2: + R = numpy.array([[0, 1], [-1, 0]]) + nhat = R @ that + else: + nface = fiat_cell.compute_scaled_normal(face) + nface /= numpy.linalg.norm(nface) + nhat = numpy.cross(that, nface) + Jn = J @ Literal(nhat) Jt = J @ Literal(that) alpha = Jn @ Jt diff --git a/finat/bell.py b/finat/bell.py index 98d4aecc6..46aaac2a5 100644 --- a/finat/bell.py +++ b/finat/bell.py @@ -12,7 +12,7 @@ class Bell(PhysicallyMappedElement, ScalarFiatElement): def __init__(self, cell, degree=5): if Citations is not None: Citations().register("Bell1969") - super().__init__(FIAT.Bell(cell)) + super().__init__(FIAT.Bell(cell, degree=degree)) reduced_dofs = deepcopy(self._element.entity_dofs()) sd = cell.get_spatial_dimension() diff --git a/finat/walkington.py b/finat/walkington.py new file mode 100644 index 000000000..527fdb880 --- /dev/null +++ b/finat/walkington.py @@ -0,0 +1,90 @@ +import FIAT +from FIAT.polynomial_set import mis +from math import comb + +from gem import ListTensor + +from finat.fiat_elements import ScalarFiatElement +from finat.physically_mapped import Citations, identity, PhysicallyMappedElement +from finat.argyris import _vertex_transform, _normal_tangential_transform +from copy import deepcopy + + +class Walkington(PhysicallyMappedElement, ScalarFiatElement): + def __init__(self, cell, degree=5): + if Citations is not None: + Citations().register("Bell1969") + super().__init__(FIAT.Walkington(cell, degree=degree)) + + reduced_dofs = deepcopy(self._element.entity_dofs()) + sd = cell.get_spatial_dimension() + for entity in reduced_dofs[sd-1]: + reduced_dofs[sd-1][entity] = reduced_dofs[sd-1][entity][:1] + self._entity_dofs = reduced_dofs + + def basis_transformation(self, coordinate_mapping): + # Jacobian at barycenter + sd = self.cell.get_spatial_dimension() + top = self.cell.get_topology() + bary, = self.cell.make_points(sd, 0, sd+1) + J = coordinate_mapping.jacobian_at(bary) + detJ = coordinate_mapping.detJ_at(bary) + + numbf = self._element.space_dimension() + ndof = self.space_dimension() + # rectangular to toss out the constraint dofs + V = identity(numbf, ndof) + + vorder = 2 + _vertex_transform(V, vorder, self.cell, coordinate_mapping) + + offset = ndof + voffset = comb(sd + vorder, vorder) + foffset = len(self._element.entity_dofs()[sd-1][0]) - len(self._entity_dofs[sd-1][0]) + + edges = self.cell.get_connectivity()[(2, 1)] + for f in sorted(top[2]): + q = len(top[0]) * voffset + f + V[q, q] *= detJ + + for j, e in enumerate(edges[f]): + s = offset + foffset * f + j + + v0id, v1id = (v * voffset for v in top[1][e]) + Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e, face=f) + + # vertex points + V[s, v1id] = 1/21 * Bnt + V[s, v0id] = -1 * V[s, v1id] + + # vertex derivatives + for i in range(sd): + V[s, v1id+1+i] = -1/42 * Bnt * Jt[i] + V[s, v0id+1+i] = V[s, v1id+1+i] + + # second derivatives + for i, alpha in enumerate(mis(sd, 2)): + ids = tuple(k for k, ak in enumerate(alpha) if ak) + a, b = ids[0], ids[-1] + tau = (1 + (a != b)) * Jt[a] * Jt[b] + V[s, v1id+sd+1+i] = 1/252 * Bnt * tau + V[s, v0id+sd+1+i] = -1 * V[s, v1id+sd+1+i] + + # Patch up conditioning + h = coordinate_mapping.cell_size() + for v in sorted(top[0]): + s = voffset * v + 1 + V[:, s:s+sd] *= 1/h[v] + V[:, s+sd:voffset*(v+1)] *= 1/(h[v]*h[v]) + + return ListTensor(V.T) + + # This wipes out the edge dofs. FIAT gives a 65 DOF element + # because we need some extra functions to help with transforming + # under the edge constraint. However, we only have an 45 DOF + # element. + def entity_dofs(self): + return self._entity_dofs + + def space_dimension(self): + return 45 diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 7841cb807..c48e75d40 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -79,10 +79,13 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): Vh, residual, *_ = np.linalg.lstsq(Phi.T, phi.T) Mh = Vh.T Mh = Mh[:num_dofs] - Mh[abs(Mh) < 1E-10] = 0.0 - M[abs(M) < 1E-10] = 0.0 - assert np.allclose(residual, 0), str(M.T - Mh.T) - assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(M.T - Mh.T) + Mh[abs(Mh) < 1E-10] = 0 + M[abs(M) < 1E-10] = 0 + + error = M - Mh + error[abs(error) < 1E-10] = 0 + assert np.allclose(residual, 0), str(error.T.tolist()) + assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(error.T.tolist()) @pytest.mark.parametrize("element", [ @@ -96,6 +99,7 @@ def test_C1_triangle(ref_to_phys, element): @pytest.mark.parametrize("element", [ finat.Morley, + finat.Walkington, ]) def test_C1_tetrahedron(ref_to_phys, element): check_zany_mapping(element, ref_to_phys[3]) From 27d688abc433b9a1ec805660aa5b7f31aee0c0ae Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 22 Jul 2025 15:41:31 +0100 Subject: [PATCH 07/25] Correct basis_transformation for Morley --- finat/morley.py | 64 ++++++++++++++++++++++++--------- test/finat/test_zany_mapping.py | 13 ++++--- 2 files changed, 56 insertions(+), 21 deletions(-) diff --git a/finat/morley.py b/finat/morley.py index 9e90032f3..c8e4efd95 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -1,6 +1,7 @@ import FIAT +import numpy -from gem import ListTensor, partial_indexed +from gem import ListTensor, partial_indexed, Literal from finat.fiat_elements import ScalarFiatElement from finat.physically_mapped import Citations, identity, PhysicallyMappedElement @@ -15,12 +16,11 @@ def __init__(self, cell, degree=2): def basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() - conn = self.cell.get_connectivity() - # Jacobians at edge midpoints + edges = self.cell.get_connectivity()[(sd-1, sd-2)] + # Jacobians at barycenter bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) - - rns = coordinate_mapping.reference_normals() + detJ = coordinate_mapping.detJ_at(bary) pns = coordinate_mapping.physical_normals() pts = coordinate_mapping.physical_tangents() pel = coordinate_mapping.physical_edge_lengths() @@ -28,19 +28,49 @@ def basis_transformation(self, coordinate_mapping): V = identity(self.space_dimension()) offset = len(top[sd-2]) - for i in top[sd-1]: - nhat = partial_indexed(rns, (i,)) - n = partial_indexed(pns, (i,)) - t = partial_indexed(pts, (i,)) - Bn = J @ nhat - Bnn = n @ Bn - Btn = t @ Bn - s = i + offset - c = list(conn[(sd-1, sd-2)][i]) - V[s, s] = Bnn - V[s, c] = Btn / pel[i] - V[s, c[0]] *= -1 + if sd == 2: + for i in top[sd-1]: + e = list(edges[i]) + s = offset + i + n = partial_indexed(pns, (i,)) + t = partial_indexed(pts, (i,)) + nhat = self.cell.compute_normal(i) + Jn = J @ Literal(nhat) + Bnn = n @ Jn + Bnt = t @ Jn + V[s, s] = Bnn + V[s, e] = Bnt / pel[i] + V[s, e[0]] *= -1 + + else: + for face in top[sd-1]: + s = offset + face + n = partial_indexed(pns, (face,)) + + nhat = self.cell.compute_scaled_normal(face) + ahat = numpy.linalg.norm(nhat) + nhat /= ahat + Jn = J @ Literal(nhat) + Bnn = n @ Jn + # Why 3/4? + V[s, s] = (3/4) * Bnn + + # Not sure where this comes from + factor = Bnn / (2 * ahat * detJ) + + xf, = self.cell.make_points(2, face, 3) + for edge in edges[face]: + t = partial_indexed(pts, (edge,)) + + xe, = self.cell.make_points(1, edge, 2) + out = numpy.array(xe) - numpy.array(xf) + that = self.cell.compute_edge_tangent(edge) + sign = numpy.sign(numpy.dot(numpy.cross(nhat, that), out)) + + # Bnt = dot(cross(n, t), J * nhat) + Bnt = sum((n[i]*t[j] - n[j]*t[i]) * Jn[k] for i, j, k in [(1, 2, 0), (2, 0, 1), (0, 1, 2)]) + V[s, edge] = Bnt * pel[edge] * sign * factor # diagonal post-scaling to patch up conditioning h = coordinate_mapping.cell_size() diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 7841cb807..94cef6469 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -79,10 +79,15 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): Vh, residual, *_ = np.linalg.lstsq(Phi.T, phi.T) Mh = Vh.T Mh = Mh[:num_dofs] - Mh[abs(Mh) < 1E-10] = 0.0 - M[abs(M) < 1E-10] = 0.0 - assert np.allclose(residual, 0), str(M.T - Mh.T) - assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(M.T - Mh.T) + Mh[abs(Mh) < 1E-10] = 0 + M[abs(M) < 1E-10] = 0 + + error = Mh.T - M.T + error[abs(error) < 1E-10] = 0 + error = error[np.ix_(*map(np.unique, np.nonzero(error)))] + + assert np.allclose(residual, 0), str(error) + assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(error) @pytest.mark.parametrize("element", [ From 051254c5869514524caca877acd5bb9098c0e494 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 22 Jul 2025 15:54:55 +0100 Subject: [PATCH 08/25] Register Morley 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 efebedb59..d0ee50513 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -115,7 +115,7 @@ def show_elements(): 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",)) -register_element("Morley", "MOR", 0, H2, "custom", (2, 2), ("triangle",)) +register_element("Morley", "MOR", 0, H2, "custom", (2, 2), simplices[1:]) # Macro elements register_element("QuadraticPowellSabin6", "PS6", 0, H2, "custom", (2, 2), ("triangle",)) From 3cf0128306a5623110ad73e166f83a449588580b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 23 Jul 2025 18:39:02 +0100 Subject: [PATCH 09/25] Fix factor of 3 --- finat/morley.py | 11 +++++++---- test/finat/test_zany_mapping.py | 5 ++++- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/finat/morley.py b/finat/morley.py index c8e4efd95..fe0d07af3 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -51,13 +51,14 @@ def basis_transformation(self, coordinate_mapping): nhat = self.cell.compute_scaled_normal(face) ahat = numpy.linalg.norm(nhat) nhat /= ahat + Jn = J @ Literal(nhat) Bnn = n @ Jn - # Why 3/4? - V[s, s] = (3/4) * Bnn + # Why 1/2? + V[s, s] = 0.5*Bnn # Not sure where this comes from - factor = Bnn / (2 * ahat * detJ) + factor = 0.5*Bnn / (ahat * detJ) xf, = self.cell.make_points(2, face, 3) for edge in edges[face]: @@ -76,6 +77,8 @@ def basis_transformation(self, coordinate_mapping): h = coordinate_mapping.cell_size() for i in top[sd-1]: s = i + offset - V[:, s] *= 2 / sum(h[v] for v in top[sd-1][i]) + verts = top[sd-1][i] + havg = sum(h[v] for v in verts) / len(verts) + V[:, s] *= havg ** -(sd-1) return ListTensor(V.T) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 94cef6469..dca4dc933 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -82,9 +82,12 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): Mh[abs(Mh) < 1E-10] = 0 M[abs(M) < 1E-10] = 0 - error = Mh.T - M.T + with np.errstate(divide='ignore', invalid='ignore'): + error = M.T / Mh.T - 1 + error[error != error] = 0 error[abs(error) < 1E-10] = 0 error = error[np.ix_(*map(np.unique, np.nonzero(error)))] + error[error != 0] += 1 assert np.allclose(residual, 0), str(error) assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(error) From da6a3c72f78cebc7c3434023726717d5d3b6917a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 23 Jul 2025 18:50:00 +0100 Subject: [PATCH 10/25] Fix for unimplemented power --- finat/morley.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/finat/morley.py b/finat/morley.py index fe0d07af3..b588acc9b 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -79,6 +79,7 @@ def basis_transformation(self, coordinate_mapping): s = i + offset verts = top[sd-1][i] havg = sum(h[v] for v in verts) / len(verts) - V[:, s] *= havg ** -(sd-1) + for i in range(sd-1): + V[:, s] *= 1/havg return ListTensor(V.T) From 932ef8a39185598fa832cafdbe7105d14dca53ca Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 24 Jul 2025 13:25:01 +0100 Subject: [PATCH 11/25] Index-free zany map --- FIAT/morley.py | 6 +---- finat/morley.py | 65 ++++++++++++++++++++++++++----------------------- 2 files changed, 36 insertions(+), 35 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index c46accdef..ad4dd92b2 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -32,12 +32,8 @@ def __init__(self, ref_el, degree): if dim == sd-1: # codim=1 dof -- average of normal derivative at each facet ell = functional.IntegralMomentOfNormalDerivative(ref_el, entity, Q_ref, scale) - elif dim == 0: - # codim=2 vertex dof -- point evaluation - pt, = ref_el.make_points(dim, entity, degree) - ell = functional.PointEvaluation(ref_el, pt) else: - # codim=2 edge dof -- integral average + # codim=2 dof -- integral average Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) ell = functional.IntegralMoment(ref_el, Q, scale / Q.jacobian_determinant()) diff --git a/finat/morley.py b/finat/morley.py index b588acc9b..5cd00a856 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -1,7 +1,7 @@ import FIAT import numpy -from gem import ListTensor, partial_indexed, Literal +from gem import ListTensor, partial_indexed, Literal, Power, Zero from finat.fiat_elements import ScalarFiatElement from finat.physically_mapped import Citations, identity, PhysicallyMappedElement @@ -16,7 +16,6 @@ def __init__(self, cell, degree=2): def basis_transformation(self, coordinate_mapping): sd = self.cell.get_spatial_dimension() top = self.cell.get_topology() - edges = self.cell.get_connectivity()[(sd-1, sd-2)] # Jacobians at barycenter bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) @@ -30,48 +29,54 @@ def basis_transformation(self, coordinate_mapping): offset = len(top[sd-2]) if sd == 2: - for i in top[sd-1]: - e = list(edges[i]) - s = offset + i - n = partial_indexed(pns, (i,)) - t = partial_indexed(pts, (i,)) - nhat = self.cell.compute_normal(i) + for e in top[sd-1]: + s = offset + e + n = partial_indexed(pns, (e,)) + t = partial_indexed(pts, (e,)) + nhat = self.cell.compute_normal(e) Jn = J @ Literal(nhat) - Bnn = n @ Jn - Bnt = t @ Jn + Bnn = Jn @ n + Bnt = Jn @ t V[s, s] = Bnn - V[s, e] = Bnt / pel[i] - V[s, e[0]] *= -1 + v = list(top[sd-1][e]) + V[s, v] = Bnt / pel[e] + V[s, v[0]] *= -1 else: + R = ListTensor([[Zero(), Literal(1)], [Literal(-1), Zero()]]) + edges = self.cell.get_connectivity()[(sd-1, sd-2)] for face in top[sd-1]: s = offset + face - n = partial_indexed(pns, (face,)) - nhat = self.cell.compute_scaled_normal(face) + te = numpy.array(list(map(self.cell.compute_edge_tangent, edges[face]))) + nhat = -numpy.cross(*te[:2]) ahat = numpy.linalg.norm(nhat) - nhat /= ahat + # Reciprocal basis + nhat /= numpy.dot(nhat, nhat) + thats = numpy.array([numpy.cross(te[0], nhat), numpy.cross(te[1], nhat)]) + + Jt = J @ Literal(thats.T) Jn = J @ Literal(nhat) - Bnn = n @ Jn - # Why 1/2? - V[s, s] = 0.5*Bnn + Jte = J @ Literal(te.T) - # Not sure where this comes from - factor = 0.5*Bnn / (ahat * detJ) + A = Jte.T @ Jte + area = Power(A[0, 0]*A[1, 1] - A[0, 1]*A[1, 0], Literal(0.5)) - xf, = self.cell.make_points(2, face, 3) - for edge in edges[face]: - t = partial_indexed(pts, (edge,)) + detF = area / ahat + Bnn = detJ / detF + V[s, s] = Bnn - xe, = self.cell.make_points(1, edge, 2) - out = numpy.array(xe) - numpy.array(xf) - that = self.cell.compute_edge_tangent(edge) - sign = numpy.sign(numpy.dot(numpy.cross(nhat, that), out)) + Gte = Jt.T @ Jte + Gnt = Jn @ Jt + Bne = (Gnt @ R) @ Gte - # Bnt = dot(cross(n, t), J * nhat) - Bnt = sum((n[i]*t[j] - n[j]*t[i]) * Jn[k] for i, j, k in [(1, 2, 0), (2, 0, 1), (0, 1, 2)]) - V[s, edge] = Bnt * pel[edge] * sign * factor + # Not sure where this comes from + factor = 2*ahat / (detF * detF) + for i, edge in enumerate(edges[face]): + # UFC convention alternates signs + sign = (-1) ** i + V[s, edge] = Bne[i] * sign * factor # diagonal post-scaling to patch up conditioning h = coordinate_mapping.cell_size() From 4a745792ae776ed6e593107119b29d5ca759b2ad Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 24 Jul 2025 14:56:13 +0100 Subject: [PATCH 12/25] Test --- test/finat/test_zany_mapping.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index dca4dc933..956480c8b 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -104,6 +104,7 @@ def test_C1_triangle(ref_to_phys, element): @pytest.mark.parametrize("element", [ finat.Morley, + finat.Walkington, ]) def test_C1_tetrahedron(ref_to_phys, element): check_zany_mapping(element, ref_to_phys[3]) From 5af9ae630988a8067b6a5df6ffcf690c4e5b51a3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 24 Jul 2025 16:29:18 +0100 Subject: [PATCH 13/25] Cleanup, citation --- FIAT/morley.py | 4 +++ finat/morley.py | 59 ++++++++++++++------------------------ finat/physically_mapped.py | 13 +++++++++ 3 files changed, 39 insertions(+), 37 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index ad4dd92b2..b870cf64a 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -9,6 +9,7 @@ from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature import numpy +import math class MorleyDualSet(dual_set.DualSet): @@ -26,6 +27,9 @@ def __init__(self, ref_el, degree): facet = ref_el.construct_subelement(dim) Q_ref = create_quadrature(facet, degree+codim-2) scale = numpy.ones(Q_ref.get_weights().shape) + if codim == 1 and sd > 2: + # normalized normals do not have unit norm! + scale /= math.factorial(sd-1) for entity in sorted(top[dim]): cur = len(nodes) diff --git a/finat/morley.py b/finat/morley.py index 5cd00a856..e3735c102 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -11,6 +11,7 @@ class Morley(PhysicallyMappedElement, ScalarFiatElement): def __init__(self, cell, degree=2): if Citations is not None: Citations().register("Morley1971") + Citations().register("MingXu2006") super().__init__(FIAT.Morley(cell, degree=degree)) def basis_transformation(self, coordinate_mapping): @@ -20,19 +21,17 @@ 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) - pns = coordinate_mapping.physical_normals() - pts = coordinate_mapping.physical_tangents() - pel = coordinate_mapping.physical_edge_lengths() - V = identity(self.space_dimension()) offset = len(top[sd-2]) - + R = ListTensor([[Zero(), Literal(1)], [Literal(-1), Zero()]]) if sd == 2: + pel = coordinate_mapping.physical_edge_lengths() + pts = coordinate_mapping.physical_tangents() for e in top[sd-1]: s = offset + e - n = partial_indexed(pns, (e,)) t = partial_indexed(pts, (e,)) + n = R @ t nhat = self.cell.compute_normal(e) Jn = J @ Literal(nhat) Bnn = Jn @ n @@ -43,48 +42,34 @@ def basis_transformation(self, coordinate_mapping): V[s, v[0]] *= -1 else: - R = ListTensor([[Zero(), Literal(1)], [Literal(-1), Zero()]]) edges = self.cell.get_connectivity()[(sd-1, sd-2)] for face in top[sd-1]: s = offset + face - - te = numpy.array(list(map(self.cell.compute_edge_tangent, edges[face]))) - nhat = -numpy.cross(*te[:2]) + thats = self.cell.compute_tangents(sd-1, face) + nhat = numpy.cross(*thats) ahat = numpy.linalg.norm(nhat) - # Reciprocal basis nhat /= numpy.dot(nhat, nhat) - thats = numpy.array([numpy.cross(te[0], nhat), numpy.cross(te[1], nhat)]) - - Jt = J @ Literal(thats.T) Jn = J @ Literal(nhat) - Jte = J @ Literal(te.T) - - A = Jte.T @ Jte - area = Power(A[0, 0]*A[1, 1] - A[0, 1]*A[1, 0], Literal(0.5)) - - detF = area / ahat - Bnn = detJ / detF + Jt = J @ Literal(thats.T) + Gnt = Jn.T @ Jt + Gtt = Jt.T @ Jt + detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] + area = Power(detG, Literal(0.5)) + + Bnn = detJ * (ahat / area) + Bnt = Gnt @ R @ Gtt + # Not sure where this factor comes from + Bnt *= ahat / detG V[s, s] = Bnn - - Gte = Jt.T @ Jte - Gnt = Jn @ Jt - Bne = (Gnt @ R) @ Gte - - # Not sure where this comes from - factor = 2*ahat / (detF * detF) - for i, edge in enumerate(edges[face]): - # UFC convention alternates signs - sign = (-1) ** i - V[s, edge] = Bne[i] * sign * factor + V[s, list(edges[face])] = (Bnt[0] - Bnt[1], Bnt[1], -1*Bnt[0]) # diagonal post-scaling to patch up conditioning h = coordinate_mapping.cell_size() - for i in top[sd-1]: - s = i + offset - verts = top[sd-1][i] + for face in top[sd-1]: + s = offset + face + verts = top[sd-1][face] havg = sum(h[v] for v in verts) / len(verts) - for i in range(sd-1): - V[:, s] *= 1/havg + V[:, s] *= 1/havg return ListTensor(V.T) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index 30315a59c..3ce923844 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -94,6 +94,19 @@ pages = {20-24}, doi = {10.1243/03093247V061020} } +""") + Citations().add("MingXu2006", """ +@article{MingXu2006, + doi = {10.1007/s00211-005-0662-x}, + title={{The Morley element for fourth order elliptic equations in any dimensions}}, + author={Ming, Wang and Xu, Jinchao}, + journal={Numerische Mathematik}, + volume={103}, + number={1}, + pages={155--169}, + year={2006}, + publisher={Springer} +} """) Citations().add("Mardal2002", """ @article{Mardal2002, From 79a1daebfd3dc2b8dfff08cb4774c2e56edc1ccb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 24 Jul 2025 17:12:01 +0100 Subject: [PATCH 14/25] cleanup --- FIAT/morley.py | 4 ++-- finat/morley.py | 19 ++++++++++--------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index b870cf64a..71e2b1f8c 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -52,9 +52,9 @@ class Morley(finite_element.CiarletElement): def __init__(self, ref_el, degree=2): if ref_el.get_shape() not in {TRIANGLE, TETRAHEDRON}: - raise ValueError("Morley only defined on simplices") + raise ValueError("Morley only defined on simplices of dimension >= 2") if degree != 2: - raise ValueError("{type(self).__name__} only defined for degree=2") + raise ValueError("{type(self).__name__} only defined for degree == 2") poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = MorleyDualSet(ref_el, degree) super().__init__(poly_set, dual, degree) diff --git a/finat/morley.py b/finat/morley.py index e3735c102..29b142013 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -1,7 +1,7 @@ import FIAT import numpy -from gem import ListTensor, partial_indexed, Literal, Power, Zero +from gem import ListTensor, partial_indexed, Literal, Power from finat.fiat_elements import ScalarFiatElement from finat.physically_mapped import Citations, identity, PhysicallyMappedElement @@ -24,14 +24,14 @@ def basis_transformation(self, coordinate_mapping): V = identity(self.space_dimension()) offset = len(top[sd-2]) - R = ListTensor([[Zero(), Literal(1)], [Literal(-1), Zero()]]) if sd == 2: pel = coordinate_mapping.physical_edge_lengths() pts = coordinate_mapping.physical_tangents() + pns = coordinate_mapping.physical_normals() for e in top[sd-1]: s = offset + e t = partial_indexed(pts, (e,)) - n = R @ t + n = partial_indexed(pns, (e,)) nhat = self.cell.compute_normal(e) Jn = J @ Literal(nhat) Bnn = Jn @ n @@ -42,14 +42,15 @@ def basis_transformation(self, coordinate_mapping): V[s, v[0]] *= -1 else: + adjugate = lambda A: ListTensor([[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]]) edges = self.cell.get_connectivity()[(sd-1, sd-2)] for face in top[sd-1]: s = offset + face thats = self.cell.compute_tangents(sd-1, face) nhat = numpy.cross(*thats) ahat = numpy.linalg.norm(nhat) - nhat /= numpy.dot(nhat, nhat) + Jn = J @ Literal(nhat) Jt = J @ Literal(thats.T) Gnt = Jn.T @ Jt @@ -57,12 +58,12 @@ def basis_transformation(self, coordinate_mapping): detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] area = Power(detG, Literal(0.5)) - Bnn = detJ * (ahat / area) - Bnt = Gnt @ R @ Gtt - # Not sure where this factor comes from - Bnt *= ahat / detG + Bnn = detJ / area + Bnt = Gnt @ adjugate(Gtt) / detG + Bnn *= ahat + Bnt *= ahat V[s, s] = Bnn - V[s, list(edges[face])] = (Bnt[0] - Bnt[1], Bnt[1], -1*Bnt[0]) + V[s, list(edges[face])] = (-1*(Bnt[0] + Bnt[1]), Bnt[0], Bnt[1]) # diagonal post-scaling to patch up conditioning h = coordinate_mapping.cell_size() From 9c32f89eb0815c052f8f52915256c4df1c0def7d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 25 Jul 2025 23:44:18 +0100 Subject: [PATCH 15/25] Remove if statements --- FIAT/morley.py | 42 +++++++++++++++++++++++------------------- 1 file changed, 23 insertions(+), 19 deletions(-) diff --git a/FIAT/morley.py b/FIAT/morley.py index 71e2b1f8c..621836210 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -22,27 +22,31 @@ def __init__(self, ref_el, degree): sd = ref_el.get_spatial_dimension() entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} nodes = [] - for codim in (2, 1): - dim = sd - codim + + def duals(ref_el, dim, degree): facet = ref_el.construct_subelement(dim) - Q_ref = create_quadrature(facet, degree+codim-2) + Q_ref = create_quadrature(facet, degree) scale = numpy.ones(Q_ref.get_weights().shape) - if codim == 1 and sd > 2: - # normalized normals do not have unit norm! - scale /= math.factorial(sd-1) - - for entity in sorted(top[dim]): - cur = len(nodes) - if dim == sd-1: - # codim=1 dof -- average of normal derivative at each facet - ell = functional.IntegralMomentOfNormalDerivative(ref_el, entity, Q_ref, scale) - else: - # codim=2 dof -- integral average - Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) - ell = functional.IntegralMoment(ref_el, Q, scale / Q.jacobian_determinant()) - - nodes.append(ell) - entity_ids[dim][entity].extend(list(range(cur, len(nodes)))) + return Q_ref, scale + + # codim=2 dof -- integral average + dim = sd - 2 + Q_ref, scale = duals(ref_el, dim, degree) + for entity in sorted(top[dim]): + cur = len(nodes) + Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) + nodes.append(functional.IntegralMoment(ref_el, Q, scale / Q.jacobian_determinant())) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + # codim=1 dof -- average of normal derivative at each facet + dim = sd - 1 + Q_ref, scale = duals(ref_el, dim, degree-1) + # normalized normals do not have unit norm! + scale /= math.factorial(sd-1) + for entity in sorted(top[dim]): + cur = len(nodes) + nodes.append(functional.IntegralMomentOfNormalDerivative(ref_el, entity, Q_ref, scale)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) From ddb71b1c121e2e03cf9a7bec7694c6acb1affab9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 28 Jul 2025 07:46:24 +0100 Subject: [PATCH 16/25] error --- test/finat/test_zany_mapping.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 956480c8b..3119d57d1 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -82,15 +82,17 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): Mh[abs(Mh) < 1E-10] = 0 M[abs(M) < 1E-10] = 0 - with np.errstate(divide='ignore', invalid='ignore'): - error = M.T / Mh.T - 1 - error[error != error] = 0 + error = M.T - Mh.T error[abs(error) < 1E-10] = 0 - error = error[np.ix_(*map(np.unique, np.nonzero(error)))] - error[error != 0] += 1 + with np.errstate(divide='ignore', invalid='ignore'): + # error /= Mh.T + error[error != error] = 0 + # error[error!=0] += 1 + # error = error[:num_dofs] + error = error[np.ix_(*map(np.unique, np.nonzero(error)))] assert np.allclose(residual, 0), str(error) - assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(error) + assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(error.tolist()) @pytest.mark.parametrize("element", [ From a812d3cf0e6c5d50f8aac60834e3c2440d3d6d28 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 14 Dec 2025 16:26:21 +0000 Subject: [PATCH 17/25] cite --- finat/walkington.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/finat/walkington.py b/finat/walkington.py index afed04acb..e715d6c84 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -4,16 +4,16 @@ from gem import ListTensor +from finat.citations import cite from finat.fiat_elements import ScalarFiatElement -from finat.physically_mapped import Citations, identity, PhysicallyMappedElement +from finat.physically_mapped import identity, PhysicallyMappedElement from finat.argyris import _vertex_transform, _normal_tangential_transform from copy import deepcopy class Walkington(PhysicallyMappedElement, ScalarFiatElement): def __init__(self, cell, degree=5): - if Citations is not None: - Citations().register("Bell1969") + cite("Bell1969") super().__init__(FIAT.Walkington(cell, degree=degree)) reduced_dofs = deepcopy(self._element.entity_dofs()) From fb8e712c1522cf77e6a14089f30f8e19fbf4d17c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 29 Dec 2025 18:45:02 -0600 Subject: [PATCH 18/25] WIP --- FIAT/walkington.py | 13 ++--- finat/walkington.py | 85 +++++++++++++++++++++------ test/FIAT/unit/test_walkington.py | 96 +++++++++++++++++++++++++++++++ test/finat/test_zany_mapping.py | 23 +++++--- 4 files changed, 186 insertions(+), 31 deletions(-) create mode 100644 test/FIAT/unit/test_walkington.py diff --git a/FIAT/walkington.py b/FIAT/walkington.py index 80fd2c128..c7809fc3b 100644 --- a/FIAT/walkington.py +++ b/FIAT/walkington.py @@ -39,13 +39,12 @@ def __init__(self, ref_el, degree): # Face dofs: moments against cubic bubble ref_face = ref_el.construct_subelement(2) Q_face, phis = make_dual_bubbles(ref_face, degree-1) + phis *= ref_face.volume()**0.5 + phis[0] = 1 for face in sorted(top[2]): cur = len(nodes) - nface = ref_el.compute_scaled_normal(face) - nface /= numpy.dot(nface, nface) - - nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phis[0], n=nface)) + nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phis[0])) entity_ids[2][face].extend(range(cur, len(nodes))) # Interior dof: point evaluation at barycenter @@ -66,8 +65,8 @@ def __init__(self, ref_el, degree): for face in sorted(top[2]): cur = len(nodes) - nface = ref_el.compute_scaled_normal(face) - nface /= numpy.dot(nface, nface) + nface = numpy.cross(*ref_el.compute_tangents(sd-1, face)) + nface /= -numpy.dot(nface, nface) for i, e in enumerate(edges[face]): Q = FacetQuadratureRule(ref_face, 1, i, Q_edge) @@ -79,7 +78,7 @@ def __init__(self, ref_el, degree): nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q, leg4_at_qpts, n=nfe)) - nodes.extend(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phi, n=nface) for phi in phis[1:]) + nodes.extend(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phi) for phi in phis[1:]) entity_ids[2][face].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/finat/walkington.py b/finat/walkington.py index e715d6c84..1012e634b 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -1,14 +1,14 @@ import FIAT from FIAT.polynomial_set import mis -from math import comb -from gem import ListTensor +from gem import ListTensor, Literal, Power from finat.citations import cite from finat.fiat_elements import ScalarFiatElement from finat.physically_mapped import identity, PhysicallyMappedElement from finat.argyris import _vertex_transform, _normal_tangential_transform from copy import deepcopy +import numpy class Walkington(PhysicallyMappedElement, ScalarFiatElement): @@ -34,34 +34,50 @@ def basis_transformation(self, coordinate_mapping): ndof = self.space_dimension() # rectangular to toss out the constraint dofs V = identity(numbf, ndof) - vorder = 2 _vertex_transform(V, vorder, self.cell, coordinate_mapping) - offset = ndof - voffset = comb(sd + vorder, vorder) - foffset = len(self._element.entity_dofs()[2][0]) - len(self.entity_dofs()[2][0]) - + entity_dofs = self._element.entity_dofs() edges = self.cell.get_connectivity()[(2, 1)] - for f in sorted(top[2]): - q = len(top[0]) * voffset + f - V[q, q] *= detJ + for f in sorted(entity_dofs[2]): + fdofs = entity_dofs[2][f] + q = fdofs[0] + c = fdofs[-2:] + + Rnn, Rnt = morley_transform(self.cell, J, detJ, f) + V[q, q] = Rnn for j, e in enumerate(edges[f]): - s = offset + foffset * f + j + s = fdofs[1+j] - v0id, v1id = (v * voffset for v in top[1][e]) + v0id, v1id = (entity_dofs[0][v][0] for v in top[1][e]) Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e, face=f) # vertex points V[s, v1id] = 1/21 * Bnt V[s, v0id] = -1 * V[s, v1id] + V[q, v1id] += Rnt[j] + V[q, v0id] += Rnt[j] + + # TODO + V[c, v1id] += Rnt[j] + V[c, v0id] += Rnt[j] + # vertex derivatives for i in range(sd): V[s, v1id+1+i] = -1/42 * Bnt * Jt[i] V[s, v0id+1+i] = V[s, v1id+1+i] + R1 = 1/5 * Rnt[j] * Jt[i] + V[q, v1id+1+i] -= R1 + V[q, v0id+1+i] += R1 + + # TODO + R1 = 1/2 * Rnt[j] * Jt[i] + V[c, v1id+1+i] -= R1 + V[c, v0id+1+i] += R1 + # second derivatives for i, alpha in enumerate(mis(sd, 2)): ids = tuple(k for k, ak in enumerate(alpha) if ak) @@ -70,13 +86,24 @@ def basis_transformation(self, coordinate_mapping): V[s, v1id+sd+1+i] = 1/252 * Bnt * tau V[s, v0id+sd+1+i] = -1 * V[s, v1id+sd+1+i] + R2 = 1/60 * Rnt[j] * tau + V[q, v1id+sd+1+i] += R2 + V[q, v0id+sd+1+i] += R2 + + # TODO + R2 = 1/12 * Rnt[j] * tau + V[c, v1id+sd+1+i] += R2 + V[c, v0id+sd+1+i] += R2 + + V[c[0], :] *= (2/5)*(6/7)**0.5 + V[c[1], :] *= -(2/5)*(2/7)**0.5 + # Patch up conditioning h = coordinate_mapping.cell_size() - for v in sorted(top[0]): - s = voffset * v + 1 - V[:, s:s+sd] *= 1/h[v] - V[:, s+sd:voffset*(v+1)] *= 1/(h[v]*h[v]) - + for v in sorted(entity_dofs[0]): + vdofs = entity_dofs[0][v] + V[:, vdofs[1:1+sd]] *= 1/h[v] + V[:, vdofs[1+sd:]] *= 1/(h[v]*h[v]) return ListTensor(V.T) # This wipes out the edge dofs. FIAT gives a 65 DOF element @@ -88,3 +115,27 @@ def entity_dofs(self): def space_dimension(self): return 45 + + +def morley_transform(cell, J, detJ, face): + adjugate = lambda A: ListTensor([[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]]) + sd = cell.get_spatial_dimension() + thats = cell.compute_tangents(sd-1, face) + nhat = numpy.cross(*thats) + ahat = numpy.linalg.norm(nhat) + nhat /= numpy.dot(nhat, nhat) + + nhat = Literal(nhat) + Jn = J @ nhat + Jt = J @ Literal(thats.T) + Gnt = Jn.T @ Jt + Gtt = Jt.T @ Jt + detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] + area = Power(detG, Literal(0.5)) + + Bnn = detJ / area + Bnt = Gnt @ adjugate(Gtt) / detG + Bnn *= ahat + Bnt *= ahat + Bnt = (-1*(Bnt[0] + Bnt[1]), Bnt[0], Bnt[1]) + return Bnn, Bnt diff --git a/test/FIAT/unit/test_walkington.py b/test/FIAT/unit/test_walkington.py new file mode 100644 index 000000000..8126d55a9 --- /dev/null +++ b/test/FIAT/unit/test_walkington.py @@ -0,0 +1,96 @@ +import pytest +import numpy + +from FIAT import Walkington +from FIAT.expansions import polynomial_dimension +from FIAT.polynomial_set import ONPolynomialSet +from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature +from FIAT.reference_element import ufc_simplex + + +@pytest.fixture(params=[0, 1], ids=["K1", "K2"]) +def cell(request): + K = ufc_simplex(3) + if request.param == 1: + K.vertices = ((0, 0, 0), + (1., 0.1, -0.37), + (0.01, 0.987, -.23), + (-0.1, -0.2, 1.38)) + return K + + +def directional_derivative(direction, tab): + return sum(direction[alpha.index(1)] * tab[alpha] + for alpha in tab if sum(alpha) == 1) + + +def inner(u, v, wts): + return numpy.dot(numpy.multiply(u, wts), v.T) + + +def test_walkington_basis_functions(cell): + degree = 5 + fe = Walkington(cell, degree) + space_dim = 45 + + ref_el = fe.get_reference_element() + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + entity_ids = fe.entity_dofs() + + degree = fe.degree() + + ref_facet = ufc_simplex(sd-1) + Qref = create_quadrature(ref_facet, 2*degree-1) + + P = ONPolynomialSet(ref_facet, degree-1) + Ptab = P.tabulate(Qref.get_points()) + P_at_qpts = Ptab[(0,)*(sd-1)] + P3dim = polynomial_dimension(ref_facet, degree-2) + P4_at_qpts = P_at_qpts[P3dim:] + + for f in top[sd-1]: + n = ref_el.compute_normal(f) + Q = FacetQuadratureRule(ref_el, sd-1, f, Qref) + qpts, qwts = Q.get_points(), Q.get_weights() + + ids = entity_ids[sd-1][f] + + tab = fe.tabulate(1, qpts) + phi = tab[(0,) * sd] + phi_n = directional_derivative(n, tab) + + # Test that normal derivative moment bfs have vanishing trace + assert numpy.allclose(phi[ids[:1]], 0) + + # Test that trace moment bfs have vanishing normal derivative + P4_moments = inner(phi_n[:space_dim], P4_at_qpts, qwts) + assert numpy.allclose(P4_moments, 0) + + +def span_greater_equal(A, B): + # span(A) >= span(B) + _, residual, *_ = numpy.linalg.lstsq(A.reshape(A.shape[0], -1).T, + B.reshape(B.shape[0], -1).T) + return numpy.allclose(residual, 0) + + +def test_walkington_space(cell): + degree = 5 + P = ONPolynomialSet(cell, degree) + + fe = Walkington(cell, degree) + V = fe.get_nodal_basis() + ref_complex = V.ref_el + + pts = [] + top = ref_complex.topology + for dim in top: + for entity in top[dim]: + pts.extend(ref_complex.make_points(dim, entity, degree)) + V_tab = V.tabulate(pts) + P_tab = P.tabulate(pts) + + sd = cell.get_spatial_dimension() + assert span_greater_equal(V_tab[(0,)*sd], P_tab[(0,)*sd]) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index a9be25164..357f07da2 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -2,6 +2,8 @@ import finat import numpy as np import pytest +import pprint + from gem.interpreter import evaluate from finat.physically_mapped import PhysicallyMappedElement @@ -80,18 +82,25 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): Vh, residual, *_ = np.linalg.lstsq(Phi.T, phi.T) Mh = Vh.T Mh = Mh[:num_dofs] - Mh[abs(Mh) < 1E-10] = 0 - M[abs(M) < 1E-10] = 0 + tol = 1E-10 + Mh[abs(Mh) < tol] = 0 + M[abs(M) < tol] = 0 + delta = M.T - Mh.T + delta[abs(delta) < tol] = 0 with np.errstate(divide='ignore', invalid='ignore'): - error = M.T / Mh.T - 1 + error = delta / Mh.T error[error != error] = 0 - error[abs(error) < 1E-10] = 0 - error = error[np.ix_(*map(np.unique, np.nonzero(error)))] + error[delta == 0] = 0 + error[abs(error) < tol] = 0 + + inds = tuple(map(np.unique, np.nonzero(error))) + error = error[np.ix_(*inds)] error[error != 0] += 1 - assert np.allclose(residual, 0), str(error) - assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), str(error) + pp = pprint.PrettyPrinter(width=140, compact=True) + assert np.allclose(residual, 0), pp.pformat((np.round(error, 8).tolist(), *inds)) + assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), pp.pformat((np.round(error, 8).tolist(), *inds)) @pytest.mark.parametrize("element", [ From cf7c9104aaa07eeca07dee7e82eab9cc5ee37729 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 30 Dec 2025 15:41:54 -0600 Subject: [PATCH 19/25] WIP --- FIAT/functional.py | 34 +++++++++++++++++++ FIAT/hierarchical.py | 4 +-- FIAT/walkington.py | 31 ++++++++++------- finat/walkington.py | 55 +++++++++++++++---------------- test/FIAT/unit/test_walkington.py | 12 +++++-- test/finat/conftest.py | 5 +++ 6 files changed, 97 insertions(+), 44 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index c9d634a96..fb2445113 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -365,6 +365,40 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts, n=None): {}, dpt_dict, "IntegralMomentOfNormalDerivative") +class IntegralMomentOfBidirectionalDerivative(Functional): + """Functional giving normal derivative integrated against some function on a facet.""" + + def __init__(self, ref_el, Q, f_at_qpts, s1, s2): + self.f_at_qpts = f_at_qpts + self.Q = Q + + sd = ref_el.get_spatial_dimension() + + # map points onto facet + points = Q.get_points() + self.dpts = points + weights = numpy.multiply(f_at_qpts, Q.get_weights()) + + tau = numpy.zeros((sd*(sd+1)//2,)) + alphas = [] + cur = 0 + for i in range(sd): + for j in range(i, sd): + alpha = [0] * sd + alpha[i] += 1 + alpha[j] += 1 + alphas.append(tuple(alpha)) + tau[cur] = s1[i] * s2[j] + (i != j) * s2[i] * s1[j] + cur += 1 + + dpt_dict = {tuple(pt): [(wt*tau[i], alphas[i], tuple()) for i in range(len(alphas))] + for pt, wt in zip(points, weights)} + + + super().__init__(ref_el, tuple(), + {}, dpt_dict, "IntegralMomentOfBidirectionalDerivative") + + class FrobeniusIntegralMoment(IntegralMoment): def __init__(self, ref_el, Q, f_at_qpts, nm=None): diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index b77553623..81020b13d 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -16,7 +16,7 @@ from FIAT.P0 import P0 -def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme=None): +def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme=None, scale="orthonormal"): """Tabulate the L2-duals of the hierarchical C0 basis.""" dim = ref_el.get_spatial_dimension() if dim == 0: @@ -25,7 +25,7 @@ def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme if interpolant_deg is None: interpolant_deg = degree Q = parse_quadrature_scheme(ref_el, degree + interpolant_deg, quad_scheme) - B = make_bubbles(ref_el, degree, codim=codim, scale="orthonormal") + B = make_bubbles(ref_el, degree, codim=codim, scale=scale) P_at_qpts = B.expansion_set.tabulate(degree, Q.get_points()) M = numpy.dot(numpy.multiply(P_at_qpts, Q.get_weights()), P_at_qpts.T) phis = numpy.linalg.solve(M, P_at_qpts) diff --git a/FIAT/walkington.py b/FIAT/walkington.py index c7809fc3b..61eaba933 100644 --- a/FIAT/walkington.py +++ b/FIAT/walkington.py @@ -9,7 +9,11 @@ # bfs, but the extra 20 are used in the transformation theory. from FIAT import finite_element, polynomial_set, dual_set, macro -from FIAT.functional import PointEvaluation, PointDerivative, IntegralMomentOfNormalDerivative +from FIAT.functional import ( + PointEvaluation, PointDerivative, + IntegralMomentOfNormalDerivative, + IntegralMomentOfBidirectionalDerivative, +) from FIAT.reference_element import TETRAHEDRON from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature @@ -18,6 +22,10 @@ import numpy +def inner(u, v, wts): + return numpy.dot(numpy.multiply(u, wts), v.T) + + class WalkingtonDualSet(dual_set.DualSet): def __init__(self, ref_el, degree): top = ref_el.get_topology() @@ -36,15 +44,13 @@ def __init__(self, ref_el, degree): for i in (1, 2) for alpha in polynomial_set.mis(sd, i)) entity_ids[0][v].extend(range(cur, len(nodes))) - # Face dofs: moments against cubic bubble + # Face dofs: moments or normal derivative ref_face = ref_el.construct_subelement(2) - Q_face, phis = make_dual_bubbles(ref_face, degree-1) - phis *= ref_face.volume()**0.5 - phis[0] = 1 - + Q_face = create_quadrature(ref_face, degree-1) + f_at_qpts = numpy.ones(Q_face.get_weights().shape) for face in sorted(top[2]): cur = len(nodes) - nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phis[0])) + nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, f_at_qpts)) entity_ids[2][face].extend(range(cur, len(nodes))) # Interior dof: point evaluation at barycenter @@ -56,17 +62,19 @@ def __init__(self, ref_el, degree): # Constraint dofs # Face-edge constraint: normal derivative along edge is cubic - # Face constraint: normal derivative is cubic edges = ref_el.get_connectivity()[(2, 1)] ref_edge = ref_el.construct_subelement(1) Q_edge = create_quadrature(ref_edge, 2*(degree-1)) x = ref_edge.compute_barycentric_coordinates(Q_edge.get_points()) leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0]) + # Face constraint: normal derivative is cubic + Q_face, phis = make_dual_bubbles(ref_face, degree-2, scale=1) for face in sorted(top[2]): cur = len(nodes) - nface = numpy.cross(*ref_el.compute_tangents(sd-1, face)) - nface /= -numpy.dot(nface, nface) + thats = ref_el.compute_tangents(sd-1, face) + nface = numpy.cross(*thats) + nface *= -1/numpy.dot(nface, nface) for i, e in enumerate(edges[face]): Q = FacetQuadratureRule(ref_face, 1, i, Q_edge) @@ -78,7 +86,8 @@ def __init__(self, ref_el, degree): nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q, leg4_at_qpts, n=nfe)) - nodes.extend(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phi) for phi in phis[1:]) + Q = FacetQuadratureRule(ref_el, 2, face, Q_face) + nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, phis[0], nface, t) for t in thats) entity_ids[2][face].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/finat/walkington.py b/finat/walkington.py index 1012e634b..18349c4c6 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -39,64 +39,63 @@ def basis_transformation(self, coordinate_mapping): entity_dofs = self._element.entity_dofs() edges = self.cell.get_connectivity()[(2, 1)] - for f in sorted(entity_dofs[2]): + for f in entity_dofs[2]: fdofs = entity_dofs[2][f] q = fdofs[0] c = fdofs[-2:] - Rnn, Rnt = morley_transform(self.cell, J, detJ, f) + Rnn, Rnt, Jn = morley_transform(self.cell, J, detJ, f) V[q, q] = Rnn for j, e in enumerate(edges[f]): + vid0, vid1 = (entity_dofs[0][v][0] for v in top[1][e]) s = fdofs[1+j] - - v0id, v1id = (entity_dofs[0][v][0] for v in top[1][e]) Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e, face=f) # vertex points - V[s, v1id] = 1/21 * Bnt - V[s, v0id] = -1 * V[s, v1id] + V[s, vid1] = 1/21 * Bnt + V[s, vid0] = -1 * V[s, vid1] - V[q, v1id] += Rnt[j] - V[q, v0id] += Rnt[j] + V[q, vid1] += Rnt[j] + V[q, vid0] += Rnt[j] # TODO - V[c, v1id] += Rnt[j] - V[c, v0id] += Rnt[j] + c = fdofs[-2:] + Qnt = (Jn @ Jt) / (Jt @ Jt) + V[c, vid1] += Qnt + V[c, vid0] += Qnt + c = [] # vertex derivatives for i in range(sd): - V[s, v1id+1+i] = -1/42 * Bnt * Jt[i] - V[s, v0id+1+i] = V[s, v1id+1+i] + V[s, vid1+i+1] = -1/42 * Bnt * Jt[i] + V[s, vid0+i+1] = V[s, vid1+1+i] R1 = 1/5 * Rnt[j] * Jt[i] - V[q, v1id+1+i] -= R1 - V[q, v0id+1+i] += R1 + V[q, vid1+i+1] -= R1 + V[q, vid0+i+1] += R1 # TODO - R1 = 1/2 * Rnt[j] * Jt[i] - V[c, v1id+1+i] -= R1 - V[c, v0id+1+i] += R1 + V[c, vid1+i+1] += Jt[i] * Qnt + V[c, vid0+i+1] -= Jt[i] * Qnt + c = [] # second derivatives - for i, alpha in enumerate(mis(sd, 2)): + for i, alpha in enumerate(mis(sd, 2), start=sd+1): ids = tuple(k for k, ak in enumerate(alpha) if ak) a, b = ids[0], ids[-1] tau = (1 + (a != b)) * Jt[a] * Jt[b] - V[s, v1id+sd+1+i] = 1/252 * Bnt * tau - V[s, v0id+sd+1+i] = -1 * V[s, v1id+sd+1+i] + V[s, vid1+i] = 1/252 * Bnt * tau + V[s, vid0+i] = -1 * V[s, vid1+i] R2 = 1/60 * Rnt[j] * tau - V[q, v1id+sd+1+i] += R2 - V[q, v0id+sd+1+i] += R2 + V[q, vid1+i] += R2 + V[q, vid0+i] += R2 # TODO - R2 = 1/12 * Rnt[j] * tau - V[c, v1id+sd+1+i] += R2 - V[c, v0id+sd+1+i] += R2 + V[c, vid1+i] += Qnt * tau + V[c, vid0+i] += Qnt * tau - V[c[0], :] *= (2/5)*(6/7)**0.5 - V[c[1], :] *= -(2/5)*(2/7)**0.5 # Patch up conditioning h = coordinate_mapping.cell_size() @@ -138,4 +137,4 @@ def morley_transform(cell, J, detJ, face): Bnn *= ahat Bnt *= ahat Bnt = (-1*(Bnt[0] + Bnt[1]), Bnt[0], Bnt[1]) - return Bnn, Bnt + return Bnn, Bnt, Jn diff --git a/test/FIAT/unit/test_walkington.py b/test/FIAT/unit/test_walkington.py index 8126d55a9..7bec435c6 100644 --- a/test/FIAT/unit/test_walkington.py +++ b/test/FIAT/unit/test_walkington.py @@ -78,19 +78,25 @@ def span_greater_equal(A, B): def test_walkington_space(cell): degree = 5 - P = ONPolynomialSet(cell, degree) - fe = Walkington(cell, degree) + space_dim = 45 + V = fe.get_nodal_basis() ref_complex = V.ref_el - pts = [] top = ref_complex.topology for dim in top: for entity in top[dim]: pts.extend(ref_complex.make_points(dim, entity, degree)) V_tab = V.tabulate(pts) + + P = ONPolynomialSet(cell, degree) P_tab = P.tabulate(pts) + # Test that the augmented space includes all quintics sd = cell.get_spatial_dimension() assert span_greater_equal(V_tab[(0,)*sd], P_tab[(0,)*sd]) + + # Test that the reduced space includes all quartics + dimP4 = polynomial_dimension(cell, degree-1) + assert span_greater_equal(V_tab[(0,)*sd][:space_dim], P_tab[(0,)*sd][:dimP4]) diff --git a/test/finat/conftest.py b/test/finat/conftest.py index c31ef7350..51f4b3106 100644 --- a/test/finat/conftest.py +++ b/test/finat/conftest.py @@ -114,6 +114,11 @@ def phys_el(): (1., 0.1, -0.37), (0.01, 0.987, -.23), (-0.1, -0.2, 1.38)) + + K[3].vertices = ((0.0, 0.0, 0.0), + (1.0, 0.0, 0.0), + (0.0, 1.0, 0.0), + (0.0, 0.0, 2.0)) return K From 1fc8f398dddc4720e074476081ff3397e29880bb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 3 Jan 2026 13:52:42 -0600 Subject: [PATCH 20/25] Zany magic --- FIAT/functional.py | 3 +- FIAT/hierarchical.py | 4 +- FIAT/walkington.py | 29 ++++++++--- finat/citations.py | 10 ++++ finat/morley.py | 48 ++++++++++------- finat/walkington.py | 115 ++++++++++++++++++----------------------- test/finat/conftest.py | 4 -- 7 files changed, 113 insertions(+), 100 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index fb2445113..adca434c7 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -366,7 +366,7 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts, n=None): class IntegralMomentOfBidirectionalDerivative(Functional): - """Functional giving normal derivative integrated against some function on a facet.""" + """Functional giving second derivative integrated against some function on a facet.""" def __init__(self, ref_el, Q, f_at_qpts, s1, s2): self.f_at_qpts = f_at_qpts @@ -394,7 +394,6 @@ def __init__(self, ref_el, Q, f_at_qpts, s1, s2): dpt_dict = {tuple(pt): [(wt*tau[i], alphas[i], tuple()) for i in range(len(alphas))] for pt, wt in zip(points, weights)} - super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfBidirectionalDerivative") diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index 81020b13d..b77553623 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -16,7 +16,7 @@ from FIAT.P0 import P0 -def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme=None, scale="orthonormal"): +def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme=None): """Tabulate the L2-duals of the hierarchical C0 basis.""" dim = ref_el.get_spatial_dimension() if dim == 0: @@ -25,7 +25,7 @@ def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None, quad_scheme if interpolant_deg is None: interpolant_deg = degree Q = parse_quadrature_scheme(ref_el, degree + interpolant_deg, quad_scheme) - B = make_bubbles(ref_el, degree, codim=codim, scale=scale) + B = make_bubbles(ref_el, degree, codim=codim, scale="orthonormal") P_at_qpts = B.expansion_set.tabulate(degree, Q.get_points()) M = numpy.dot(numpy.multiply(P_at_qpts, Q.get_weights()), P_at_qpts.T) phis = numpy.linalg.solve(M, P_at_qpts) diff --git a/FIAT/walkington.py b/FIAT/walkington.py index 61eaba933..3a4b24231 100644 --- a/FIAT/walkington.py +++ b/FIAT/walkington.py @@ -18,7 +18,6 @@ from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature from FIAT.jacobi import eval_jacobi -from FIAT.hierarchical import make_dual_bubbles import numpy @@ -68,13 +67,13 @@ def __init__(self, ref_el, degree): x = ref_edge.compute_barycentric_coordinates(Q_edge.get_points()) leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0]) # Face constraint: normal derivative is cubic - Q_face, phis = make_dual_bubbles(ref_face, degree-2, scale=1) + Q_face, phi = face_constraints(ref_face, degree) for face in sorted(top[2]): cur = len(nodes) thats = ref_el.compute_tangents(sd-1, face) - nface = numpy.cross(*thats) - nface *= -1/numpy.dot(nface, nface) + nface = -numpy.cross(*thats) + nface /= numpy.linalg.norm(nface) for i, e in enumerate(edges[face]): Q = FacetQuadratureRule(ref_face, 1, i, Q_edge) @@ -83,13 +82,15 @@ def __init__(self, ref_el, degree): nfe = numpy.cross(te, nface) nfe /= numpy.linalg.norm(nfe) nfe /= Q.jacobian_determinant() - nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q, leg4_at_qpts, n=nfe)) Q = FacetQuadratureRule(ref_el, 2, face, Q_face) - nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, phis[0], nface, t) for t in thats) + nface /= Q.jacobian_determinant() + nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, phi, nface, t) for t in thats) entity_ids[2][face].extend(range(cur, len(nodes))) + self.Q_face = Q_face + self.phi = phi super().__init__(nodes, ref_el, entity_ids) @@ -106,3 +107,19 @@ def __init__(self, ref_el, degree=5): ref_complex = macro.AlfeldSplit(ref_el) poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=4, variant="bubble") super().__init__(poly_set, dual, degree) + + +def face_constraints(ref_face, degree): + from FIAT.expansions import polynomial_dimension + sd = ref_face.get_spatial_dimension() + Q_face = create_quadrature(ref_face, 2*(degree-2)) + dimPkm1 = polynomial_dimension(ref_face, degree-3) + + pts = list(Q_face.get_points()[:3]) + pts.append(Q_face.get_points()[-1]) + P = polynomial_set.ONPolynomialSet(ref_face, degree-2, scale=1) + Pk = P.tabulate(pts)[(0,)*sd][dimPkm1:] + c = numpy.linalg.solve(Pk.T, [0, 0, 0, 1/Q_face.get_weights()[-1]]) + Pk = P.tabulate(Q_face.get_points())[(0,)*sd][dimPkm1:] + phi = numpy.dot(c, Pk) + return Q_face, phi diff --git a/finat/citations.py b/finat/citations.py index c80227e66..4f0baea1a 100644 --- a/finat/citations.py +++ b/finat/citations.py @@ -259,3 +259,13 @@ def cite(*args, **kwargs): doi={https://doi.org/10.1023/A:1004420829610}, } """) + petsctools.add_citation("Walkington2010", """ +author = {Walkington, Noel J.}, +title = {{A \\$C^1\\$ Tetrahedral Finite Element without Edge Degrees of Freedom}}, +journal = {SIAM Journal on Numerical Analysis}, +volume = {52}, +number = {1}, +pages = {330-342}, +year = {2014}, +doi = {10.1137/130912013}, +}""") diff --git a/finat/morley.py b/finat/morley.py index 64c442943..ed58ed63e 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -8,6 +8,30 @@ from finat.physically_mapped import identity, PhysicallyMappedElement +def morley_transform(cell, J, detJ, face): + adjugate = lambda A: ListTensor([[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]]) + sd = cell.get_spatial_dimension() + thats = cell.compute_tangents(sd-1, face) + nhat = numpy.cross(*thats) + ahat = numpy.linalg.norm(nhat) + nhat /= numpy.dot(nhat, nhat) + + nhat = Literal(nhat) + Jn = J @ nhat + Jt = J @ Literal(thats.T) + Gnt = Jn.T @ Jt + Gtt = Jt.T @ Jt + detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] + area = Power(detG, Literal(0.5)) + + Bnn = detJ / area + Bnt = Gnt @ adjugate(Gtt) / detG + Bnn *= ahat + Bnt *= ahat + Bnt = (-1*(Bnt[0] + Bnt[1]), Bnt[0], Bnt[1]) + return Bnn, Bnt + + class Morley(PhysicallyMappedElement, ScalarFiatElement): def __init__(self, cell, degree=2): cite("Morley1971") @@ -42,28 +66,12 @@ def basis_transformation(self, coordinate_mapping): V[s, v[0]] *= -1 else: - adjugate = lambda A: ListTensor([[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]]) edges = self.cell.get_connectivity()[(sd-1, sd-2)] for face in top[sd-1]: - s = offset + face - thats = self.cell.compute_tangents(sd-1, face) - nhat = numpy.cross(*thats) - ahat = numpy.linalg.norm(nhat) - nhat /= numpy.dot(nhat, nhat) - - Jn = J @ Literal(nhat) - Jt = J @ Literal(thats.T) - Gnt = Jn.T @ Jt - Gtt = Jt.T @ Jt - detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] - area = Power(detG, Literal(0.5)) - - Bnn = detJ / area - Bnt = Gnt @ adjugate(Gtt) / detG - Bnn *= ahat - Bnt *= ahat - V[s, s] = Bnn - V[s, list(edges[face])] = (-1*(Bnt[0] + Bnt[1]), Bnt[0], Bnt[1]) + Bnn, Bnt = morley_transform(self.cell, J, detJ, face) + fid = offset + face + V[fid, fid] = Bnn + V[fid, list(edges[face])] = Bnt # diagonal post-scaling to patch up conditioning h = coordinate_mapping.cell_size() diff --git a/finat/walkington.py b/finat/walkington.py index 18349c4c6..01d70cca4 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -1,19 +1,21 @@ import FIAT -from FIAT.polynomial_set import mis +import numpy -from gem import ListTensor, Literal, Power +from FIAT.polynomial_set import mis +from gem import ListTensor, Zero from finat.citations import cite from finat.fiat_elements import ScalarFiatElement from finat.physically_mapped import identity, PhysicallyMappedElement from finat.argyris import _vertex_transform, _normal_tangential_transform +from finat.morley import morley_transform from copy import deepcopy -import numpy +from itertools import chain class Walkington(PhysicallyMappedElement, ScalarFiatElement): def __init__(self, cell, degree=5): - cite("Bell1969") + cite("Walkington2010") super().__init__(FIAT.Walkington(cell, degree=degree)) reduced_dofs = deepcopy(self._element.entity_dofs()) @@ -39,63 +41,68 @@ def basis_transformation(self, coordinate_mapping): entity_dofs = self._element.entity_dofs() edges = self.cell.get_connectivity()[(2, 1)] + + Q_face = self._element.dual.Q_face + phi = self._element.dual.phi + wts = numpy.multiply(Q_face.get_weights(), phi) + pts = Q_face.get_points() + pts = pts[abs(wts) > 1E-10] + wts = wts[abs(wts) > 1E-10] + for f in entity_dofs[2]: + Rnn, Rnt = morley_transform(self.cell, J, detJ, f) fdofs = entity_dofs[2][f] - q = fdofs[0] - c = fdofs[-2:] - - Rnn, Rnt, Jn = morley_transform(self.cell, J, detJ, f) - V[q, q] = Rnn + fid = fdofs[0] + V[fid, fid] = Rnn for j, e in enumerate(edges[f]): - vid0, vid1 = (entity_dofs[0][v][0] for v in top[1][e]) - s = fdofs[1+j] Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e, face=f) + vid0, vid1 = (entity_dofs[0][v][0] for v in top[1][e]) + eid = fdofs[1+j] # vertex points - V[s, vid1] = 1/21 * Bnt - V[s, vid0] = -1 * V[s, vid1] - - V[q, vid1] += Rnt[j] - V[q, vid0] += Rnt[j] - - # TODO - c = fdofs[-2:] - Qnt = (Jn @ Jt) / (Jt @ Jt) - V[c, vid1] += Qnt - V[c, vid0] += Qnt - c = [] + V[fid, vid1] += Rnt[j] + V[fid, vid0] += Rnt[j] + V[eid, vid1] = 1/21 * Bnt + V[eid, vid0] = -1 * V[eid, vid1] # vertex derivatives for i in range(sd): - V[s, vid1+i+1] = -1/42 * Bnt * Jt[i] - V[s, vid0+i+1] = V[s, vid1+1+i] - R1 = 1/5 * Rnt[j] * Jt[i] - V[q, vid1+i+1] -= R1 - V[q, vid0+i+1] += R1 - - # TODO - V[c, vid1+i+1] += Jt[i] * Qnt - V[c, vid0+i+1] -= Jt[i] * Qnt - c = [] + V[fid, vid1+i+1] -= R1 + V[fid, vid0+i+1] += R1 + V[eid, vid1+i+1] = -1/42 * Bnt * Jt[i] + V[eid, vid0+i+1] = V[eid, vid1+1+i] # second derivatives for i, alpha in enumerate(mis(sd, 2), start=sd+1): ids = tuple(k for k, ak in enumerate(alpha) if ak) a, b = ids[0], ids[-1] tau = (1 + (a != b)) * Jt[a] * Jt[b] - V[s, vid1+i] = 1/252 * Bnt * tau - V[s, vid0+i] = -1 * V[s, vid1+i] R2 = 1/60 * Rnt[j] * tau - V[q, vid1+i] += R2 - V[q, vid0+i] += R2 - - # TODO - V[c, vid1+i] += Qnt * tau - V[c, vid0+i] += Qnt * tau - + V[fid, vid1+i] += R2 + V[fid, vid0+i] += R2 + V[eid, vid1+i] = 1/252 * Bnt * tau + V[eid, vid0+i] = -1 * V[eid, vid1+i] + + # Evaluate nodal completion of the face constraints + tab = self._element.tabulate(2, pts, entity=(2, f)) + thats = self.cell.compute_tangents(2, f) + tt = [numpy.outer(thats[i], thats[j]) for i in range(2) for j in range(i, 2)] + T = [[a[i, j] + (i != j) * a[j, i] for i in range(sd) for j in range(i, sd)] for a in tt] + C = numpy.array([tab[alpha] @ wts for alpha in tab if sum(alpha) == 2]) + C = numpy.dot(T, C) + C = C.astype(object) + C[abs(C) < 1e-10] = Zero() + + # Recombine the physical basis functions to satify the constraints + Gnt = numpy.asarray(Rnt[1:]) + vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) + CV = C @ V[:, vids] + c0, c1 = fdofs[-2:] + V[c0, vids] = -1 * Gnt @ CV[[0, 1]] + V[c1, vids] = -1 * Gnt @ CV[[1, 2]] # Patch up conditioning h = coordinate_mapping.cell_size() @@ -114,27 +121,3 @@ def entity_dofs(self): def space_dimension(self): return 45 - - -def morley_transform(cell, J, detJ, face): - adjugate = lambda A: ListTensor([[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]]) - sd = cell.get_spatial_dimension() - thats = cell.compute_tangents(sd-1, face) - nhat = numpy.cross(*thats) - ahat = numpy.linalg.norm(nhat) - nhat /= numpy.dot(nhat, nhat) - - nhat = Literal(nhat) - Jn = J @ nhat - Jt = J @ Literal(thats.T) - Gnt = Jn.T @ Jt - Gtt = Jt.T @ Jt - detG = Gtt[0, 0]*Gtt[1, 1] - Gtt[0, 1]*Gtt[1, 0] - area = Power(detG, Literal(0.5)) - - Bnn = detJ / area - Bnt = Gnt @ adjugate(Gtt) / detG - Bnn *= ahat - Bnt *= ahat - Bnt = (-1*(Bnt[0] + Bnt[1]), Bnt[0], Bnt[1]) - return Bnn, Bnt, Jn diff --git a/test/finat/conftest.py b/test/finat/conftest.py index 51f4b3106..56a32c8ea 100644 --- a/test/finat/conftest.py +++ b/test/finat/conftest.py @@ -115,10 +115,6 @@ def phys_el(): (0.01, 0.987, -.23), (-0.1, -0.2, 1.38)) - K[3].vertices = ((0.0, 0.0, 0.0), - (1.0, 0.0, 0.0), - (0.0, 1.0, 0.0), - (0.0, 0.0, 2.0)) return K From a0339cec9071dc93b94dee7081be28376ae85064 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 3 Jan 2026 14:47:42 -0600 Subject: [PATCH 21/25] register element --- FIAT/walkington.py | 21 ++++++++------------- finat/element_factory.py | 1 + finat/ufl/elementlist.py | 1 + finat/walkington.py | 16 +++++++++------- 4 files changed, 19 insertions(+), 20 deletions(-) diff --git a/FIAT/walkington.py b/FIAT/walkington.py index 3a4b24231..e9d2cf6c9 100644 --- a/FIAT/walkington.py +++ b/FIAT/walkington.py @@ -9,6 +9,7 @@ # bfs, but the extra 20 are used in the transformation theory. from FIAT import finite_element, polynomial_set, dual_set, macro +from FIAT.expansions import polynomial_dimension from FIAT.functional import ( PointEvaluation, PointDerivative, IntegralMomentOfNormalDerivative, @@ -21,10 +22,6 @@ import numpy -def inner(u, v, wts): - return numpy.dot(numpy.multiply(u, wts), v.T) - - class WalkingtonDualSet(dual_set.DualSet): def __init__(self, ref_el, degree): top = ref_el.get_topology() @@ -67,7 +64,7 @@ def __init__(self, ref_el, degree): x = ref_edge.compute_barycentric_coordinates(Q_edge.get_points()) leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0]) # Face constraint: normal derivative is cubic - Q_face, phi = face_constraints(ref_face, degree) + Q_face, phi = face_constraint(ref_face) for face in sorted(top[2]): cur = len(nodes) @@ -89,8 +86,6 @@ def __init__(self, ref_el, degree): nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, phi, nface, t) for t in thats) entity_ids[2][face].extend(range(cur, len(nodes))) - self.Q_face = Q_face - self.phi = phi super().__init__(nodes, ref_el, entity_ids) @@ -109,17 +104,17 @@ def __init__(self, ref_el, degree=5): super().__init__(poly_set, dual, degree) -def face_constraints(ref_face, degree): - from FIAT.expansions import polynomial_dimension +def face_constraint(ref_face): + k = 3 sd = ref_face.get_spatial_dimension() - Q_face = create_quadrature(ref_face, 2*(degree-2)) - dimPkm1 = polynomial_dimension(ref_face, degree-3) + Q_face = create_quadrature(ref_face, 2*k) + dimPkm1 = polynomial_dimension(ref_face, k-1) pts = list(Q_face.get_points()[:3]) pts.append(Q_face.get_points()[-1]) - P = polynomial_set.ONPolynomialSet(ref_face, degree-2, scale=1) + P = polynomial_set.ONPolynomialSet(ref_face, k) Pk = P.tabulate(pts)[(0,)*sd][dimPkm1:] - c = numpy.linalg.solve(Pk.T, [0, 0, 0, 1/Q_face.get_weights()[-1]]) + c = numpy.linalg.solve(Pk.T, [0, 0, 0, 1]) Pk = P.tabulate(Q_face.get_points())[(0,)*sd][dimPkm1:] phi = numpy.dot(c, Pk) return Q_face, phi diff --git a/finat/element_factory.py b/finat/element_factory.py index b688ba09e..be49975dc 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -90,6 +90,7 @@ "Nonconforming Arnold-Winther": finat.ArnoldWintherNC, "Hu-Zhang": finat.HuZhang, "Mardal-Tai-Winther": finat.MardalTaiWinther, + "Walkington": finat.Walkington, # These require special treatment "Q": None, "DQ": None, diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index 6293b39e6..8af76008b 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -123,6 +123,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("Walkington", "WALK", 0, H2, "custom", (5, 5), ("tetrahedron",)) register_element("Arnold-Qin", "AQ", 1, H1, "identity", (2, 2), ("triangle",)) register_element("Reduced-Arnold-Qin", "AQ-red", 1, H1, "contravariant Piola", (2, 2), ("triangle",)) diff --git a/finat/walkington.py b/finat/walkington.py index 01d70cca4..d98c5c63b 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -1,6 +1,7 @@ import FIAT import numpy +from FIAT.walkington import face_constraint from FIAT.polynomial_set import mis from gem import ListTensor, Zero @@ -42,8 +43,7 @@ def basis_transformation(self, coordinate_mapping): entity_dofs = self._element.entity_dofs() edges = self.cell.get_connectivity()[(2, 1)] - Q_face = self._element.dual.Q_face - phi = self._element.dual.phi + Q_face, phi = face_constraint(self.cell.construct_subelement(2)) wts = numpy.multiply(Q_face.get_weights(), phi) pts = Q_face.get_points() pts = pts[abs(wts) > 1E-10] @@ -53,7 +53,6 @@ def basis_transformation(self, coordinate_mapping): Rnn, Rnt = morley_transform(self.cell, J, detJ, f) fdofs = entity_dofs[2][f] fid = fdofs[0] - V[fid, fid] = Rnn for j, e in enumerate(edges[f]): Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e, face=f) @@ -86,6 +85,8 @@ def basis_transformation(self, coordinate_mapping): V[eid, vid1+i] = 1/252 * Bnt * tau V[eid, vid0+i] = -1 * V[eid, vid1+i] + vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) + # Evaluate nodal completion of the face constraints tab = self._element.tabulate(2, pts, entity=(2, f)) thats = self.cell.compute_tangents(2, f) @@ -93,13 +94,14 @@ def basis_transformation(self, coordinate_mapping): T = [[a[i, j] + (i != j) * a[j, i] for i in range(sd) for j in range(i, sd)] for a in tt] C = numpy.array([tab[alpha] @ wts for alpha in tab if sum(alpha) == 2]) C = numpy.dot(T, C) + C[abs(C) < 1e-10] = 0 + supp = numpy.unique(numpy.nonzero(C)[1]) C = C.astype(object) - C[abs(C) < 1e-10] = Zero() + C[C == 0] = Zero() + CV = C[:, supp] @ V[numpy.ix_(supp, vids)] - # Recombine the physical basis functions to satify the constraints + # Add completion basis functions to satisfy the physical constraints Gnt = numpy.asarray(Rnt[1:]) - vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) - CV = C @ V[:, vids] c0, c1 = fdofs[-2:] V[c0, vids] = -1 * Gnt @ CV[[0, 1]] V[c1, vids] = -1 * Gnt @ CV[[1, 2]] From 36fdd4b5daa4cd333ce3af3cc4f011e1349a39ac Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 3 Jan 2026 23:47:29 -0600 Subject: [PATCH 22/25] cleanup --- finat/walkington.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/finat/walkington.py b/finat/walkington.py index d98c5c63b..b0ee80a79 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -85,22 +85,20 @@ def basis_transformation(self, coordinate_mapping): V[eid, vid1+i] = 1/252 * Bnt * tau V[eid, vid0+i] = -1 * V[eid, vid1+i] - vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) - # Evaluate nodal completion of the face constraints tab = self._element.tabulate(2, pts, entity=(2, f)) thats = self.cell.compute_tangents(2, f) tt = [numpy.outer(thats[i], thats[j]) for i in range(2) for j in range(i, 2)] T = [[a[i, j] + (i != j) * a[j, i] for i in range(sd) for j in range(i, sd)] for a in tt] - C = numpy.array([tab[alpha] @ wts for alpha in tab if sum(alpha) == 2]) - C = numpy.dot(T, C) + C = numpy.dot(T, [tab[alpha] @ wts for alpha in tab if sum(alpha) == 2]) C[abs(C) < 1e-10] = 0 supp = numpy.unique(numpy.nonzero(C)[1]) C = C.astype(object) C[C == 0] = Zero() - CV = C[:, supp] @ V[numpy.ix_(supp, vids)] - # Add completion basis functions to satisfy the physical constraints + # Recombine with nodal completion to satisfy the physical constraints + vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) + CV = C[:, supp] @ V[numpy.ix_(supp, vids)] Gnt = numpy.asarray(Rnt[1:]) c0, c1 = fdofs[-2:] V[c0, vids] = -1 * Gnt @ CV[[0, 1]] From 64907c0dc86c1331633c86fe0c47951ee3b12de1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 7 Jan 2026 20:56:55 -0600 Subject: [PATCH 23/25] add mass conditioning test --- test/finat/conftest.py | 1 - test/finat/test_mass_conditioning.py | 26 +++++++++++++------------- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/test/finat/conftest.py b/test/finat/conftest.py index 56a32c8ea..c31ef7350 100644 --- a/test/finat/conftest.py +++ b/test/finat/conftest.py @@ -114,7 +114,6 @@ def phys_el(): (1., 0.1, -0.37), (0.01, 0.987, -.23), (-0.1, -0.2, 1.38)) - return K diff --git a/test/finat/test_mass_conditioning.py b/test/finat/test_mass_conditioning.py index c6de89cc9..60b270e30 100644 --- a/test/finat/test_mass_conditioning.py +++ b/test/finat/test_mass_conditioning.py @@ -4,20 +4,20 @@ from gem.interpreter import evaluate -@pytest.mark.parametrize("element, degree, variant", [ - (finat.Hermite, 3, None), - (finat.QuadraticPowellSabin6, 2, None), - (finat.QuadraticPowellSabin12, 2, None), - (finat.ReducedHsiehCloughTocher, 3, None), - (finat.HsiehCloughTocher, 3, None), - (finat.HsiehCloughTocher, 4, None), - (finat.Bell, 5, None), - (finat.Argyris, 5, "point"), - (finat.Argyris, 5, None), - (finat.Argyris, 6, None), +@pytest.mark.parametrize("sd,element,degree,variant", [ + (2, finat.Hermite, 3, None), + (2, finat.QuadraticPowellSabin6, 2, None), + (2, finat.QuadraticPowellSabin12, 2, None), + (2, finat.ReducedHsiehCloughTocher, 3, None), + (2, finat.HsiehCloughTocher, 3, None), + (2, finat.HsiehCloughTocher, 4, None), + (2, finat.Bell, 5, None), + (2, finat.Argyris, 5, "point"), + (2, finat.Argyris, 5, None), + (2, finat.Argyris, 6, None), + (3, finat.Walkington, 5, None), ]) -def test_mass_scaling(scaled_ref_to_phys, element, degree, variant): - sd = 2 +def test_mass_scaling(scaled_ref_to_phys, sd, element, degree, variant): ref_cell = scaled_ref_to_phys[sd][0].ref_cell if variant is not None: ref_element = element(ref_cell, degree, variant=variant) From b966b4963334629e4a6a13a20a32b2c82aa58e18 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 8 Jan 2026 16:16:38 -0600 Subject: [PATCH 24/25] style --- FIAT/walkington.py | 36 ++++++++++++++++++++++++++---------- finat/walkington.py | 24 +++++++++--------------- 2 files changed, 35 insertions(+), 25 deletions(-) diff --git a/FIAT/walkington.py b/FIAT/walkington.py index e9d2cf6c9..4d2a5c0a7 100644 --- a/FIAT/walkington.py +++ b/FIAT/walkington.py @@ -8,7 +8,8 @@ # functionals. The first 45 basis functions are the reference element # bfs, but the extra 20 are used in the transformation theory. -from FIAT import finite_element, polynomial_set, dual_set, macro +from FIAT import finite_element, polynomial_set, macro +from FIAT.dual_set import DualSet from FIAT.expansions import polynomial_dimension from FIAT.functional import ( PointEvaluation, PointDerivative, @@ -16,13 +17,13 @@ IntegralMomentOfBidirectionalDerivative, ) from FIAT.reference_element import TETRAHEDRON -from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature import FacetQuadratureRule, QuadratureRule from FIAT.quadrature_schemes import create_quadrature from FIAT.jacobi import eval_jacobi import numpy -class WalkingtonDualSet(dual_set.DualSet): +class WalkingtonDualSet(DualSet): def __init__(self, ref_el, degree): top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() @@ -66,6 +67,9 @@ def __init__(self, ref_el, degree): # Face constraint: normal derivative is cubic Q_face, phi = face_constraint(ref_face) + extra_entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + extra_nodes = [] + for face in sorted(top[2]): cur = len(nodes) thats = ref_el.compute_tangents(sd-1, face) @@ -82,10 +86,16 @@ def __init__(self, ref_el, degree): nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q, leg4_at_qpts, n=nfe)) Q = FacetQuadratureRule(ref_el, 2, face, Q_face) - nface /= Q.jacobian_determinant() - nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, phi, nface, t) for t in thats) + f = phi / Q.jacobian_determinant() + nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, f, nface, t) for t in thats) entity_ids[2][face].extend(range(cur, len(nodes))) + cur = len(extra_nodes) + extra_nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, f, thats[i], thats[j]) + for i in range(2) for j in range(i, 2)) + extra_entity_ids[2][face].extend(range(cur, len(extra_nodes))) + + self.nodal_completion = DualSet(extra_nodes, ref_el, extra_entity_ids) super().__init__(nodes, ref_el, entity_ids) @@ -107,14 +117,20 @@ def __init__(self, ref_el, degree=5): def face_constraint(ref_face): k = 3 sd = ref_face.get_spatial_dimension() - Q_face = create_quadrature(ref_face, 2*k) + Q = create_quadrature(ref_face, 2*k) dimPkm1 = polynomial_dimension(ref_face, k-1) - pts = list(Q_face.get_points()[:3]) - pts.append(Q_face.get_points()[-1]) + pts = list(Q.get_points()[:3]) + pts.append(Q.get_points()[-1]) P = polynomial_set.ONPolynomialSet(ref_face, k) Pk = P.tabulate(pts)[(0,)*sd][dimPkm1:] c = numpy.linalg.solve(Pk.T, [0, 0, 0, 1]) - Pk = P.tabulate(Q_face.get_points())[(0,)*sd][dimPkm1:] + Pk = P.tabulate(Q.get_points())[(0,)*sd][dimPkm1:] phi = numpy.dot(c, Pk) - return Q_face, phi + + supp = abs(phi) > 1E-12 + pts = Q.get_points()[supp] + wts = Q.get_weights()[supp] + Q = QuadratureRule(ref_face, pts, wts) + phi = phi[supp] + return Q, phi diff --git a/finat/walkington.py b/finat/walkington.py index b0ee80a79..b9d11ae38 100644 --- a/finat/walkington.py +++ b/finat/walkington.py @@ -1,7 +1,6 @@ import FIAT import numpy -from FIAT.walkington import face_constraint from FIAT.polynomial_set import mis from gem import ListTensor, Zero @@ -43,11 +42,12 @@ def basis_transformation(self, coordinate_mapping): entity_dofs = self._element.entity_dofs() edges = self.cell.get_connectivity()[(2, 1)] - Q_face, phi = face_constraint(self.cell.construct_subelement(2)) - wts = numpy.multiply(Q_face.get_weights(), phi) - pts = Q_face.get_points() - pts = pts[abs(wts) > 1E-10] - wts = wts[abs(wts) > 1E-10] + # Evaluate nodal completion of the face constraints + P = self._element.poly_set + L = self._element.dual.nodal_completion + coeffs = P.get_coeffs() + tangential_dofs = numpy.dot(L.to_riesz(P), coeffs.T) + tangential_dofs[abs(tangential_dofs) < 1e-10] = 0 for f in entity_dofs[2]: Rnn, Rnt = morley_transform(self.cell, J, detJ, f) @@ -85,19 +85,13 @@ def basis_transformation(self, coordinate_mapping): V[eid, vid1+i] = 1/252 * Bnt * tau V[eid, vid0+i] = -1 * V[eid, vid1+i] - # Evaluate nodal completion of the face constraints - tab = self._element.tabulate(2, pts, entity=(2, f)) - thats = self.cell.compute_tangents(2, f) - tt = [numpy.outer(thats[i], thats[j]) for i in range(2) for j in range(i, 2)] - T = [[a[i, j] + (i != j) * a[j, i] for i in range(sd) for j in range(i, sd)] for a in tt] - C = numpy.dot(T, [tab[alpha] @ wts for alpha in tab if sum(alpha) == 2]) - C[abs(C) < 1e-10] = 0 + vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) + # Recombine with nodal completion to satisfy the physical constraints + C = tangential_dofs[L.entity_ids[2][f]] supp = numpy.unique(numpy.nonzero(C)[1]) C = C.astype(object) C[C == 0] = Zero() - # Recombine with nodal completion to satisfy the physical constraints - vids = list(chain.from_iterable(entity_dofs[0][v] for v in top[2][f])) CV = C[:, supp] @ V[numpy.ix_(supp, vids)] Gnt = numpy.asarray(Rnt[1:]) c0, c1 = fdofs[-2:] From 95adf367308c0e24fffdb67368172aba67b4fac3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 12 Jan 2026 08:56:06 -0600 Subject: [PATCH 25/25] tidy --- finat/morley.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/finat/morley.py b/finat/morley.py index ed58ed63e..ca3e89e88 100644 --- a/finat/morley.py +++ b/finat/morley.py @@ -16,8 +16,7 @@ def morley_transform(cell, J, detJ, face): ahat = numpy.linalg.norm(nhat) nhat /= numpy.dot(nhat, nhat) - nhat = Literal(nhat) - Jn = J @ nhat + Jn = J @ Literal(nhat) Jt = J @ Literal(thats.T) Gnt = Jn.T @ Jt Gtt = Jt.T @ Jt