diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 3583fac18..f222a70b3 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -6,7 +6,7 @@ # Import finite element classes from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401 -from FIAT.argyris import Argyris, QuinticArgyris +from FIAT.argyris import Argyris from FIAT.bernstein import Bernstein from FIAT.bell import Bell from FIAT.hct import HsiehCloughTocher @@ -102,5 +102,4 @@ "Mardal-Tai-Winther": MardalTaiWinther} # List of extra elements -extra_elements = {"P0": P0, - "Quintic Argyris": QuinticArgyris} +extra_elements = {"P0": P0} diff --git a/FIAT/argyris.py b/FIAT/argyris.py index ce2d0ec36..ba35b61d7 100644 --- a/FIAT/argyris.py +++ b/FIAT/argyris.py @@ -4,142 +4,107 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import TRIANGLE +from FIAT import finite_element, polynomial_set, dual_set +from FIAT.check_format_variant import check_format_variant +from FIAT.functional import (PointEvaluation, PointDerivative, PointNormalDerivative, + IntegralMoment, + IntegralMomentOfNormalDerivative) +from FIAT.jacobi import eval_jacobi_batch, eval_jacobi_deriv_batch +from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature +from FIAT.reference_element import TRIANGLE, ufc_simplex class ArgyrisDualSet(dual_set.DualSet): - def __init__(self, ref_el, degree): - entity_ids = {} - nodes = [] - cur = 0 - - top = ref_el.get_topology() - verts = ref_el.get_vertices() - sd = ref_el.get_spatial_dimension() - + def __init__(self, ref_el, degree, variant, interpolant_deg): if ref_el.get_shape() != TRIANGLE: raise ValueError("Argyris only defined on triangles") - pe = functional.PointEvaluation - pd = functional.PointDerivative - pnd = functional.PointNormalDerivative - - # get jet at each vertex - - entity_ids[0] = {} - for v in sorted(top[0]): - nodes.append(pe(ref_el, verts[v])) - - # first derivatives - for i in range(sd): - alpha = [0] * sd - alpha[i] = 1 - nodes.append(pd(ref_el, verts[v], alpha)) - - # second derivatives - alphas = [[2, 0], [1, 1], [0, 2]] - for alpha in alphas: - nodes.append(pd(ref_el, verts[v], alpha)) - - entity_ids[0][v] = list(range(cur, cur + 6)) - cur += 6 - - # edge dof - entity_ids[1] = {} - for e in sorted(top[1]): - # normal derivatives at degree - 4 points on each edge - ndpts = ref_el.make_points(1, e, degree - 3) - ndnds = [pnd(ref_el, e, pt) for pt in ndpts] - nodes.extend(ndnds) - entity_ids[1][e] = list(range(cur, cur + len(ndpts))) - cur += len(ndpts) - - # point value at degree-5 points on each edge - if degree > 5: - ptvalpts = ref_el.make_points(1, e, degree - 4) - ptvalnds = [pe(ref_el, pt) for pt in ptvalpts] - nodes.extend(ptvalnds) - entity_ids[1][e] += list(range(cur, cur + len(ptvalpts))) - cur += len(ptvalpts) - - # internal dof - entity_ids[2] = {} - if degree > 5: - internalpts = ref_el.make_points(2, 0, degree - 3) - internalnds = [pe(ref_el, pt) for pt in internalpts] - nodes.extend(internalnds) - entity_ids[2][0] = list(range(cur, cur + len(internalpts))) - cur += len(internalpts) - else: - entity_ids[2] = {0: []} - - super(ArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids) - - -class QuinticArgyrisDualSet(dual_set.DualSet): - 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() - if ref_el.get_shape() != TRIANGLE: - raise ValueError("Argyris only defined on triangles") - - pd = functional.PointDerivative - - # get jet at each vertex + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + nodes = [] - entity_ids[0] = {} + # get second order jet at each vertex + verts = ref_el.get_vertices() + alphas = [(1, 0), (0, 1), (2, 0), (1, 1), (0, 2)] for v in sorted(top[0]): - nodes.append(functional.PointEvaluation(ref_el, verts[v])) - - # first derivatives - for i in range(sd): - alpha = [0] * sd - alpha[i] = 1 - nodes.append(pd(ref_el, verts[v], alpha)) + cur = len(nodes) + nodes.append(PointEvaluation(ref_el, verts[v])) + nodes.extend(PointDerivative(ref_el, verts[v], alpha) for alpha in alphas) + entity_ids[0][v] = list(range(cur, len(nodes))) + + if variant == "integral": + # edge dofs + k = degree - 5 + rline = ufc_simplex(1) + Q = create_quadrature(rline, interpolant_deg+k-1) + x = 2.0 * Q.get_points() - 1.0 + phis = eval_jacobi_batch(2, 2, k, x) + dphis = eval_jacobi_deriv_batch(2, 2, k, x) + for e in sorted(top[1]): + Q_mapped = FacetQuadratureRule(ref_el, 1, e, Q) + scale = 2 / Q_mapped.jacobian_determinant() + cur = len(nodes) + nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis) + nodes.extend(IntegralMoment(ref_el, Q_mapped, dphi * scale) for dphi in dphis[1:]) + entity_ids[1][e].extend(range(cur, len(nodes))) + + # 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))) + + elif variant == "point": + # edge dofs + for e in sorted(top[1]): + cur = len(nodes) + # normal derivatives at degree - 4 points on each edge + ndpts = ref_el.make_points(1, e, degree - 3) + nodes.extend(PointNormalDerivative(ref_el, e, pt) for pt in ndpts) + + # point value at degree - 5 points on each edge + ptvalpts = ref_el.make_points(1, e, degree - 4) + nodes.extend(PointEvaluation(ref_el, pt) for pt in ptvalpts) + entity_ids[1][e] = list(range(cur, len(nodes))) - # second derivatives - alphas = [[2, 0], [1, 1], [0, 2]] - for alpha in alphas: - nodes.append(pd(ref_el, verts[v], alpha)) + # 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))) + else: + raise ValueError("Invalid variant for Argyris") + super(ArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids) - entity_ids[0][v] = list(range(cur, cur + 6)) - cur += 6 - # edge dof -- normal at each edge midpoint - entity_ids[1] = {} - for e in sorted(top[1]): - pt = ref_el.make_points(1, e, 2)[0] - n = functional.PointNormalDerivative(ref_el, e, pt) - nodes.append(n) - entity_ids[1][e] = [cur] - cur += 1 +class Argyris(finite_element.CiarletElement): + """ + The Argyris finite element. - entity_ids[2] = {0: []} + :arg ref_el: The reference element. + :arg degree: The degree. + :arg variant: optional variant specifying the types of nodes. - super(QuinticArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids) + variant can be chosen from ["point", "integral", "integral(q)"] + "point" -> dofs are evaluated by point evaluation. + "integral" -> dofs are evaluated by quadrature rules with the minimum + degree required for unisolvence. + "integral(q)" -> dofs are evaluated by quadrature rules with the minimum + degree required for unisolvence plus q. + """ + def __init__(self, ref_el, degree=5, variant=None): -class Argyris(finite_element.CiarletElement): - """The Argyris finite element.""" + variant, interpolant_deg = check_format_variant(variant, degree) - def __init__(self, ref_el, degree): - poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) - dual = ArgyrisDualSet(ref_el, degree) + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble") + dual = ArgyrisDualSet(ref_el, degree, variant, interpolant_deg) super(Argyris, self).__init__(poly_set, dual, degree) - - -class QuinticArgyris(finite_element.CiarletElement): - """The Argyris finite element.""" - - def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(ref_el, 5) - dual = QuinticArgyrisDualSet(ref_el) - super().__init__(poly_set, dual, 5) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 26c092081..bf38f9795 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -439,6 +439,39 @@ def tabulate_normal_jumps(self, n, ref_pts, facet, order=0): results[r][indices] += vr return results + def tabulate_jumps(self, n, points, order=0): + """Tabulates derivative jumps on given points. + + :arg n: the polynomial degree. + :arg points: an iterable of points on the cell complex. + :kwarg order: the order of differentiation. + + :returns: a dictionary of tabulations of derivative jumps across interior facets. + """ + + from FIAT.polynomial_set import mis + num_members = self.get_num_members(n) + cell_node_map = self.get_cell_node_map(n) + + top = self.ref_el.get_topology() + sd = self.ref_el.get_spatial_dimension() + interior_facets = self.ref_el.get_interior_facets(sd-1) + derivs = {cell: self._tabulate_on_cell(n, points, order=order, cell=cell) + for cell in top[sd]} + jumps = {} + for r in range(order+1): + cur = 0 + alphas = mis(sd, r) + jumps[r] = numpy.zeros((num_members, len(alphas)*len(interior_facets), len(points))) + for facet in interior_facets: + facet_verts = set(top[sd-1][facet]) + c0, c1 = tuple(k for k in top[sd] if facet_verts < set(top[sd][k])) + for alpha in alphas: + jumps[r][cell_node_map[c1], cur] += derivs[c1][alpha] + jumps[r][cell_node_map[c0], cur] -= derivs[c0][alpha] + cur = cur + 1 + return jumps + def get_dmats(self, degree, cell=0): """Returns a numpy array with the expansion coefficients dmat[k, j, i] of the gradient of each member of the expansion set: diff --git a/FIAT/hct.py b/FIAT/hct.py index 00dde6791..8b2b5102c 100644 --- a/FIAT/hct.py +++ b/FIAT/hct.py @@ -1,15 +1,28 @@ +# Copyright (C) 2024 Pablo D. Brubeck +# +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 + from FIAT.functional import (PointEvaluation, PointDerivative, + IntegralMoment, IntegralMomentOfNormalDerivative) from FIAT import finite_element, dual_set, macro, polynomial_set from FIAT.reference_element import TRIANGLE, ufc_simplex +from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -from FIAT.jacobi import eval_jacobi +from FIAT.jacobi import eval_jacobi, eval_jacobi_batch, eval_jacobi_deriv_batch class HCTDualSet(dual_set.DualSet): - def __init__(self, ref_el, degree, reduced=False): - if degree != 3: - raise ValueError("HCT only defined for degree=3") + def __init__(self, ref_complex, degree, reduced=False): + if reduced and degree != 3: + raise ValueError("Reduced HCT only defined for degree = 3") + if degree < 3: + raise ValueError("HCT only defined for degree >= 3") + ref_el = ref_complex.get_parent() if ref_el.get_shape() != TRIANGLE: raise ValueError("HCT only defined on triangles") top = ref_el.get_topology() @@ -27,24 +40,49 @@ def __init__(self, ref_el, degree, reduced=False): nodes.extend(PointDerivative(ref_el, pt, alpha) for alpha in alphas) entity_ids[0][v].extend(range(cur, len(nodes))) + k = 2 if reduced else degree - 3 rline = ufc_simplex(1) - k = 2 if reduced else 0 Q = create_quadrature(rline, degree-1+k) qpts = Q.get_points() - x, = qpts.T - f_at_qpts = eval_jacobi(0, 0, k, 2.0*x - 1) - for e in sorted(top[1]): - cur = len(nodes) - nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q, f_at_qpts)) - entity_ids[1][e].extend(range(cur, len(nodes))) + if reduced: + x, = qpts.T + f_at_qpts = eval_jacobi(0, 0, k, 2.0*x - 1) + for e in sorted(top[1]): + cur = len(nodes) + nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q, f_at_qpts)) + entity_ids[1][e].extend(range(cur, len(nodes))) + else: + x = 2.0*qpts - 1 + phis = eval_jacobi_batch(2, 2, k, x) + dphis = eval_jacobi_deriv_batch(2, 2, k, x) + for e in sorted(top[1]): + Q_mapped = FacetQuadratureRule(ref_el, 1, e, Q) + scale = 2 / Q_mapped.jacobian_determinant() + cur = len(nodes) + nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis) + nodes.extend(IntegralMoment(ref_el, Q_mapped, dphi * scale) for dphi in dphis[1:]) + entity_ids[1][e].extend(range(cur, len(nodes))) + + q = degree - 4 + if q >= 0: + Q = create_quadrature(ref_complex, degree + q) + Pq = polynomial_set.ONPolynomialSet(ref_el, q, scale=1) + phis = Pq.tabulate(Q.get_points())[(0,) * sd] + scale = 1 / 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))) super(HCTDualSet, self).__init__(nodes, ref_el, entity_ids) class HsiehCloughTocher(finite_element.CiarletElement): - """The HCT finite element.""" - + """The HCT macroelement. For degree higher than 3, we implement the + super-smooth C^1 space from Groselj and Knez (2022) on a barycentric split, + although there the basis functions are positive on an incenter split. + """ def __init__(self, ref_el, degree=3, reduced=False): - dual = HCTDualSet(ref_el, degree, reduced=reduced) - poly_set = macro.CkPolynomialSet(macro.AlfeldSplit(ref_el), degree, variant=None) + ref_complex = macro.AlfeldSplit(ref_el) + dual = HCTDualSet(ref_complex, degree, reduced=reduced) + poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=degree-1, variant="bubble") super(HsiehCloughTocher, self).__init__(poly_set, dual, degree) diff --git a/FIAT/macro.py b/FIAT/macro.py index 84df2cede..968c316f2 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -360,11 +360,12 @@ class CkPolynomialSet(polynomial_set.PolynomialSet): :arg ref_el: The simplicial complex. :arg degree: The polynomial degree. :kwarg order: The order of continuity across subcells. + :kwarg vorder: The order of super-smoothness at interior vertices. :kwarg shape: The value shape. :kwarg variant: The variant for the underlying ExpansionSet. :kwarg scale: The scale for the underlying ExpansionSet. """ - def __init__(self, ref_el, degree, order=1, shape=(), **kwargs): + def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): from FIAT.quadrature_schemes import create_quadrature expansion_set = expansions.ExpansionSet(ref_el, **kwargs) k = 1 if expansion_set.continuity == "C0" else 0 @@ -390,11 +391,25 @@ def __init__(self, ref_el, degree, order=1, shape=(), **kwargs): if len(rows) > 0: dual_mat = numpy.vstack(rows) _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) + tol = sig[0] * 1E-10 + num_sv = len([s for s in sig if abs(s) > tol]) coeffs = vt[num_sv:] else: coeffs = numpy.eye(expansion_set.get_num_members(degree)) + if vorder > order + 1: + # Impose C^vorder super-smoothness at interior vertices + # C^order automatically gives C^{order+1} at the interior vertex + verts = ref_el.get_vertices() + points = [verts[i] for i in ref_el.get_interior_facets(0)] + jumps = expansion_set.tabulate_jumps(degree, points, order=vorder) + for r in range(order+2, vorder+1): + dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T) + _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) + tol = sig[0] * 1E-10 + num_sv = len([s for s in sig if abs(s) > tol]) + coeffs = numpy.dot(vt[num_sv:], coeffs) + if shape != tuple(): m, n = coeffs.shape coeffs = coeffs.reshape((m,) + (1,)*len(shape) + (n,)) diff --git a/test/unit/test_argyris.py b/test/unit/test_argyris.py new file mode 100644 index 000000000..e41192179 --- /dev/null +++ b/test/unit/test_argyris.py @@ -0,0 +1,72 @@ +import pytest +import numpy + +from FIAT import Argyris +from FIAT.jacobi import eval_jacobi_batch, eval_jacobi_deriv_batch +from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature +from FIAT.reference_element import ufc_simplex + + +@pytest.fixture +def cell(): + return ufc_simplex(2) + + +@pytest.mark.parametrize("degree", range(6, 10)) +def test_argyris_basis_functions(cell, degree): + fe = Argyris(cell, degree) + + entity_ids = fe.entity_dofs() + sd = cell.get_spatial_dimension() + top = cell.get_topology() + rline = ufc_simplex(1) + Qref = create_quadrature(rline, 2*degree) + + q = degree - 5 + xref = 2.0*Qref.get_points()-1 + ell_at_qpts = eval_jacobi_batch(0, 0, degree, xref) + P_at_qpts = eval_jacobi_batch(2, 2, degree, xref) + DP_at_qpts = eval_jacobi_deriv_batch(2, 2, degree, xref) + + for e in top[1]: + n = cell.compute_normal(e) + t = cell.compute_edge_tangent(e) + Q = FacetQuadratureRule(cell, 1, e, Qref) + qpts, qwts = Q.get_points(), Q.get_weights() + + ids = entity_ids[1][e] + ids1 = ids[1:-q] + ids0 = ids[-q:] + + tab = fe.tabulate(1, qpts) + # Test that normal derivative moment bfs have vanishing trace + assert numpy.allclose(tab[(0,) * sd][ids[:-q]], 0) + + phi = tab[(0,) * sd][ids0] + + phi_n = sum(n[alpha.index(1)] * tab[alpha][ids1] + for alpha in tab if sum(alpha) == 1) + + phi_t = sum(t[alpha.index(1)] * tab[alpha][ids0] + for alpha in tab if sum(alpha) == 1) + + # Test that facet bfs are hierarchical + coeffs = numpy.dot(numpy.multiply(phi, qwts), ell_at_qpts[6:].T) + assert numpy.allclose(numpy.triu(coeffs, k=1), 0) + + # Test duality of normal derivarive moments + coeffs = numpy.dot(numpy.multiply(phi_n, qwts), P_at_qpts[1:].T) + assert numpy.allclose(coeffs[:, q:], 0) + assert numpy.allclose(coeffs[:, :q], numpy.diag(numpy.diag(coeffs[:, :q]))) + + # Test duality of trace moments + coeffs = numpy.dot(numpy.multiply(phi, qwts), DP_at_qpts[1:].T) + assert numpy.allclose(coeffs[:, q:], 0) + assert numpy.allclose(coeffs[:, :q], numpy.diag(numpy.diag(coeffs[:, :q]))) + + # Test the integration by parts property arising from the choice + # of normal derivative and trace moments DOFs. + # The normal derivative of the normal derviative bfs must be equal + # to minus the tangential derivative of the trace moment bfs + assert numpy.allclose(numpy.dot((phi_n + phi_t)**2, qwts), 0) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 65abfaef9..bed985f5d 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -42,7 +42,7 @@ from FIAT.tensor_product import TensorProductElement # noqa: F401 from FIAT.tensor_product import FlattenedDimensions # noqa: F401 from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 -from FIAT.argyris import Argyris, QuinticArgyris # noqa: F401 +from FIAT.argyris import Argyris # noqa: F401 from FIAT.hermite import CubicHermite # noqa: F401 from FIAT.morley import Morley # noqa: F401 from FIAT.bubble import Bubble @@ -291,8 +291,9 @@ def __init__(self, a, b): " Regge(S, 1)," " RestrictedElement(Regge(S, 2), restriction_domain='interior')" ")", - "Argyris(T, 5)", - "QuinticArgyris(T)", + "Argyris(T, 5, 'point')", + "Argyris(T, 5, 'integral')", + "Argyris(T, 6, 'integral')", "CubicHermite(I)", "CubicHermite(T)", "CubicHermite(S)", diff --git a/test/unit/test_pointwise_dual.py b/test/unit/test_pointwise_dual.py index 4d2698408..6e2325571 100644 --- a/test/unit/test_pointwise_dual.py +++ b/test/unit/test_pointwise_dual.py @@ -8,7 +8,7 @@ import numpy from FIAT import ( - BrezziDouglasMarini, Morley, QuinticArgyris, CubicHermite) + BrezziDouglasMarini, Morley, Argyris, CubicHermite) from FIAT.reference_element import ( UFCTriangle, @@ -22,7 +22,7 @@ @pytest.mark.parametrize("element", [CubicHermite(T), Morley(T), - QuinticArgyris(T), + Argyris(T), BrezziDouglasMarini(T, 1, variant="integral")]) def test_pw_dual(element): deg = element.degree()