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 0014b5277..e462e5c20 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -9,7 +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 class BDMDualSet(dual_set.DualSet): @@ -18,54 +17,34 @@ 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] = [] - 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) 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.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]: - 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.append(functional.PointDirectionalEvaluationBlock(ref_el, sd-1, f, direction="normal", degree=degree+sd)) # internal nodes 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) - 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))) - super().__init__(nodes, ref_el, entity_ids) + 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()) + Phi = Nedel.get_nodal_basis() + for entity in top[sd]: + nodes.append(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Phi, mapping=mapping)) + + super().__init__(nodes, ref_el) class BrezziDouglasMarini(finite_element.CiarletElement): @@ -90,7 +69,9 @@ class BrezziDouglasMarini(finite_element.CiarletElement): def __init__(self, ref_el, degree, variant=None): - 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") diff --git a/FIAT/check_format_variant.py b/FIAT/check_format_variant.py index 1431876de..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): @@ -61,7 +63,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 +83,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: 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/dual_set.py b/FIAT/dual_set.py index c09db974f..4432bfe40 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -12,7 +12,20 @@ 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 entity_ids is None: + # flatten nodes + top = ref_el.get_topology() + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + blocks = nodes + 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 self.ref_el = ref_el diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 69d0d4b5a..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() @@ -371,7 +378,7 @@ 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 = {} @@ -599,7 +606,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 +630,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/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/functional.py b/FIAT/functional.py index 18de29942..962c8f621 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,172 @@ 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, 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 + + +class FunctionalBlockView(FunctionalBlock): + def __init__(self, block, subset): + self.block = block + nodes = [block.nodes[i] for i in subset] + super().__init__(block.ref_el, block.entity_dim, block.entity_id, nodes) + + def completion(self, entity_dim, entity_id): + return self.block + + +class RelabeledFunctionalBlock(FunctionalBlock): + def __init__(self, block, entity_dim, entity_id): + super().__init__(block.ref_el, entity_dim, entity_id, block.nodes) + + +class PointEvaluationBlock(FunctionalBlock): + def __init__(self, ref_el, entity_dim, entity_id, order=0, degree=0, variant=None, comp=(), shape=None): + from FIAT.polynomial_set import mis + pts = ref_el.make_points(entity_dim, entity_id, degree, variant=variant) + if order == 0: + 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 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) + if direction == "normal": + 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) + + +class PointDirectionalDerivativeBlock(PointEvaluationBlock): + 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, 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, 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, degree=degree, variant=variant) + + +class PointNormalDerivativeView(FunctionalBlockView): + 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]) + + +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() + 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 = pullback(mapping, Q.jacobian(), Q.jacobian_determinant(), 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) + + def mapping(self, Phis): + return Phis + + +class FacetDirectionalIntegralMomentBlock(FacetIntegralMomentBlock): + 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 mapping(self, Phis): + udir = self.direction + 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, 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): 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/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py index 656cee6bb..6fe608864 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.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 @@ -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,7 +55,14 @@ class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement): the weak symmetry constraint. """ - def __init__(self, ref_el, degree): + def __init__(self, ref_el, degree, variant=None): + + splitting, variant, interpolant_deg = check_format_variant(variant, degree) + assert variant == "integral" + + 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() @@ -68,7 +71,7 @@ def __init__(self, ref_el, degree): 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. @@ -77,12 +80,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..52697f0a5 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/FIAT/hellan_herrmann_johnson.py @@ -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) - entity_ids[sd][0].extend(range(cur, len(nodes))) + # 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][entity].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) @@ -93,7 +100,10 @@ def __init__(self, ref_el, degree=0, variant=None): if degree < 0: raise ValueError(f"{type(self).__name__} only defined for degree >= 0") - variant, qdegree = check_format_variant(variant, 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() 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/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/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/macro.py b/FIAT/macro.py index 07a809acb..7d2abb90b 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) @@ -485,6 +485,67 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): super().__init__(ref_el, degree, degree, expansion_set, coeffs) +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() + 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) + ax = tuple(range(1, weights.ndim)) + + rows = [] + for facet in ref_el.get_interior_facets(sd-1): + 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): + njump = numpy.dot(ncoeffs, jumps[r]) + rows.append(numpy.tensordot(weights, njump, 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. + :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): + sd = ref_el.get_spatial_dimension() + U = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), **kwargs) + coeffs = hdiv_conforming_coefficients(U, order=order) + super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) + + class HDivSymPolynomialSet(polynomial_set.PolynomialSet): """Constructs a symmetric tensor-valued PolynomialSet with continuous normal components on a simplicial complex. @@ -496,36 +557,90 @@ 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 + coeffs = hdiv_conforming_coefficients(U, order=order) + super().__init__(ref_el, degree, degree, U.expansion_set, coeffs) + +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,), + "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}") + + 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 = J / Jdet + + for i, k in enumerate(formdegree): + if k == 0: + continue + F = F1 if k == 1 else F2 + + 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. + + :arg ref_el: The SimplicialComplex. + :arg element: The CiarletElement. + """ + def __init__(self, ref_el, element): + top = ref_el.get_topology() 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=(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) + mapping, = set(element.mapping()) + base_ref_el = element.get_reference_element() + base_entity_ids = element.entity_dofs() + n = element.degree() - 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))) + base_expansion_set = element.get_nodal_basis().get_expansion_set() + expansion_set = base_expansion_set.reconstruct(ref_el=ref_el) - 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)) + shp = element.value_shape() + 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() - super().__init__(ref_el, degree, degree, expansion_set, 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[sd]): + cell_verts = ref_el.get_vertices_of_subcomplex(top[sd][cell]) + 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(base_coeffs, mapping, J=A) + + super().__init__(ref_el, element.degree(), element.degree(), expansion_set, coeffs) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index 6cf122772..31ad36d8c 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -6,12 +6,11 @@ # 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 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): @@ -104,13 +103,6 @@ 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] = [] - if variant == "integral": # edge nodes are \int_F v\cdot t p ds where p \in P_{q-1}(edge) # degree is q - 1 @@ -124,54 +116,29 @@ 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] - 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) - entity_ids[dim][entity] = list(range(cur, len(nodes))) + nodes.append(functional.FacetIntegralMomentBlock(ref_el, dim, entity, Q_ref, Pqmd, mapping=mapping)) 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))) - - 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))) + for dim in range(1, sd): + 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) - 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))) - super().__init__(nodes, ref_el, entity_ids) + cell = ref_el.construct_subelement(sd) + Q_ref = create_quadrature(cell, interpolant_deg + phi_deg) + Pqmd = polynomial_set.ONPolynomialSet(cell, phi_deg) + for entity in top[sd]: + 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) class Nedelec(finite_element.CiarletElement): @@ -195,9 +162,14 @@ class Nedelec(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - - variant, interpolant_deg = check_format_variant(variant, degree) - if ref_el.get_spatial_dimension() == 3: + 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) + 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..a2e93eca2 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -5,15 +5,11 @@ # # 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, PointDirectionalEvaluationBlock 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 @@ -46,126 +42,58 @@ 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 - (dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg) - # Call init of super-class - super().__init__(dofs, cell, ids) - - def generate_degrees_of_freedom(self, cell, degree, variant, interpolant_deg): - "Generate dofs and geometry-to-dof maps (ids)." - - dofs = [] - ids = {} - # 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 (d+1 of these) - ids[0] = {i: [] for i in range(d + 1)} - - # (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) - - # 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) - - # 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) - - return (dofs, ids) - - def _generate_edge_dofs(self, cell, degree, offset, variant, interpolant_deg): - """Generate degrees of freedom (dofs) for entities of - codimension 1 (edges).""" + sd = cell.get_spatial_dimension() + assert (sd in (2, 3)), "Second kind Nedelecs only implemented in 2/3D." - if variant == "integral": - return self._generate_facet_dofs(1, cell, degree, offset, variant, interpolant_deg) + nodes = [] + for dim in range(1, sd+1): + nodes.extend(self._generate_facet_dofs(dim, cell, degree, variant, interpolant_deg)) - # (degree+1) tangential component point evaluation degrees of - # freedom per entity of codimension 1 (edges) - dofs = [] - ids = {} - if variant == "point": - for edge in range(len(cell.get_topology()[1])): + super().__init__(nodes, cell) - # 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) - - # Associate these dofs with this edge - i = len(points) * edge - ids[edge] = list(range(offset + i, offset + i + len(points))) - - return (dofs, ids) - - def _generate_facet_dofs(self, codim, 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()[codim]) - dofs = [] - ids = {i: [] for i in range(num_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 - codim + 1 + rt_degree = degree - dim + 1 if rt_degree < 1: - return (dofs, ids) + return nodes if interpolant_deg is None: 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,)) + + # Construct Raviart-Thomas on the reference facet + if dim == 1: + 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() - - # 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 + mapping, = set(RT.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, 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)) - + 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(IntegralMoment(cell, Q_facet, phi) for phi in phis) + nodes.append(FacetIntegralMomentBlock(cell, dim, entity, Q_ref, Phi, mapping=mapping)) - # Assign identifiers (num RTs per face + previous edge dofs) - ids[facet].extend(range(cur, cur + len(phis))) - cur += len(phis) - - return (dofs, ids) + return nodes class NedelecSecondKind(CiarletElement): @@ -191,21 +119,21 @@ class NedelecSecondKind(CiarletElement): interpolation. """ - def __init__(self, cell, degree, variant=None): - - variant, interpolant_deg = check_format_variant(variant, degree) + def __init__(self, ref_el, degree, variant=None): + 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!" # 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, )) + 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..234a2a9d8 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -6,12 +6,11 @@ # 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.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule def RTSpace(ref_el, degree): @@ -63,60 +62,35 @@ 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] = [] - 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) 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.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: - 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) + for entity in top[sd]: + 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 + # Normal component evaluation at lattice points on 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.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: - 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))) + 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, entity_ids) + super().__init__(nodes, ref_el) class RaviartThomas(finite_element.CiarletElement): @@ -140,10 +114,15 @@ class RaviartThomas(finite_element.CiarletElement): """ def __init__(self, ref_el, degree, variant=None): - - variant, interpolant_deg = check_format_variant(variant, degree) - - poly_set = RTSpace(ref_el, 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) + 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/reference_element.py b/FIAT/reference_element.py index 547da2c04..72dc140ba 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 @@ -617,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 31dd6bb44..2f0f4bc7d 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -65,7 +65,10 @@ def __init__(self, ref_el, degree=0, variant=None): if degree < 0: raise ValueError(f"{type(self).__name__} only defined for degree >= 0") - variant, qdegree = check_format_variant(variant, 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) 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 8e367bf20..a92eeae62 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -1,8 +1,11 @@ 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, 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 from FIAT.expansions import polynomial_entity_ids, polynomial_cell_node_map @@ -438,3 +441,76 @@ 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), (HellanHerrmannJohnson, 0), + (GopalakrishnanLedererSchoberlSecondKind, 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 = 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 for pt in pts_cur) + + P = MacroPolynomialSet(A, fe) + phis = P.tabulate(pts)[(0,)*dim] + shape = phis.shape + shp = shape[1:-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(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, 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): + 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) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index d5136e651..e31074552 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -162,5 +162,14 @@ 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]) +@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)