From b901896f28c9506a6650e66e67fd2a346b331903 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 2 Nov 2024 14:23:12 +0000 Subject: [PATCH 1/2] IntegralMoment factory --- FIAT/brezzi_douglas_marini.py | 32 ++++++--------------- FIAT/functional.py | 39 +++++++++++++++++++++++++- FIAT/nedelec.py | 31 ++++++-------------- FIAT/nedelec_second_kind.py | 53 ++++++++++------------------------- FIAT/raviart_thomas.py | 35 +++++++---------------- 5 files changed, 81 insertions(+), 109 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 0014b5277..0b2bf96f1 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -9,38 +9,26 @@ polynomial_set, nedelec) from FIAT.check_format_variant import check_format_variant from FIAT.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule 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] = [] + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] if variant == "integral": facet = ref_el.get_facet_element() - # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} + # Facet nodes are \int_F v.n p ds where p \in P_{q} # degree is q - Q_ref = create_quadrature(facet, interpolant_deg + degree) + Q = 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 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) - entity_ids[sd - 1][f] = list(range(cur, len(nodes))) + nodes.extend(functional.NormalMoments(ref_el, Q, Pq, f)) + entity_ids[sd - 1][f].extend(range(cur, len(nodes))) elif variant == "point": # Define each functional for the dual set @@ -50,7 +38,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): 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))) + entity_ids[sd - 1][f].extend(range(cur, len(nodes))) # internal nodes if degree > 1: @@ -60,10 +48,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q = create_quadrature(ref_el, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(ref_el, 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))) + nodes.extend(functional.FrobeniusIntegralMoments(ref_el, Q, Nedfs)) + entity_ids[sd][0].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/functional.py b/FIAT/functional.py index 1284cd5e2..c829d86ac 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -16,7 +16,7 @@ import numpy import sympy -from FIAT import polynomial_set, jacobi +from FIAT import polynomial_set, quadrature, jacobi from FIAT.quadrature import GaussLegendreQuadratureLineRule from FIAT.reference_element import UFCInterval as interval @@ -754,3 +754,40 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): pt_dict = {tuple(pt): [(wt*nnT[idx], idx) for idx in index_iterator(shp)] for pt, wt in zip(points, weights)} super().__init__(ref_el, shp, pt_dict, {}, "IntegralMomentOfNormalNormalEvaluation") + + +def IntegralMoments(ref_el, Q, P, **kwargs): + sd = ref_el.get_spatial_dimension() + phis = P.tabulate(Q.get_points())[(0,)*sd] + for phi in phis: + yield IntegralMoment(ref_el, Q, phi, **kwargs) + + +def FrobeniusIntegralMoments(ref_el, Q, P, **kwargs): + sd = ref_el.get_spatial_dimension() + phis = P.tabulate(Q.get_points())[(0,)*sd] + for phi in phis: + yield FrobeniusIntegralMoment(ref_el, Q, phi, **kwargs) + + +def FacetIntegralMoments(ref_el, Q, P, dim, entity_id, transform=None): + if transform is None: + transform = lambda f: f + phis = P.tabulate(Q.get_points())[(0,)*dim] + Q_mapped = quadrature.FacetQuadratureRule(ref_el, dim, entity_id, Q) + phis /= Q_mapped.jacobian_determinant() + for phi in phis: + yield FrobeniusIntegralMoment(ref_el, Q_mapped, transform(phi)) + + +def NormalMoments(ref_el, Q, P, face_id): + dim = ref_el.get_spatial_dimension() - 1 + n = ref_el.compute_scaled_normal(face_id) + transform = lambda f: numpy.multiply(n[..., None], f) + yield from FacetIntegralMoments(ref_el, Q, P, dim, face_id, transform) + + +def TangentialMoments(ref_el, Q, P, dim, entity_id): + ts = ref_el.compute_tangents(dim, entity_id) + transform = lambda f: numpy.dot(ts.T, f) + yield from FacetIntegralMoments(ref_el, Q, P, dim, entity_id, transform) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c92fa06be..0ab615a6a 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -11,7 +11,6 @@ import numpy 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,21 +121,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): phi_deg = degree - dim if phi_deg >= 0: facet = ref_el.construct_subelement(dim) - Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) + Q = 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)) - 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) - entity_ids[dim][entity] = list(range(cur, len(nodes))) + nodes.extend(functional.TangentialMoments(ref_el, Q, Pqmd, dim, entity)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) elif variant == "point": for i in top[1]: @@ -145,7 +135,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): 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))) + entity_ids[1][i].extend(range(cur, len(nodes))) if sd > 2 and degree > 1: # face tangents for i in top[2]: # loop over faces @@ -155,21 +145,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): 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))) + entity_ids[2][i].extend(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))) + Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg, shape=(sd,)) + nodes.extend(functional.FrobeniusIntegralMoments(ref_el, Q, Pqmd)) + entity_ids[sd][0].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 481dcf64e..f36761da0 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 TangentialMoments 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 @@ -109,16 +106,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]) + top = cell.get_topology() dofs = [] - ids = {i: [] for i in range(num_facets)} + ids = {i: [] for i in top[dim]} # Return empty info if not applicable - rt_degree = degree - codim + 1 + rt_degree = degree - dim + 1 if rt_degree < 1: return (dofs, ids) @@ -126,44 +123,24 @@ 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) - Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree) - if codim == 1: - Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,)) + facet = cell.construct_subelement(dim) + Q = create_quadrature(facet, interpolant_deg + rt_degree) + if dim == 1: + P = ONPolynomialSet(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] - # 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, codim, 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)) + RT = RaviartThomas(facet, rt_degree, variant) + P = RT.get_nodal_basis() + for entity in sorted(top[dim]): + cur = 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(TangentialMoments(cell, Q, P, 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, len(dofs))) return (dofs, ids) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index db0058d12..76532edb5 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): @@ -59,44 +58,30 @@ 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() - entity_ids = {} - # set to empty - for dim in top: - entity_ids[dim] = {} - for entity in top[dim]: - entity_ids[dim][entity] = [] + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] if variant == "integral": facet = ref_el.get_facet_element() # 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) + Q = 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)] 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) - entity_ids[sd - 1][f] = list(range(cur, len(nodes))) + nodes.extend(functional.NormalMoments(ref_el, Q, Pq, f)) + entity_ids[sd - 1][f].extend(range(cur, len(nodes))) # 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))) + Pqm1 = polynomial_set.ONPolynomialSet(ref_el, q - 1, shape=(sd,)) + nodes.extend(functional.FrobeniusIntegralMoments(ref_el, Q, Pqm1)) + entity_ids[sd][0].extend(range(cur, len(nodes))) elif variant == "point": # codimension 1 facets @@ -105,7 +90,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): 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))) + entity_ids[sd - 1][i].extend(range(cur, len(nodes))) # internal nodes. Let's just use points at a lattice if degree > 1: @@ -114,7 +99,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): 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))) + entity_ids[sd][0].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) From b5b15bac78b60304ddf81034c0d07169d84b39b3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 2 Nov 2024 15:52:01 +0000 Subject: [PATCH 2/2] FunctionalGenerator classes to reuse tabulations on facets --- FIAT/brezzi_douglas_marini.py | 6 ++- FIAT/functional.py | 78 +++++++++++++++++++++++------------ FIAT/nedelec.py | 6 ++- FIAT/nedelec_second_kind.py | 3 +- FIAT/raviart_thomas.py | 6 ++- 5 files changed, 65 insertions(+), 34 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 0b2bf96f1..446962aea 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -25,9 +25,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # degree is q Q = create_quadrature(facet, interpolant_deg + degree) Pq = polynomial_set.ONPolynomialSet(facet, degree) + ells = functional.NormalMoments(ref_el, Q, Pq) for f in top[sd - 1]: cur = len(nodes) - nodes.extend(functional.NormalMoments(ref_el, Q, Pq, f)) + nodes.extend(ells.generate(sd-1, f)) entity_ids[sd - 1][f].extend(range(cur, len(nodes))) elif variant == "point": @@ -48,7 +49,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): Q = create_quadrature(ref_el, interpolant_deg + degree - 1) Nedel = nedelec.Nedelec(ref_el, degree - 1, variant) Nedfs = Nedel.get_nodal_basis() - nodes.extend(functional.FrobeniusIntegralMoments(ref_el, Q, Nedfs)) + ells = functional.FrobeniusIntegralMoments(ref_el, Q, Nedfs) + nodes.extend(ells.generate(sd, 0)) entity_ids[sd][0].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/functional.py b/FIAT/functional.py index c829d86ac..acb1583f0 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -756,38 +756,62 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): super().__init__(ref_el, shp, pt_dict, {}, "IntegralMomentOfNormalNormalEvaluation") -def IntegralMoments(ref_el, Q, P, **kwargs): - sd = ref_el.get_spatial_dimension() - phis = P.tabulate(Q.get_points())[(0,)*sd] - for phi in phis: - yield IntegralMoment(ref_el, Q, phi, **kwargs) +class FunctionalGenerator(): + def __init__(self, ref_el): + self.ref_el = ref_el + + def generate(self, dim, entity): + raise NotImplementedError + + +class IntegralMoments(FunctionalGenerator): + def __init__(self, ref_el, Q, P, comp=tuple(), shp=None): + if shp is None: + shp = P.get_shape() + sd = P.ref_el.get_spatial_dimension() + self.phis = P.tabulate(Q.get_points())[(0,)*sd] + self.Q = Q + self.comp = comp + self.shape = shp + super().__init__(ref_el) + + def generate(self, dim, entity): + if dim == self.ref_el.get_spatial_dimension(): + assert entity == 0 + for phi in self.phis: + yield IntegralMoment(self.ref_el, self.Q, phi, comp=self.comp, shp=self.shape) + + +class FrobeniusIntegralMoments(IntegralMoments): + def __init__(self, ref_el, Q, P): + super().__init__(ref_el, Q, P, comp=slice(None)) + def generate(self, dim, entity): + if dim == self.ref_el.get_spatial_dimension(): + assert entity == 0 + for phi in self.phis: + yield FrobeniusIntegralMoment(self.ref_el, self.Q, phi) -def FrobeniusIntegralMoments(ref_el, Q, P, **kwargs): - sd = ref_el.get_spatial_dimension() - phis = P.tabulate(Q.get_points())[(0,)*sd] - for phi in phis: - yield FrobeniusIntegralMoment(ref_el, Q, phi, **kwargs) +class FacetIntegralMoments(FrobeniusIntegralMoments): + def get_entity_transform(self, dim, entity): + return lambda f: f -def FacetIntegralMoments(ref_el, Q, P, dim, entity_id, transform=None): - if transform is None: - transform = lambda f: f - phis = P.tabulate(Q.get_points())[(0,)*dim] - Q_mapped = quadrature.FacetQuadratureRule(ref_el, dim, entity_id, Q) - phis /= Q_mapped.jacobian_determinant() - for phi in phis: - yield FrobeniusIntegralMoment(ref_el, Q_mapped, transform(phi)) + def generate(self, dim, entity): + transform = self.get_entity_transform(dim, entity) + Q_mapped = quadrature.FacetQuadratureRule(self.ref_el, dim, entity, self.Q) + Jdet = Q_mapped.jacobian_determinant() + for phi in self.phis: + yield FrobeniusIntegralMoment(self.ref_el, Q_mapped, transform(phi)/Jdet) -def NormalMoments(ref_el, Q, P, face_id): - dim = ref_el.get_spatial_dimension() - 1 - n = ref_el.compute_scaled_normal(face_id) - transform = lambda f: numpy.multiply(n[..., None], f) - yield from FacetIntegralMoments(ref_el, Q, P, dim, face_id, transform) +class NormalMoments(FacetIntegralMoments): + def get_entity_transform(self, dim, entity): + n = self.ref_el.compute_scaled_normal(entity) + return lambda f: numpy.multiply(n[..., None], f) -def TangentialMoments(ref_el, Q, P, dim, entity_id): - ts = ref_el.compute_tangents(dim, entity_id) - transform = lambda f: numpy.dot(ts.T, f) - yield from FacetIntegralMoments(ref_el, Q, P, dim, entity_id, transform) +class TangentialMoments(FacetIntegralMoments): + def get_entity_transform(self, dim, entity): + ts = self.ref_el.compute_tangents(dim, entity) + return lambda f: numpy.dot(ts.T, f) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 0ab615a6a..1ee2e348b 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -123,9 +123,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): facet = ref_el.construct_subelement(dim) Q = create_quadrature(facet, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) + ells = functional.TangentialMoments(ref_el, Q, Pqmd) for entity in top[dim]: cur = len(nodes) - nodes.extend(functional.TangentialMoments(ref_el, Q, Pqmd, dim, entity)) + nodes.extend(ells.generate(dim, entity)) entity_ids[dim][entity].extend(range(cur, len(nodes))) elif variant == "point": @@ -155,7 +156,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cur = len(nodes) Q = create_quadrature(ref_el, interpolant_deg + phi_deg) Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg, shape=(sd,)) - nodes.extend(functional.FrobeniusIntegralMoments(ref_el, Q, Pqmd)) + ells = functional.FrobeniusIntegralMoments(ref_el, Q, Pqmd) + nodes.extend(ells.generate(sd, 0)) entity_ids[sd][0].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index f36761da0..1f40cb464 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -132,12 +132,13 @@ def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_d RT = RaviartThomas(facet, rt_degree, variant) P = RT.get_nodal_basis() + ells = TangentialMoments(cell, Q, P) for entity in sorted(top[dim]): cur = 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(TangentialMoments(cell, Q, P, dim, entity)) + dofs.extend(ells.generate(dim, entity)) # Assign identifiers (num RTs per face + previous edge dofs) ids[entity].extend(range(cur, len(dofs))) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 76532edb5..e7ea53264 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -70,9 +70,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): q = degree - 1 Q = create_quadrature(facet, interpolant_deg + q) Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0) + ells = functional.NormalMoments(ref_el, Q, Pq) for f in top[sd - 1]: cur = len(nodes) - nodes.extend(functional.NormalMoments(ref_el, Q, Pq, f)) + nodes.extend(ells.generate(sd-1, f)) entity_ids[sd - 1][f].extend(range(cur, len(nodes))) # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d @@ -80,7 +81,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): cur = len(nodes) Q = create_quadrature(ref_el, interpolant_deg + q - 1) Pqm1 = polynomial_set.ONPolynomialSet(ref_el, q - 1, shape=(sd,)) - nodes.extend(functional.FrobeniusIntegralMoments(ref_el, Q, Pqm1)) + ells = functional.FrobeniusIntegralMoments(ref_el, Q, Pqm1) + nodes.extend(ells.generate(sd, 0)) entity_ids[sd][0].extend(range(cur, len(nodes))) elif variant == "point":