Changes to FIAT to enable FUSE

FIAT/discontinuous_pc.py

@@ -7,29 +7,17 @@
 # Modified by David A. Ham (david.ham@imperial.ac.uk), 2018
 
 from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.reference_element import (Point,
-                                    DefaultLine,
-                                    UFCInterval,
-                                    UFCQuadrilateral,
-                                    UFCHexahedron,
-                                    UFCTriangle,
-                                    UFCTetrahedron,
-                                    make_affine_mapping,
-                                    flatten_reference_cube)
+from FIAT.reference_element import (make_affine_mapping,
+                                    flatten_reference_cube,
+                                    cell_to_simplex)
 from FIAT.P0 import P0Dual
 import numpy as np
 
-hypercube_simplex_map = {Point(): Point(),
-                         DefaultLine(): DefaultLine(),
-                         UFCInterval(): UFCInterval(),
-                         UFCQuadrilateral(): UFCTriangle(),
-                         UFCHexahedron(): UFCTetrahedron()}
-
 
 class DPC0(finite_element.CiarletElement):
     def __init__(self, ref_el):
         flat_el = flatten_reference_cube(ref_el)
-        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], 0)
+        poly_set = polynomial_set.ONPolynomialSet(cell_to_simplex(flat_el), 0)
         dual = P0Dual(ref_el)
         # Implement entity_permutations when we handle that for HigherOrderDPC.
         # Currently, orientation_tuples in P0Dual(ref_el).entity_permutations
@@ -58,7 +46,7 @@ def __init__(self, ref_el, flat_el, degree):
 
         # Change coordinates here.
         # Vertices of the simplex corresponding to the reference element.
-        v_simplex = hypercube_simplex_map[flat_el].get_vertices()
+        v_simplex = cell_to_simplex(flat_el).get_vertices()
         # Vertices of the reference element.
         v_hypercube = flat_el.get_vertices()
         # For the mapping, first two vertices are unchanged in all dimensions.
@@ -74,12 +62,12 @@ def __init__(self, ref_el, flat_el, degree):
 
         # make nodes by getting points
         # need to do this dimension-by-dimension, facet-by-facet
-        top = hypercube_simplex_map[flat_el].get_topology()
+        top = cell_to_simplex(flat_el).get_topology()
 
         cur = 0
         for dim in sorted(top):
             for entity in sorted(top[dim]):
-                pts_cur = hypercube_simplex_map[flat_el].make_points(dim, entity, degree)
+                pts_cur = cell_to_simplex(flat_el).make_points(dim, entity, degree)
                 pts_cur = [tuple(np.matmul(A, np.array(x)) + b) for x in pts_cur]
                 nodes_cur = [functional.PointEvaluation(flat_el, x)
                              for x in pts_cur]
@@ -102,7 +90,7 @@ class HigherOrderDPC(finite_element.CiarletElement):
 
     def __init__(self, ref_el, degree):
         flat_el = flatten_reference_cube(ref_el)
-        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], degree)
+        poly_set = polynomial_set.ONPolynomialSet(cell_to_simplex(flat_el), degree)
         dual = DPCDualSet(ref_el, flat_el, degree)
         formdegree = flat_el.get_spatial_dimension()  # n-form
         super().__init__(poly_set=poly_set,

FIAT/orientation_utils.py

@@ -97,7 +97,7 @@ def _make_axis_perms_tensorproduct(cells, dim):
     This is the single sources of extrinsic orientations and corresponding axis permutations.
 
     """
-    from FIAT.reference_element import UFCInterval
+    from FIAT.reference_element import LINE
 
     # Handle extrinsic orientations.
     # This is complex and we need to think to make this function more general.
@@ -116,8 +116,8 @@ def _make_axis_perms_tensorproduct(cells, dim):
         #          dim == (2, 1) ->
         #          triangle x interval (1 possible extrinsic orientation).
         axis_perms = (tuple(range(nprod)), )  # Identity: no permutations
-    elif len(set(cells)) == 1 and isinstance(cells[0], UFCInterval):
-        # Tensor product of intervals.
+    elif len(set(cells)) == 1 and cells[0].get_shape() == LINE:
+        # Tensor product of intervals (any 1D reference cell implementation)
         # Example: interval x interval x interval x interval
         #          dim == (0, 1, 1, 1) ->
         #          point x interval x interval x interval  (1! * 3! possible extrinsic orientations).

FIAT/reference_element.py

@@ -133,7 +133,7 @@ class Cell:
     """Abstract class for a reference cell.  Provides accessors for
     geometry (vertex coordinates) as well as topology (orderings of
     vertices that make up edges, faces, etc."""
-    def __init__(self, shape, vertices, topology):
+    def __init__(self, shape, vertices, topology, sub_entities=None):
         """The constructor takes a shape code, the physical vertices expressed
         as a list of tuples of numbers, and the topology of a cell.
 
@@ -145,24 +145,25 @@ def __init__(self, shape, vertices, topology):
         self.vertices = vertices
         self.topology = topology
 
-        # Given the topology, work out for each entity in the cell,
-        # which other entities it contains.
-        self.sub_entities = {}
-        for dim, entities in topology.items():
-            self.sub_entities[dim] = {}
-
-            for e, v in entities.items():
-                vertices = frozenset(v)
-                sub_entities = []
-
-                for dim_, entities_ in topology.items():
-                    for e_, vertices_ in entities_.items():
-                        if vertices.issuperset(vertices_):
-                            sub_entities.append((dim_, e_))
-
-                # Sort for the sake of determinism and by UFC conventions
-                self.sub_entities[dim][e] = sorted(sub_entities)
-
+        if sub_entities is not None:
+            self.sub_entities = sub_entities
+        else:
+            # If sub entity list not provided
+            # Given the topology, work out for each entity in the cell,
+            # which other entities it contains.
+            self.sub_entities = {}
+            for dim, entities in topology.items():
+                self.sub_entities[dim] = {}
+
+                for e, v in entities.items():
+                    vertices = frozenset(v)
+                    sub_entities = []
+                    for dim_, entities_ in topology.items():
+                        for e_, vertices_ in entities_.items():
+                            if vertices.issuperset(vertices_):
+                                sub_entities.append((dim_, e_))
+
+                    self.sub_entities[dim][e] = sorted(list(sub_entities))
         # Build super-entity dictionary by inverting the sub-entity dictionary
         self.super_entities = {dim: {entity: [] for entity in topology[dim]} for dim in topology}
         for dim0 in topology:
@@ -183,7 +184,6 @@ def __init__(self, shape, vertices, topology):
                     neighbors = children if dim1 < dim0 else parents
                     d01_entities = tuple(e for d, e in neighbors if d == dim1)
                     self.connectivity[(dim0, dim1)].append(d01_entities)
-
         # Dictionary with derived cells
         self._split_cache = {}
 
@@ -387,14 +387,14 @@ class SimplicialComplex(Cell):
 
     This consists of list of vertex locations and a topology map defining facets.
     """
-    def __init__(self, shape, vertices, topology):
+    def __init__(self, shape, vertices, topology, sub_entities=None):
         # Make sure that every facet has the right number of vertices to be
         # a simplex.
         for dim in topology:
             for entity in topology[dim]:
                 assert len(topology[dim][entity]) == dim + 1
 
-        super().__init__(shape, vertices, topology)
+        super().__init__(shape, vertices, topology, sub_entities)
 
     def compute_normal(self, facet_i, cell=None):
         """Returns the unit normal vector to facet i of codimension 1."""
@@ -528,7 +528,7 @@ def make_points(self, dim, entity_id, order, variant=None, interior=1):
         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]],)
         elif 0 < dim <= self.get_spatial_dimension():
             entity_verts = \
                 self.get_vertices_of_subcomplex(
@@ -1155,7 +1155,7 @@ def compute_normal(self, i):
 class TensorProductCell(Cell):
     """A cell that is the product of FIAT cells."""
 
-    def __init__(self, *cells):
+    def __init__(self, *cells, sub_entities=None):
         # Vertices
         vertices = tuple(tuple(chain(*coords))
                          for coords in product(*[cell.get_vertices()
@@ -1178,7 +1178,7 @@ def __init__(self, *cells):
             topology[dim] = dict(enumerate(topology[dim][key]
                                            for key in sorted(topology[dim])))
 
-        super().__init__(TENSORPRODUCT, vertices, topology)
+        super().__init__(TENSORPRODUCT, vertices, topology, sub_entities)
         self.cells = tuple(cells)
 
     def __repr__(self):
@@ -1421,15 +1421,15 @@ def is_macrocell(self):
 class Hypercube(Cell):
     """Abstract class for a reference hypercube"""
 
-    def __init__(self, dimension, product):
+    def __init__(self, dimension, product, sub_entities=None):
         self.dimension = dimension
         self.shape = hypercube_shapes[dimension]
 
         pt = product.get_topology()
         verts = product.get_vertices()
         topology = flatten_entities(pt)
 
-        super().__init__(self.shape, verts, topology)
+        super().__init__(self.shape, verts, topology, sub_entities)
 
         self.product = product
         self.unflattening_map = compute_unflattening_map(pt)
@@ -1872,3 +1872,10 @@ def max_complex(complexes):
         return max_cell
     else:
         raise ValueError("Cannot find the maximal complex")
+
+
+def cell_to_simplex(cell):
+    if cell.is_simplex():
+        return cell
+    else:
+        return ufc_simplex(cell.get_dimension())

finat/__init__.py

@@ -8,7 +8,8 @@
                             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 .spectral import (GaussLobattoLegendre, GaussLegendre, KongMulderVeldhuizen,            # noqa: F401
                        Legendre, IntegratedLegendre,                                         # noqa: F401
                        FDMLagrange, FDMQuadrature, FDMDiscontinuousLagrange,                 # noqa: F401

finat/element_factory.py

@@ -109,15 +109,24 @@
 element is supported, but must be handled specially because it doesn't
 have a direct FInAT equivalent."""
 
+hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
+quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
+
 
 @cache
 def as_fiat_cell(cell):
     """Convert a ufl cell to a FIAT cell.
 
     :arg cell: the :class:`ufl.Cell` to convert."""
+    if isinstance(cell, str):
+        cell = finat.ufl.as_cell(cell)
     if not isinstance(cell, ufl.AbstractCell):
         raise ValueError("Expecting a UFL Cell")
-    return ufc_cell(cell)
+
+    if hasattr(cell, "to_fiat"):
+        return cell.to_fiat()
+    else:
+        return ufc_cell(cell)
 
 
 @singledispatch
@@ -325,23 +334,28 @@ def convert_restrictedelement(element, **kwargs):
     return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps
 
 
-hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
-quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
+@convert.register(finat.ufl.FuseElement)
+def convert_fuse_element(element, **kwargs):
+    return finat.fiat_elements.FuseElement(element.triple), set()
+
+
 _cache = weakref.WeakKeyDictionary()
 
 
-def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None):
+def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None, cell_backend=finat.ufl.CellBackend.FIAT):
     """Create a FInAT element (suitable for tabulating with) given a UFL element.
 
     :arg ufl_element: The UFL element to create a FInAT element from.
     :arg shape_innermost: Vector/tensor indices come after basis function indices
     :arg restriction: cell restriction in interior facet integrals
                       (only for runtime tabulated elements)
+    :arg cell_backend: Enum determining the default cell type to use.
     """
     finat_element, deps = _create_element(ufl_element,
                                           shape_innermost=shape_innermost,
                                           shift_axes=shift_axes,
-                                          restriction=restriction)
+                                          restriction=restriction,
+                                          cell_backend=cell_backend)
     return finat_element
 
 

finat/fiat_elements.py

@@ -439,3 +439,8 @@ def __init__(self, cell, degree, **kwargs):
 class NedelecSecondKind(VectorFiatElement):
     def __init__(self, cell, degree, **kwargs):
         super().__init__(FIAT.NedelecSecondKind(cell, degree, **kwargs))
+
+
+class FuseElement(FiatElement):
+    def __init__(self, triple):
+        super(FuseElement, self).__init__(triple.to_fiat())

finat/ufl/__init__.py

@@ -14,8 +14,9 @@
 from finat.ufl.brokenelement import BrokenElement  # noqa: F401
 from finat.ufl.enrichedelement import EnrichedElement, NodalEnrichedElement  # noqa: F401
 from finat.ufl.finiteelement import FiniteElement  # noqa: F401
-from finat.ufl.finiteelementbase import FiniteElementBase  # noqa: F401
+from finat.ufl.finiteelementbase import FiniteElementBase, as_cell, CellBackend  # noqa: F401
 from finat.ufl.hdivcurl import HCurlElement, HDivElement, WithMapping, HDiv, HCurl  # noqa: F401
 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/elementlist.py

@@ -18,6 +18,7 @@
 # Modified by Pablo Brubeck, 2024
 
 import warnings
+import inspect
 
 from numpy import asarray
 
@@ -416,6 +417,10 @@ def canonical_element_description(family, cell, order, form_degree):
             raise ValueError("Need dimension to handle element aliases.")
         (family, order) = aliases[family](family, tdim, order, form_degree)
 
+    # Check that we don't have a raw FUSE object
+    if not isinstance(family, str) and 'fuse' in inspect.getmodule(family).__name__:
+        raise ValueError("Received unconverted FUSE triple - did you forget to call to_ufl()?")
+
     # Check that the element family exists
     if family not in ufl_elements:
         raise ValueError(f"Unknown finite element '{family}'.")

finat/ufl/finiteelement.py

@@ -11,9 +11,9 @@
 # Modified by Massimiliano Leoni, 2016
 # Modified by Matthew Scroggs, 2023
 
-from ufl.cell import TensorProductCell, as_cell
+from ufl.cell import TensorProductCell
 from finat.ufl.elementlist import canonical_element_description, simplices
-from finat.ufl.finiteelementbase import FiniteElementBase
+from finat.ufl.finiteelementbase import FiniteElementBase, as_cell
 from ufl.utils.formatting import istr
 
 

finat/ufl/finiteelementbase.py

@@ -14,9 +14,10 @@
 from abc import abstractmethod, abstractproperty
 from hashlib import md5
 from typing import Sequence, Union
+from enum import Enum
 
 from ufl import pullback
-from ufl.cell import AbstractCell, as_cell
+from ufl.cell import AbstractCell, as_cell as as_cell_ufl
 from ufl.finiteelement import AbstractFiniteElement
 from ufl.utils.sequences import product
 
@@ -291,3 +292,21 @@ def pullback(self):
             return supported_pullbacks[self.mapping()]
         except KeyError:
             raise ValueError(f"Unsupported mapping: {self.mapping()}")
+
+
+class CellBackend(Enum):
+    FIAT = 1
+    FUSE = 2
+
+
+def as_cell(cell: AbstractCell | str | tuple[AbstractCell, ...], cell_backend: CellBackend = CellBackend.FIAT) -> AbstractCell:
+    if isinstance(cell, str) and cell_backend == CellBackend.FUSE:
+        try:
+            import fuse
+        except ModuleNotFoundError as exc:
+            raise ModuleNotFoundError(
+                "FUSE cell construction requires the optional 'fuse' dependency. "
+            ) from exc
+        return fuse.constructCellComplex(cell)
+    else:
+        return as_cell_ufl(cell)