From 7aa40e677993db5dc7b0c3592e46ef70dd44bfd5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 15 Jan 2025 16:24:12 +0000 Subject: [PATCH] CrouzeixRaviart: integral variant (#118) * CrouzeixRaviart: integral variant * High-order CR variants --- FIAT/crouzeix_raviart.py | 100 +++++++++++++++--------- FIAT/hierarchical.py | 7 +- FIAT/reference_element.py | 4 +- finat/fiat_elements.py | 4 +- finat/ufl/elementlist.py | 2 +- gem/utils.py | 20 +---- test/FIAT/regression/test_regression.py | 2 +- test/FIAT/unit/test_fiat.py | 7 ++ 8 files changed, 79 insertions(+), 67 deletions(-) diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index b809f7a87..8447a1269 100644 --- a/FIAT/crouzeix_raviart.py +++ b/FIAT/crouzeix_raviart.py @@ -9,62 +9,84 @@ # # Last changed: 2010-01-28 +import numpy from FIAT import finite_element, polynomial_set, dual_set, functional - - -def _initialize_entity_ids(topology): - entity_ids = {} - for (i, entity) in list(topology.items()): - entity_ids[i] = {} - for j in entity: - entity_ids[i][j] = [] - return entity_ids +from FIAT.check_format_variant import check_format_variant +from FIAT.quadrature_schemes import create_quadrature +from FIAT.quadrature import FacetQuadratureRule class CrouzeixRaviartDualSet(dual_set.DualSet): - """Dual basis for Crouzeix-Raviart element (linears continuous at - boundary midpoints).""" - - def __init__(self, cell, degree): + def __init__(self, ref_el, degree, variant, interpolant_deg): # Get topology dictionary - d = cell.get_spatial_dimension() - topology = cell.get_topology() + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + + if degree > 1 and sd != 2: + raise NotImplementedError("High-order Crouzeix-Raviart is only implemented on triangles.") # Initialize empty nodes and entity_ids - entity_ids = _initialize_entity_ids(topology) - nodes = [None for i in list(topology[d - 1].keys())] + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] # Construct nodes and entity_ids - for i in topology[d - 1]: - - # Construct midpoint - x = cell.make_points(d - 1, i, d)[0] - - # Degree of freedom number i is evaluation at midpoint - nodes[i] = functional.PointEvaluation(cell, x) - entity_ids[d - 1][i] += [i] - - # Initialize super-class - super().__init__(nodes, cell, entity_ids) + if variant == "integral": + for dim in sorted(top): + if dim == 0 and dim != sd-1: + # Skip vertex dofs + continue + facet = ref_el.construct_subelement(dim) + if dim == 0: + Q_facet = create_quadrature(facet, degree + interpolant_deg-1) + Phis = numpy.ones((1, len(Q_facet.pts))) + else: + k = degree - 1 if dim == sd-1 else degree - (1+dim) + if k < 0: + continue + Q_facet = create_quadrature(facet, k + interpolant_deg) + poly_set = polynomial_set.ONPolynomialSet(facet, k) + Phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim] + + for i in sorted(top[dim]): + cur = len(nodes) + Q = FacetQuadratureRule(ref_el, dim, i, Q_facet) + scale = 1 / Q.jacobian_determinant() + phis = scale * Phis + nodes.extend(functional.IntegralMoment(ref_el, Q, phi) for phi in phis) + entity_ids[dim][i].extend(range(cur, len(nodes))) + else: + for dim in sorted(top): + if dim == 0 and dim != sd-1: + # Skip vertex dofs + continue + for i in sorted(top[dim]): + cur = len(nodes) + if dim == sd-1 and dim != 0: + pts = ref_el.make_points(dim, i, degree-1, variant="gl", interior=0) + else: + pts = ref_el.make_points(dim, i, degree, variant="gll") + nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts) + entity_ids[dim][i].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) class CrouzeixRaviart(finite_element.CiarletElement): """The Crouzeix-Raviart finite element: K: Triangle/Tetrahedron - Polynomial space: P_1 - Dual basis: Evaluation at facet midpoints + Polynomial space: P_k + Dual basis: Evaluation at points or integral moments """ - def __init__(self, cell, degree): + def __init__(self, ref_el, degree, variant=None): + + variant, interpolant_deg = check_format_variant(variant, degree) - # Crouzeix Raviart is only defined for polynomial degree == 1 - if not (degree == 1): - raise Exception("Crouzeix-Raviart only defined for degree 1") + if degree % 2 != 1: + raise ValueError("Crouzeix-Raviart only defined for odd degree") - # Construct polynomial spaces, dual basis and initialize - # FiniteElement - space = polynomial_set.ONPolynomialSet(cell, 1) - dual = CrouzeixRaviartDualSet(cell, 1) - super().__init__(space, dual, 1) + space = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble") + dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg) + super().__init__(space, dual, degree) diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index e8666572b..17c415710 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -17,12 +17,13 @@ from FIAT.check_format_variant import parse_lagrange_variant -def make_dual_bubbles(ref_el, degree, codim=0): +def make_dual_bubbles(ref_el, degree, codim=0, interpolant_deg=None): """Tabulate the L2-duals of the hierarchical C0 basis.""" - Q = create_quadrature(ref_el, 2 * degree) + if interpolant_deg is None: + interpolant_deg = degree + Q = create_quadrature(ref_el, degree + interpolant_deg) B = make_bubbles(ref_el, degree, codim=codim, scale="orthonormal") P_at_qpts = B.expansion_set.tabulate(degree, Q.get_points()) - M = numpy.dot(numpy.multiply(P_at_qpts, Q.get_weights()), P_at_qpts.T) phis = numpy.linalg.solve(M, P_at_qpts) phis = numpy.dot(B.get_coeffs(), phis) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 7d9ff2968..b4cd2c63e 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -510,7 +510,7 @@ def compute_face_edge_tangents(self, dim, entity_id): v1.append(dest) return vert_coords[v1] - vert_coords[v0] - def make_points(self, dim, entity_id, order, variant=None): + def make_points(self, dim, entity_id, order, variant=None, interior=1): """Constructs a lattice of points on the entity_id:th facet of dimension dim. Order indicates how many points to include in each direction.""" @@ -520,7 +520,7 @@ def make_points(self, dim, entity_id, order, variant=None): entity_verts = \ self.get_vertices_of_subcomplex( self.get_topology()[dim][entity_id]) - return make_lattice(entity_verts, order, 1, variant=variant) + return make_lattice(entity_verts, order, interior=interior, variant=variant) else: raise ValueError("illegal dimension") diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 75c6596a0..0efea5e00 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -356,8 +356,8 @@ def __init__(self, cell, degree): class CrouzeixRaviart(ScalarFiatElement): - def __init__(self, cell, degree): - super().__init__(FIAT.CrouzeixRaviart(cell, degree)) + def __init__(self, cell, degree, variant=None): + super().__init__(FIAT.CrouzeixRaviart(cell, degree, variant=variant)) class Lagrange(ScalarFiatElement): diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index d487f83e0..efebedb59 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -94,7 +94,7 @@ def show_elements(): # Elements not in the periodic table # TODO: Implement generic Tear operator for elements instead of this: register_element("Brezzi-Douglas-Fortin-Marini", "BDFM", 1, HDiv, "contravariant Piola", (1, None), simplices[1:]) -register_element("Crouzeix-Raviart", "CR", 0, L2, "identity", (1, 1), simplices[1:]) +register_element("Crouzeix-Raviart", "CR", 0, L2, "identity", (1, None), simplices[1:]) register_element("Discontinuous Raviart-Thomas", "DRT", 1, L2, "contravariant Piola", (1, None), simplices[1:]) register_element("Kong-Mulder-Veldhuizen", "KMV", 0, H1, "identity", (1, None), simplices[1:]) diff --git a/gem/utils.py b/gem/utils.py index 12e0e0f6c..7e855f768 100644 --- a/gem/utils.py +++ b/gem/utils.py @@ -1,23 +1,5 @@ import collections - - -# This is copied from PyOP2, and it is here to be available for both -# FInAT and TSFC without depending on PyOP2. -class cached_property(object): - """A read-only @property that is only evaluated once. The value is cached - on the object itself rather than the function or class; this should prevent - memory leakage.""" - def __init__(self, fget, doc=None): - self.fget = fget - self.__doc__ = doc or fget.__doc__ - self.__name__ = fget.__name__ - self.__module__ = fget.__module__ - - def __get__(self, obj, cls): - if obj is None: - return self - obj.__dict__[self.__name__] = result = self.fget(obj) - return result +from functools import cached_property # noqa: F401 def groupby(iterable, key=None): diff --git a/test/FIAT/regression/test_regression.py b/test/FIAT/regression/test_regression.py index a4101d0cb..fc5ee01e5 100644 --- a/test/FIAT/regression/test_regression.py +++ b/test/FIAT/regression/test_regression.py @@ -281,7 +281,7 @@ def create_data(family, dim, degree): '''Create the reference data. ''' kwargs = {} - if family in {"Regge", "Hellan-Herrmann-Johnson"}: + if family in {"Crouzeix-Raviart", "Regge", "Hellan-Herrmann-Johnson"}: kwargs["variant"] = "point" # Get domain and element class domain = ufc_simplex(dim) diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index 5a1dfb21b..da97387f2 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -146,6 +146,13 @@ def __init__(self, a, b): "CrouzeixRaviart(I, 1)", "CrouzeixRaviart(T, 1)", "CrouzeixRaviart(S, 1)", + "CrouzeixRaviart(T, 3)", + "CrouzeixRaviart(T, 5)", + "CrouzeixRaviart(I, 1, variant='point')", + "CrouzeixRaviart(T, 1, variant='point')", + "CrouzeixRaviart(S, 1, variant='point')", + "CrouzeixRaviart(T, 3, variant='point')", + "CrouzeixRaviart(T, 5, variant='point')", "RaviartThomas(I, 1)", "RaviartThomas(I, 2)", "RaviartThomas(I, 3)",