From dbc1d560537e335c8facbdf6d514969c5c7143a6 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 24 Jul 2026 12:04:58 +0100 Subject: [PATCH 01/10] Minimal changes to FIAT for simplex fuse to work --- FIAT/discontinuous_pc.py | 28 +++------- FIAT/orientation_utils.py | 6 +-- FIAT/reference_element.py | 87 ++++++++++++++++++++----------- finat/__init__.py | 3 +- finat/element_factory.py | 23 ++++++-- finat/fiat_elements.py | 5 ++ finat/ufl/__init__.py | 3 +- finat/ufl/elementlist.py | 2 + finat/ufl/finiteelement.py | 4 +- finat/ufl/finiteelementbase.py | 13 ++++- finat/ufl/fuseelement.py | 44 ++++++++++++++++ finat/ufl/mixedelement.py | 4 +- finat/ufl/tensorproductelement.py | 5 +- pyproject.toml | 1 + 14 files changed, 162 insertions(+), 66 deletions(-) create mode 100644 finat/ufl/fuseelement.py diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index dbb566906..8973a58f6 100644 --- a/FIAT/discontinuous_pc.py +++ b/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, diff --git a/FIAT/orientation_utils.py b/FIAT/orientation_utils.py index 2852f33f4..92cc68db6 100644 --- a/FIAT/orientation_utils.py +++ b/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). diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 87794b506..321b3cad1 100644 --- a/FIAT/reference_element.py +++ b/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,42 @@ 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: + 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 = {} + self.sub_entities_old = {} + for dim, entities in topology.items(): + self.sub_entities[dim] = {} + self.sub_entities_old[dim] = {} + + for e, v in entities.items(): + vertices = frozenset(v) + sub_entities = [] + sub_entities_old = [] + for dim_, entities_ in topology.items(): + for e_, vertices_ in entities_.items(): + if vertices.issuperset(vertices_): + sub_entities_old.append((dim_, e_)) + + # in order to maintain ordering, extract subentities from vertex numbering + entities_of_dim_ = list(entities_.values()) + + from itertools import permutations + # generate all possible sub entities + sub_list = permutations(v, len(entities_of_dim_[0])) + for s in sub_list: + # add the sub entities in the same order as in topology + for i, val in entities_.items(): + if set(s) == set(val) and (dim_, i) not in sub_entities: + sub_entities.append((dim_, i)) + + self.sub_entities[dim][e] = list(sub_entities) + self.sub_entities_old[dim][e] = list(sub_entities_old) + self.sub_entities = self.sub_entities_old # 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 +201,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 +404,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_ents=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_ents) def compute_normal(self, facet_i, cell=None): """Returns the unit normal vector to facet i of codimension 1.""" @@ -528,7 +545,8 @@ 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]],) + # return (self.get_vertices()[entity_id], ) elif 0 < dim <= self.get_spatial_dimension(): entity_verts = \ self.get_vertices_of_subcomplex( @@ -1144,7 +1162,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() @@ -1167,7 +1185,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): @@ -1317,10 +1335,14 @@ def compare(self, op, other): This is done dimension by dimension.""" if hasattr(other, "product"): other = other.product - if isinstance(other, type(self)): + if isinstance(other, TensorProductCell): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: - return op(self, other) + if op == operator.gt or op == operator.lt or operator.ne: + return not op(other, self) + if op == operator.ge or op == operator.le: + return not op(other, self) or operator.eq(self, other) + raise ValueError("Unknown operator in cell comparison") def __gt__(self, other): return self.compare(operator.gt, other) @@ -1410,7 +1432,7 @@ 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] @@ -1418,7 +1440,7 @@ def __init__(self, dimension, product): 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) @@ -1845,3 +1867,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()) diff --git a/finat/__init__.py b/finat/__init__.py index bde435a1f..34ea55b88 100644 --- a/finat/__init__.py +++ b/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 diff --git a/finat/element_factory.py b/finat/element_factory.py index 8563aa60f..e8947d80e 100644 --- a/finat/element_factory.py +++ b/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,12 +334,15 @@ 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, use_fuse=False): """Create a FInAT element (suitable for tabulating with) given a UFL element. :arg ufl_element: The UFL element to create a FInAT element from. @@ -341,7 +353,8 @@ def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction= finat_element, deps = _create_element(ufl_element, shape_innermost=shape_innermost, shift_axes=shift_axes, - restriction=restriction) + restriction=restriction, + use_fuse=use_fuse) return finat_element diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 22408cf07..f52c55ca0 100644 --- a/finat/fiat_elements.py +++ b/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()) diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py index 21a7d13d1..1cb7f3dca 100644 --- a/finat/ufl/__init__.py +++ b/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 # 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 diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index 7bd949a3f..b7f37feb5 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -418,6 +418,8 @@ def canonical_element_description(family, cell, order, form_degree): # Check that the element family exists if family not in ufl_elements: + if repr(family).startswith("FuseTriple"): + raise ValueError("Recieved unconverted FUSE triple - did you forget to call to_ufl()?") raise ValueError(f"Unknown finite element '{family}'.") # Check that element data is valid (and also get common family diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py index dade48b06..bc2a0bdb0 100644 --- a/finat/ufl/finiteelement.py +++ b/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 diff --git a/finat/ufl/finiteelementbase.py b/finat/ufl/finiteelementbase.py index 41602b025..2dfee9aa1 100644 --- a/finat/ufl/finiteelementbase.py +++ b/finat/ufl/finiteelementbase.py @@ -15,7 +15,7 @@ from hashlib import md5 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 @@ -272,3 +272,14 @@ def pullback(self): return supported_pullbacks[self.mapping()] except KeyError: raise ValueError(f"Unsupported mapping: {self.mapping()}") + + +def as_cell(cell: AbstractCell | str | tuple[AbstractCell, ...], use_fuse: bool = False) -> AbstractCell: + if isinstance(cell, str) and use_fuse: + try: + import fuse + except ModuleNotFoundError: + raise ModuleNotFoundError("Cannot create FUSE cell without FUSE") + return fuse.constructCellComplex(cell) + else: + return as_cell_ufl(cell) diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py new file mode 100644 index 000000000..d5b259851 --- /dev/null +++ b/finat/ufl/fuseelement.py @@ -0,0 +1,44 @@ +"""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() + + degree = self.triple.degree + self.sobolev_space = self.triple.spaces[1].to_ufl() + super(FuseElement, self).__init__("IT", cell, degree, None, triple.get_value_shape()) + + def __repr__(self): + return repr(self.triple) + + def __str__(self): + return "" % (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) diff --git a/finat/ufl/mixedelement.py b/finat/ufl/mixedelement.py index 0fda85303..75c7cae8c 100644 --- a/finat/ufl/mixedelement.py +++ b/finat/ufl/mixedelement.py @@ -13,10 +13,10 @@ import numpy as np -from ufl.cell import CellSequence, as_cell +from ufl.cell import CellSequence from ufl.domain import MeshSequence from finat.ufl.finiteelement import FiniteElement -from finat.ufl.finiteelementbase import FiniteElementBase +from finat.ufl.finiteelementbase import FiniteElementBase, as_cell from ufl.permutation import compute_indices from ufl.pullback import MixedPullback, SymmetricPullback from ufl.utils.indexflattening import flatten_multiindex, shape_to_strides, unflatten_index diff --git a/finat/ufl/tensorproductelement.py b/finat/ufl/tensorproductelement.py index a5e81a35c..063822248 100644 --- a/finat/ufl/tensorproductelement.py +++ b/finat/ufl/tensorproductelement.py @@ -13,8 +13,9 @@ from itertools import chain -from ufl.cell import TensorProductCell, as_cell -from finat.ufl.finiteelementbase import FiniteElementBase +from ufl.cell import TensorProductCell +from finat.ufl.finiteelementbase import FiniteElementBase, as_cell + from ufl.sobolevspace import DirectionalSobolevSpace diff --git a/pyproject.toml b/pyproject.toml index e2668d681..20d3d0766 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,6 +14,7 @@ dependencies = [ "scipy", "symengine", "sympy", + "fuse-element @ git+https://github.com/firedrakeproject/fuse.git", ] requires-python = ">=3.11" authors = [ From 68ac357d982e1eed1d2b9494b14e9e8237768f21 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 24 Jul 2026 14:49:43 +0100 Subject: [PATCH 02/10] remove old subentity logic --- FIAT/reference_element.py | 23 +++-------------------- 1 file changed, 3 insertions(+), 20 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 321b3cad1..73d56abd7 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -152,35 +152,18 @@ def __init__(self, shape, vertices, topology, sub_entities=None): # Given the topology, work out for each entity in the cell, # which other entities it contains. self.sub_entities = {} - self.sub_entities_old = {} for dim, entities in topology.items(): self.sub_entities[dim] = {} - self.sub_entities_old[dim] = {} for e, v in entities.items(): vertices = frozenset(v) sub_entities = [] - sub_entities_old = [] for dim_, entities_ in topology.items(): for e_, vertices_ in entities_.items(): if vertices.issuperset(vertices_): - sub_entities_old.append((dim_, e_)) - - # in order to maintain ordering, extract subentities from vertex numbering - entities_of_dim_ = list(entities_.values()) - - from itertools import permutations - # generate all possible sub entities - sub_list = permutations(v, len(entities_of_dim_[0])) - for s in sub_list: - # add the sub entities in the same order as in topology - for i, val in entities_.items(): - if set(s) == set(val) and (dim_, i) not in sub_entities: - sub_entities.append((dim_, i)) - - self.sub_entities[dim][e] = list(sub_entities) - self.sub_entities_old[dim][e] = list(sub_entities_old) - self.sub_entities = self.sub_entities_old + 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: From 6000936334233203968550b68d363bc364ac4af3 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 31 Jul 2026 10:42:00 +0100 Subject: [PATCH 03/10] move fuse to optional dependencies --- pyproject.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 20d3d0766..96f46ad79 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,6 @@ dependencies = [ "scipy", "symengine", "sympy", - "fuse-element @ git+https://github.com/firedrakeproject/fuse.git", ] requires-python = ">=3.11" authors = [ @@ -35,6 +34,9 @@ doc = [ "sphinx", ] test = ["pytest"] +fuse = [ + "fuse-element @ git+https://github.com/firedrakeproject/fuse.git", +] [tool.setuptools] packages = ["FIAT", "finat", "finat.ufl", "gem"] From 22b59881616c76bbe6d4e742d99f23785f2a390c Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 3 Aug 2026 10:27:52 +0100 Subject: [PATCH 04/10] remove comparison change --- FIAT/reference_element.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index efe386b26..7d59b6614 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -529,7 +529,6 @@ def make_points(self, dim, entity_id, order, variant=None, interior=1): include in each direction.""" if dim == 0: 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( @@ -1332,10 +1331,10 @@ def compare(self, op, other): if isinstance(other, TensorProductCell): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: - if op == operator.gt or op == operator.lt or operator.ne: - return not op(other, self) - if op == operator.ge or op == operator.le: - return not op(other, self) or operator.eq(self, other) + #if op == operator.gt or op == operator.lt or operator.ne: + # return not op(other, self) + #if op == operator.ge or op == operator.le: + # return not op(other, self) or operator.eq(self, other) raise ValueError("Unknown operator in cell comparison") def __gt__(self, other): From d014892e97e92e3f59efb5903e10276f92e94f57 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 3 Aug 2026 16:06:24 +0100 Subject: [PATCH 05/10] remove commented out code --- FIAT/reference_element.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 7d59b6614..61a94a9ec 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1331,10 +1331,6 @@ def compare(self, op, other): if isinstance(other, TensorProductCell): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: - #if op == operator.gt or op == operator.lt or operator.ne: - # return not op(other, self) - #if op == operator.ge or op == operator.le: - # return not op(other, self) or operator.eq(self, other) raise ValueError("Unknown operator in cell comparison") def __gt__(self, other): From 5a353b20a5b980a47620ab23e44009058b706e69 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 6 Aug 2026 09:57:29 +0100 Subject: [PATCH 06/10] rename subents --- FIAT/reference_element.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index efe386b26..b757c7ed6 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -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, sub_ents=None): + 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, sub_ents) + 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.""" From c73600cc49a945c304bcc46a000eca4c4ae7615f Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 6 Aug 2026 11:54:45 +0100 Subject: [PATCH 07/10] add tests for subentities --- FIAT/reference_element.py | 4 +-- finat/element_factory.py | 1 + finat/ufl/elementlist.py | 6 ++-- finat/ufl/finiteelementbase.py | 6 ++-- finat/ufl/fuseelement.py | 13 +++++-- test/FIAT/unit/test_reference_element.py | 46 ++++++++++++++++++++++++ 6 files changed, 68 insertions(+), 8 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index b757c7ed6..13dae2ce4 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -145,7 +145,7 @@ def __init__(self, shape, vertices, topology, sub_entities=None): self.vertices = vertices self.topology = topology - if sub_entities: + if sub_entities is not None: self.sub_entities = sub_entities else: # If sub entity list not provided @@ -1332,7 +1332,7 @@ def compare(self, op, other): if isinstance(other, TensorProductCell): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: - if op == operator.gt or op == operator.lt or operator.ne: + if op == operator.gt or op == operator.lt or op == operator.ne: return not op(other, self) if op == operator.ge or op == operator.le: return not op(other, self) or operator.eq(self, other) diff --git a/finat/element_factory.py b/finat/element_factory.py index e8947d80e..49a27ed78 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -349,6 +349,7 @@ def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction= :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 use_fuse: When creating cells, ensure they are fuse compatible. inherited from the mesh. """ finat_element, deps = _create_element(ufl_element, shape_innermost=shape_innermost, diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index b7f37feb5..49db23e0b 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -416,10 +416,12 @@ 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 repr(family).startswith("FuseTriple"): + raise ValueError("Recieved unconverted FUSE triple - did you forget to call to_ufl()?") + # Check that the element family exists if family not in ufl_elements: - if repr(family).startswith("FuseTriple"): - raise ValueError("Recieved unconverted FUSE triple - did you forget to call to_ufl()?") raise ValueError(f"Unknown finite element '{family}'.") # Check that element data is valid (and also get common family diff --git a/finat/ufl/finiteelementbase.py b/finat/ufl/finiteelementbase.py index 2dfee9aa1..a2892c9e5 100644 --- a/finat/ufl/finiteelementbase.py +++ b/finat/ufl/finiteelementbase.py @@ -278,8 +278,10 @@ def as_cell(cell: AbstractCell | str | tuple[AbstractCell, ...], use_fuse: bool if isinstance(cell, str) and use_fuse: try: import fuse - except ModuleNotFoundError: - raise ModuleNotFoundError("Cannot create FUSE cell without 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) diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py index d5b259851..eb42cfe62 100644 --- a/finat/ufl/fuseelement.py +++ b/finat/ufl/fuseelement.py @@ -11,17 +11,26 @@ class FuseElement(FiniteElementBase): """ A finite element defined using FUSE. - TODO: Need to deal with cases where value shape and reference value shape are different + :arg triple: An ElementTriple object defined with FUSE + :arg cell: Optional (defaults to triple.cell) The cell the element is defined on + """ def __init__(self, triple, cell=None): + try: + import fuse + except ModuleNotFoundError as exc: + raise ModuleNotFoundError( + "FUSE element creation requires the optional 'fuse' dependency. " + ) from exc + assert isinstance(triple, fuse.ElementTriple) self.triple = triple if not cell: cell = self.triple.cell.to_ufl() degree = self.triple.degree self.sobolev_space = self.triple.spaces[1].to_ufl() - super(FuseElement, self).__init__("IT", cell, degree, None, triple.get_value_shape()) + super().__init__("IT", cell, degree, None, triple.get_value_shape()) def __repr__(self): return repr(self.triple) diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index d166ba6ef..0ce630025 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -22,6 +22,7 @@ from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron from FIAT.reference_element import is_ufc, is_hypercube, default_simplex, flatten_reference_cube, Hypercube +from FIAT.reference_element import Cell point = Point() interval = UFCInterval() @@ -86,6 +87,51 @@ def test_ufc_connectivity_Dx(cell): assert connectivity[0] == tuple(range(len(connectivity[0]))) +def test_explicit_sub_entities_reproduces_default(): + """An explicit sub_entities mapping equal to the one Cell would + compute automatically must produce identical super_entities and + connectivity.""" + topology = {0: {0: (0,), 1: (1,)}, 1: {0: (0, 1)}} + vertices = ((0.0,), (1.0,)) + + reference = Cell("interval", vertices, topology) + explicit = Cell("interval", vertices, topology, + sub_entities=reference.sub_entities) + + assert explicit.sub_entities == reference.sub_entities + assert explicit.super_entities == reference.super_entities + assert explicit.connectivity == reference.connectivity + + +def test_explicit_sub_entities_custom_mapping(): + """A deliberately different (but self-consistent) explicit + sub_entities mapping should be used as-is, and super_entities / + connectivity should be derived from it rather than recomputed + from the topology.""" + topology = {0: {0: (0,), 1: (1,)}, 1: {0: (0, 1)}} + vertices = ((0.0,), (1.0,)) + + # Deliberately omit the "each entity is a sub-entity of itself" + # entries that the automatic computation would include. + custom_sub_entities = { + 0: {0: [], 1: []}, + 1: {0: [(0, 0), (0, 1)]}, + } + + cell = Cell("interval", vertices, topology, + sub_entities=custom_sub_entities) + + assert cell.sub_entities == custom_sub_entities + assert cell.super_entities == { + 0: {0: [(1, 0)], 1: [(1, 0)]}, + 1: {0: []}, + } + assert cell.connectivity[(1, 0)] == [(0, 1)] + assert cell.connectivity[(0, 1)] == [(0,), (0,)] + assert cell.connectivity[(1, 1)] == [()] + assert cell.connectivity[(0, 0)] == [(), ()] + + @pytest.mark.parametrize(('cell', 'volume'), [pytest.param(point, 1, marks=pytest.mark.xfail(conditional=sys.version_info < (3, 6))), (interval, 1), From d0e5e7729c946ecadb374d085237582772db9a79 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 6 Aug 2026 17:05:48 +0100 Subject: [PATCH 08/10] remove sobolev space, use pullback --- finat/ufl/fuseelement.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py index eb42cfe62..38c80473d 100644 --- a/finat/ufl/fuseelement.py +++ b/finat/ufl/fuseelement.py @@ -29,7 +29,6 @@ def __init__(self, triple, cell=None): cell = self.triple.cell.to_ufl() degree = self.triple.degree - self.sobolev_space = self.triple.spaces[1].to_ufl() super().__init__("IT", cell, degree, None, triple.get_value_shape()) def __repr__(self): @@ -38,16 +37,13 @@ def __repr__(self): def __str__(self): return "" % (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" - + @property def sobolev_space(self): - return self.triple.spaces[1] + return self.triple.spaces[2].ufl_sobolev_space(self.triple.form_degree, + self.triple.cell.dim()) + + def mapping(self): + return self.triple.spaces[2].mapping() def reconstruct(self, family=None, cell=None, degree=None, quad_scheme=None, variant=None): return FuseElement(self.triple, cell=cell) From 4fb0438b157a62ad8bcfdd2291d5eafb189d16cd Mon Sep 17 00:00:00 2001 From: India Marsden Date: Wed, 12 Aug 2026 09:39:00 +0100 Subject: [PATCH 09/10] Change to enum, remove out of date bits --- FIAT/reference_element.py | 4 ++-- finat/element_factory.py | 6 +++--- finat/ufl/__init__.py | 2 +- finat/ufl/elementlist.py | 5 +++-- finat/ufl/finiteelementbase.py | 8 ++++++-- finat/ufl/fuseelement.py | 4 ++-- 6 files changed, 17 insertions(+), 12 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 7b26ff1b7..aae357bdb 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1328,10 +1328,10 @@ def compare(self, op, other): This is done dimension by dimension.""" if hasattr(other, "product"): other = other.product - if isinstance(other, TensorProductCell): + if isinstance(other, type(self)): return all(op(a, b) for a, b in zip(self.cells, other.cells)) else: - raise ValueError("Unknown operator in cell comparison") + return op(self, other) def __gt__(self, other): return self.compare(operator.gt, other) diff --git a/finat/element_factory.py b/finat/element_factory.py index 49a27ed78..b5d5715fc 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -342,20 +342,20 @@ def convert_fuse_element(element, **kwargs): _cache = weakref.WeakKeyDictionary() -def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None, use_fuse=False): +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 use_fuse: When creating cells, ensure they are fuse compatible. inherited from the mesh. + :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, - use_fuse=use_fuse) + cell_backend=cell_backend) return finat_element diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py index 1cb7f3dca..9da7f2dbc 100644 --- a/finat/ufl/__init__.py +++ b/finat/ufl/__init__.py @@ -14,7 +14,7 @@ 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, as_cell # 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 diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index 49db23e0b..865ad00ff 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -18,6 +18,7 @@ # Modified by Pablo Brubeck, 2024 import warnings +import inspect from numpy import asarray @@ -417,8 +418,8 @@ def canonical_element_description(family, cell, order, form_degree): (family, order) = aliases[family](family, tdim, order, form_degree) # Check that we don't have a raw FUSE object - if repr(family).startswith("FuseTriple"): - raise ValueError("Recieved unconverted FUSE triple - did you forget to call to_ufl()?") + 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: diff --git a/finat/ufl/finiteelementbase.py b/finat/ufl/finiteelementbase.py index 26601c858..fa995cc12 100644 --- a/finat/ufl/finiteelementbase.py +++ b/finat/ufl/finiteelementbase.py @@ -14,6 +14,7 @@ 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 as as_cell_ufl @@ -292,9 +293,12 @@ def pullback(self): except KeyError: raise ValueError(f"Unsupported mapping: {self.mapping()}") +class CellBackend(Enum): + FIAT = 1 + FUSE = 2 -def as_cell(cell: AbstractCell | str | tuple[AbstractCell, ...], use_fuse: bool = False) -> AbstractCell: - if isinstance(cell, str) and use_fuse: +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: diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py index eb42cfe62..b36f68a8a 100644 --- a/finat/ufl/fuseelement.py +++ b/finat/ufl/fuseelement.py @@ -30,13 +30,13 @@ def __init__(self, triple, cell=None): degree = self.triple.degree self.sobolev_space = self.triple.spaces[1].to_ufl() - super().__init__("IT", cell, degree, None, triple.get_value_shape()) + super().__init__("FUSE", cell, degree, None, triple.get_value_shape()) def __repr__(self): return repr(self.triple) def __str__(self): - return "" % (self.triple.spaces[0], self.triple.cell) + return f"" def mapping(self): if str(self.sobolev_space) == "HCurl": From 0f46f319f734b3a74127a84891c9aa23ac88e5ed Mon Sep 17 00:00:00 2001 From: India Marsden Date: Thu, 13 Aug 2026 17:32:58 +0100 Subject: [PATCH 10/10] lint --- finat/ufl/finiteelementbase.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/finat/ufl/finiteelementbase.py b/finat/ufl/finiteelementbase.py index fa995cc12..1dcac4891 100644 --- a/finat/ufl/finiteelementbase.py +++ b/finat/ufl/finiteelementbase.py @@ -293,10 +293,12 @@ def pullback(self): 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: