merge release

Author: pbrubeck

.github/workflows/pypi.yml

@@ -1,38 +1,28 @@
 name: Publish on PyPI
 
-# By default this action does not push to test or production PyPI. The wheels
-# are available as an artifact that can be downloaded and tested locally.
+# By default this action does not push to PyPI. The wheels are available
+# as an artifact that can be downloaded and tested locally.
 
 on:
   workflow_dispatch:
     inputs:
       ref:
         description: Git ref to publish
-        default: master
         type: string
       pypi:
         description: Publish to PyPI
         default: false
         type: boolean
-      testpypi:
-        description: Publish to TestPyPI
-        default: false
-        type: boolean
 
   workflow_call:
     inputs:
       ref:
         description: Git ref to publish
-        default: master
         type: string
       pypi:
         description: Publish to PyPI
         default: false
         type: boolean
-      testpypi:
-        description: Publish to TestPyPI
-        default: false
-        type: boolean
 
 jobs:
   build:
@@ -65,8 +55,6 @@ jobs:
     runs-on: ubuntu-latest
     environment:
       name: pypi
-    permissions:
-      id-token: write
     steps:
       - uses: actions/download-artifact@v4
         with:
@@ -75,23 +63,5 @@ jobs:
 
       - name: Push to PyPI
         uses: pypa/gh-action-pypi-publish@release/v1
-
-  upload_test_pypi:
-    name: Upload to TestPyPI (optional)
-    if: inputs.testpypi
-    needs: build
-    runs-on: ubuntu-latest
-    environment:
-      name: testpypi
-    permissions:
-      id-token: write
-    steps:
-      - uses: actions/download-artifact@v4
-        with:
-          name: dist
-          path: dist
-
-      - name: Push to TestPyPI
-        uses: pypa/gh-action-pypi-publish@release/v1
         with:
-          repository-url: https://test.pypi.org/legacy/
+          password: ${{ secrets.PYPI_API_TOKEN }}

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)

FIAT/brezzi_douglas_marini.py

@@ -10,6 +10,7 @@
 from FIAT.check_format_variant import check_format_variant
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.quadrature import FacetQuadratureRule
+import numpy
 
 
 class BDMDualSet(dual_set.DualSet):
@@ -26,7 +27,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 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)
@@ -56,14 +57,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
         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)))
+
+            cell = ref_el.construct_subelement(sd)
+            Q_ref = create_quadrature(cell, interpolant_deg + degree - 1)
+            Nedel = nedelec.Nedelec(cell, degree - 1, variant)
+            Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd]
+            for entity in top[sd]:
+                Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
+                Jinv = numpy.linalg.inv(Q.jacobian())
+                phis = numpy.tensordot(Jinv.T, Ned_at_qpts, (1, 1)).transpose((1, 0, 2))
+                cur = len(nodes)
+                nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis)
+                entity_ids[sd][entity] = list(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
@@ -89,6 +94,11 @@ class BrezziDouglasMarini(finite_element.CiarletElement):
     """
 
     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)
+
         if degree < 1:
             raise ValueError(f"{type(self).__name__} elements only valid for k >= 1")
 
@@ -99,7 +109,6 @@ def __init__(self, ref_el, degree, variant=None):
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "HDiv", type(self))
         else:
-            variant, interpolant_deg = check_format_variant(variant, degree)
             dual = BDMDualSet(ref_el, degree, variant, interpolant_deg)
         formdegree = sd - 1  # (n-1)-form
         super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")

FIAT/bubble.py

@@ -7,15 +7,20 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 from FIAT.lagrange import Lagrange
+from FIAT.hierarchical import IntegratedLegendre
 from FIAT.restricted import RestrictedElement
 from itertools import chain
 
 
 class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
-    def __init__(self, ref_el, degree, codim):
-        element = Lagrange(ref_el, degree)
+    def __init__(self, ref_el, degree, codim, variant=None):
+        if variant and variant.startswith("integral"):
+            CG = IntegratedLegendre
+        else:
+            CG = Lagrange
+        element = CG(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
         assert cell_dim == max(element.entity_dofs().keys())
@@ -29,12 +34,12 @@ def __init__(self, ref_el, degree, codim):
 class Bubble(CodimBubble):
     """The bubble finite element: the dofs of the Lagrange FE in the interior of the cell"""
 
-    def __init__(self, ref_el, degree):
-        super().__init__(ref_el, degree, codim=0)
+    def __init__(self, ref_el, degree, variant=None):
+        super().__init__(ref_el, degree, codim=0, variant=variant)
 
 
 class FacetBubble(CodimBubble):
     """The facet bubble finite element: the dofs of the Lagrange FE in the interior of the facets"""
 
-    def __init__(self, ref_el, degree):
-        super().__init__(ref_el, degree, codim=1)
+    def __init__(self, ref_el, degree, variant=None):
+        super().__init__(ref_el, degree, codim=1, variant=variant)

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"
 
@@ -43,7 +45,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):
@@ -60,7 +62,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False):
 
     default = "integral" if integral else "spectral"
     if integral:
-        supported_point_variants = {"integral": None, "fdm": "fdm", "demkowicz": "demkowicz", "demkowiczmass": "demkowiczmass"}
+        supported_point_variants = {"integral": None, "point": "point", "fdm": "fdm", "demkowicz": "demkowicz", "demkowiczmass": "demkowiczmass"}
     elif discontinuous:
         supported_point_variants = supported_dg_variants
     else:
@@ -79,6 +81,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:

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)

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

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.")