From bdbe08e0ef4184d98ab9dd73b25bbe4bb5d8ef92 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 7 Jul 2025 13:57:44 +0100 Subject: [PATCH 01/17] Macro H(div) elements --- FIAT/brezzi_douglas_marini.py | 36 +++++++++------ FIAT/macro.py | 87 ++++++++++++++++++++++------------- 2 files changed, 79 insertions(+), 44 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 0014b5277..5d46c01d0 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -6,10 +6,11 @@ # SPDX-License-Identifier: LGPL-3.0-or-later from FIAT import (finite_element, functional, dual_set, - polynomial_set, nedelec) -from FIAT.check_format_variant import check_format_variant + polynomial_set, nedelec, macro) +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule +import numpy class BDMDualSet(dual_set.DualSet): @@ -26,7 +27,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): entity_ids[dim][entity] = [] if variant == "integral": - facet = ref_el.get_facet_element() + facet = ref_el.construct_subelement(sd-1) # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) @@ -56,14 +57,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): if degree > 1: if interpolant_deg is None: interpolant_deg = degree - cur = len(nodes) - Q = create_quadrature(ref_el, interpolant_deg + degree - 1) - Nedel = nedelec.Nedelec(ref_el, degree - 1, variant) + + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, interpolant_deg + degree - 1) + Nedel = nedelec.Nedelec(cell, degree - 1, variant) Nedfs = Nedel.get_nodal_basis() - Ned_at_qpts = Nedfs.tabulate(Q.get_points())[(0,) * sd] - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in Ned_at_qpts) - entity_ids[sd][0] = list(range(cur, len(nodes))) + Ned_at_qpts = Nedfs.tabulate(Q_ref.get_points())[(0,) * sd] + for entity in top[sd]: + Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) + JinvT = numpy.linalg.inv(Q.jacobian().T) + cur = len(nodes) + nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, JinvT @ phi) for phi in Ned_at_qpts) + entity_ids[sd][entity] = list(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -89,14 +94,19 @@ class BrezziDouglasMarini(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) variant, interpolant_deg = check_format_variant(variant, degree) - if degree < 1: raise Exception("BDM_k elements only valid for k >= 1") sd = ref_el.get_spatial_dimension() - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) + if ref_el.is_macrocell(): + poly_set = macro.HDivPolynomialSet(ref_el, degree) + else: + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) dual = BDMDualSet(ref_el, degree, variant, interpolant_deg) formdegree = sd - 1 # (n-1)-form super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/macro.py b/FIAT/macro.py index 07a809acb..08b633eb2 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -485,8 +485,46 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): super().__init__(ref_el, degree, degree, expansion_set, coeffs) -class HDivSymPolynomialSet(polynomial_set.PolynomialSet): - """Constructs a symmetric tensor-valued PolynomialSet with continuous +def hdiv_conforming(U, order=0): + from FIAT.quadrature_schemes import create_quadrature + degree = U.degree + ref_el = U.get_reference_element() + coeffs = U.get_coeffs() + shape = U.get_shape() + expansion_set = U.get_expansion_set() + k = 1 if expansion_set.continuity == "C0" else 0 + + sd = ref_el.get_spatial_dimension() + facet_el = ref_el.construct_subelement(sd-1) + + phi_deg = 0 if sd == 1 else degree - k + phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg, shape=shape[1:]) + Q = create_quadrature(facet_el, 2 * phi_deg) + qpts, qwts = Q.get_points(), Q.get_weights() + phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] + weights = numpy.multiply(phi_at_qpts, qwts) + + rows = [] + for facet in ref_el.get_interior_facets(sd-1): + normal = ref_el.compute_normal(facet) + normal = normal[(None,) * len(shape) + (slice(None), None)] + jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) + for r in range(k, order+1): + jump = numpy.dot(coeffs, jumps[r]) + # num_wt = 1 if sd == 1 else expansions.polynomial_dimension(facet_el, degree-r) + wn = weights[..., None, :] * normal + ax = tuple(range(1, len(wn.shape))) + rows.append(numpy.tensordot(wn, jump, axes=(ax, ax))) + + if len(rows) > 0: + dual_mat = numpy.vstack(rows) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0)) + return coeffs + + +class HDivPolynomialSet(polynomial_set.PolynomialSet): + """Constructs a vector-valued PolynomialSet with continuous normal components on a simplicial complex. :arg ref_el: The simplicial complex. @@ -496,36 +534,23 @@ class HDivSymPolynomialSet(polynomial_set.PolynomialSet): :kwarg scale: The scale for the underlying ExpansionSet. """ def __init__(self, ref_el, degree, order=0, **kwargs): - from FIAT.quadrature_schemes import create_quadrature - U = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree, **kwargs) - coeffs = U.get_coeffs() - expansion_set = U.get_expansion_set() - k = 1 if expansion_set.continuity == "C0" else 0 - sd = ref_el.get_spatial_dimension() - facet_el = ref_el.construct_subelement(sd-1) + U = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), **kwargs) + coeffs = hdiv_conforming(U, order=order) + super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) - phi_deg = 0 if sd == 1 else degree - k - phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg, shape=(sd,)) - Q = create_quadrature(facet_el, 2 * phi_deg) - qpts, qwts = Q.get_points(), Q.get_weights() - phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] - weights = numpy.multiply(phi_at_qpts, qwts) - - rows = [] - for facet in ref_el.get_interior_facets(sd-1): - normal = ref_el.compute_normal(facet) - jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) - for r in range(k, order+1): - jump = numpy.dot(coeffs, jumps[r]) - # num_wt = 1 if sd == 1 else expansions.polynomial_dimension(facet_el, degree-r) - wn = weights[:, :, None, :] * normal[None, None, :, None] - ax = tuple(range(1, len(wn.shape))) - rows.append(numpy.tensordot(wn, jump, axes=(ax, ax))) - if len(rows) > 0: - dual_mat = numpy.vstack(rows) - nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0)) +class HDivSymPolynomialSet(polynomial_set.PolynomialSet): + """Constructs a symmetric tensor-valued PolynomialSet with continuous + normal components on a simplicial complex. - super().__init__(ref_el, degree, degree, expansion_set, coeffs) + :arg ref_el: The simplicial complex. + :arg degree: The polynomial degree. + :kwarg order: The order of continuity across subcells. + :kwarg variant: The variant for the underlying ExpansionSet. + :kwarg scale: The scale for the underlying ExpansionSet. + """ + def __init__(self, ref_el, degree, order=0, **kwargs): + U = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree, **kwargs) + coeffs = hdiv_conforming(U, order=order) + super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) From c5a82e0bbd46d647fec9999b617b15056e52968f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Jul 2025 11:28:18 +0100 Subject: [PATCH 02/17] MacroExpansionSet --- FIAT/expansions.py | 17 +++++++--- FIAT/macro.py | 85 +++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 96 insertions(+), 6 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 69d0d4b5a..fbc4658a7 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -371,18 +371,20 @@ def _tabulate(self, n, pts, order=0): phi[alpha] /= mult[None, ipts] # Insert subcell tabulations into the corresponding submatrices - idx = lambda *args: args if args[1] is Ellipsis else numpy.ix_(*args) + idx = lambda *args: args if args[-1] is Ellipsis else numpy.ix_(*args) num_phis = self.get_num_members(n) cell_node_map = self.get_cell_node_map(n) result = {} base_phi = tuple(phis.values())[0] for alpha in base_phi: dtype = base_phi[alpha].dtype - result[alpha] = numpy.zeros((num_phis, *pts.shape[:-1]), dtype=dtype) + phi_shape = (num_phis,) + base_phi[alpha].shape[1:-1] + shape_idx = tuple(range(i) for i in phi_shape[1:]) + result[alpha] = numpy.zeros(phi_shape + pts.shape[:-1], dtype=dtype) for cell in cell_point_map: ibfs = cell_node_map[cell] ipts = cell_point_map[cell] - result[alpha][idx(ibfs, ipts)] += phis[cell][alpha] + result[alpha][idx(ibfs, *shape_idx, ipts)] += phis[cell][alpha] return result def tabulate_normal_jumps(self, n, ref_pts, facet, order=0): @@ -599,7 +601,10 @@ def polynomial_dimension(ref_el, n, continuity=None): raise ValueError("Only degree zero polynomials supported on point elements.") return 1 top = ref_el.get_topology() - if continuity == "C0": + + if isinstance(continuity, dict): + space_dimension = sum(len(continuity[dim][0]) * len(top[dim]) for dim in top) + elif continuity == "C0": space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top) else: dim = ref_el.get_spatial_dimension() @@ -620,7 +625,9 @@ def polynomial_entity_ids(ref_el, n, continuity=None): entity_ids = {} cur = 0 for dim in sorted(top): - if continuity == "C0": + if isinstance(continuity, dict): + dofs, = set(len(continuity[dim][entity]) for entity in continuity[dim]) + elif continuity == "C0": dofs = math.comb(n - 1, dim) else: # DG numbering diff --git a/FIAT/macro.py b/FIAT/macro.py index 08b633eb2..6018feb1c 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -2,7 +2,7 @@ import numpy -from FIAT import expansions, polynomial_set +from FIAT import expansions, polynomial_set, reference_element from FIAT.quadrature import FacetQuadratureRule, QuadratureRule from FIAT.reference_element import (TRIANGLE, SimplicialComplex, lattice_iter, make_lattice) @@ -554,3 +554,86 @@ def __init__(self, ref_el, degree, order=0, **kwargs): U = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree, **kwargs) coeffs = hdiv_conforming(U, order=order) super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) + + +class MacroExpansionSet(expansions.ExpansionSet): + + def __init__(self, ref_el, element): + self.element = element + self.ref_el = ref_el + self.variant = None + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + base_ref_el = element.ref_el + base_verts = base_ref_el.get_vertices() + + self.affine_mappings = [reference_element.make_affine_mapping( + ref_el.get_vertices_of_subcomplex(top[sd][cell]), + base_verts) for cell in top[sd]] + self.scale = 1 + self.continuity = element.entity_dofs() + self.recurrence_order = element.degree + self._dmats_cache = {} + self._cell_node_map_cache = {} + + def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): + A, b = self.affine_mappings[cell] + ref_pts = numpy.add(numpy.dot(pts, A.T), b) + Jinv = A if direction is None else numpy.dot(A, direction)[:, None] + phis = self.element.tabulate(order, ref_pts) + for alpha in phis: + phis[alpha] = self.pullback(Jinv, phis[alpha]) + return phis + + def pullback(self, Jinv, phi): + mapping = self.element.mapping()[0] + if mapping == "covariant piola": + phi = numpy.tensordot(Jinv.T, phi, (1, 1)).transpose((1, 0, 2)) + elif mapping == "contravariant piola": + J = numpy.linalg.inv(Jinv) + Jdet = abs(numpy.linalg.det(J)) + phi = numpy.tensordot(J/Jdet, phi, (1, 1)).transpose((1, 0, 2)) + elif mapping != "identity": + raise ValueError(f"Unrecognized mapping {mapping}") + + return phi + + +if __name__ == "__main__": + from FIAT import Nedelec, RaviartThomas, ufc_simplex + from FIAT.reference_element import symmetric_simplex + + dim = 3 + K = symmetric_simplex(dim) + # A = AlfeldSplit(K) + A = IsoSplit(K) + + degree = 1 + fe = RaviartThomas(K, degree) + # fe = Nedelec(K, degree) + + U = MacroExpansionSet(A, fe) + + mapping = fe.mapping()[0] + fdim = fe.formdegree + comps = [] + pts = [] + for entity in A.topology[fdim]: + pts_cur = A.make_points(fdim, entity, degree+fdim) + pts.extend(pts_cur) + if mapping == "covariant piola": + comp = A.compute_edge_tangent(entity) + elif mapping == "contravariant piola": + comp = A.compute_scaled_normal(entity) + comps.extend([comp] * len(pts_cur)) + + phis = U._tabulate(degree, pts) + for alpha in phis: + shape = phis[alpha].shape + result = numpy.zeros((shape[0], shape[-1])) + for i, comp in enumerate(comps): + result[..., i] = numpy.dot(phis[alpha][..., i], comp) + result[abs(result) < 1E-14] = 0 + phis[alpha] = result + + print(phis) From 5a0fb483575ce7b60f11d3e36c84fe6a2bc878f5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Jul 2025 15:07:19 +0100 Subject: [PATCH 03/17] BDM: use consistent normals --- FIAT/brezzi_douglas_marini.py | 27 +++++++++++++++------------ FIAT/macro.py | 3 +-- test/finat/test_zany_mapping.py | 8 ++++++++ 3 files changed, 24 insertions(+), 14 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 5d46c01d0..c0035bc8d 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -18,26 +18,28 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): nodes = [] sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - entity_ids = {} # set to empty - for dim in top: + for dim in sorted(top): entity_ids[dim] = {} - for entity in top[dim]: + for entity in sorted(top[dim]): entity_ids[dim][entity] = [] if variant == "integral": + R = numpy.array([[0, 1], [-1, 0]]) facet = ref_el.construct_subelement(sd-1) # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) - Pq = polynomial_set.ONPolynomialSet(facet, degree) + Pq = polynomial_set.ONPolynomialSet(facet, degree, scale=1) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] - for f in top[sd - 1]: + for f in sorted(top[sd - 1]): cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet + thats = ref_el.compute_tangents(sd-1, f) + nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats) + n = nhat / Jdet phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) @@ -46,7 +48,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): elif variant == "point": # Define each functional for the dual set # codimension 1 facets - for f in top[sd - 1]: + for f in sorted(top[sd - 1]): cur = len(nodes) pts_cur = ref_el.make_points(sd - 1, f, sd + degree) nodes.extend(functional.PointScaledNormalEvaluation(ref_el, f, pt) @@ -61,13 +63,13 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cell = ref_el.construct_subelement(sd) Q_ref = create_quadrature(cell, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(cell, degree - 1, variant) - Nedfs = Nedel.get_nodal_basis() - Ned_at_qpts = Nedfs.tabulate(Q_ref.get_points())[(0,) * sd] - for entity in top[sd]: + Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd] + for entity in sorted(top[sd]): Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - JinvT = numpy.linalg.inv(Q.jacobian().T) + F = numpy.linalg.inv(Q.jacobian().T) + phis = numpy.tensordot(F, Ned_at_qpts, (1, 1)).transpose((1, 0, 2)) cur = len(nodes) - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, JinvT @ phi) for phi in Ned_at_qpts) + nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) entity_ids[sd][entity] = list(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -98,6 +100,7 @@ def __init__(self, ref_el, degree, variant=None): splitting, variant = parse_lagrange_variant(variant, integral=True) if splitting is not None: ref_el = splitting(ref_el) + variant, interpolant_deg = check_format_variant(variant, degree) if degree < 1: raise Exception("BDM_k elements only valid for k >= 1") diff --git a/FIAT/macro.py b/FIAT/macro.py index 6018feb1c..3dc7b5b5a 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -600,12 +600,11 @@ def pullback(self, Jinv, phi): if __name__ == "__main__": - from FIAT import Nedelec, RaviartThomas, ufc_simplex + from FIAT import RaviartThomas from FIAT.reference_element import symmetric_simplex dim = 3 K = symmetric_simplex(dim) - # A = AlfeldSplit(K) A = IsoSplit(K) degree = 1 diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index d5136e651..290329ba7 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -164,3 +164,11 @@ def test_piola_triangle_high_order(ref_to_phys, element, degree, variant): @pytest.mark.parametrize("dimension", [2, 3]) def test_affine(ref_to_phys, element, degree, dimension): check_zany_mapping(element, ref_to_phys[dimension], degree) + + +@pytest.mark.parametrize("element", [finat.BrezziDouglasMarini]) +@pytest.mark.parametrize("degree", [1, 2]) +@pytest.mark.parametrize("dimension", [2, 3]) +@pytest.mark.parametrize("variant", [None, "iso"]) +def test_macro_piola(ref_to_phys, element, degree, variant, dimension): + check_zany_mapping(element, ref_to_phys[dimension], degree, variant=variant) From bd9c33b0a5f690a26ca9117c83ab3ca28c05f228 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Jul 2025 16:27:59 +0100 Subject: [PATCH 04/17] Fix variants --- FIAT/brezzi_douglas_marini.py | 13 ++++++------- FIAT/check_format_variant.py | 4 +++- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index c0035bc8d..9c972197c 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -31,14 +31,14 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) - Pq = polynomial_set.ONPolynomialSet(facet, degree, scale=1) + Pq = polynomial_set.ONPolynomialSet(facet, degree) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in sorted(top[sd - 1]): cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) Jdet = Q.jacobian_determinant() thats = ref_el.compute_tangents(sd-1, f) - nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats) + nhat = numpy.dot(R, *thats) if sd == 2 else -numpy.cross(*thats) n = nhat / Jdet phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) @@ -66,8 +66,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd] for entity in sorted(top[sd]): Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - F = numpy.linalg.inv(Q.jacobian().T) - phis = numpy.tensordot(F, Ned_at_qpts, (1, 1)).transpose((1, 0, 2)) + Jinv = numpy.linalg.inv(Q.jacobian()) + phis = numpy.tensordot(Jinv.T, Ned_at_qpts, (1, 1)).transpose((1, 0, 2)) cur = len(nodes) nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) entity_ids[sd][entity] = list(range(cur, len(nodes))) @@ -96,14 +96,13 @@ class BrezziDouglasMarini(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): + if degree < 1: + raise Exception("BDM_k elements only valid for k >= 1") if variant is not None: splitting, variant = parse_lagrange_variant(variant, integral=True) if splitting is not None: ref_el = splitting(ref_el) - variant, interpolant_deg = check_format_variant(variant, degree) - if degree < 1: - raise Exception("BDM_k elements only valid for k >= 1") sd = ref_el.get_spatial_dimension() if ref_el.is_macrocell(): diff --git a/FIAT/check_format_variant.py b/FIAT/check_format_variant.py index 1431876de..9a50f3262 100644 --- a/FIAT/check_format_variant.py +++ b/FIAT/check_format_variant.py @@ -61,7 +61,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False): default = "integral" if integral else "spectral" if integral: - supported_point_variants = {"integral": None} + supported_point_variants = {"integral": None, "point": "point"} elif discontinuous: supported_point_variants = supported_dg_variants else: @@ -81,6 +81,8 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False): k, = match.groups() call_split = IsoSplit splitting_args = (int(k),) + elif opt.startswith("integral"): + point_variant = opt elif opt in supported_point_variants: point_variant = supported_point_variants[opt] else: From 013c5f8cb41f4620d996bf525162c1a4624508da Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Jul 2025 11:48:42 +0100 Subject: [PATCH 05/17] Always return UFC-oriented normals in compute_scaled_normal --- FIAT/brezzi_douglas_marini.py | 16 +++++++--------- FIAT/reference_element.py | 6 ++++++ 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 9c972197c..9ad9fa82f 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -18,28 +18,26 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): nodes = [] sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() + entity_ids = {} # set to empty - for dim in sorted(top): + for dim in top: entity_ids[dim] = {} - for entity in sorted(top[dim]): + for entity in top[dim]: entity_ids[dim][entity] = [] if variant == "integral": - R = numpy.array([[0, 1], [-1, 0]]) facet = ref_el.construct_subelement(sd-1) # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) Pq = polynomial_set.ONPolynomialSet(facet, degree) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] - for f in sorted(top[sd - 1]): + for f in top[sd - 1]: cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) Jdet = Q.jacobian_determinant() - thats = ref_el.compute_tangents(sd-1, f) - nhat = numpy.dot(R, *thats) if sd == 2 else -numpy.cross(*thats) - n = nhat / Jdet + n = ref_el.compute_scaled_normal(f) / Jdet phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) @@ -48,7 +46,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): elif variant == "point": # Define each functional for the dual set # codimension 1 facets - for f in sorted(top[sd - 1]): + for f in top[sd - 1]: cur = len(nodes) pts_cur = ref_el.make_points(sd - 1, f, sd + degree) nodes.extend(functional.PointScaledNormalEvaluation(ref_el, f, pt) @@ -64,7 +62,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(cell, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(cell, degree - 1, variant) Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd] - for entity in sorted(top[sd]): + for entity in top[sd]: Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) Jinv = numpy.linalg.inv(Q.jacobian()) phis = numpy.tensordot(Jinv.T, Ned_at_qpts, (1, 1)).transpose((1, 0, 2)) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 547da2c04..3baca43e9 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -548,6 +548,12 @@ def compute_scaled_normal(self, facet_i): """Returns the unit normal to facet_i of scaled by the volume of that facet.""" dim = self.get_spatial_dimension() + if dim == 2: + n, = self.compute_tangents(dim-1, facet_i) + n[0], n[1] = n[1], -n[0] + return n + elif dim == 3: + return -numpy.cross(*self.compute_tangents(dim-1, facet_i)) v = self.volume_of_subcomplex(dim - 1, facet_i) return self.compute_normal(facet_i) * v From 74ccf80f07ad50a07d8d47543cb65cb299aeb088 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 10 Jul 2025 16:08:02 +0100 Subject: [PATCH 06/17] MacroPolynomialSet --- FIAT/brezzi_douglas_marini.py | 2 + FIAT/expansions.py | 13 +++- FIAT/macro.py | 136 ++++++++++++++++++++++------------ 3 files changed, 98 insertions(+), 53 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 9ad9fa82f..c6c12869f 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -105,6 +105,8 @@ def __init__(self, ref_el, degree, variant=None): sd = ref_el.get_spatial_dimension() if ref_el.is_macrocell(): poly_set = macro.HDivPolynomialSet(ref_el, degree) + # base_element = BrezziDouglasMarini(ref_el.get_parent(), degree) + # poly_set = macro.MacroPolynomialSet(ref_el, base_element) else: poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) dual = BDMDualSet(ref_el, degree, variant, interpolant_deg) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index fbc4658a7..b7d846f26 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -274,11 +274,18 @@ def __init__(self, ref_el, scale=None, variant=None): if scale is None: scale = math.sqrt(1.0 / base_ref_el.volume()) self.scale = scale + self.variant = variant self.continuity = "C0" if variant == "bubble" else None self.recurrence_order = 2 self._dmats_cache = {} self._cell_node_map_cache = {} + def reconstruct(self, ref_el=None, scale=None, variant=None): + """Reconstructs this ExpansionSet with modified arguments.""" + return ExpansionSet(ref_el or self.ref_el, + scale=scale or self.scale, + variant=variant or self.variant) + def get_scale(self, n, cell=0): scale = self.scale sd = self.ref_el.get_spatial_dimension() @@ -378,13 +385,11 @@ def _tabulate(self, n, pts, order=0): base_phi = tuple(phis.values())[0] for alpha in base_phi: dtype = base_phi[alpha].dtype - phi_shape = (num_phis,) + base_phi[alpha].shape[1:-1] - shape_idx = tuple(range(i) for i in phi_shape[1:]) - result[alpha] = numpy.zeros(phi_shape + pts.shape[:-1], dtype=dtype) + result[alpha] = numpy.zeros((num_phis, *pts.shape[:-1]), dtype=dtype) for cell in cell_point_map: ibfs = cell_node_map[cell] ipts = cell_point_map[cell] - result[alpha][idx(ibfs, *shape_idx, ipts)] += phis[cell][alpha] + result[alpha][idx(ibfs, ipts)] += phis[cell][alpha] return result def tabulate_normal_jumps(self, n, ref_pts, facet, order=0): diff --git a/FIAT/macro.py b/FIAT/macro.py index 3dc7b5b5a..741e69d9c 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -556,82 +556,120 @@ def __init__(self, ref_el, degree, order=0, **kwargs): super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) -class MacroExpansionSet(expansions.ExpansionSet): - +def pullback(mapping, Jinv, phi): + try: + formdegree = { + "affine": (0,), + "covariant piola": (1,), + "contravariant piola": (2,), + "double covariant piola": (1, 1), + "double contravariant piola": (2, 2), + "covariant contravariant piola": (1, 2), + "contravariant covariant piola": (2, 1), }[mapping] + except KeyError: + raise ValueError(f"Unrecognized mapping {mapping}") + + F1 = Jinv.T + F2 = numpy.linalg.inv(Jinv) * numpy.linalg.det(Jinv) + + perm = [*range(1, phi.ndim-1), 0, -1] + phi = phi.transpose(perm) + + for i, k in enumerate(formdegree): + if k == 0: + continue + F = F1 if k == 1 else F2 + phi = numpy.tensordot(F, phi, (1, i)) + + perm = [-2, *range(phi.ndim-2), -1] + phi = phi.transpose(perm) + return phi + + +class MacroPolynomialSet(polynomial_set.PolynomialSet): + """Constructs a PolynomialSet by tiling a CiarletElement on a + SimplicialComplex. + + :arg ref_el: The SimplicialComplex. + :arg element: The CiarletElement. + """ def __init__(self, ref_el, element): - self.element = element - self.ref_el = ref_el - self.variant = None - sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - base_ref_el = element.ref_el - base_verts = base_ref_el.get_vertices() - - self.affine_mappings = [reference_element.make_affine_mapping( - ref_el.get_vertices_of_subcomplex(top[sd][cell]), - base_verts) for cell in top[sd]] - self.scale = 1 - self.continuity = element.entity_dofs() - self.recurrence_order = element.degree - self._dmats_cache = {} - self._cell_node_map_cache = {} - - def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): - A, b = self.affine_mappings[cell] - ref_pts = numpy.add(numpy.dot(pts, A.T), b) - Jinv = A if direction is None else numpy.dot(A, direction)[:, None] - phis = self.element.tabulate(order, ref_pts) - for alpha in phis: - phis[alpha] = self.pullback(Jinv, phis[alpha]) - return phis - - def pullback(self, Jinv, phi): - mapping = self.element.mapping()[0] - if mapping == "covariant piola": - phi = numpy.tensordot(Jinv.T, phi, (1, 1)).transpose((1, 0, 2)) - elif mapping == "contravariant piola": - J = numpy.linalg.inv(Jinv) - Jdet = abs(numpy.linalg.det(J)) - phi = numpy.tensordot(J/Jdet, phi, (1, 1)).transpose((1, 0, 2)) - elif mapping != "identity": - raise ValueError(f"Unrecognized mapping {mapping}") - - return phi + sd = ref_el.get_spatial_dimension() + + mapping, = set(element.mapping()) + base_ref_el = element.get_reference_element() + base_entity_ids = element.entity_dofs() + n = element.degree() + + base_expansion_set = element.get_nodal_basis().get_expansion_set() + expansion_set = base_expansion_set.reconstruct(ref_el=ref_el) + + base_entity_ids = expansions.polynomial_entity_ids(base_ref_el, n, base_entity_ids) + entity_ids = expansions.polynomial_entity_ids(ref_el, n, base_entity_ids) + + shp = element.value_shape() + num_bfs = expansions.polynomial_dimension(ref_el, n, entity_ids) + num_members = expansion_set.get_num_members(n) + coeffs = numpy.zeros((num_bfs, *shp, num_members)) + base_coeffs = element.get_coeffs() + + rmap = expansions.polynomial_cell_node_map(ref_el, n, base_entity_ids) + cmap = expansion_set.get_cell_node_map(n) + for cell in sorted(top): + cell_verts = ref_el.get_vertices_of_subcomplex(top[sd][cell]) + A, b = reference_element.make_affine_mapping(cell_verts, base_ref_el.vertices) + + indices = numpy.ix_(rmap[cell], *map(range, shp), cmap[cell]) + coeffs[indices] = pullback(mapping, A, base_coeffs) + + super().__init__(ref_el, element.degree(), element.degree(), expansion_set, coeffs) if __name__ == "__main__": - from FIAT import RaviartThomas + from FIAT import Lagrange, Nedelec, RaviartThomas, Regge from FIAT.reference_element import symmetric_simplex + degree = 1 dim = 3 K = symmetric_simplex(dim) - A = IsoSplit(K) + A = AlfeldSplit(K) - degree = 1 + fe = Lagrange(K, degree) + fe = Nedelec(K, degree) fe = RaviartThomas(K, degree) - # fe = Nedelec(K, degree) - - U = MacroExpansionSet(A, fe) + fe = Regge(K, 0) mapping = fe.mapping()[0] fdim = fe.formdegree + if isinstance(fdim, tuple): + fdim, = set(fdim) comps = [] pts = [] for entity in A.topology[fdim]: pts_cur = A.make_points(fdim, entity, degree+fdim) pts.extend(pts_cur) - if mapping == "covariant piola": + if mapping == "affine": + comp = numpy.ones(()) + elif mapping.endswith("covariant piola"): comp = A.compute_edge_tangent(entity) - elif mapping == "contravariant piola": + elif mapping.endswith("contravariant piola"): comp = A.compute_scaled_normal(entity) + if mapping.startswith("double"): + comp = numpy.outer(comp, comp) comps.extend([comp] * len(pts_cur)) - phis = U._tabulate(degree, pts) + P = MacroPolynomialSet(A, fe) + phis = P.tabulate(pts) + for alpha in phis: shape = phis[alpha].shape + shp = shape[1:-1] result = numpy.zeros((shape[0], shape[-1])) + + ax = (tuple(range(-len(shp), 0)), )*2 for i, comp in enumerate(comps): - result[..., i] = numpy.dot(phis[alpha][..., i], comp) + result[..., i] = numpy.tensordot(phis[alpha][..., i], comp, ax) result[abs(result) < 1E-14] = 0 phis[alpha] = result From 8449d7dbbe71f8f85f83b1995dc124e82a2e0769 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 11 Jul 2025 10:24:55 +0100 Subject: [PATCH 07/17] More macro variants --- FIAT/brezzi_douglas_marini.py | 11 ++-- FIAT/gopalakrishnan_lederer_schoberl.py | 69 ++++++++++++++----------- FIAT/macro.py | 59 ++------------------- FIAT/nedelec.py | 40 +++++++++----- FIAT/nedelec_second_kind.py | 40 ++++++++------ FIAT/raviart_thomas.py | 37 ++++++++----- FIAT/regge.py | 15 ++++-- test/FIAT/unit/test_macro.py | 59 ++++++++++++++++++++- test/finat/test_zany_mapping.py | 2 +- 9 files changed, 192 insertions(+), 140 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index c6c12869f..c9c944e75 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -94,19 +94,20 @@ class BrezziDouglasMarini(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - if degree < 1: - raise Exception("BDM_k elements only valid for k >= 1") + if variant is not None: splitting, variant = parse_lagrange_variant(variant, integral=True) if splitting is not None: ref_el = splitting(ref_el) variant, interpolant_deg = check_format_variant(variant, degree) + if degree < 1: + raise Exception("BDM_k elements only valid for k >= 1") + sd = ref_el.get_spatial_dimension() if ref_el.is_macrocell(): - poly_set = macro.HDivPolynomialSet(ref_el, degree) - # base_element = BrezziDouglasMarini(ref_el.get_parent(), degree) - # poly_set = macro.MacroPolynomialSet(ref_el, base_element) + base_element = BrezziDouglasMarini(ref_el.get_parent(), degree) + poly_set = macro.MacroPolynomialSet(ref_el, base_element) else: poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) dual = BDMDualSet(ref_el, degree, variant, interpolant_deg) diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py index 656cee6bb..d25637728 100644 --- a/FIAT/gopalakrishnan_lederer_schoberl.py +++ b/FIAT/gopalakrishnan_lederer_schoberl.py @@ -1,4 +1,5 @@ -from FIAT import finite_element, dual_set, polynomial_set, expansions +from FIAT import finite_element, dual_set, polynomial_set, expansions, macro +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant from FIAT.functional import TensorBidirectionalIntegralMoment as BidirectionalMoment from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule @@ -12,34 +13,29 @@ def __init__(self, ref_el, degree): nodes = [] entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} - # Face dofs: moments of normal-tangential components against a basis for Pk - ref_facet = ref_el.construct_subelement(sd-1) - Qref = create_quadrature(ref_facet, 2*degree) - P = polynomial_set.ONPolynomialSet(ref_facet, degree) - phis = P.tabulate(Qref.get_points())[(0,) * (sd-1)] - for f in sorted(top[sd-1]): - cur = len(nodes) - Q = FacetQuadratureRule(ref_el, sd-1, f, Qref) - Jdet = Q.jacobian_determinant() - normal = ref_el.compute_scaled_normal(f) - tangents = ref_el.compute_tangents(sd-1, f) - n = normal / Jdet - nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi) - for phi in phis for t in tangents) - entity_ids[sd-1][f].extend(range(cur, len(nodes))) - + # Face dofs: moments of normal-tangential components against a basis for P_k # Interior dofs: moments of normal-tangential components against a basis for P_{k-1} - if degree > 0: - cur = len(nodes) - Q = create_quadrature(ref_el, 2*degree-1) - P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola") - phis = P.tabulate(Q.get_points())[(0,) * sd] - for f in sorted(top[sd-1]): - n = ref_el.compute_scaled_normal(f) - tangents = ref_el.compute_tangents(sd-1, f) - nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi) - for phi in phis for t in tangents) - entity_ids[sd][0].extend(range(cur, len(nodes))) + for dim in (sd-1, sd): + q = degree + sd-1 - dim + if q < 0: + continue + + ref_facet = ref_el.construct_subelement(dim) + Q_ref = create_quadrature(ref_facet, degree + q) + P = polynomial_set.ONPolynomialSet(ref_facet, q, scale=1) + phis = P.tabulate(Q_ref.get_points())[(0,) * dim] + + for entity in sorted(top[dim]): + cur = len(nodes) + Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) + Jdet = Q.jacobian_determinant() + for f in ref_el.get_connectivity()[(dim, sd-1)][entity]: + normal = ref_el.compute_scaled_normal(f) + tangents = ref_el.compute_tangents(sd-1, f) + n = normal / Jdet + nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi) + for phi in phis for t in tangents) + entity_ids[dim][entity].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -59,8 +55,21 @@ class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement): the weak symmetry constraint. """ - def __init__(self, ref_el, degree): - poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree) + def __init__(self, ref_el, degree, variant=None): + + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) + variant, interpolant_deg = check_format_variant(variant, degree) + assert variant == "integral" + + if ref_el.is_macrocell(): + base_element = GopalakrishnanLedererSchoberlSecondKind(ref_el.get_parent(), degree) + poly_set = macro.MacroPolynomialSet(ref_el, base_element) + else: + poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree) + dual = GLSDual(ref_el, degree) sd = ref_el.get_spatial_dimension() formdegree = (1, sd-1) diff --git a/FIAT/macro.py b/FIAT/macro.py index 741e69d9c..6327b6724 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -581,7 +581,7 @@ def pullback(mapping, Jinv, phi): F = F1 if k == 1 else F2 phi = numpy.tensordot(F, phi, (1, i)) - perm = [-2, *range(phi.ndim-2), -1] + perm = [-2, *reversed(range(phi.ndim-2)), -1] phi = phi.transpose(perm) return phi @@ -605,18 +605,15 @@ def __init__(self, ref_el, element): base_expansion_set = element.get_nodal_basis().get_expansion_set() expansion_set = base_expansion_set.reconstruct(ref_el=ref_el) - base_entity_ids = expansions.polynomial_entity_ids(base_ref_el, n, base_entity_ids) - entity_ids = expansions.polynomial_entity_ids(ref_el, n, base_entity_ids) - shp = element.value_shape() - num_bfs = expansions.polynomial_dimension(ref_el, n, entity_ids) + num_bfs = expansions.polynomial_dimension(ref_el, n, base_entity_ids) num_members = expansion_set.get_num_members(n) coeffs = numpy.zeros((num_bfs, *shp, num_members)) base_coeffs = element.get_coeffs() rmap = expansions.polynomial_cell_node_map(ref_el, n, base_entity_ids) cmap = expansion_set.get_cell_node_map(n) - for cell in sorted(top): + for cell in sorted(top[sd]): cell_verts = ref_el.get_vertices_of_subcomplex(top[sd][cell]) A, b = reference_element.make_affine_mapping(cell_verts, base_ref_el.vertices) @@ -624,53 +621,3 @@ def __init__(self, ref_el, element): coeffs[indices] = pullback(mapping, A, base_coeffs) super().__init__(ref_el, element.degree(), element.degree(), expansion_set, coeffs) - - -if __name__ == "__main__": - from FIAT import Lagrange, Nedelec, RaviartThomas, Regge - from FIAT.reference_element import symmetric_simplex - - degree = 1 - dim = 3 - K = symmetric_simplex(dim) - A = AlfeldSplit(K) - - fe = Lagrange(K, degree) - fe = Nedelec(K, degree) - fe = RaviartThomas(K, degree) - fe = Regge(K, 0) - - mapping = fe.mapping()[0] - fdim = fe.formdegree - if isinstance(fdim, tuple): - fdim, = set(fdim) - comps = [] - pts = [] - for entity in A.topology[fdim]: - pts_cur = A.make_points(fdim, entity, degree+fdim) - pts.extend(pts_cur) - if mapping == "affine": - comp = numpy.ones(()) - elif mapping.endswith("covariant piola"): - comp = A.compute_edge_tangent(entity) - elif mapping.endswith("contravariant piola"): - comp = A.compute_scaled_normal(entity) - if mapping.startswith("double"): - comp = numpy.outer(comp, comp) - comps.extend([comp] * len(pts_cur)) - - P = MacroPolynomialSet(A, fe) - phis = P.tabulate(pts) - - for alpha in phis: - shape = phis[alpha].shape - shp = shape[1:-1] - result = numpy.zeros((shape[0], shape[-1])) - - ax = (tuple(range(-len(shp), 0)), )*2 - for i, comp in enumerate(comps): - result[..., i] = numpy.tensordot(phis[alpha][..., i], comp, ax) - result[abs(result) < 1E-14] = 0 - phis[alpha] = result - - print(phis) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 6cf122772..6c8257baf 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -6,10 +6,10 @@ # SPDX-License-Identifier: LGPL-3.0-or-later from FIAT import (polynomial_set, expansions, dual_set, - finite_element, functional) -from itertools import chain + finite_element, functional, macro) import numpy -from FIAT.check_format_variant import check_format_variant +from itertools import chain +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule @@ -158,18 +158,23 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): entity_ids[2][i] = list(range(cur, len(nodes))) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T) - dim = sd - phi_deg = degree - dim + phi_deg = degree - sd if phi_deg >= 0: if interpolant_deg is None: interpolant_deg = degree - cur = len(nodes) - Q = create_quadrature(ref_el, interpolant_deg + phi_deg) - Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg) - Phis = Pqmd.tabulate(Q.get_points())[(0,) * dim] - nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (dim,)) - for d in range(dim) for phi in Phis) - entity_ids[dim][0] = list(range(cur, len(nodes))) + + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, interpolant_deg + phi_deg) + Pqmd = polynomial_set.ONPolynomialSet(cell, phi_deg) + Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * sd] + + for entity in top[sd]: + Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) + cur = len(nodes) + nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) + for d in range(sd) + for phi in Phis) + entity_ids[sd][entity] = list(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -195,9 +200,16 @@ class Nedelec(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) variant, interpolant_deg = check_format_variant(variant, degree) - if ref_el.get_spatial_dimension() == 3: + + if ref_el.is_macrocell(): + base_element = Nedelec(ref_el.get_parent(), degree) + poly_set = macro.MacroPolynomialSet(ref_el, base_element) + elif ref_el.get_spatial_dimension() == 3: poly_set = NedelecSpace3D(ref_el, degree) elif ref_el.get_spatial_dimension() == 2: poly_set = NedelecSpace2D(ref_el, degree) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 481dcf64e..622c4ef34 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -15,7 +15,8 @@ from FIAT.raviart_thomas import RaviartThomas from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -from FIAT.check_format_variant import check_format_variant +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant +from FIAT import macro class NedelecSecondKindDual(DualSet): @@ -64,8 +65,8 @@ def generate_degrees_of_freedom(self, cell, degree, variant, interpolant_deg): d = cell.get_spatial_dimension() assert (d in (2, 3)), "Second kind Nedelecs only implemented in 2/3D." - # Zero vertex-based degrees of freedom (d+1 of these) - ids[0] = {i: [] for i in range(d + 1)} + # Zero vertex-based degrees of freedom + ids[0] = {i: [] for i in sorted(cell.topology[0])} # (degree+1) degrees of freedom per entity of codimension 1 (edges) (edge_dofs, ids[1]) = self._generate_edge_dofs(cell, degree, 0, variant, interpolant_deg) @@ -109,16 +110,16 @@ def _generate_edge_dofs(self, cell, degree, offset, variant, interpolant_deg): return (dofs, ids) - def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant_deg): + def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_deg): """Generate degrees of freedom (dofs) for facets.""" # Initialize empty dofs and identifiers (ids) - num_facets = len(cell.get_topology()[codim]) + num_facets = len(cell.get_topology()[dim]) dofs = [] ids = {i: [] for i in range(num_facets)} # Return empty info if not applicable - rt_degree = degree - codim + 1 + rt_degree = degree - dim + 1 if rt_degree < 1: return (dofs, ids) @@ -126,17 +127,17 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant interpolant_deg = degree # Construct quadrature scheme for the reference facet - ref_facet = cell.construct_subelement(codim) + ref_facet = cell.construct_subelement(dim) Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree) - if codim == 1: - Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,)) + if dim == 1: + Phi = ONPolynomialSet(ref_facet, rt_degree, (dim,)) else: # Construct Raviart-Thomas on the reference facet RT = RaviartThomas(ref_facet, rt_degree, variant) Phi = RT.get_nodal_basis() # Evaluate basis functions at reference quadrature points - Phis = Phi.tabulate(Q_ref.get_points())[(0,) * codim] + Phis = Phi.tabulate(Q_ref.get_points())[(0,) * dim] # Note: Phis has dimensions: # num_basis_functions x num_components x num_quad_points Phis = numpy.transpose(Phis, (0, 2, 1)) @@ -147,7 +148,7 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant cur = offset for facet in range(num_facets): # Get the quadrature and Jacobian on this facet - Q_facet = FacetQuadratureRule(cell, codim, facet, Q_ref) + Q_facet = FacetQuadratureRule(cell, dim, facet, Q_ref) J = Q_facet.jacobian() detJ = Q_facet.jacobian_determinant() @@ -191,21 +192,28 @@ class NedelecSecondKind(CiarletElement): interpolation. """ - def __init__(self, cell, degree, variant=None): + def __init__(self, ref_el, degree, variant=None): + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) variant, interpolant_deg = check_format_variant(variant, degree) # Check degree assert degree >= 1, "Second kind Nedelecs start at 1!" # Get dimension - d = cell.get_spatial_dimension() - + d = ref_el.get_spatial_dimension() # Construct polynomial basis for d-vector fields - Ps = ONPolynomialSet(cell, degree, (d, )) + if ref_el.is_macrocell(): + base_element = NedelecSecondKind(ref_el.get_parent(), degree) + Ps = macro.MacroPolynomialSet(ref_el, base_element) + else: + Ps = ONPolynomialSet(ref_el, degree, (d, )) # Construct dual space - Ls = NedelecSecondKindDual(cell, degree, variant, interpolant_deg) + Ls = NedelecSecondKindDual(ref_el, degree, variant, interpolant_deg) # Set form degree formdegree = 1 # 1-form diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 2c7d72c83..329a06489 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -6,10 +6,10 @@ # SPDX-License-Identifier: LGPL-3.0-or-later from FIAT import (expansions, polynomial_set, dual_set, - finite_element, functional) + finite_element, functional, macro) import numpy from itertools import chain -from FIAT.check_format_variant import check_format_variant +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule @@ -71,7 +71,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): entity_ids[dim][entity] = [] if variant == "integral": - facet = ref_el.get_facet_element() + facet = ref_el.construct_subelement(sd-1) # Facet nodes are \int_F v\cdot n p ds where p \in P_q q = degree - 1 Q_ref = create_quadrature(facet, interpolant_deg + q) @@ -89,14 +89,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d if q > 0: - cur = len(nodes) - Q = create_quadrature(ref_el, interpolant_deg + q - 1) - Pqm1 = polynomial_set.ONPolynomialSet(ref_el, q - 1) - Pqm1_at_qpts = Pqm1.tabulate(Q.get_points())[(0,) * sd] - nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) - for phi in Pqm1_at_qpts) - entity_ids[sd][0] = list(range(cur, len(nodes))) + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, interpolant_deg + q - 1) + Pqm1 = polynomial_set.ONPolynomialSet(cell, q - 1) + Pqm1_at_qpts = Pqm1.tabulate(Q_ref.get_points())[(0,) * sd] + + for entity in top[sd]: + Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) + cur = len(nodes) + nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) + for d in range(sd) + for phi in Pqm1_at_qpts) + entity_ids[sd][entity] = list(range(cur, len(nodes))) elif variant == "point": # codimension 1 facets @@ -140,10 +144,17 @@ class RaviartThomas(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) variant, interpolant_deg = check_format_variant(variant, degree) - poly_set = RTSpace(ref_el, degree) + if ref_el.is_macrocell(): + base_element = RaviartThomas(ref_el.get_parent(), degree) + poly_set = macro.MacroPolynomialSet(ref_el, base_element) + else: + poly_set = RTSpace(ref_el, degree) dual = RTDualSet(ref_el, degree, variant, interpolant_deg) formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/regge.py b/FIAT/regge.py index 31dd6bb44..5c4e758c3 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -8,8 +8,8 @@ # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT import dual_set, finite_element, polynomial_set -from FIAT.check_format_variant import check_format_variant +from FIAT import dual_set, finite_element, polynomial_set, macro +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant from FIAT.functional import (PointwiseInnerProductEvaluation, TensorBidirectionalIntegralMoment as BidirectionalMoment) from FIAT.quadrature import FacetQuadratureRule @@ -65,8 +65,17 @@ def __init__(self, ref_el, degree=0, variant=None): if degree < 0: raise ValueError(f"{type(self).__name__} only defined for degree >= 0") + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) variant, qdegree = check_format_variant(variant, degree) - poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + if ref_el.is_macrocell(): + base_element = Regge(ref_el.get_parent(), degree) + poly_set = macro.MacroPolynomialSet(ref_el, base_element) + else: + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + dual = ReggeDual(ref_el, degree, variant, qdegree) formdegree = (1, 1) mapping = "double covariant piola" diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py index 8e367bf20..fde89b7e2 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -1,8 +1,10 @@ import math import numpy import pytest -from FIAT import DiscontinuousLagrange, Lagrange, Legendre, P0 -from FIAT.macro import AlfeldSplit, IsoSplit, PowellSabinSplit, CkPolynomialSet +from FIAT import (Lagrange, Nedelec, RaviartThomas, DiscontinuousLagrange, Legendre, P0, + NedelecSecondKind, BrezziDouglasMarini, + Regge, GopalakrishnanLedererSchoberlSecondKind) +from FIAT.macro import AlfeldSplit, IsoSplit, PowellSabinSplit, CkPolynomialSet, MacroPolynomialSet from FIAT.quadrature_schemes import create_quadrature from FIAT.reference_element import ufc_simplex from FIAT.expansions import polynomial_entity_ids, polynomial_cell_node_map @@ -438,3 +440,56 @@ def test_macro_sympy(cell, element): results = phis(*numpy.transpose(pts)) tab_numpy = fe.tabulate(0, pts)[(0,) * dim] assert numpy.allclose(results, tab_numpy) + + +@pytest.mark.parametrize("element,degree", [ + (Lagrange, 1), (Nedelec, 1), (RaviartThomas, 1), (DiscontinuousLagrange, 0), + (Regge, 0)]) +@pytest.mark.parametrize("dim", (2, 3)) +def test_macro_polynomial_set(dim, element, degree): + K = ufc_simplex(dim) + A = IsoSplit(K) + + fe = element(K, degree) + mapping = fe.mapping()[0] + fdim = fe.formdegree + if isinstance(fdim, tuple): + fdim, = set(fdim) + comps = [] + pts = [] + for entity in A.topology[fdim]: + pts_cur = A.make_points(fdim, entity, 1+fdim) + pts.extend(pts_cur) + if mapping == "affine": + comp = numpy.ones(()) + elif mapping.endswith("covariant piola"): + comp = A.compute_edge_tangent(entity) + elif mapping.endswith("contravariant piola"): + comp = A.compute_scaled_normal(entity) + if mapping.startswith("double"): + comp = numpy.outer(comp, comp) + comps.extend([comp] * len(pts_cur)) + + P = MacroPolynomialSet(A, fe) + phis = P.tabulate(pts) + + alpha = (0,) * dim + shape = phis[alpha].shape + shp = shape[1:-1] + result = numpy.zeros((shape[0], shape[-1])) + + ax = (tuple(range(-len(shp), 0)), )*2 + for i, comp in enumerate(comps): + result[..., i] = numpy.tensordot(phis[alpha][..., i], comp, ax) + result[abs(result) < 1E-14] = 0 + + assert numpy.allclose(numpy.diag(numpy.diag(result)), result) + + +@pytest.mark.parametrize("element,degree", [ + (Lagrange, 4), (Nedelec, 3), (RaviartThomas, 2), (DiscontinuousLagrange, 1), + (NedelecSecondKind, 3), (BrezziDouglasMarini, 2), (Regge, 1), (GopalakrishnanLedererSchoberlSecondKind, 1)]) +@pytest.mark.parametrize("dim", (2, 3)) +@pytest.mark.parametrize("variant", ("alfeld", "iso")) +def test_macro_variants(dim, element, degree, variant): + element(ufc_simplex(dim), degree, variant=variant) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 290329ba7..b56e6969b 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -166,7 +166,7 @@ def test_affine(ref_to_phys, element, degree, dimension): check_zany_mapping(element, ref_to_phys[dimension], degree) -@pytest.mark.parametrize("element", [finat.BrezziDouglasMarini]) +@pytest.mark.parametrize("element", [finat.BrezziDouglasMarini, finat.NedelecSecondKind]) @pytest.mark.parametrize("degree", [1, 2]) @pytest.mark.parametrize("dimension", [2, 3]) @pytest.mark.parametrize("variant", [None, "iso"]) From 9f75c4c7cde3948ced02a8684ecbd3c499a44c7f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 12 Jul 2025 22:32:18 +0100 Subject: [PATCH 08/17] Macro HHJ --- FIAT/gopalakrishnan_lederer_schoberl.py | 14 +++++--- FIAT/hellan_herrmann_johnson.py | 44 +++++++++++++++++-------- finat/fiat_elements.py | 8 ++--- test/FIAT/unit/test_macro.py | 35 +++++++++++++------- test/finat/test_zany_mapping.py | 5 +-- 5 files changed, 70 insertions(+), 36 deletions(-) diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py index d25637728..15625d4de 100644 --- a/FIAT/gopalakrishnan_lederer_schoberl.py +++ b/FIAT/gopalakrishnan_lederer_schoberl.py @@ -77,7 +77,7 @@ def __init__(self, ref_el, degree, variant=None): super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) -def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree): +def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree, variant=None): """The GLS element used for the Mass-Conserving mixed Stress (MCS) formulation for Stokes flow. @@ -86,12 +86,16 @@ def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree): Reference: https://doi.org/10.1093/imanum/drz022 """ - fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree) + fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree, variant=variant) entity_dofs = fe.entity_dofs() sd = ref_el.get_spatial_dimension() - dimPkm1 = (sd-1)*expansions.polynomial_dimension(ref_el.construct_subelement(sd-1), degree-1) + facet = ref_el.construct_subelement(sd-1) + dimPkm1 = (sd-1)*expansions.polynomial_dimension(facet, degree-1) + indices = [] - for f in entity_dofs[sd-1]: + for f in sorted(entity_dofs[sd-1]): indices.extend(entity_dofs[sd-1][f][:dimPkm1]) - indices.extend(entity_dofs[sd][0]) + for cell in sorted(entity_dofs[sd]): + indices.extend(entity_dofs[sd][cell]) + return RestrictedElement(fe, indices=indices) diff --git a/FIAT/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py index d90a9c19f..8e916ddf5 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/FIAT/hellan_herrmann_johnson.py @@ -8,8 +8,8 @@ # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT import dual_set, finite_element, polynomial_set -from FIAT.check_format_variant import check_format_variant +from FIAT import dual_set, finite_element, polynomial_set, macro +from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant from FIAT.functional import (PointwiseInnerProductEvaluation, ComponentPointEvaluation, TensorBidirectionalIntegralMoment as BidirectionalMoment) @@ -66,20 +66,27 @@ def __init__(self, ref_el, degree, variant, qdegree): nodes.extend(PointwiseInnerProductEvaluation(ref_el, n[i+1], n[i+2], pt) for pt in pts for i in range((sd-1)*(sd-2))) else: - Q = create_quadrature(ref_el, qdegree + degree) - P = polynomial_set.ONPolynomialSet(ref_el, degree) - phis = P.tabulate(Q.get_points())[(0,)*sd] - phis /= ref_el.volume() + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, qdegree + degree) + P = polynomial_set.ONPolynomialSet(cell, degree) + Phis = P.tabulate(Q_ref.get_points())[(0,) * sd] dimPkm1 = P.expansion_set.get_num_members(degree-1) - # n[f]^T u n[f] integrated against a basis for P_{k-1} - nodes.extend(BidirectionalMoment(ref_el, n[f], n[f], Q, phi) - for phi in phis[:dimPkm1] for f in sorted(top[sd-1])) - # n[i+1]^T u n[i+2] integrated against a basis for Pk - nodes.extend(BidirectionalMoment(ref_el, n[i+1], n[i+2], Q, phi) - for phi in phis for i in range((sd-1)*(sd-2))) + for entity in sorted(top[sd]): + cur = len(nodes) + faces = ref_el.get_connectivity()[(sd, sd-1)][entity] + Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) + phis = Phis / Q.jacobian_determinant() + + # n[f]^T u n[f] integrated against a basis for P_{k-1} + nodes.extend(BidirectionalMoment(ref_el, n[f], n[f], Q, phi) + for phi in phis[:dimPkm1] for f in faces) + + # n[i+1]^T u n[i+2] integrated against a basis for Pk + nodes.extend(BidirectionalMoment(ref_el, n[faces[i+1]], n[faces[i+2]], Q, phi) + for phi in phis for i in range((sd-1)*(sd-2))) - entity_ids[sd][0].extend(range(cur, len(nodes))) + entity_ids[sd][entity].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -93,8 +100,17 @@ def __init__(self, ref_el, degree=0, variant=None): if degree < 0: raise ValueError(f"{type(self).__name__} only defined for degree >= 0") + if variant is not None: + splitting, variant = parse_lagrange_variant(variant, integral=True) + if splitting is not None: + ref_el = splitting(ref_el) variant, qdegree = check_format_variant(variant, degree) - poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + if ref_el.is_macrocell(): + base_element = HellanHerrmannJohnson(ref_el.get_parent(), degree) + poly_set = macro.MacroPolynomialSet(ref_el, base_element) + else: + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + dual = HellanHerrmannJohnsonDual(ref_el, degree, variant, qdegree) sd = ref_el.get_spatial_dimension() formdegree = (sd-1, sd-1) diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 3677ef965..a5d61b3b5 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -326,13 +326,13 @@ def __init__(self, cell, degree, variant=None): class GopalakrishnanLedererSchoberlFirstKind(FiatElement): # traceless matrix valued - def __init__(self, cell, degree): - super().__init__(FIAT.GopalakrishnanLedererSchoberlFirstKind(cell, degree)) + def __init__(self, cell, degree, variant=None): + super().__init__(FIAT.GopalakrishnanLedererSchoberlFirstKind(cell, degree, variant=variant)) class GopalakrishnanLedererSchoberlSecondKind(FiatElement): # traceless matrix valued - def __init__(self, cell, degree): - super().__init__(FIAT.GopalakrishnanLedererSchoberlSecondKind(cell, degree)) + def __init__(self, cell, degree, variant=None): + super().__init__(FIAT.GopalakrishnanLedererSchoberlSecondKind(cell, degree, variant=variant)) class ScalarFiatElement(FiatElement): diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py index fde89b7e2..b3c49033a 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -3,7 +3,7 @@ import pytest from FIAT import (Lagrange, Nedelec, RaviartThomas, DiscontinuousLagrange, Legendre, P0, NedelecSecondKind, BrezziDouglasMarini, - Regge, GopalakrishnanLedererSchoberlSecondKind) + Regge, HellanHerrmannJohnson, GopalakrishnanLedererSchoberlSecondKind) from FIAT.macro import AlfeldSplit, IsoSplit, PowellSabinSplit, CkPolynomialSet, MacroPolynomialSet from FIAT.quadrature_schemes import create_quadrature from FIAT.reference_element import ufc_simplex @@ -444,7 +444,8 @@ def test_macro_sympy(cell, element): @pytest.mark.parametrize("element,degree", [ (Lagrange, 1), (Nedelec, 1), (RaviartThomas, 1), (DiscontinuousLagrange, 0), - (Regge, 0)]) + (Regge, 0), (HellanHerrmannJohnson, 0), + (GopalakrishnanLedererSchoberlSecondKind, 0)]) @pytest.mark.parametrize("dim", (2, 3)) def test_macro_polynomial_set(dim, element, degree): K = ufc_simplex(dim) @@ -454,41 +455,53 @@ def test_macro_polynomial_set(dim, element, degree): mapping = fe.mapping()[0] fdim = fe.formdegree if isinstance(fdim, tuple): - fdim, = set(fdim) + fdim = max(fdim) comps = [] pts = [] for entity in A.topology[fdim]: pts_cur = A.make_points(fdim, entity, 1+fdim) pts.extend(pts_cur) + if mapping == "affine": comp = numpy.ones(()) + elif mapping == "covariant contravariant piola": + ts = A.compute_tangents(fdim, entity) + n = A.compute_scaled_normal(entity) + comp = ts[..., None] * n[None, None, :] + elif mapping.endswith("covariant piola"): comp = A.compute_edge_tangent(entity) elif mapping.endswith("contravariant piola"): comp = A.compute_scaled_normal(entity) if mapping.startswith("double"): comp = numpy.outer(comp, comp) - comps.extend([comp] * len(pts_cur)) - P = MacroPolynomialSet(A, fe) - phis = P.tabulate(pts) + comps.extend(comp for pt in pts_cur) - alpha = (0,) * dim - shape = phis[alpha].shape + P = MacroPolynomialSet(A, fe) + phis = P.tabulate(pts)[(0,)*dim] + shape = phis.shape shp = shape[1:-1] - result = numpy.zeros((shape[0], shape[-1])) + ncomp = comp.shape[:-len(shp)] + result = numpy.zeros((shape[0], *ncomp, shape[-1])) ax = (tuple(range(-len(shp), 0)), )*2 for i, comp in enumerate(comps): - result[..., i] = numpy.tensordot(phis[alpha][..., i], comp, ax) + result[..., i] = numpy.tensordot(comp, phis[..., i], ax).T result[abs(result) < 1E-14] = 0 + result = result[:len(pts)*numpy.prod(ncomp, dtype=int)] + if ncomp: + result = result.transpose((*range(1, len(ncomp)+1), 0, -1)) + result = result.reshape((shape[0], -1)) assert numpy.allclose(numpy.diag(numpy.diag(result)), result) @pytest.mark.parametrize("element,degree", [ (Lagrange, 4), (Nedelec, 3), (RaviartThomas, 2), (DiscontinuousLagrange, 1), - (NedelecSecondKind, 3), (BrezziDouglasMarini, 2), (Regge, 1), (GopalakrishnanLedererSchoberlSecondKind, 1)]) + (NedelecSecondKind, 3), (BrezziDouglasMarini, 2), + (Regge, 3), (HellanHerrmannJohnson, 1), + (GopalakrishnanLedererSchoberlSecondKind, 1)]) @pytest.mark.parametrize("dim", (2, 3)) @pytest.mark.parametrize("variant", ("alfeld", "iso")) def test_macro_variants(dim, element, degree, variant): diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index b56e6969b..e31074552 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -162,8 +162,9 @@ def test_piola_triangle_high_order(ref_to_phys, element, degree, variant): *((finat.GopalakrishnanLedererSchoberlSecondKind, k) for k in range(0, 3)), ]) @pytest.mark.parametrize("dimension", [2, 3]) -def test_affine(ref_to_phys, element, degree, dimension): - check_zany_mapping(element, ref_to_phys[dimension], degree) +@pytest.mark.parametrize("variant", [None, "alfeld"]) +def test_affine(ref_to_phys, element, degree, variant, dimension): + check_zany_mapping(element, ref_to_phys[dimension], degree, variant=variant) @pytest.mark.parametrize("element", [finat.BrezziDouglasMarini, finat.NedelecSecondKind]) From 312c447dc500cd30413b83535fec89c4f4f72564 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 15 Jul 2025 14:57:26 -0400 Subject: [PATCH 09/17] Macroize all the non-zany elements --- FIAT/argyris.py | 29 ++++++++++------- FIAT/brezzi_douglas_marini.py | 18 ++++------- FIAT/check_format_variant.py | 4 ++- FIAT/crouzeix_raviart.py | 11 ++++--- FIAT/finite_element.py | 6 ++++ FIAT/gopalakrishnan_lederer_schoberl.py | 18 ++++------- FIAT/hellan_herrmann_johnson.py | 18 ++++------- FIAT/hu_zhang.py | 43 ++++++++++++++++--------- FIAT/nedelec.py | 10 +++--- FIAT/nedelec_second_kind.py | 18 +++-------- FIAT/raviart_thomas.py | 10 +++--- FIAT/reference_element.py | 2 +- FIAT/regge.py | 18 ++++------- test/FIAT/unit/test_macro.py | 10 +++++- 14 files changed, 107 insertions(+), 108 deletions(-) diff --git a/FIAT/argyris.py b/FIAT/argyris.py index d98dbbe3b..e3eaeb349 100644 --- a/FIAT/argyris.py +++ b/FIAT/argyris.py @@ -53,13 +53,17 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # interior dofs q = degree - 6 if q >= 0: - Q = create_quadrature(ref_el, interpolant_deg + q) - Pq = polynomial_set.ONPolynomialSet(ref_el, q, scale=1) - phis = Pq.tabulate(Q.get_points())[(0,) * sd] - scale = ref_el.volume() - cur = len(nodes) - nodes.extend(IntegralMoment(ref_el, Q, phi/scale) for phi in phis) - entity_ids[sd][0] = list(range(cur, len(nodes))) + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, interpolant_deg + q) + Pq = polynomial_set.ONPolynomialSet(cell, q, scale=1) + Phis = Pq.tabulate(Q_ref.get_points())[(0,) * sd] + for entity in sorted(top[sd]): + Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) + scale = 1 / Q.jacobian_determinant() + phis = scale * Phis + cur = len(nodes) + nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis) + entity_ids[sd][entity] = list(range(cur, len(nodes))) elif variant == "point": # edge dofs @@ -77,9 +81,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # interior dofs if degree > 5: cur = len(nodes) - internalpts = ref_el.make_points(2, 0, degree - 3) - nodes.extend(PointEvaluation(ref_el, pt) for pt in internalpts) - entity_ids[2][0] = list(range(cur, len(nodes))) + for entity in sorted(top[sd]): + internalpts = ref_el.make_points(sd, entity, degree - 3) + nodes.extend(PointEvaluation(ref_el, pt) for pt in internalpts) + entity_ids[sd][entity] = list(range(cur, len(nodes))) else: raise ValueError("Invalid variant for Argyris") super().__init__(nodes, ref_el, entity_ids) @@ -103,7 +108,9 @@ class Argyris(finite_element.CiarletElement): def __init__(self, ref_el, degree=5, variant=None): - variant, interpolant_deg = check_format_variant(variant, degree) + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + if splitting is not None: + raise NotImplementedError(f"{type(self).__name__} is not implemented as a macroelement.") poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble") dual = ArgyrisDualSet(ref_el, degree, variant, interpolant_deg) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index c9c944e75..13a2687e6 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -6,8 +6,8 @@ # SPDX-License-Identifier: LGPL-3.0-or-later from FIAT import (finite_element, functional, dual_set, - polynomial_set, nedelec, macro) -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant + polynomial_set, nedelec) +from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule import numpy @@ -95,21 +95,15 @@ class BrezziDouglasMarini(finite_element.CiarletElement): def __init__(self, ref_el, degree, variant=None): - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, interpolant_deg = check_format_variant(variant, degree) + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) if degree < 1: raise Exception("BDM_k elements only valid for k >= 1") sd = ref_el.get_spatial_dimension() - if ref_el.is_macrocell(): - base_element = BrezziDouglasMarini(ref_el.get_parent(), degree) - poly_set = macro.MacroPolynomialSet(ref_el, base_element) - else: - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) dual = BDMDualSet(ref_el, degree, variant, interpolant_deg) formdegree = sd - 1 # (n-1)-form super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/check_format_variant.py b/FIAT/check_format_variant.py index 9a50f3262..65f1892c4 100644 --- a/FIAT/check_format_variant.py +++ b/FIAT/check_format_variant.py @@ -27,6 +27,8 @@ def check_format_variant(variant, degree): + splitting, variant = parse_lagrange_variant(variant, integral=True) + if variant is None: variant = "integral" @@ -44,7 +46,7 @@ def check_format_variant(variant, degree): raise ValueError('Choose either variant="point" or variant="integral"' 'or variant="integral(q)"') - return variant, interpolant_degree + return splitting, variant, interpolant_degree def parse_lagrange_variant(variant, discontinuous=False, integral=False): diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index 8447a1269..7b834f40d 100644 --- a/FIAT/crouzeix_raviart.py +++ b/FIAT/crouzeix_raviart.py @@ -81,12 +81,13 @@ class CrouzeixRaviart(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - - variant, interpolant_deg = check_format_variant(variant, degree) - if degree % 2 != 1: raise ValueError("Crouzeix-Raviart only defined for odd degree") - space = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble") + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) + + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg) - super().__init__(space, dual, degree) + super().__init__(poly_set, dual, degree) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 96cf26286..a93554680 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -15,6 +15,7 @@ from FIAT.dual_set import DualSet from FIAT.polynomial_set import PolynomialSet from FIAT.quadrature_schemes import create_quadrature +from FIAT.macro import MacroPolynomialSet class FiniteElement(object): @@ -134,6 +135,11 @@ def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref ref_complex = ref_complex or poly_set.get_reference_element() super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex) + # Tile the poly_set + if ref_complex.is_macrocell() and len(poly_set) > len(dual): + base_element = type(self)(ref_el, order) + poly_set = MacroPolynomialSet(ref_complex, base_element) + if len(poly_set) != len(dual): raise ValueError(f"Dimension of function space is {len(poly_set)}, but got {len(dual)} nodes.") diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py index 15625d4de..6fe608864 100644 --- a/FIAT/gopalakrishnan_lederer_schoberl.py +++ b/FIAT/gopalakrishnan_lederer_schoberl.py @@ -1,5 +1,5 @@ -from FIAT import finite_element, dual_set, polynomial_set, expansions, macro -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant +from FIAT import finite_element, dual_set, polynomial_set, expansions +from FIAT.check_format_variant import check_format_variant from FIAT.functional import TensorBidirectionalIntegralMoment as BidirectionalMoment from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule @@ -57,19 +57,13 @@ class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, interpolant_deg = check_format_variant(variant, degree) + splitting, variant, interpolant_deg = check_format_variant(variant, degree) assert variant == "integral" - if ref_el.is_macrocell(): - base_element = GopalakrishnanLedererSchoberlSecondKind(ref_el.get_parent(), degree) - poly_set = macro.MacroPolynomialSet(ref_el, base_element) - else: - poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree) + if splitting is not None: + ref_el = splitting(ref_el) + poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree) dual = GLSDual(ref_el, degree) sd = ref_el.get_spatial_dimension() formdegree = (1, sd-1) diff --git a/FIAT/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py index 8e916ddf5..52697f0a5 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/FIAT/hellan_herrmann_johnson.py @@ -8,8 +8,8 @@ # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT import dual_set, finite_element, polynomial_set, macro -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant +from FIAT import dual_set, finite_element, polynomial_set +from FIAT.check_format_variant import check_format_variant from FIAT.functional import (PointwiseInnerProductEvaluation, ComponentPointEvaluation, TensorBidirectionalIntegralMoment as BidirectionalMoment) @@ -100,17 +100,11 @@ def __init__(self, ref_el, degree=0, variant=None): if degree < 0: raise ValueError(f"{type(self).__name__} only defined for degree >= 0") - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, qdegree = check_format_variant(variant, degree) - if ref_el.is_macrocell(): - base_element = HellanHerrmannJohnson(ref_el.get_parent(), degree) - poly_set = macro.MacroPolynomialSet(ref_el, base_element) - else: - poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + splitting, variant, qdegree = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) dual = HellanHerrmannJohnsonDual(ref_el, degree, variant, qdegree) sd = ref_el.get_spatial_dimension() formdegree = (sd-1, sd-1) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 5068a6b39..074c457b7 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -11,6 +11,7 @@ from FIAT import finite_element, polynomial_set, dual_set from FIAT.check_format_variant import check_format_variant from FIAT.reference_element import TRIANGLE +from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import (ComponentPointEvaluation, PointwiseInnerProductEvaluation, @@ -54,23 +55,30 @@ def __init__(self, ref_el, degree, variant, qdegree): entity_ids[1][entity].extend(range(cur, len(nodes))) # interior dofs - cur = len(nodes) - if variant == "point": - # unique components evaluated at interior points - pts = ref_el.make_points(sd, 0, degree+1) - nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt) - for pt in pts for i in range(sd) for j in range(i, sd)) + if variant == "integral": + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, 2*degree-2) + P = polynomial_set.ONPolynomialSet(cell, degree-2, scale=1) + Phis = P.tabulate(Q_ref.get_points())[(0,)*sd] + + for entity in sorted(top[sd]): + cur = len(nodes) + if variant == "point": + # unique components evaluated at interior points + pts = ref_el.make_points(sd, entity, degree+1) + nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt) + for pt in pts for i in range(sd) for j in range(i, sd)) - elif variant == "integral": - # Moments of unique components against a basis for P_{k-2} - n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) - Q = create_quadrature(ref_el, 2*degree-2) - P = polynomial_set.ONPolynomialSet(ref_el, degree-2, scale="L2 piola") - phis = P.tabulate(Q.get_points())[(0,)*sd] - nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi) - for phi in phis for i in range(sd) for j in range(i, sd)) + elif variant == "integral": + # Moments of unique components against a basis for P_{k-2} + faces = ref_el.get_connectivity()[(sd, sd-1)][entity] + n = list(map(ref_el.compute_scaled_normal, faces)) + Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) + phis = Phis / Q.jacobian_determinant() + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi) + for phi in phis for i in range(sd) for j in range(i, sd)) - entity_ids[2][0].extend(range(cur, len(nodes))) + entity_ids[sd][entity].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -81,7 +89,10 @@ def __init__(self, ref_el, degree=3, variant=None): raise ValueError(f"{type(self).__name__} only defined for degree >= 3") if ref_el.shape != TRIANGLE: raise ValueError(f"{type(self).__name__} only defined on triangles") - variant, qdegree = check_format_variant(variant, degree) + splitting, variant, qdegree = check_format_variant(variant, degree) + if splitting is not None: + raise NotImplementedError(f"{type(self).__name__} is not implemented as a macroelement.") + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) dual = HuZhangDual(ref_el, degree, variant, qdegree) formdegree = ref_el.get_spatial_dimension() - 1 diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 6c8257baf..8208f7918 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -9,7 +9,7 @@ finite_element, functional, macro) import numpy from itertools import chain -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant +from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule @@ -200,11 +200,9 @@ class Nedelec(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, interpolant_deg = check_format_variant(variant, degree) + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) if ref_el.is_macrocell(): base_element = Nedelec(ref_el.get_parent(), degree) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 622c4ef34..dc1670574 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -15,8 +15,7 @@ from FIAT.raviart_thomas import RaviartThomas from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant -from FIAT import macro +from FIAT.check_format_variant import check_format_variant class NedelecSecondKindDual(DualSet): @@ -193,12 +192,9 @@ class NedelecSecondKind(CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, interpolant_deg = check_format_variant(variant, degree) + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) # Check degree assert degree >= 1, "Second kind Nedelecs start at 1!" @@ -206,11 +202,7 @@ def __init__(self, ref_el, degree, variant=None): # Get dimension d = ref_el.get_spatial_dimension() # Construct polynomial basis for d-vector fields - if ref_el.is_macrocell(): - base_element = NedelecSecondKind(ref_el.get_parent(), degree) - Ps = macro.MacroPolynomialSet(ref_el, base_element) - else: - Ps = ONPolynomialSet(ref_el, degree, (d, )) + Ps = ONPolynomialSet(ref_el, degree, (d, )) # Construct dual space Ls = NedelecSecondKindDual(ref_el, degree, variant, interpolant_deg) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 329a06489..220833609 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -9,7 +9,7 @@ finite_element, functional, macro) import numpy from itertools import chain -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant +from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature from FIAT.quadrature import FacetQuadratureRule @@ -144,11 +144,9 @@ class RaviartThomas(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, interpolant_deg = check_format_variant(variant, degree) + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) if ref_el.is_macrocell(): base_element = RaviartThomas(ref_el.get_parent(), degree) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 3baca43e9..72dc140ba 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -623,7 +623,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): indices = slice(None) subcomplex = top[entity_dim][entity_id] if entity_dim != sd: - cell_id = self.connectivity[(entity_dim, sd)][0][0] + cell_id = self.connectivity[(entity_dim, sd)][entity_id][0] indices = [i for i, v in enumerate(top[sd][cell_id]) if v in subcomplex] subcomplex = top[sd][cell_id] diff --git a/FIAT/regge.py b/FIAT/regge.py index 5c4e758c3..2f0f4bc7d 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -8,8 +8,8 @@ # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT import dual_set, finite_element, polynomial_set, macro -from FIAT.check_format_variant import check_format_variant, parse_lagrange_variant +from FIAT import dual_set, finite_element, polynomial_set +from FIAT.check_format_variant import check_format_variant from FIAT.functional import (PointwiseInnerProductEvaluation, TensorBidirectionalIntegralMoment as BidirectionalMoment) from FIAT.quadrature import FacetQuadratureRule @@ -65,17 +65,11 @@ def __init__(self, ref_el, degree=0, variant=None): if degree < 0: raise ValueError(f"{type(self).__name__} only defined for degree >= 0") - if variant is not None: - splitting, variant = parse_lagrange_variant(variant, integral=True) - if splitting is not None: - ref_el = splitting(ref_el) - variant, qdegree = check_format_variant(variant, degree) - if ref_el.is_macrocell(): - base_element = Regge(ref_el.get_parent(), degree) - poly_set = macro.MacroPolynomialSet(ref_el, base_element) - else: - poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + splitting, variant, qdegree = check_format_variant(variant, degree) + if splitting is not None: + ref_el = splitting(ref_el) + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) dual = ReggeDual(ref_el, degree, variant, qdegree) formdegree = (1, 1) mapping = "double covariant piola" diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py index b3c49033a..a92eeae62 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -3,7 +3,8 @@ import pytest from FIAT import (Lagrange, Nedelec, RaviartThomas, DiscontinuousLagrange, Legendre, P0, NedelecSecondKind, BrezziDouglasMarini, - Regge, HellanHerrmannJohnson, GopalakrishnanLedererSchoberlSecondKind) + Regge, HellanHerrmannJohnson, GopalakrishnanLedererSchoberlSecondKind, + CrouzeixRaviart) from FIAT.macro import AlfeldSplit, IsoSplit, PowellSabinSplit, CkPolynomialSet, MacroPolynomialSet from FIAT.quadrature_schemes import create_quadrature from FIAT.reference_element import ufc_simplex @@ -506,3 +507,10 @@ def test_macro_polynomial_set(dim, element, degree): @pytest.mark.parametrize("variant", ("alfeld", "iso")) def test_macro_variants(dim, element, degree, variant): element(ufc_simplex(dim), degree, variant=variant) + + +@pytest.mark.parametrize("element,degree", [ + (CrouzeixRaviart, 3)]) +@pytest.mark.parametrize("variant", ("alfeld", "iso")) +def test_macro_variants_triangle(element, degree, variant): + element(ufc_simplex(2), degree, variant=variant) From 3d15fee77af83e2532e9c8da1621c9799d1a7269 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 17 Jul 2025 10:47:11 -0400 Subject: [PATCH 10/17] FunctionalBlock --- FIAT/brezzi_douglas_marini.py | 15 +--- FIAT/functional.py | 124 +++++++++++++++++++++++++++++++++- FIAT/nedelec.py | 11 +-- FIAT/nedelec_second_kind.py | 36 +++------- FIAT/raviart_thomas.py | 7 +- 5 files changed, 138 insertions(+), 55 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 13a2687e6..e16b93a25 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -9,8 +9,6 @@ polynomial_set, nedelec) from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule -import numpy class BDMDualSet(dual_set.DualSet): @@ -35,12 +33,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in top[sd - 1]: cur = len(nodes) - Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet - phis = n[None, :, None] * Pq_at_qpts[:, None, :] - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in phis) + nodes.extend(functional.FacetNormalIntegralMomentBlock(ref_el, f, Q_ref, Pq_at_qpts)) entity_ids[sd - 1][f] = list(range(cur, len(nodes))) elif variant == "point": @@ -61,13 +54,11 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cell = ref_el.construct_subelement(sd) Q_ref = create_quadrature(cell, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(cell, degree - 1, variant) + mapping, = set(Nedel.mapping()) Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd] for entity in top[sd]: - Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - Jinv = numpy.linalg.inv(Q.jacobian()) - phis = numpy.tensordot(Jinv.T, Ned_at_qpts, (1, 1)).transpose((1, 0, 2)) cur = len(nodes) - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) + nodes.extend(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Ned_at_qpts, mapping=mapping)) entity_ids[sd][entity] = list(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/functional.py b/FIAT/functional.py index 18de29942..f3636931d 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -15,7 +15,7 @@ from itertools import chain import numpy -from FIAT import polynomial_set, jacobi, quadrature_schemes +from FIAT import polynomial_set, jacobi, quadrature, quadrature_schemes def index_iterator(shp): @@ -25,7 +25,127 @@ def index_iterator(shp): return numpy.ndindex(shp) -class Functional(object): +def pullback(mapping, J, detJ, phi): + try: + formdegree = { + "identity": (0,), + "L2 piola": (3,), + "covariant piola": (1,), + "contravariant piola": (2,), + "double covariant piola": (1, 1), + "double contravariant piola": (2, 2), + "covariant contravariant piola": (1, 2), + "contravariant covariant piola": (2, 1), }[mapping] + except KeyError: + raise ValueError(f"Unrecognized mapping {mapping}") + + F1 = numpy.linalg.pinv(J).T + F2 = J / detJ + + perm = [*range(1, phi.ndim-1), 0, -1] + phi = phi.transpose(perm) + + for i, k in enumerate(formdegree): + if k == 0: + continue + elif k == 3: + phi = phi / detJ + else: + F = F1 if k == 1 else F2 + phi = numpy.tensordot(F, phi, (1, i)) + + perm = [-2, *reversed(range(phi.ndim-2)), -1] + phi = phi.transpose(perm) + return phi + + +class FunctionalBlock: + def __init__(self, nodes): + self.nodes = nodes + + def __iter__(self): + for ell in self.nodes: + yield ell + + def __len__(self): + return len(self.nodes) + + def completion(self): + return self + + +class FunctionalBlockView(FunctionalBlock): + def __init__(self, block, index): + self.block = block + nodes = [block.nodes[index]] + super().__init__(nodes) + + def completion(self): + return self.block + + +class PointFunctionalBlock(FunctionalBlock): + def __init__(self, ref_el, x, order): + from FIAT.polynomial_set import mis + sd = ref_el.get_spatial_dimension() + nodes = [PointDerivative(ref_el, x, alpha) for alpha in mis(sd, order)] + super().__init__(nodes) + + +class PointGradient(PointFunctionalBlock): + def __init__(self, ref_el, x): + super().__init__(ref_el, x, 1) + + +class PointHessian(PointFunctionalBlock): + def __init__(self, ref_el, x): + super().__init__(ref_el, x, 2) + + +class PointDirectionalDerivativeBlock(PointFunctionalBlock): + def __init__(self, ref_el, x, basis): + nodes = [PointDirectionalDerivative(ref_el, x, e) for e in basis] + super().__init__(nodes) + + +class PointNormalTangentialDerivativeBlock(PointDirectionalDerivativeBlock): + def __init__(self, ref_el, x, facet_id): + sd = ref_el.get_spatial_dimension() + basis = numpy.zeros((sd, sd)) + basis[0] = ref_el.compute_scaled_normal(facet_id) + basis[1:] = ref_el.compute_tangents(sd-1, facet_id) + super().__init__(self, x, basis) + + +class PointNormalDerivativeView(FunctionalBlockView): + def __init__(self, ref_el, x, facet_id): + block = PointNormalTangentialDerivativeBlock(ref_el, x, facet_id) + super().__init__(block, 0) + + +class FacetIntegralMomentBlock(FunctionalBlock): + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis, mapping="L2 piola"): + Q = quadrature.FacetQuadratureRule(ref_el, entity_dim, entity_id, Q_ref) + phis = pullback(mapping, Q.jacobian(), Q.jacobian_determinant(), Phis) + nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis] + super().__init__(nodes) + + +class FacetDirectionalIntegralMomentBlock(FacetIntegralMomentBlock): + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis, direction): + direction = numpy.asarray(direction) + Phis = Phis[(slice(None), *(None for _ in range(direction.ndim)), slice(None))] * direction[(*(None for _ in range(Phis.ndim-1)), Ellipsis, None)] + super().__init__(ref_el, entity_dim, entity_id, Q_ref, Phis) + + +class FacetNormalIntegralMomentBlock(FacetDirectionalIntegralMomentBlock): + def __init__(self, ref_el, facet_id, Q_ref, Phis): + direction = ref_el.compute_scaled_normal(facet_id) + sd = ref_el.get_spatial_dimension() + super().__init__(ref_el, sd-1, facet_id, Q_ref, Phis, direction) + + +class Functional(FunctionalBlock): r"""Abstract class representing a linear functional. All FIAT functionals are discrete in the sense that they are written as a weighted sum of (derivatives of components of) their diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 8208f7918..771283b5f 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -125,17 +125,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim] - Phis = numpy.transpose(Phis, (0, 2, 1)) - + mapping = "contravariant piola" for entity in top[dim]: cur = len(nodes) - Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) - Jdet = Q.jacobian_determinant() - R = numpy.array(ref_el.compute_tangents(dim, entity)) - phis = numpy.dot(Phis, R / Jdet) - phis = numpy.transpose(phis, (0, 2, 1)) - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in phis) + nodes.extend(functional.FacetIntegralMomentBlock(ref_el, dim, entity, Q_ref, Phis, mapping=mapping)) entity_ids[dim][entity] = list(range(cur, len(nodes))) elif variant == "point": diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index dc1670574..6e970b4e4 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -5,15 +5,12 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later -import numpy - from FIAT.finite_element import CiarletElement from FIAT.dual_set import DualSet from FIAT.polynomial_set import ONPolynomialSet from FIAT.functional import PointEdgeTangentEvaluation as Tangent -from FIAT.functional import FrobeniusIntegralMoment as IntegralMoment +from FIAT.functional import FacetIntegralMomentBlock from FIAT.raviart_thomas import RaviartThomas -from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature from FIAT.check_format_variant import check_format_variant @@ -128,42 +125,29 @@ def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_d # Construct quadrature scheme for the reference facet ref_facet = cell.construct_subelement(dim) Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree) + + # Construct Raviart-Thomas on the reference facet if dim == 1: - Phi = ONPolynomialSet(ref_facet, rt_degree, (dim,)) + Phi = ONPolynomialSet(ref_facet, rt_degree, shape=(dim,)) + mapping = "contravariant piola" else: - # Construct Raviart-Thomas on the reference facet RT = RaviartThomas(ref_facet, rt_degree, variant) Phi = RT.get_nodal_basis() + mapping, = set(RT.mapping()) # Evaluate basis functions at reference quadrature points Phis = Phi.tabulate(Q_ref.get_points())[(0,) * dim] - # Note: Phis has dimensions: - # num_basis_functions x num_components x num_quad_points - Phis = numpy.transpose(Phis, (0, 2, 1)) - # Note: Phis has dimensions: - # num_basis_functions x num_quad_points x num_components # Iterate over the facets - cur = offset - for facet in range(num_facets): - # Get the quadrature and Jacobian on this facet - Q_facet = FacetQuadratureRule(cell, dim, facet, Q_ref) - J = Q_facet.jacobian() - detJ = Q_facet.jacobian_determinant() - - # Map Phis -> phis (reference values to physical values) - piola_map = J / detJ - phis = numpy.dot(Phis, piola_map.T) - phis = numpy.transpose(phis, (0, 2, 1)) - + for entity in range(num_facets): + cur = offset + len(dofs) # Construct degrees of freedom as integral moments on this cell, # using the face quadrature weighted against the values # of the (physical) Raviart--Thomas'es on the face - dofs.extend(IntegralMoment(cell, Q_facet, phi) for phi in phis) + dofs.extend(FacetIntegralMomentBlock(cell, dim, entity, Q_ref, Phis, mapping=mapping)) # Assign identifiers (num RTs per face + previous edge dofs) - ids[facet].extend(range(cur, cur + len(phis))) - cur += len(phis) + ids[entity].extend(range(cur, offset + len(dofs))) return (dofs, ids) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 220833609..e9aa9738e 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -79,12 +79,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in top[sd - 1]: cur = len(nodes) - Q = FacetQuadratureRule(ref_el, sd-1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet - phis = n[None, :, None] * Pq_at_qpts[:, None, :] - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in phis) + nodes.extend(functional.FacetNormalIntegralMomentBlock(ref_el, f, Q_ref, Pq_at_qpts)) entity_ids[sd - 1][f] = list(range(cur, len(nodes))) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d From c477bd02d5b2fffe5f785a461f94b84d76eb1cf5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 17 Jul 2025 10:47:11 -0400 Subject: [PATCH 11/17] FunctionalBlock: generate nodes on multiple entities --- FIAT/brezzi_douglas_marini.py | 17 ++-- FIAT/functional.py | 159 +++++++++++++++++++++++++++++++++- FIAT/nedelec.py | 12 +-- FIAT/nedelec_second_kind.py | 37 +++----- FIAT/raviart_thomas.py | 8 +- 5 files changed, 178 insertions(+), 55 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 13a2687e6..8be412c98 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -9,8 +9,6 @@ polynomial_set, nedelec) from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule -import numpy class BDMDualSet(dual_set.DualSet): @@ -33,14 +31,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(facet, interpolant_deg + degree) Pq = polynomial_set.ONPolynomialSet(facet, degree) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] + ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq_at_qpts) for f in top[sd - 1]: cur = len(nodes) - Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet - phis = n[None, :, None] * Pq_at_qpts[:, None, :] - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in phis) + nodes.extend(ells.nodes(sd-1, f)) entity_ids[sd - 1][f] = list(range(cur, len(nodes))) elif variant == "point": @@ -61,13 +55,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cell = ref_el.construct_subelement(sd) Q_ref = create_quadrature(cell, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(cell, degree - 1, variant) + mapping, = set(Nedel.mapping()) Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd] + ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Ned_at_qpts, mapping=mapping) for entity in top[sd]: - Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - Jinv = numpy.linalg.inv(Q.jacobian()) - phis = numpy.tensordot(Jinv.T, Ned_at_qpts, (1, 1)).transpose((1, 0, 2)) cur = len(nodes) - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) + nodes.extend(ells.nodes(sd, entity)) entity_ids[sd][entity] = list(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/functional.py b/FIAT/functional.py index 18de29942..01174f7c5 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -15,7 +15,7 @@ from itertools import chain import numpy -from FIAT import polynomial_set, jacobi, quadrature_schemes +from FIAT import polynomial_set, jacobi, quadrature, quadrature_schemes def index_iterator(shp): @@ -25,7 +25,162 @@ def index_iterator(shp): return numpy.ndindex(shp) -class Functional(object): +def pullback(mapping, J, detJ, phi): + try: + formdegree = { + "identity": (0,), + "L2 piola": (3,), + "covariant piola": (1,), + "contravariant piola": (2,), + "double covariant piola": (1, 1), + "double contravariant piola": (2, 2), + "covariant contravariant piola": (1, 2), + "contravariant covariant piola": (2, 1), }[mapping] + except KeyError: + raise ValueError(f"Unrecognized mapping {mapping}") + + F1 = numpy.linalg.pinv(J).T + F2 = J / detJ + + perm = [*range(1, phi.ndim-1), 0, -1] + phi = phi.transpose(perm) + + for i, k in enumerate(formdegree): + if k == 0: + continue + elif k == 3: + phi = phi / detJ + else: + F = F1 if k == 1 else F2 + phi = numpy.tensordot(F, phi, (1, i)) + + perm = [-2, *reversed(range(phi.ndim-2)), -1] + phi = phi.transpose(perm) + return phi + + +class FunctionalBlock: + def __init__(self, ref_el): + self.ref_el = ref_el + + def completion(self): + return self + + def nodes(self, entity_dim, entity_id): + raise NotImplementedError + + +class FunctionalBlockView(FunctionalBlock): + def __init__(self, block, subset): + self.subset = subset + self.block = block + super().__init__(block.ref_el) + + def completion(self): + return self.block + + def nodes(self, entity_dim, entity_id): + nodes = list(self.block.nodes(entity_dim, entity_id)) + for i in self.subset: + yield nodes[i] + + +class PointFunctionalBlock(FunctionalBlock): + def __init__(self, ref_el, order): + self.order = order + super().__init__(ref_el) + + def nodes(self, entity_dim, entity_id): + from FIAT.polynomial_set import mis + sd = self.ref_el.get_spatial_dimension() + for pt in self.ref_el.make_points(entity_dim, entity_id, 0): + for alpha in mis(sd, self.order): + yield PointDerivative(self.ref_el, pt, alpha) + + +class PointGradient(PointFunctionalBlock): + def __init__(self, ref_el, x): + super().__init__(ref_el, 1) + + +class PointHessian(PointFunctionalBlock): + def __init__(self, ref_el, x): + super().__init__(ref_el, 2) + + +class PointDirectionalDerivativeBlock(PointFunctionalBlock): + def __init__(self, ref_el, basis): + self.basis = basis + super().__init__(ref_el) + + def get_basis(self, entity_dim, entity_id): + return self.basis + + def nodes(self, entity_dim, entity_id): + for pt in self.ref_el.make_points(entity_dim, entity_id, 0): + for e in self.get_basis(entity_dim, entity_id): + yield PointDirectionalDerivative(self.ref_el, pt, e) + + +class PointNormalTangentialDerivativeBlock(PointDirectionalDerivativeBlock): + def __init__(self, ref_el, facet_id): + super().__init__(ref_el, None) + + def get_basis(self, entity_dim, entity_id): + sd = self.ref_el.get_spatial_dimension() + basis = numpy.zeros((sd, sd)) + basis[0] = self.ref_el.compute_scaled_normal(entity_id) + basis[1:] = self.ref_el.compute_tangents(sd-1, entity_id) + return basis + + +class PointNormalDerivativeView(FunctionalBlockView): + def __init__(self, ref_el): + block = PointNormalTangentialDerivativeBlock(ref_el) + super().__init__(block, [0]) + + +class FacetIntegralMomentBlock(FunctionalBlock): + def __init__(self, ref_el, Q_ref, Phis, mapping="L2 piola"): + self.Q_ref = Q_ref + self.Phis = Phis + self.mapping = mapping + super().__init__(ref_el) + + def get_functions(self, entity_dim, entity_id): + return self.Phis + + def nodes(self, entity_dim, entity_id): + Q = quadrature.FacetQuadratureRule(self.ref_el, entity_dim, entity_id, self.Q_ref) + Phis = self.get_functions(entity_dim, entity_id) + phis = pullback(self.mapping, Q.jacobian(), Q.jacobian_determinant(), Phis) + for phi in phis: + yield FrobeniusIntegralMoment(self.ref_el, Q, phi) + + +class FacetDirectionalIntegralMomentBlock(FacetIntegralMomentBlock): + def __init__(self, ref_el, Q_ref, Phis, direction): + self.direction = direction + super().__init__(ref_el, Q_ref, Phis) + + def get_direction(self, entity_dim, entity_id): + return self.direction + + def get_functions(self, entity_dim, entity_id): + udir = self.get_direction(entity_dim, entity_id) + Phis = self.Phis + return Phis[(slice(None), *(None for _ in range(udir.ndim)), slice(None))] * udir[(*(None for _ in range(Phis.ndim-1)), Ellipsis, None)] + + +class FacetNormalIntegralMomentBlock(FacetDirectionalIntegralMomentBlock): + def __init__(self, ref_el, Q_ref, Phis): + super().__init__(ref_el, Q_ref, Phis, "normal") + + def get_direction(self, entity_dim, entity_id): + return self.ref_el.compute_scaled_normal(entity_id) + + +class Functional(FunctionalBlock): r"""Abstract class representing a linear functional. All FIAT functionals are discrete in the sense that they are written as a weighted sum of (derivatives of components of) their diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 8208f7918..1310bb760 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -125,17 +125,11 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim] - Phis = numpy.transpose(Phis, (0, 2, 1)) - + mapping = "contravariant piola" + ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Phis, mapping=mapping) for entity in top[dim]: cur = len(nodes) - Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref) - Jdet = Q.jacobian_determinant() - R = numpy.array(ref_el.compute_tangents(dim, entity)) - phis = numpy.dot(Phis, R / Jdet) - phis = numpy.transpose(phis, (0, 2, 1)) - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in phis) + nodes.extend(ells.nodes(dim, entity)) entity_ids[dim][entity] = list(range(cur, len(nodes))) elif variant == "point": diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index dc1670574..d3257da57 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -5,15 +5,12 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later -import numpy - from FIAT.finite_element import CiarletElement from FIAT.dual_set import DualSet from FIAT.polynomial_set import ONPolynomialSet from FIAT.functional import PointEdgeTangentEvaluation as Tangent -from FIAT.functional import FrobeniusIntegralMoment as IntegralMoment +from FIAT.functional import FacetIntegralMomentBlock from FIAT.raviart_thomas import RaviartThomas -from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature from FIAT.check_format_variant import check_format_variant @@ -128,42 +125,30 @@ def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_d # Construct quadrature scheme for the reference facet ref_facet = cell.construct_subelement(dim) Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree) + + # Construct Raviart-Thomas on the reference facet if dim == 1: - Phi = ONPolynomialSet(ref_facet, rt_degree, (dim,)) + Phi = ONPolynomialSet(ref_facet, rt_degree, shape=(dim,)) + mapping = "contravariant piola" else: - # Construct Raviart-Thomas on the reference facet RT = RaviartThomas(ref_facet, rt_degree, variant) Phi = RT.get_nodal_basis() + mapping, = set(RT.mapping()) # Evaluate basis functions at reference quadrature points Phis = Phi.tabulate(Q_ref.get_points())[(0,) * dim] - # Note: Phis has dimensions: - # num_basis_functions x num_components x num_quad_points - Phis = numpy.transpose(Phis, (0, 2, 1)) - # Note: Phis has dimensions: - # num_basis_functions x num_quad_points x num_components + ells = FacetIntegralMomentBlock(cell, Q_ref, Phis, mapping=mapping) # Iterate over the facets - cur = offset - for facet in range(num_facets): - # Get the quadrature and Jacobian on this facet - Q_facet = FacetQuadratureRule(cell, dim, facet, Q_ref) - J = Q_facet.jacobian() - detJ = Q_facet.jacobian_determinant() - - # Map Phis -> phis (reference values to physical values) - piola_map = J / detJ - phis = numpy.dot(Phis, piola_map.T) - phis = numpy.transpose(phis, (0, 2, 1)) - + for entity in range(num_facets): + cur = offset + len(dofs) # Construct degrees of freedom as integral moments on this cell, # using the face quadrature weighted against the values # of the (physical) Raviart--Thomas'es on the face - dofs.extend(IntegralMoment(cell, Q_facet, phi) for phi in phis) + dofs.extend(ells.nodes(dim, entity)) # Assign identifiers (num RTs per face + previous edge dofs) - ids[facet].extend(range(cur, cur + len(phis))) - cur += len(phis) + ids[entity].extend(range(cur, offset + len(dofs))) return (dofs, ids) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 220833609..4216bc701 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -77,14 +77,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(facet, interpolant_deg + q) Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] + ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq_at_qpts) for f in top[sd - 1]: cur = len(nodes) - Q = FacetQuadratureRule(ref_el, sd-1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet - phis = n[None, :, None] * Pq_at_qpts[:, None, :] - nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) - for phi in phis) + nodes.extend(ells.nodes(sd-1, f)) entity_ids[sd - 1][f] = list(range(cur, len(nodes))) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d From ac856a0434a2e6cf6fee9440e20b17f90b53c257 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 17 Jul 2025 14:24:22 -0400 Subject: [PATCH 12/17] Allow nodes to be passed as a dict --- FIAT/brezzi_douglas_marini.py | 30 +++++----------- FIAT/dual_set.py | 17 ++++++++- FIAT/functional.py | 8 +++-- FIAT/nedelec.py | 36 ++++++------------- FIAT/nedelec_second_kind.py | 66 ++++++++++++++--------------------- FIAT/raviart_thomas.py | 35 ++++++------------- 6 files changed, 77 insertions(+), 115 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 8be412c98..c3b59b901 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -13,16 +13,11 @@ class BDMDualSet(dual_set.DualSet): def __init__(self, ref_el, degree, variant, interpolant_deg): - nodes = [] sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - entity_ids = {} # set to empty - for dim in top: - entity_ids[dim] = {} - for entity in top[dim]: - entity_ids[dim][entity] = [] + nodes = {dim: {entity: [] for entity in top[dim]} for dim in top} if variant == "integral": facet = ref_el.construct_subelement(sd-1) @@ -30,22 +25,17 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) Pq = polynomial_set.ONPolynomialSet(facet, degree) - Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] - ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq_at_qpts) + ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq) for f in top[sd - 1]: - cur = len(nodes) - nodes.extend(ells.nodes(sd-1, f)) - entity_ids[sd - 1][f] = list(range(cur, len(nodes))) + nodes[sd-1][f] = [ells] elif variant == "point": # Define each functional for the dual set # codimension 1 facets for f in top[sd - 1]: - cur = len(nodes) pts_cur = ref_el.make_points(sd - 1, f, sd + degree) - nodes.extend(functional.PointScaledNormalEvaluation(ref_el, f, pt) - for pt in pts_cur) - entity_ids[sd - 1][f] = list(range(cur, len(nodes))) + nodes[sd-1][f].extend(functional.PointScaledNormalEvaluation(ref_el, f, pt) + for pt in pts_cur) # internal nodes if degree > 1: @@ -56,14 +46,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(cell, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(cell, degree - 1, variant) mapping, = set(Nedel.mapping()) - Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd] - ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Ned_at_qpts, mapping=mapping) + Phi = Nedel.get_nodal_basis() + ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Phi, mapping=mapping) for entity in top[sd]: - cur = len(nodes) - nodes.extend(ells.nodes(sd, entity)) - entity_ids[sd][entity] = list(range(cur, len(nodes))) + nodes[sd][entity] = [ells] - super().__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el) class BrezziDouglasMarini(finite_element.CiarletElement): diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index c09db974f..033ab77f5 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -12,7 +12,22 @@ class DualSet(object): - def __init__(self, nodes, ref_el, entity_ids, entity_permutations=None): + def __init__(self, nodes, ref_el, entity_ids=None, entity_permutations=None): + if isinstance(nodes, dict): + assert entity_ids is None + # flatten nodes + top = ref_el.get_topology() + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + entity_nodes = nodes + nodes = [] + for dim in entity_nodes: + for entity in entity_nodes[dim]: + cur = len(nodes) + ells = entity_nodes[dim][entity] + for ell in ells: + nodes.extend(ell.nodes(dim, entity)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + nodes, ref_el, entity_ids, entity_permutations = merge_entities(nodes, ref_el, entity_ids, entity_permutations) self.nodes = nodes self.ref_el = ref_el diff --git a/FIAT/functional.py b/FIAT/functional.py index 01174f7c5..cb40c94ca 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -141,9 +141,10 @@ def __init__(self, ref_el): class FacetIntegralMomentBlock(FunctionalBlock): - def __init__(self, ref_el, Q_ref, Phis, mapping="L2 piola"): + def __init__(self, ref_el, Q_ref, P, mapping="L2 piola"): + dim = P.ref_el.get_spatial_dimension() self.Q_ref = Q_ref - self.Phis = Phis + self.Phis = P.tabulate(Q_ref.get_points())[(0,) * dim] self.mapping = mapping super().__init__(ref_el) @@ -216,6 +217,9 @@ def __init__(self, ref_el, target_shape, pt_dict, deriv_dict, else: self.max_deriv_order = 0 + def nodes(self, entity_dim, entity_id): + yield self + def evaluate(self, f): """Obsolete and broken functional evaluation. diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 1310bb760..373269359 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -104,12 +104,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - entity_ids = {} # set to empty - for dim in top: - entity_ids[dim] = {} - for entity in top[dim]: - entity_ids[dim][entity] = [] + nodes = {dim: {entity: [] for entity in top[dim]} for dim in top} if variant == "integral": # edge nodes are \int_F v\cdot t p ds where p \in P_{q-1}(edge) @@ -124,32 +120,23 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): facet = ref_el.construct_subelement(dim) Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) - Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim] mapping = "contravariant piola" - ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Phis, mapping=mapping) + ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Pqmd, mapping=mapping) for entity in top[dim]: - cur = len(nodes) - nodes.extend(ells.nodes(dim, entity)) - entity_ids[dim][entity] = list(range(cur, len(nodes))) + nodes[dim][entity].append(ells) elif variant == "point": for i in top[1]: - cur = len(nodes) # points to specify P_k on each edge pts_cur = ref_el.make_points(1, i, degree + 1) - nodes.extend(functional.PointEdgeTangentEvaluation(ref_el, i, pt) - for pt in pts_cur) - entity_ids[1][i] = list(range(cur, len(nodes))) + nodes[1][i].extend(functional.PointEdgeTangentEvaluation(ref_el, i, pt) + for pt in pts_cur) if sd > 2 and degree > 1: # face tangents for i in top[2]: # loop over faces - cur = len(nodes) pts_cur = ref_el.make_points(2, i, degree + 1) - nodes.extend(functional.PointFaceTangentEvaluation(ref_el, i, k, pt) - for k in range(2) # loop over tangents - for pt in pts_cur # loop over points - ) - entity_ids[2][i] = list(range(cur, len(nodes))) + nodes[2][1].extend(functional.PointFaceTangentEvaluation(ref_el, i, k, pt) + for k in range(2) for pt in pts_cur) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T) phi_deg = degree - sd @@ -164,13 +151,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for entity in top[sd]: Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - cur = len(nodes) - nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) - for phi in Phis) - entity_ids[sd][entity] = list(range(cur, len(nodes))) + nodes[sd][entity].extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) + for d in range(sd) for phi in Phis) - super().__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el) class Nedelec(finite_element.CiarletElement): diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index d3257da57..8fb735e49 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -47,77 +47,68 @@ class NedelecSecondKindDual(DualSet): def __init__(self, cell, degree, variant, interpolant_deg): # Define degrees of freedom - (dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg) + nodes = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg) # Call init of super-class - super().__init__(dofs, cell, ids) + super().__init__(nodes, cell) def generate_degrees_of_freedom(self, cell, degree, variant, interpolant_deg): - "Generate dofs and geometry-to-dof maps (ids)." + """Generate degrees of freedom for each entity.""" - dofs = [] - ids = {} + nodes = {} # Extract spatial dimension and topology d = cell.get_spatial_dimension() assert (d in (2, 3)), "Second kind Nedelecs only implemented in 2/3D." # Zero vertex-based degrees of freedom - ids[0] = {i: [] for i in sorted(cell.topology[0])} + nodes[0] = {i: [] for i in sorted(cell.topology[0])} # (degree+1) degrees of freedom per entity of codimension 1 (edges) - (edge_dofs, ids[1]) = self._generate_edge_dofs(cell, degree, 0, variant, interpolant_deg) - dofs.extend(edge_dofs) + nodes[1] = self._generate_edge_dofs(cell, degree, variant, interpolant_deg) # Include face degrees of freedom if 3D if d == 3: - face_dofs, ids[d-1] = self._generate_facet_dofs(d-1, cell, degree, - len(dofs), variant, interpolant_deg) - dofs.extend(face_dofs) + nodes[d-1] = self._generate_facet_dofs(d-1, cell, degree, + variant, interpolant_deg) # Varying degrees of freedom (possibly zero) per cell - cell_dofs, ids[d] = self._generate_facet_dofs(d, cell, degree, len(dofs), variant, interpolant_deg) - dofs.extend(cell_dofs) + nodes[d] = self._generate_facet_dofs(d, cell, degree, variant, interpolant_deg) - return (dofs, ids) + return nodes - def _generate_edge_dofs(self, cell, degree, offset, variant, interpolant_deg): + def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg): """Generate degrees of freedom (dofs) for entities of codimension 1 (edges).""" if variant == "integral": - return self._generate_facet_dofs(1, cell, degree, offset, variant, interpolant_deg) + return self._generate_facet_dofs(1, cell, degree, variant, interpolant_deg) # (degree+1) tangential component point evaluation degrees of # freedom per entity of codimension 1 (edges) - dofs = [] - ids = {} + top = cell.get_topology() + nodes = {entity: [] for entity in top[1]} if variant == "point": - for edge in range(len(cell.get_topology()[1])): + for edge in top[1]: # Create points for evaluation of tangential components points = cell.make_points(1, edge, degree + 2) # A tangential component evaluation for each point - dofs.extend(Tangent(cell, edge, point) for point in points) + nodes[edge].extend(Tangent(cell, edge, point) for point in points) - # Associate these dofs with this edge - i = len(points) * edge - ids[edge] = list(range(offset + i, offset + i + len(points))) + return nodes - return (dofs, ids) - - def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_deg): + def _generate_facet_dofs(self, dim, cell, degree, variant, interpolant_deg): """Generate degrees of freedom (dofs) for facets.""" - # Initialize empty dofs and identifiers (ids) - num_facets = len(cell.get_topology()[dim]) - dofs = [] - ids = {i: [] for i in range(num_facets)} + # Initialize empty dofs + top = cell.get_topology() + nodes = {entity: [] for entity in top[dim]} # Return empty info if not applicable rt_degree = degree - dim + 1 if rt_degree < 1: - return (dofs, ids) + return nodes if interpolant_deg is None: interpolant_deg = degree @@ -136,21 +127,16 @@ def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_d mapping, = set(RT.mapping()) # Evaluate basis functions at reference quadrature points - Phis = Phi.tabulate(Q_ref.get_points())[(0,) * dim] - ells = FacetIntegralMomentBlock(cell, Q_ref, Phis, mapping=mapping) + ells = FacetIntegralMomentBlock(cell, Q_ref, Phi, mapping=mapping) # Iterate over the facets - for entity in range(num_facets): - cur = offset + len(dofs) + for entity in top[dim]: # Construct degrees of freedom as integral moments on this cell, # using the face quadrature weighted against the values # of the (physical) Raviart--Thomas'es on the face - dofs.extend(ells.nodes(dim, entity)) - - # Assign identifiers (num RTs per face + previous edge dofs) - ids[entity].extend(range(cur, offset + len(dofs))) + nodes[entity].append(ells) - return (dofs, ids) + return nodes class NedelecSecondKind(CiarletElement): diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 4216bc701..25f603d12 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -63,12 +63,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - entity_ids = {} # set to empty - for dim in top: - entity_ids[dim] = {} - for entity in top[dim]: - entity_ids[dim][entity] = [] + nodes = {dim: {entity: [] for entity in top[dim]} for dim in top} if variant == "integral": facet = ref_el.construct_subelement(sd-1) @@ -76,12 +72,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): q = degree - 1 Q_ref = create_quadrature(facet, interpolant_deg + q) Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0) - Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] - ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq_at_qpts) + ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq) for f in top[sd - 1]: - cur = len(nodes) - nodes.extend(ells.nodes(sd-1, f)) - entity_ids[sd - 1][f] = list(range(cur, len(nodes))) + nodes[sd - 1][f].append(ells) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d if q > 0: @@ -92,31 +85,23 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for entity in top[sd]: Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - cur = len(nodes) - nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) - for phi in Pqm1_at_qpts) - entity_ids[sd][entity] = list(range(cur, len(nodes))) + nodes[sd][entity].extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) + for d in range(sd) for phi in Pqm1_at_qpts) elif variant == "point": # codimension 1 facets for i in top[sd - 1]: - cur = len(nodes) pts_cur = ref_el.make_points(sd - 1, i, sd + degree - 1) - nodes.extend(functional.PointScaledNormalEvaluation(ref_el, i, pt) - for pt in pts_cur) - entity_ids[sd - 1][i] = list(range(cur, len(nodes))) + nodes[sd - 1][i].extend(functional.PointScaledNormalEvaluation(ref_el, i, pt) + for pt in pts_cur) # internal nodes. Let's just use points at a lattice if degree > 1: - cur = len(nodes) pts = ref_el.make_points(sd, 0, sd + degree - 1) - nodes.extend(functional.ComponentPointEvaluation(ref_el, d, (sd,), pt) - for d in range(sd) - for pt in pts) - entity_ids[sd][0] = list(range(cur, len(nodes))) + nodes[sd][0].extend(functional.ComponentPointEvaluation(ref_el, d, (sd,), pt) + for d in range(sd) for pt in pts) - super().__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el) class RaviartThomas(finite_element.CiarletElement): From 601f6d48562ef338ec6b3506e449170cb023a6c4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 17 Jul 2025 17:46:42 -0400 Subject: [PATCH 13/17] Attach entity to FunctionalBlock --- FIAT/brezzi_douglas_marini.py | 13 ++-- FIAT/dual_set.py | 18 +++-- FIAT/functional.py | 120 +++++++++++++++------------------- FIAT/hermite.py | 42 ++---------- FIAT/lagrange.py | 15 ++--- FIAT/nedelec.py | 21 +++--- FIAT/nedelec_second_kind.py | 36 ++++------ FIAT/raviart_thomas.py | 21 +++--- 8 files changed, 108 insertions(+), 178 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index c3b59b901..4ff66cdb8 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -17,7 +17,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): top = ref_el.get_topology() # set to empty - nodes = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] if variant == "integral": facet = ref_el.construct_subelement(sd-1) @@ -25,17 +25,17 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) Pq = polynomial_set.ONPolynomialSet(facet, degree) - ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq) for f in top[sd - 1]: - nodes[sd-1][f] = [ells] + nodes.append(functional.FacetNormalIntegralMomentBlock(ref_el, sd-1, f, Q_ref, Pq)) elif variant == "point": # Define each functional for the dual set # codimension 1 facets for f in top[sd - 1]: pts_cur = ref_el.make_points(sd - 1, f, sd + degree) - nodes[sd-1][f].extend(functional.PointScaledNormalEvaluation(ref_el, f, pt) - for pt in pts_cur) + nodes.append(functional.FunctionalBlock(ref_el, sd-1, f, [ + functional.PointScaledNormalEvaluation(ref_el, f, pt) + for pt in pts_cur])) # internal nodes if degree > 1: @@ -47,9 +47,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Nedel = nedelec.Nedelec(cell, degree - 1, variant) mapping, = set(Nedel.mapping()) Phi = Nedel.get_nodal_basis() - ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Phi, mapping=mapping) for entity in top[sd]: - nodes[sd][entity] = [ells] + nodes.append(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Phi, mapping=mapping)) super().__init__(nodes, ref_el) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 033ab77f5..4432bfe40 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -13,20 +13,18 @@ class DualSet(object): def __init__(self, nodes, ref_el, entity_ids=None, entity_permutations=None): - if isinstance(nodes, dict): - assert entity_ids is None + if entity_ids is None: # flatten nodes top = ref_el.get_topology() entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} - entity_nodes = nodes + blocks = nodes nodes = [] - for dim in entity_nodes: - for entity in entity_nodes[dim]: - cur = len(nodes) - ells = entity_nodes[dim][entity] - for ell in ells: - nodes.extend(ell.nodes(dim, entity)) - entity_ids[dim][entity].extend(range(cur, len(nodes))) + for ell in blocks: + cur = len(nodes) + nodes.extend(ell.nodes) + dim = ell.entity_dim + entity = ell.entity_id + entity_ids[dim][entity].extend(range(cur, len(nodes))) nodes, ref_el, entity_ids, entity_permutations = merge_entities(nodes, ref_el, entity_ids, entity_permutations) self.nodes = nodes diff --git a/FIAT/functional.py b/FIAT/functional.py index cb40c94ca..fe5d98071 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -60,125 +60,107 @@ def pullback(mapping, J, detJ, phi): class FunctionalBlock: - def __init__(self, ref_el): + def __init__(self, ref_el, entity_dim, entity_id, nodes): self.ref_el = ref_el + self.entity_dim = entity_dim + self.entity_id = entity_id + self.nodes = nodes def completion(self): return self - def nodes(self, entity_dim, entity_id): - raise NotImplementedError - class FunctionalBlockView(FunctionalBlock): def __init__(self, block, subset): - self.subset = subset self.block = block - super().__init__(block.ref_el) + nodes = [block.nodes[i] for i in subset] + super().__init__(block.ref_el, block.entity_dim, block.entity_id, nodes) - def completion(self): + def completion(self, entity_dim, entity_id): return self.block - def nodes(self, entity_dim, entity_id): - nodes = list(self.block.nodes(entity_dim, entity_id)) - for i in self.subset: - yield nodes[i] +class RelabeledFunctionalBlock(FunctionalBlock): + def __init__(self, block, entity_dim, entity_id): + super().__init__(block.ref_el, entity_dim, entity_id, block.nodes) -class PointFunctionalBlock(FunctionalBlock): - def __init__(self, ref_el, order): - self.order = order - super().__init__(ref_el) - def nodes(self, entity_dim, entity_id): +class PointEvaluationBlock(FunctionalBlock): + def __init__(self, ref_el, entity_dim, entity_id, order=0, degree=0, variant=None): from FIAT.polynomial_set import mis - sd = self.ref_el.get_spatial_dimension() - for pt in self.ref_el.make_points(entity_dim, entity_id, 0): - for alpha in mis(sd, self.order): - yield PointDerivative(self.ref_el, pt, alpha) - - -class PointGradient(PointFunctionalBlock): - def __init__(self, ref_el, x): - super().__init__(ref_el, 1) + sd = ref_el.get_spatial_dimension() + pts = ref_el.make_points(entity_dim, entity_id, degree, variant=variant) + if order == 0: + nodes = [PointEvaluation(ref_el, pt) for pt in pts] + else: + nodes = [PointDerivative(ref_el, pt, alpha) for pt in pts for alpha in mis(sd, order)] + super().__init__(ref_el, entity_dim, entity_id, nodes) -class PointHessian(PointFunctionalBlock): - def __init__(self, ref_el, x): - super().__init__(ref_el, 2) +class PointGradientBlock(PointEvaluationBlock): + def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): + super().__init__(ref_el, entity_dim, entity_id, order=1, degree=degree, variant=variant) -class PointDirectionalDerivativeBlock(PointFunctionalBlock): - def __init__(self, ref_el, basis): - self.basis = basis - super().__init__(ref_el) +class PointHessianBlock(PointEvaluationBlock): + def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): + super().__init__(ref_el, entity_dim, entity_id, order=2, degree=degree, variant=None) - def get_basis(self, entity_dim, entity_id): - return self.basis - def nodes(self, entity_dim, entity_id): - for pt in self.ref_el.make_points(entity_dim, entity_id, 0): - for e in self.get_basis(entity_dim, entity_id): - yield PointDirectionalDerivative(self.ref_el, pt, e) +class PointDirectionalDerivativeBlock(PointEvaluationBlock): + def __init__(self, ref_el, entity_dim, entity_id, basis): + nodes = [PointDirectionalDerivative(self.ref_el, pt, e) + for pt in self.ref_el.make_points(entity_dim, entity_id, 0) + for e in basis] + super().__init__(ref_el, entity_dim, entity_id, nodes) class PointNormalTangentialDerivativeBlock(PointDirectionalDerivativeBlock): - def __init__(self, ref_el, facet_id): - super().__init__(ref_el, None) - - def get_basis(self, entity_dim, entity_id): + def __init__(self, ref_el, entity_dim, entity_id): sd = self.ref_el.get_spatial_dimension() basis = numpy.zeros((sd, sd)) basis[0] = self.ref_el.compute_scaled_normal(entity_id) basis[1:] = self.ref_el.compute_tangents(sd-1, entity_id) - return basis + super().__init__(ref_el, entity_dim, entity_id, basis) class PointNormalDerivativeView(FunctionalBlockView): - def __init__(self, ref_el): - block = PointNormalTangentialDerivativeBlock(ref_el) + def __init__(self, ref_el, entity_dim, entity_id): + block = PointNormalTangentialDerivativeBlock(ref_el, entity_dim, entity_id) super().__init__(block, [0]) class FacetIntegralMomentBlock(FunctionalBlock): - def __init__(self, ref_el, Q_ref, P, mapping="L2 piola"): + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, P, mapping="L2 piola"): dim = P.ref_el.get_spatial_dimension() - self.Q_ref = Q_ref self.Phis = P.tabulate(Q_ref.get_points())[(0,) * dim] - self.mapping = mapping - super().__init__(ref_el) - def get_functions(self, entity_dim, entity_id): - return self.Phis + Q = quadrature.FacetQuadratureRule(ref_el, entity_dim, entity_id, Q_ref) + Phis = self.get_functions() + phis = pullback(mapping, Q.jacobian(), Q.jacobian_determinant(), Phis) + nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis] - def nodes(self, entity_dim, entity_id): - Q = quadrature.FacetQuadratureRule(self.ref_el, entity_dim, entity_id, self.Q_ref) - Phis = self.get_functions(entity_dim, entity_id) - phis = pullback(self.mapping, Q.jacobian(), Q.jacobian_determinant(), Phis) - for phi in phis: - yield FrobeniusIntegralMoment(self.ref_el, Q, phi) + super().__init__(ref_el, entity_dim, entity_id, nodes) + + def get_functions(self): + return self.Phis class FacetDirectionalIntegralMomentBlock(FacetIntegralMomentBlock): - def __init__(self, ref_el, Q_ref, Phis, direction): + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis, direction): self.direction = direction - super().__init__(ref_el, Q_ref, Phis) + super().__init__(ref_el, entity_dim, entity_id, Q_ref, Phis) - def get_direction(self, entity_dim, entity_id): - return self.direction - - def get_functions(self, entity_dim, entity_id): - udir = self.get_direction(entity_dim, entity_id) + def get_functions(self): + udir = self.direction Phis = self.Phis return Phis[(slice(None), *(None for _ in range(udir.ndim)), slice(None))] * udir[(*(None for _ in range(Phis.ndim-1)), Ellipsis, None)] class FacetNormalIntegralMomentBlock(FacetDirectionalIntegralMomentBlock): - def __init__(self, ref_el, Q_ref, Phis): - super().__init__(ref_el, Q_ref, Phis, "normal") - - def get_direction(self, entity_dim, entity_id): - return self.ref_el.compute_scaled_normal(entity_id) + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis): + normal = ref_el.compute_scaled_normal(entity_id) + super().__init__(ref_el, entity_dim, entity_id, Q_ref, Phis, normal) class Functional(FunctionalBlock): diff --git a/FIAT/hermite.py b/FIAT/hermite.py index 79c83ff97..564a352d5 100644 --- a/FIAT/hermite.py +++ b/FIAT/hermite.py @@ -14,57 +14,25 @@ class CubicHermiteDualSet(dual_set.DualSet): equispaced points.""" 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 top = ref_el.get_topology() - verts = ref_el.get_vertices() sd = ref_el.get_spatial_dimension() # get jet at each vertex - - entity_ids[0] = {} - for v in sorted(top[0]): - nodes.append(functional.PointEvaluation(ref_el, verts[v])) - pd = functional.PointDerivative - for i in range(sd): - alpha = [0] * sd - alpha[i] = 1 - - nodes.append(pd(ref_el, verts[v], alpha)) - - entity_ids[0][v] = list(range(cur, cur + 1 + sd)) - cur += sd + 1 + nodes = [functional.PointEvaluationBlock(ref_el, 0, v, order=order) + for v in sorted(top[0]) for order in (0, 1)] # now only have dofs at the barycenter, which is the # maximal dimension # no edge dof - - entity_ids[1] = {} - for i in top[1]: - entity_ids - entity_ids[1][i] = [] - if sd > 1: # face dof # point evaluation at barycenter - entity_ids[2] = {} - for f in sorted(top[2]): - pt = ref_el.make_points(2, f, 3)[0] - n = functional.PointEvaluation(ref_el, pt) - nodes.append(n) - entity_ids[2][f] = list(range(cur, cur + 1)) - cur += 1 - - for dim in range(3, sd + 1): - entity_ids[dim] = {} - for facet in top[dim]: - entity_ids[dim][facet] = [] + nodes.extend(functional.PointEvaluationBlock(ref_el, 2, f, degree=3) + for f in sorted(top[2])) - super().__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el) class CubicHermite(finite_element.CiarletElement): diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index eead0cfa6..b598bee9e 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -26,12 +26,10 @@ class LagrangeDualSet(dual_set.DualSet): and then lexicographically by lattice multiindex. """ def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=False): - nodes = [] - entity_ids = {} - entity_permutations = {} top = ref_el.get_topology() + + entity_permutations = {} for dim in sorted(top): - entity_ids[dim] = {} entity_permutations[dim] = {} perms = {0: [0]} if dim == 0 else make_entity_permutations_simplex(dim, degree - dim) for entity in sorted(top[dim]): @@ -45,12 +43,9 @@ def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=Fal # make nodes by getting points # need to do this entity-by-entity - for dim, entity in entities: - cur = len(nodes) - pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant) - nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts_cur) - entity_ids[dim][entity] = list(range(cur, len(nodes))) - super().__init__(nodes, ref_el, entity_ids, entity_permutations) + nodes = [functional.PointEvaluationBlock(ref_el, dim, entity, degree=degree, variant=point_variant) + for dim, entity in entities] + super().__init__(nodes, ref_el, entity_permutations=entity_permutations) class Lagrange(finite_element.CiarletElement): diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 373269359..4760e21e8 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -104,9 +104,6 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - # set to empty - nodes = {dim: {entity: [] for entity in top[dim]} for dim in top} - if variant == "integral": # edge nodes are \int_F v\cdot t p ds where p \in P_{q-1}(edge) # degree is q - 1 @@ -121,22 +118,23 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) mapping = "contravariant piola" - ells = functional.FacetIntegralMomentBlock(ref_el, Q_ref, Pqmd, mapping=mapping) for entity in top[dim]: - nodes[dim][entity].append(ells) + nodes.append(functional.FacetIntegralMomentBlock(ref_el, dim, entity, Q_ref, Pqmd, mapping=mapping)) elif variant == "point": for i in top[1]: # points to specify P_k on each edge pts_cur = ref_el.make_points(1, i, degree + 1) - nodes[1][i].extend(functional.PointEdgeTangentEvaluation(ref_el, i, pt) - for pt in pts_cur) + nodes.append(functional.FunctionalBlock(ref_el, 1, i, [ + functional.PointEdgeTangentEvaluation(ref_el, i, pt) + for pt in pts_cur])) if sd > 2 and degree > 1: # face tangents for i in top[2]: # loop over faces pts_cur = ref_el.make_points(2, i, degree + 1) - nodes[2][1].extend(functional.PointFaceTangentEvaluation(ref_el, i, k, pt) - for k in range(2) for pt in pts_cur) + nodes.append(functional.FunctionalBlock(ref_el, 2, i, [ + functional.PointFaceTangentEvaluation(ref_el, i, k, pt) + for k in range(2) for pt in pts_cur])) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T) phi_deg = degree - sd @@ -151,8 +149,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for entity in top[sd]: Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - nodes[sd][entity].extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) for phi in Phis) + nodes.append(functional.FunctionalBlock(ref_el, sd, entity, [ + functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) + for d in range(sd) for phi in Phis])) super().__init__(nodes, ref_el) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 8fb735e49..6994b6f74 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -9,7 +9,7 @@ from FIAT.dual_set import DualSet from FIAT.polynomial_set import ONPolynomialSet from FIAT.functional import PointEdgeTangentEvaluation as Tangent -from FIAT.functional import FacetIntegralMomentBlock +from FIAT.functional import FacetIntegralMomentBlock, FunctionalBlock from FIAT.raviart_thomas import RaviartThomas from FIAT.quadrature_schemes import create_quadrature from FIAT.check_format_variant import check_format_variant @@ -43,38 +43,27 @@ class NedelecSecondKindDual(DualSet): Higher spatial dimensions are not yet implemented. (For d = 1, these elements coincide with the CG_k elements.) """ - def __init__(self, cell, degree, variant, interpolant_deg): - - # Define degrees of freedom - nodes = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg) - # Call init of super-class - super().__init__(nodes, cell) - - def generate_degrees_of_freedom(self, cell, degree, variant, interpolant_deg): - """Generate degrees of freedom for each entity.""" - - nodes = {} - # Extract spatial dimension and topology d = cell.get_spatial_dimension() assert (d in (2, 3)), "Second kind Nedelecs only implemented in 2/3D." # Zero vertex-based degrees of freedom - nodes[0] = {i: [] for i in sorted(cell.topology[0])} + nodes = [] # (degree+1) degrees of freedom per entity of codimension 1 (edges) - nodes[1] = self._generate_edge_dofs(cell, degree, variant, interpolant_deg) + nodes.extend(self._generate_edge_dofs(cell, degree, variant, interpolant_deg)) # Include face degrees of freedom if 3D if d == 3: - nodes[d-1] = self._generate_facet_dofs(d-1, cell, degree, - variant, interpolant_deg) + nodes.extend(self._generate_facet_dofs(d-1, cell, degree, + variant, interpolant_deg)) # Varying degrees of freedom (possibly zero) per cell - nodes[d] = self._generate_facet_dofs(d, cell, degree, variant, interpolant_deg) + nodes.extend(self._generate_facet_dofs(d, cell, degree, variant, interpolant_deg)) - return nodes + # Call init of super-class + super().__init__(nodes, cell) def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg): """Generate degrees of freedom (dofs) for entities of @@ -86,7 +75,7 @@ def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg): # (degree+1) tangential component point evaluation degrees of # freedom per entity of codimension 1 (edges) top = cell.get_topology() - nodes = {entity: [] for entity in top[1]} + nodes = [] if variant == "point": for edge in top[1]: @@ -94,7 +83,7 @@ def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg): points = cell.make_points(1, edge, degree + 2) # A tangential component evaluation for each point - nodes[edge].extend(Tangent(cell, edge, point) for point in points) + nodes.append(FunctionalBlock(cell, 1, edge, [Tangent(cell, edge, point) for point in points])) return nodes @@ -103,7 +92,7 @@ def _generate_facet_dofs(self, dim, cell, degree, variant, interpolant_deg): # Initialize empty dofs top = cell.get_topology() - nodes = {entity: [] for entity in top[dim]} + nodes = [] # Return empty info if not applicable rt_degree = degree - dim + 1 @@ -127,14 +116,13 @@ def _generate_facet_dofs(self, dim, cell, degree, variant, interpolant_deg): mapping, = set(RT.mapping()) # Evaluate basis functions at reference quadrature points - ells = FacetIntegralMomentBlock(cell, Q_ref, Phi, mapping=mapping) # Iterate over the facets for entity in top[dim]: # Construct degrees of freedom as integral moments on this cell, # using the face quadrature weighted against the values # of the (physical) Raviart--Thomas'es on the face - nodes[entity].append(ells) + nodes.append(FacetIntegralMomentBlock(cell, dim, entity, Q_ref, Phi, mapping=mapping)) return nodes diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 25f603d12..df9597ffb 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -59,12 +59,11 @@ class RTDualSet(dual_set.DualSet): moments against polynomials""" def __init__(self, ref_el, degree, variant, interpolant_deg): - nodes = [] sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() # set to empty - nodes = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] if variant == "integral": facet = ref_el.construct_subelement(sd-1) @@ -72,9 +71,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): q = degree - 1 Q_ref = create_quadrature(facet, interpolant_deg + q) Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0) - ells = functional.FacetNormalIntegralMomentBlock(ref_el, Q_ref, Pq) for f in top[sd - 1]: - nodes[sd - 1][f].append(ells) + nodes.append(functional.FacetNormalIntegralMomentBlock(ref_el, sd-1, f, Q_ref, Pq)) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d if q > 0: @@ -85,21 +83,24 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for entity in top[sd]: Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - nodes[sd][entity].extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) for phi in Pqm1_at_qpts) + nodes.append(functional.FunctionalBlock(ref_el, sd, entity, [ + functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) + for d in range(sd) for phi in Pqm1_at_qpts])) elif variant == "point": # codimension 1 facets for i in top[sd - 1]: pts_cur = ref_el.make_points(sd - 1, i, sd + degree - 1) - nodes[sd - 1][i].extend(functional.PointScaledNormalEvaluation(ref_el, i, pt) - for pt in pts_cur) + nodes.append(functional.FunctionalBlock(ref_el, sd-1, i, [ + functional.PointScaledNormalEvaluation(ref_el, i, pt) + for pt in pts_cur])) # internal nodes. Let's just use points at a lattice if degree > 1: pts = ref_el.make_points(sd, 0, sd + degree - 1) - nodes[sd][0].extend(functional.ComponentPointEvaluation(ref_el, d, (sd,), pt) - for d in range(sd) for pt in pts) + nodes.append(functional.FunctionalBlock(ref_el, sd, 0, [ + functional.ComponentPointEvaluation(ref_el, d, (sd,), pt) + for d in range(sd) for pt in pts])) super().__init__(nodes, ref_el) From 198df5413542ae20fff09ec6e57ab157b80b085c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 17 Jul 2025 22:31:11 -0400 Subject: [PATCH 14/17] More FunctionalBlocks --- FIAT/functional.py | 32 +++++++++++++++++++++++++++----- FIAT/nedelec.py | 25 +++++-------------------- FIAT/nedelec_second_kind.py | 14 ++++---------- FIAT/raviart_thomas.py | 20 +++++--------------- 4 files changed, 41 insertions(+), 50 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index fe5d98071..b3d0890b1 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -86,17 +86,35 @@ def __init__(self, block, entity_dim, entity_id): class PointEvaluationBlock(FunctionalBlock): - def __init__(self, ref_el, entity_dim, entity_id, order=0, degree=0, variant=None): + def __init__(self, ref_el, entity_dim, entity_id, order=0, degree=0, variant=None, comp=(), shape=None): from FIAT.polynomial_set import mis - sd = ref_el.get_spatial_dimension() pts = ref_el.make_points(entity_dim, entity_id, degree, variant=variant) if order == 0: - nodes = [PointEvaluation(ref_el, pt) for pt in pts] + if shape is None: + nodes = [PointEvaluation(ref_el, pt) for pt in pts] + else: + nodes = [ComponentPointEvaluation(ref_el, comp, shape, pt) for pt in pts] else: + sd = ref_el.get_spatial_dimension() nodes = [PointDerivative(ref_el, pt, alpha) for pt in pts for alpha in mis(sd, order)] super().__init__(ref_el, entity_dim, entity_id, nodes) +class PointDirectionalEvaluationBlock(FunctionalBlock): + def __init__(self, ref_el, entity_dim, entity_id, direction=None, degree=0, variant=None): + pts = ref_el.make_points(entity_dim, entity_id, degree, variant=variant) + if direction == "normal": + nodes = [PointScaledNormalEvaluation(ref_el, entity_id, pt) for pt in pts] + elif direction == "tangential": + if entity_dim == 1: + nodes = [PointEdgeTangentEvaluation(ref_el, entity_id, pt) for pt in pts] + else: + nodes = [PointFaceTangentEvaluation(ref_el, entity_id, i, pt) + for i in range(entity_dim) for pt in pts] + + super().__init__(ref_el, entity_dim, entity_id, nodes) + + class PointGradientBlock(PointEvaluationBlock): def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): super().__init__(ref_el, entity_dim, entity_id, order=1, degree=degree, variant=variant) @@ -131,14 +149,18 @@ def __init__(self, ref_el, entity_dim, entity_id): class FacetIntegralMomentBlock(FunctionalBlock): - def __init__(self, ref_el, entity_dim, entity_id, Q_ref, P, mapping="L2 piola"): + def __init__(self, ref_el, entity_dim, entity_id, Q_ref, P, mapping="L2 piola", comp=(), shape=None): dim = P.ref_el.get_spatial_dimension() self.Phis = P.tabulate(Q_ref.get_points())[(0,) * dim] Q = quadrature.FacetQuadratureRule(ref_el, entity_dim, entity_id, Q_ref) Phis = self.get_functions() phis = pullback(mapping, Q.jacobian(), Q.jacobian_determinant(), Phis) - nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis] + + if shape is not None: + nodes = [IntegralMoment(ref_el, Q, phi, comp=comp, shp=shape) for phi in phis] + else: + nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis] super().__init__(ref_el, entity_dim, entity_id, nodes) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 4760e21e8..ea2aee0ba 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -11,7 +11,6 @@ from itertools import chain from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule def NedelecSpace2D(ref_el, degree): @@ -122,19 +121,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): nodes.append(functional.FacetIntegralMomentBlock(ref_el, dim, entity, Q_ref, Pqmd, mapping=mapping)) elif variant == "point": - for i in top[1]: - # points to specify P_k on each edge - pts_cur = ref_el.make_points(1, i, degree + 1) - nodes.append(functional.FunctionalBlock(ref_el, 1, i, [ - functional.PointEdgeTangentEvaluation(ref_el, i, pt) - for pt in pts_cur])) - - if sd > 2 and degree > 1: # face tangents - for i in top[2]: # loop over faces - pts_cur = ref_el.make_points(2, i, degree + 1) - nodes.append(functional.FunctionalBlock(ref_el, 2, i, [ - functional.PointFaceTangentEvaluation(ref_el, i, k, pt) - for k in range(2) for pt in pts_cur])) + for dim in range(1, sd): + for i in top[2]: + nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, dim, i, direction="tangential", degree=degree+1)) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T) phi_deg = degree - sd @@ -145,13 +134,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cell = ref_el.construct_subelement(sd) Q_ref = create_quadrature(cell, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(cell, phi_deg) - Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * sd] - for entity in top[sd]: - Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - nodes.append(functional.FunctionalBlock(ref_el, sd, entity, [ - functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) for phi in Phis])) + nodes.extend(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Pqmd, comp=(i,), shape=(sd,)) + for i in range(sd)) super().__init__(nodes, ref_el) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 6994b6f74..96d6d27d5 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -8,8 +8,7 @@ from FIAT.finite_element import CiarletElement from FIAT.dual_set import DualSet from FIAT.polynomial_set import ONPolynomialSet -from FIAT.functional import PointEdgeTangentEvaluation as Tangent -from FIAT.functional import FacetIntegralMomentBlock, FunctionalBlock +from FIAT.functional import FacetIntegralMomentBlock, PointDirectionalEvaluationBlock from FIAT.raviart_thomas import RaviartThomas from FIAT.quadrature_schemes import create_quadrature from FIAT.check_format_variant import check_format_variant @@ -76,14 +75,9 @@ def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg): # freedom per entity of codimension 1 (edges) top = cell.get_topology() nodes = [] - if variant == "point": - for edge in top[1]: - - # Create points for evaluation of tangential components - points = cell.make_points(1, edge, degree + 2) - - # A tangential component evaluation for each point - nodes.append(FunctionalBlock(cell, 1, edge, [Tangent(cell, edge, point) for point in points])) + for edge in top[1]: + # A tangential component evaluation for each point + nodes.append(PointDirectionalEvaluationBlock(cell, 1, edge, direction="tangential", degree=degree+2)) return nodes diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index df9597ffb..6581b4dfc 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -11,7 +11,6 @@ from itertools import chain from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule def RTSpace(ref_el, degree): @@ -79,28 +78,19 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cell = ref_el.construct_subelement(sd) Q_ref = create_quadrature(cell, interpolant_deg + q - 1) Pqm1 = polynomial_set.ONPolynomialSet(cell, q - 1) - Pqm1_at_qpts = Pqm1.tabulate(Q_ref.get_points())[(0,) * sd] - for entity in top[sd]: - Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref) - nodes.append(functional.FunctionalBlock(ref_el, sd, entity, [ - functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) - for d in range(sd) for phi in Pqm1_at_qpts])) + nodes.extend(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Pqm1, comp=(i,), shape=(sd,)) + for i in range(sd)) elif variant == "point": # codimension 1 facets for i in top[sd - 1]: - pts_cur = ref_el.make_points(sd - 1, i, sd + degree - 1) - nodes.append(functional.FunctionalBlock(ref_el, sd-1, i, [ - functional.PointScaledNormalEvaluation(ref_el, i, pt) - for pt in pts_cur])) + nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, sd-1, i, degree=degree+sd-1, direction="normal")) # internal nodes. Let's just use points at a lattice if degree > 1: - pts = ref_el.make_points(sd, 0, sd + degree - 1) - nodes.append(functional.FunctionalBlock(ref_el, sd, 0, [ - functional.ComponentPointEvaluation(ref_el, d, (sd,), pt) - for d in range(sd) for pt in pts])) + nodes.extend(functional.PointEvaluationBlock(ref_el, sd, 0, degree=degree+sd-1, comp=(d,), shape=(sd,)) + for d in range(sd)) super().__init__(nodes, ref_el) From e38875d4f413e4d9a2cec3221d44d448d9f01086 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 17 Jul 2025 23:29:46 -0400 Subject: [PATCH 15/17] cleanup --- FIAT/brezzi_douglas_marini.py | 12 +++------- FIAT/functional.py | 22 ++++++++++-------- FIAT/nedelec.py | 2 +- FIAT/nedelec_second_kind.py | 43 ++++++++--------------------------- FIAT/raviart_thomas.py | 12 ++++------ 5 files changed, 30 insertions(+), 61 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 4ff66cdb8..e462e5c20 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -13,12 +13,10 @@ class BDMDualSet(dual_set.DualSet): def __init__(self, ref_el, degree, variant, interpolant_deg): + nodes = [] sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - # set to empty - nodes = [] - if variant == "integral": facet = ref_el.construct_subelement(sd-1) # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} @@ -29,13 +27,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): nodes.append(functional.FacetNormalIntegralMomentBlock(ref_el, sd-1, f, Q_ref, Pq)) elif variant == "point": - # Define each functional for the dual set - # codimension 1 facets + # Normal component evaluation at lattice points on facets for f in top[sd - 1]: - pts_cur = ref_el.make_points(sd - 1, f, sd + degree) - nodes.append(functional.FunctionalBlock(ref_el, sd-1, f, [ - functional.PointScaledNormalEvaluation(ref_el, f, pt) - for pt in pts_cur])) + nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, sd-1, f, direction="normal", degree=degree+sd)) # internal nodes if degree > 1: diff --git a/FIAT/functional.py b/FIAT/functional.py index b3d0890b1..f28b25c73 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -104,13 +104,19 @@ class PointDirectionalEvaluationBlock(FunctionalBlock): def __init__(self, ref_el, entity_dim, entity_id, direction=None, degree=0, variant=None): pts = ref_el.make_points(entity_dim, entity_id, degree, variant=variant) if direction == "normal": - nodes = [PointScaledNormalEvaluation(ref_el, entity_id, pt) for pt in pts] + sd = ref_el.get_spatial_dimension() + if entity_dim == sd-1: + nodes = [PointScaledNormalEvaluation(ref_el, entity_id, pt) for pt in pts] + else: + raise ValueError("Normals are only defined on codim=1 facets.") elif direction == "tangential": if entity_dim == 1: nodes = [PointEdgeTangentEvaluation(ref_el, entity_id, pt) for pt in pts] else: nodes = [PointFaceTangentEvaluation(ref_el, entity_id, i, pt) for i in range(entity_dim) for pt in pts] + else: + raise ValueError(f"Unrecognized direction {direction}.") super().__init__(ref_el, entity_dim, entity_id, nodes) @@ -151,10 +157,10 @@ def __init__(self, ref_el, entity_dim, entity_id): class FacetIntegralMomentBlock(FunctionalBlock): def __init__(self, ref_el, entity_dim, entity_id, Q_ref, P, mapping="L2 piola", comp=(), shape=None): dim = P.ref_el.get_spatial_dimension() - self.Phis = P.tabulate(Q_ref.get_points())[(0,) * dim] + Phis = P.tabulate(Q_ref.get_points())[(0,) * dim] + Phis = self.mapping(Phis) Q = quadrature.FacetQuadratureRule(ref_el, entity_dim, entity_id, Q_ref) - Phis = self.get_functions() phis = pullback(mapping, Q.jacobian(), Q.jacobian_determinant(), Phis) if shape is not None: @@ -164,8 +170,8 @@ def __init__(self, ref_el, entity_dim, entity_id, Q_ref, P, mapping="L2 piola", super().__init__(ref_el, entity_dim, entity_id, nodes) - def get_functions(self): - return self.Phis + def mapping(self, Phis): + return Phis class FacetDirectionalIntegralMomentBlock(FacetIntegralMomentBlock): @@ -173,9 +179,8 @@ def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis, direction): self.direction = direction super().__init__(ref_el, entity_dim, entity_id, Q_ref, Phis) - def get_functions(self): + def mapping(self, Phis): udir = self.direction - Phis = self.Phis return Phis[(slice(None), *(None for _ in range(udir.ndim)), slice(None))] * udir[(*(None for _ in range(Phis.ndim-1)), Ellipsis, None)] @@ -221,9 +226,6 @@ def __init__(self, ref_el, target_shape, pt_dict, deriv_dict, else: self.max_deriv_order = 0 - def nodes(self, entity_dim, entity_id): - yield self - def evaluate(self, f): """Obsolete and broken functional evaluation. diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index ea2aee0ba..31ad36d8c 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -122,7 +122,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): elif variant == "point": for dim in range(1, sd): - for i in top[2]: + for i in top[dim]: nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, dim, i, direction="tangential", degree=degree+1)) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 96d6d27d5..a2e93eca2 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -44,49 +44,26 @@ class NedelecSecondKindDual(DualSet): """ def __init__(self, cell, degree, variant, interpolant_deg): # Extract spatial dimension and topology - d = cell.get_spatial_dimension() - assert (d in (2, 3)), "Second kind Nedelecs only implemented in 2/3D." + sd = cell.get_spatial_dimension() + assert (sd in (2, 3)), "Second kind Nedelecs only implemented in 2/3D." - # Zero vertex-based degrees of freedom nodes = [] + for dim in range(1, sd+1): + nodes.extend(self._generate_facet_dofs(dim, cell, degree, variant, interpolant_deg)) - # (degree+1) degrees of freedom per entity of codimension 1 (edges) - nodes.extend(self._generate_edge_dofs(cell, degree, variant, interpolant_deg)) - - # Include face degrees of freedom if 3D - if d == 3: - nodes.extend(self._generate_facet_dofs(d-1, cell, degree, - variant, interpolant_deg)) - - # Varying degrees of freedom (possibly zero) per cell - nodes.extend(self._generate_facet_dofs(d, cell, degree, variant, interpolant_deg)) - - # Call init of super-class super().__init__(nodes, cell) - def _generate_edge_dofs(self, cell, degree, variant, interpolant_deg): - """Generate degrees of freedom (dofs) for entities of - codimension 1 (edges).""" - - if variant == "integral": - return self._generate_facet_dofs(1, cell, degree, variant, interpolant_deg) - - # (degree+1) tangential component point evaluation degrees of - # freedom per entity of codimension 1 (edges) - top = cell.get_topology() - nodes = [] - for edge in top[1]: - # A tangential component evaluation for each point - nodes.append(PointDirectionalEvaluationBlock(cell, 1, edge, direction="tangential", degree=degree+2)) - - return nodes - def _generate_facet_dofs(self, dim, cell, degree, variant, interpolant_deg): """Generate degrees of freedom (dofs) for facets.""" # Initialize empty dofs top = cell.get_topology() nodes = [] + if dim == 1 and variant == "point": + for entity in top[dim]: + # A tangential component evaluation for each point + nodes.append(PointDirectionalEvaluationBlock(cell, dim, entity, direction="tangential", degree=degree+2)) + return nodes # Return empty info if not applicable rt_degree = degree - dim + 1 @@ -109,8 +86,6 @@ def _generate_facet_dofs(self, dim, cell, degree, variant, interpolant_deg): Phi = RT.get_nodal_basis() mapping, = set(RT.mapping()) - # Evaluate basis functions at reference quadrature points - # Iterate over the facets for entity in top[dim]: # Construct degrees of freedom as integral moments on this cell, diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 6581b4dfc..234a2a9d8 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -58,12 +58,10 @@ class RTDualSet(dual_set.DualSet): moments against polynomials""" def __init__(self, ref_el, degree, variant, interpolant_deg): + nodes = [] sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - # set to empty - nodes = [] - if variant == "integral": facet = ref_el.construct_subelement(sd-1) # Facet nodes are \int_F v\cdot n p ds where p \in P_q @@ -83,13 +81,13 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for i in range(sd)) elif variant == "point": - # codimension 1 facets + # Normal component evaluation at lattice points on facets for i in top[sd - 1]: - nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, sd-1, i, degree=degree+sd-1, direction="normal")) + nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, sd-1, i, direction="normal", degree=degree+sd-1)) # internal nodes. Let's just use points at a lattice - if degree > 1: - nodes.extend(functional.PointEvaluationBlock(ref_el, sd, 0, degree=degree+sd-1, comp=(d,), shape=(sd,)) + for i in top[sd]: + nodes.extend(functional.PointEvaluationBlock(ref_el, sd, i, degree=degree+sd-1, comp=(d,), shape=(sd,)) for d in range(sd)) super().__init__(nodes, ref_el) From 8db851a12d45c17ebd2825d7316ff61b1d02ac03 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 19 Jul 2025 14:59:28 -0400 Subject: [PATCH 16/17] Cleanup --- FIAT/functional.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index f28b25c73..962c8f621 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -100,6 +100,16 @@ def __init__(self, ref_el, entity_dim, entity_id, order=0, degree=0, variant=Non super().__init__(ref_el, entity_dim, entity_id, nodes) +class PointGradientBlock(PointEvaluationBlock): + def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): + super().__init__(ref_el, entity_dim, entity_id, order=1, degree=degree, variant=variant) + + +class PointHessianBlock(PointEvaluationBlock): + def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): + super().__init__(ref_el, entity_dim, entity_id, order=2, degree=degree, variant=variant) + + class PointDirectionalEvaluationBlock(FunctionalBlock): def __init__(self, ref_el, entity_dim, entity_id, direction=None, degree=0, variant=None): pts = ref_el.make_points(entity_dim, entity_id, degree, variant=variant) @@ -121,36 +131,26 @@ def __init__(self, ref_el, entity_dim, entity_id, direction=None, degree=0, vari super().__init__(ref_el, entity_dim, entity_id, nodes) -class PointGradientBlock(PointEvaluationBlock): - def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): - super().__init__(ref_el, entity_dim, entity_id, order=1, degree=degree, variant=variant) - - -class PointHessianBlock(PointEvaluationBlock): - def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): - super().__init__(ref_el, entity_dim, entity_id, order=2, degree=degree, variant=None) - - class PointDirectionalDerivativeBlock(PointEvaluationBlock): - def __init__(self, ref_el, entity_dim, entity_id, basis): + def __init__(self, ref_el, entity_dim, entity_id, basis, degree=0, variant=None): nodes = [PointDirectionalDerivative(self.ref_el, pt, e) - for pt in self.ref_el.make_points(entity_dim, entity_id, 0) + for pt in self.ref_el.make_points(entity_dim, entity_id, degree, variant=variant) for e in basis] super().__init__(ref_el, entity_dim, entity_id, nodes) class PointNormalTangentialDerivativeBlock(PointDirectionalDerivativeBlock): - def __init__(self, ref_el, entity_dim, entity_id): + def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): sd = self.ref_el.get_spatial_dimension() basis = numpy.zeros((sd, sd)) basis[0] = self.ref_el.compute_scaled_normal(entity_id) basis[1:] = self.ref_el.compute_tangents(sd-1, entity_id) - super().__init__(ref_el, entity_dim, entity_id, basis) + super().__init__(ref_el, entity_dim, entity_id, basis, degree=degree, variant=variant) class PointNormalDerivativeView(FunctionalBlockView): - def __init__(self, ref_el, entity_dim, entity_id): - block = PointNormalTangentialDerivativeBlock(ref_el, entity_dim, entity_id) + def __init__(self, ref_el, entity_dim, entity_id, degree=0, variant=None): + block = PointNormalTangentialDerivativeBlock(ref_el, entity_dim, entity_id, degree=0, variant=None) super().__init__(block, [0]) From a7854d636af2beda1efac38771c43356154906d9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 21 Jul 2025 11:07:05 +0100 Subject: [PATCH 17/17] Docstrings, cleanup --- FIAT/macro.py | 65 ++++++++++++++++++++++++++++++++++----------------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/FIAT/macro.py b/FIAT/macro.py index 6327b6724..7d2abb90b 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -485,7 +485,15 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): super().__init__(ref_el, degree, degree, expansion_set, coeffs) -def hdiv_conforming(U, order=0): +def hdiv_conforming_coefficients(U, order=0): + """Constrains a PolynomialSet to have vanishing normal jumps on the + interior facets of a simplicial complex. + + :arg U: A PolynomialSet, it may be vector- or tensor-valued. + :arg order: The maximum order of differentiation on the normal jump constraints. + + :returns: The coefficients of the constrained PolynomialSet. + """ from FIAT.quadrature_schemes import create_quadrature degree = U.degree ref_el = U.get_reference_element() @@ -503,18 +511,16 @@ def hdiv_conforming(U, order=0): qpts, qwts = Q.get_points(), Q.get_weights() phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] weights = numpy.multiply(phi_at_qpts, qwts) + ax = tuple(range(1, weights.ndim)) rows = [] for facet in ref_el.get_interior_facets(sd-1): - normal = ref_el.compute_normal(facet) - normal = normal[(None,) * len(shape) + (slice(None), None)] + normal = ref_el.compute_scaled_normal(facet) + ncoeffs = numpy.tensordot(coeffs, normal, axes=(len(shape), 0)) jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) for r in range(k, order+1): - jump = numpy.dot(coeffs, jumps[r]) - # num_wt = 1 if sd == 1 else expansions.polynomial_dimension(facet_el, degree-r) - wn = weights[..., None, :] * normal - ax = tuple(range(1, len(wn.shape))) - rows.append(numpy.tensordot(wn, jump, axes=(ax, ax))) + njump = numpy.dot(ncoeffs, jumps[r]) + rows.append(numpy.tensordot(weights, njump, axes=(ax, ax))) if len(rows) > 0: dual_mat = numpy.vstack(rows) @@ -536,7 +542,7 @@ class HDivPolynomialSet(polynomial_set.PolynomialSet): def __init__(self, ref_el, degree, order=0, **kwargs): sd = ref_el.get_spatial_dimension() U = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), **kwargs) - coeffs = hdiv_conforming(U, order=order) + coeffs = hdiv_conforming_coefficients(U, order=order) super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) @@ -552,11 +558,21 @@ class HDivSymPolynomialSet(polynomial_set.PolynomialSet): """ def __init__(self, ref_el, degree, order=0, **kwargs): U = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree, **kwargs) - coeffs = hdiv_conforming(U, order=order) + coeffs = hdiv_conforming_coefficients(U, order=order) super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) -def pullback(mapping, Jinv, phi): +def pullback(phi, mapping, J=None, Jinv=None, Jdet=None): + """Transforms a reference tabulation into physical space. + + :arg phi: The tabulation on reference space. + :arg mapping: A string indicating the pullback type. + :arg J: The Jacobian of the transformation from reference to physical space. + :arg Jinv: The inverse of the Jacobian. + :arg Jdet: The determinant of the Jacobian. + + :returns: The tabulation on physical space. + """ try: formdegree = { "affine": (0,), @@ -569,26 +585,33 @@ def pullback(mapping, Jinv, phi): except KeyError: raise ValueError(f"Unrecognized mapping {mapping}") + if J is None: + J = numpy.linalg.pinv(Jinv) + if Jinv is None: + Jinv = numpy.linalg.pinv(J) + if Jdet is None: + Jdet = numpy.linalg.det(J) F1 = Jinv.T - F2 = numpy.linalg.inv(Jinv) * numpy.linalg.det(Jinv) - - perm = [*range(1, phi.ndim-1), 0, -1] - phi = phi.transpose(perm) + F2 = J / Jdet for i, k in enumerate(formdegree): if k == 0: continue F = F1 if k == 1 else F2 - phi = numpy.tensordot(F, phi, (1, i)) - perm = [-2, *reversed(range(phi.ndim-2)), -1] - phi = phi.transpose(perm) + perm = list(range(phi.ndim)) + perm[i+1], perm[-1] = perm[-1], perm[i+1] + + phi = phi.transpose(perm) + phi = phi.dot(F.T) + phi = phi.transpose(perm) + return phi class MacroPolynomialSet(polynomial_set.PolynomialSet): """Constructs a PolynomialSet by tiling a CiarletElement on a - SimplicialComplex. + SimplicialComplex. :arg ref_el: The SimplicialComplex. :arg element: The CiarletElement. @@ -615,9 +638,9 @@ def __init__(self, ref_el, element): cmap = expansion_set.get_cell_node_map(n) for cell in sorted(top[sd]): cell_verts = ref_el.get_vertices_of_subcomplex(top[sd][cell]) - A, b = reference_element.make_affine_mapping(cell_verts, base_ref_el.vertices) + A, b = reference_element.make_affine_mapping(base_ref_el.vertices, cell_verts) indices = numpy.ix_(rmap[cell], *map(range, shp), cmap[cell]) - coeffs[indices] = pullback(mapping, A, base_coeffs) + coeffs[indices] = pullback(base_coeffs, mapping, J=A) super().__init__(ref_el, element.degree(), element.degree(), expansion_set, coeffs)