Merge branch 'master' into pbrubeck/fix/gll-ordering

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}

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)

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:

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 ([email protected]), 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)