diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 0014b5277..446962aea 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -9,38 +9,27 @@ 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)] + ells = functional.NormalMoments(ref_el, Q, Pq) 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(ells.generate(sd-1, f)) + entity_ids[sd - 1][f].extend(range(cur, len(nodes))) elif variant == "point": # Define each functional for the dual set @@ -50,7 +39,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 +49,9 @@ 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))) + 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 1284cd5e2..acb1583f0 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,64 @@ 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") + + +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) + + +class FacetIntegralMoments(FrobeniusIntegralMoments): + def get_entity_transform(self, dim, entity): + return lambda f: f + + 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) + + +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) + + +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 c92fa06be..1ee2e348b 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,13 @@ 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)) - + ells = functional.TangentialMoments(ref_el, Q, Pqmd) 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(ells.generate(dim, entity)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) elif variant == "point": for i in top[1]: @@ -145,7 +136,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 +146,19 @@ 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,)) + 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 481dcf64e..1f40cb464 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,25 @@ 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() + 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(IntegralMoment(cell, Q_facet, phi) for phi in phis) + dofs.extend(ells.generate(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..e7ea53264 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,32 @@ 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)] + ells = functional.NormalMoments(ref_el, Q, Pq) 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(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 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,)) + 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": # codimension 1 facets @@ -105,7 +92,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 +101,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)