Add fuse element infrastructure

FIAT/lagrange.py

@@ -40,7 +40,7 @@ def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=Fal
         entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])]
         if sort_entities:
             # sort the entities by support vertex ids
-            support = [top[dim][entity] for dim, entity in entities]
+            support = [sorted(top[dim][entity]) for dim, entity in entities]
             entities = [entity for verts, entity in sorted(zip(support, entities))]
 
         # make nodes by getting points

FIAT/nedelec.py

@@ -22,7 +22,7 @@ def NedelecSpace2D(ref_el, degree):
         raise Exception("NedelecSpace2D requires 2d reference element")
 
     k = degree - 1
-    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,))
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal")
 
     dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1)
     dimPk = expansions.polynomial_dimension(ref_el, k)
@@ -32,7 +32,7 @@ def NedelecSpace2D(ref_el, degree):
                                   for i in range(sd))))
     vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
 
-    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1)
+    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal")
     PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
 
     Q = create_quadrature(ref_el, 2 * (k + 1))
@@ -43,8 +43,8 @@ def NedelecSpace2D(ref_el, degree):
 
     CrossX = numpy.dot(numpy.array([[0.0, 1.0], [-1.0, 0.0]]), Qpts.T)
     PkHCrossX_at_Qpts = PkH_at_Qpts[:, None, :] * CrossX[None, :, :]
-    PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), Pkp1_at_Qpts.T)
-
+    PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts),
+                                 Pkp1_at_Qpts.T)
     PkHcrossX = polynomial_set.PolynomialSet(ref_el,
                                              k + 1,
                                              k + 1,

FIAT/polynomial_set.py

@@ -176,15 +176,49 @@ def polynomial_set_union_normalized(A, B):
     not contain any of the same members of the set, as we construct a
     span via SVD.
     """
-    new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0)
+    assert A.get_reference_element() == B.get_reference_element()
+    new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B)
+
+    deg = max(A.get_degree(), B.get_degree())
+    em_deg = max(A.get_embedded_degree(), B.get_embedded_degree())
     coeffs = spanning_basis(new_coeffs)
     return PolynomialSet(A.get_reference_element(),
-                         A.get_degree(),
-                         A.get_embedded_degree(),
+                         deg,
+                         em_deg,
                          A.get_expansion_set(),
                          coeffs)
 
 
+def construct_new_coeffs(ref_el, A, B):
+    """
+    Constructs new coefficients for the set union of A and B
+    If A and B are discontinuous and do not have the same degree the smaller one
+    is extended to match the larger.
+
+    This does not handle the case that A and B have continuity but not the same degree.
+    """
+
+    if A.get_expansion_set().continuity != B.get_expansion_set().continuity:
+        raise ValueError("Continuity of expansion sets does not match.")
+
+    if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None:
+        higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B
+        lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A
+
+        diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1]
+
+        # pad only the 0th axis with the difference in size
+        padding = [(0, 0) for i in range(len(lower.coeffs.shape) - 1)] + [(0, diff)]
+        embedded_coeffs = numpy.pad(lower.coeffs, padding)
+
+        new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0)
+    elif A.get_embedded_degree() == B.get_embedded_degree():
+        new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0)
+    else:
+        raise NotImplementedError("Extending of coefficients is not implemented for PolynomialSets with continuity and different degrees")
+    return new_coeffs
+
+
 class ONSymTensorPolynomialSet(PolynomialSet):
     """Constructs an orthonormal basis for symmetric-tensor-valued
     polynomials on a reference element.

FIAT/raviart_thomas.py

@@ -20,7 +20,7 @@ def RTSpace(ref_el, degree):
     sd = ref_el.get_spatial_dimension()
 
     k = degree - 1
-    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,))
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal")
 
     dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1)
     dimPk = expansions.polynomial_dimension(ref_el, k)
@@ -30,7 +30,7 @@ def RTSpace(ref_el, degree):
                                   for i in range(sd))))
     vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
 
-    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1)
+    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal")
     PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
 
     Q = create_quadrature(ref_el, 2 * (k + 1))
@@ -39,12 +39,12 @@ def RTSpace(ref_el, degree):
     # have to work on this through "tabulate" interface
     # first, tabulate PkH at quadrature points
     PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd]
+
     Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd]
 
     x = Qpts.T
     PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :]
     PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T)
-
     PkHx = polynomial_set.PolynomialSet(ref_el,
                                         k,
                                         k + 1,

FIAT/reference_element.py

@@ -515,7 +515,8 @@ def make_points(self, dim, entity_id, order, variant=None):
         facet of dimension dim.  Order indicates how many points to
         include in each direction."""
         if dim == 0:
-            return (self.get_vertices()[entity_id], )
+            return (self.get_vertices()[self.get_topology()[dim][entity_id][0]],)
+            # return (self.get_vertices()[entity_id], )
         elif 0 < dim <= self.get_spatial_dimension():
             entity_verts = \
                 self.get_vertices_of_subcomplex(

finat/__init__.py

@@ -7,7 +7,8 @@
                             KongMulderVeldhuizen, Lagrange, Real, Serendipity,               # noqa: F401
                             TrimmedSerendipityCurl, TrimmedSerendipityDiv,                   # noqa: F401
                             TrimmedSerendipityEdge, TrimmedSerendipityFace,                  # noqa: F401
-                            Nedelec, NedelecSecondKind, RaviartThomas, Regge)                # noqa: F401
+                            Nedelec, NedelecSecondKind, RaviartThomas, Regge,                # noqa: F401
+                            FuseElement)                                                     # noqa: F401
 
 from .argyris import Argyris                                       # noqa: F401
 from .aw import ArnoldWinther, ArnoldWintherNC                     # noqa: F401

finat/element_factory.py

@@ -112,7 +112,11 @@ def as_fiat_cell(cell):
     :arg cell: the :class:`ufl.Cell` to convert."""
     if not isinstance(cell, ufl.AbstractCell):
         raise ValueError("Expecting a UFL Cell")
-    return ufc_cell(cell)
+
+    try:
+        return cell.to_fiat()
+    except AttributeError:
+        return ufc_cell(cell)
 
 
 @singledispatch
@@ -302,6 +306,11 @@ def convert_restrictedelement(element, **kwargs):
     return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps
 
 
+@convert.register(finat.ufl.FuseElement)
+def convert_fuse_element(element, **kwargs):
+    return finat.fiat_elements.FuseElement(element.triple), set()
+
+
 hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
 quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
 _cache = weakref.WeakKeyDictionary()

finat/fiat_elements.py

@@ -480,3 +480,8 @@ def __init__(self, cell, degree, variant=None):
 class NedelecSecondKind(VectorFiatElement):
     def __init__(self, cell, degree, variant=None):
         super().__init__(FIAT.NedelecSecondKind(cell, degree, variant=variant))
+
+
+class FuseElement(FiatElement):
+    def __init__(self, triple):
+        super(FuseElement, self).__init__(triple.to_fiat())

finat/ufl/__init__.py

@@ -19,3 +19,4 @@
 from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement  # noqa: F401
 from finat.ufl.restrictedelement import RestrictedElement  # noqa: F401
 from finat.ufl.tensorproductelement import TensorProductElement  # noqa: F401
+from finat.ufl.fuseelement import FuseElement  # noqa: F401

finat/ufl/fuseelement.py

@@ -0,0 +1,45 @@
+"""Element."""
+# -*- coding: utf-8 -*-
+# Copyright (C) 2025 India Marsden
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+from finat.ufl.finiteelementbase import FiniteElementBase
+
+
+class FuseElement(FiniteElementBase):
+    """
+    A finite element defined using FUSE.
+
+    TODO: Need to deal with cases where value shape and reference value shape are different
+    """
+
+    def __init__(self, triple, cell=None):
+        self.triple = triple
+        if not cell:
+            cell = self.triple.cell.to_ufl()
+
+        # this isn't really correct
+        degree = self.triple.spaces[0].degree()
+        super(FuseElement, self).__init__("IT", cell, degree, None, triple.get_value_shape())
+
+    def __repr__(self):
+        return "FiniteElement(%s, %s, (%s, %s, %s), %s)" % (
+               repr(self.triple.DOFGenerator), repr(self.triple.cell), repr(self.triple.spaces[0]), repr(self.triple.spaces[1]), repr(self.triple.spaces[2]), "X")
+
+    def __str__(self):
+        return "<Fuse%sElem on %s>" % (self.triple.spaces[0], self.triple.cell)
+
+    def mapping(self):
+        if str(self.sobolev_space) == "HCurl":
+            return "covariant Piola"
+        elif str(self.sobolev_space) == "HDiv":
+            return "contravariant Piola"
+        else:
+            return "identity"
+
+    def sobolev_space(self):
+        return self.triple.spaces[1]
+
+    def reconstruct(self, family=None, cell=None, degree=None, quad_scheme=None, variant=None):
+        return FuseElement(self.triple, cell=cell)

test/FIAT/unit/test_polynomial.py

@@ -21,6 +21,7 @@
 
 from FIAT import expansions, polynomial_set, reference_element
 from FIAT.quadrature_schemes import create_quadrature
+from itertools import chain
 
 
 @pytest.fixture(params=(1, 2, 3))
@@ -107,3 +108,45 @@ def test_bubble_duality(cell, degree):
     results = scale * numpy.dot(numpy.multiply(phi_dual, qwts), phi.T)
     assert numpy.allclose(results, numpy.diag(numpy.diag(results)))
     assert numpy.allclose(numpy.diag(results), 1.0)
+
+
+@pytest.mark.parametrize("degree", [10])
+def test_union_of_polysets(cell, degree):
+    """ demonstrates that polysets don't need to have the same degree for union
+    using RT space as an example"""
+
+    sd = cell.get_spatial_dimension()
+    k = degree
+    vecPk = polynomial_set.ONPolynomialSet(cell, degree, (sd,))
+
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, (sd,), scale="orthonormal")
+
+    dimPkp1 = expansions.polynomial_dimension(cell, k + 1)
+    dimPk = expansions.polynomial_dimension(cell, k)
+    dimPkm1 = expansions.polynomial_dimension(cell, k - 1)
+
+    vec_Pk_indices = list(chain(*(range(i * dimPkp1, i * dimPkp1 + dimPk)
+                                  for i in range(sd))))
+    vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
+
+    Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, scale="orthonormal")
+    PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
+
+    Q = create_quadrature(cell, 2 * (k + 1))
+    Qpts, Qwts = Q.get_points(), Q.get_weights()
+
+    PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd]
+    Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd]
+    x = Qpts.T
+    PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :]
+    PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T)
+    PkHx = polynomial_set.PolynomialSet(cell, k, k + 1, vec_Pkp1.get_expansion_set(), PkHx_coeffs)
+
+    same_deg = polynomial_set.polynomial_set_union_normalized(vec_Pk_from_Pkp1, PkHx)
+    different_deg = polynomial_set.polynomial_set_union_normalized(vecPk, PkHx)
+
+    Q = create_quadrature(cell, 2*(degree))
+    Qpts, _ = Q.get_points(), Q.get_weights()
+    same_vals = same_deg.tabulate(Qpts)[(0,) * sd]
+    diff_vals = different_deg.tabulate(Qpts)[(0,) * sd]
+    assert numpy.allclose(same_vals - diff_vals, 0)