From 770b0a4c8bdea9fbc968198b4ba05424a08d23c0 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 12:13:28 +0000 Subject: [PATCH 1/8] Restrict PhysicallyMappedElement --- finat/physically_mapped.py | 15 +++++++------ finat/restricted.py | 39 ++++++++++++++++++++++++++++++---- test/finat/test_restriction.py | 22 ++++++++++++++++++- 3 files changed, 65 insertions(+), 11 deletions(-) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index 978bbaef1..c5a6a8727 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -19,13 +19,16 @@ class MappedTabulation(Mapping): """A lazy tabulation dict that applies the basis transformation only on the requested derivatives.""" - def __init__(self, M, ref_tabulation): + def __init__(self, M, ref_tabulation, indices=None): self.M = M self.ref_tabulation = ref_tabulation + if indices is None: + indices = list(range(M.shape[0])) + self.indices = indices # we expect M to be sparse with O(1) nonzeros per row # for each row, get the column index of each nonzero entry csr = [[j for j in range(M.shape[1]) if not isinstance(M.array[i, j], gem.Zero)] - for i in range(M.shape[0])] + for i in indices] self.csr = csr self._tabulation_cache = {} @@ -35,8 +38,7 @@ def matvec(self, table): phi = [gem.Indexed(table, (j, *ii)) for j in range(self.M.shape[1])] # the sum approach is faster than calling numpy.dot or gem.IndexSum exprs = [gem.ComponentTensor(gem.Sum(*(self.M.array[i, j] * phi[j] for j in js)), ii) - for i, js in enumerate(self.csr)] - + for i, js in zip(self.indices, self.csr)] val = gem.ListTensor(exprs) # val = self.M @ table return gem.optimise.aggressive_unroll(val) @@ -63,6 +65,7 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) cite("Kirby2018zany") cite("Kirby2019zany") + self.indices = None @abstractmethod def basis_transformation(self, coordinate_mapping): @@ -74,7 +77,7 @@ def basis_transformation(self, coordinate_mapping): def map_tabulation(self, ref_tabulation, coordinate_mapping): assert coordinate_mapping is not None M = self.basis_transformation(coordinate_mapping) - return MappedTabulation(M, ref_tabulation) + return MappedTabulation(M, ref_tabulation, indices=self.indices) def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None): result = super().basis_evaluation(order, ps, entity=entity) @@ -93,7 +96,7 @@ def dual_transformation(self, Q, coordinate_mapping=None): M_dual = gem.ListTensor(inverse(M.T)) key = None - return MappedTabulation(M_dual, {key: Q})[key] + return MappedTabulation(M_dual, {key: Q}, indices=self.indices)[key] class DirectlyDefinedElement(NeedsCoordinateMappingElement): diff --git a/finat/restricted.py b/finat/restricted.py index a85e0d393..2100d8868 100644 --- a/finat/restricted.py +++ b/finat/restricted.py @@ -1,5 +1,6 @@ from functools import singledispatch from itertools import chain +from copy import deepcopy import FIAT from FIAT.polynomial_set import mis @@ -13,6 +14,24 @@ null_element = object() +class RestrictedPhysicallyMappedElement(PhysicallyMappedElement, FiatElement): + """A restricted PhysicallyMappedElement.""" + def __init__(self, element, indices, entity_dofs): + super().__init__(element._element) + self.indices = indices + self._entity_dofs = entity_dofs + self.full_element = element + + def basis_transformation(self, coordinate_mapping): + return self.full_element.basis_transformation(coordinate_mapping) + + def space_dimension(self): + return len(self.indices) + + def entity_dofs(self): + return self._entity_dofs + + @singledispatch def restrict(element, domain, take_closure): """Restrict an element to a given subentity. @@ -30,6 +49,7 @@ def restrict(element, domain, take_closure): @restrict.register(FiatElement) +@restrict.register(PhysicallyMappedElement) def restrict_fiat(element, domain, take_closure): try: re = FIAT.RestrictedElement(element._element, @@ -43,12 +63,23 @@ def restrict_fiat(element, domain, take_closure): # In case the restriction is trivial we return the original element # to avoid reconstructing the space with an undesired permutation. return element - return FiatElement(re) + if isinstance(element, PhysicallyMappedElement) and not (domain == "interior" and not take_closure): + dual = element._element.get_dual_set() + indices = dual.get_indices(domain, take_closure) + if element.space_dimension() < element._element.space_dimension(): + # Throw away constraint dofs + indices = [i for i in indices if i < element.space_dimension()] + entity_dofs = deepcopy(element.entity_dofs()) + for dim in entity_dofs: + for entity in entity_dofs[dim]: + ids = entity_dofs[dim][entity] + entity_dofs[dim][entity] = [indices.index(i) for i in ids if i in indices] + else: + entity_dofs = re.entity_dofs() + return RestrictedPhysicallyMappedElement(element, indices, entity_dofs) -@restrict.register(PhysicallyMappedElement) -def restrict_physically_mapped(element, domain, take_closure): - raise NotImplementedError("Can't restrict Physically Mapped things") + return FiatElement(re) @restrict.register(finat.FlattenedDimensions) diff --git a/test/finat/test_restriction.py b/test/finat/test_restriction.py index 443292a04..8ba91d7a2 100644 --- a/test/finat/test_restriction.py +++ b/test/finat/test_restriction.py @@ -5,10 +5,18 @@ from finat.point_set import PointSet from finat.restricted import r_to_codim from gem.interpreter import evaluate +from finat.physically_mapped import NeedsCoordinateMappingElement +from .conftest import MyMapping def tabulate(element, ps): - tabulation, = element.basis_evaluation(0, ps).values() + coordinate_mapping = None + if isinstance(element, NeedsCoordinateMappingElement): + ref_el = element.cell + sd = ref_el.get_spatial_dimension() + phys_el = FIAT.reference_element.symmetric_simplex(sd) + coordinate_mapping = MyMapping(ref_el, phys_el) + tabulation, = element.basis_evaluation(0, ps, coordinate_mapping=coordinate_mapping).values() result, = evaluate([tabulation]) # Singleton point shape = (int(numpy.prod(element.index_shape)), ) + element.value_shape @@ -142,3 +150,15 @@ def test_hdiv_restriction(hdiv_element, restriction, ps): def test_hcurl_restriction(hcurl_element, restriction, ps): run_restriction(hcurl_element, restriction, ps) + + +@pytest.fixture +def zany_element(cell): + if len(cell) == 1: + return finat.Walkington(cell[0]) + else: + pytest.skip() + + +def test_zany_restriction(zany_element, restriction, ps): + run_restriction(zany_element, restriction, ps) From c6b2842b5dae1beb226fcb023c1592ac4ded9bc0 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 13:05:56 +0000 Subject: [PATCH 2/8] fixes --- finat/restricted.py | 2 +- finat/ufl/restrictedelement.py | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/finat/restricted.py b/finat/restricted.py index 2100d8868..4f6f8ebee 100644 --- a/finat/restricted.py +++ b/finat/restricted.py @@ -59,7 +59,7 @@ def restrict_fiat(element, domain, take_closure): return null_element if element.space_dimension() == re.space_dimension(): - # FIAT.RestrictedElement wipes out entity_permuations. + # FIAT.RestrictedElement wipes out entity_permutations. # In case the restriction is trivial we return the original element # to avoid reconstructing the space with an undesired permutation. return element diff --git a/finat/ufl/restrictedelement.py b/finat/ufl/restrictedelement.py index 921acd026..c2e518a6e 100644 --- a/finat/ufl/restrictedelement.py +++ b/finat/ufl/restrictedelement.py @@ -61,10 +61,7 @@ def __repr__(self): @property def sobolev_space(self): """Doc.""" - if self._restriction_domain == "interior": - return L2 - else: - return self._element.sobolev_space + return self._element.sobolev_space def is_cellwise_constant(self): """Return whether the basis functions of this element is spatially constant over each cell.""" From 6ae8c56aa3284ab2edcdc13d892a60c6fe3970df Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 14:20:11 +0000 Subject: [PATCH 3/8] fix --- finat/ufl/restrictedelement.py | 1 - 1 file changed, 1 deletion(-) diff --git a/finat/ufl/restrictedelement.py b/finat/ufl/restrictedelement.py index c2e518a6e..225bd33c9 100644 --- a/finat/ufl/restrictedelement.py +++ b/finat/ufl/restrictedelement.py @@ -13,7 +13,6 @@ from finat.ufl.finiteelementbase import FiniteElementBase from finat.ufl.mixedelement import MixedElement, VectorElement, TensorElement -from ufl.sobolevspace import L2 valid_restriction_domains = ("interior", "facet", "ridge", "face", "edge", "vertex", "reduced") From 1a67f875d560bf73967a424445cad5409c50dc7a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Feb 2026 14:39:27 +0000 Subject: [PATCH 4/8] dual_transformation --- finat/physically_mapped.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index c5a6a8727..9bdd1f503 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -93,10 +93,13 @@ def dual_transformation(self, Q, coordinate_mapping=None): M = M.array if M.shape[0] != M.shape[1]: M = M[:, :self.space_dimension()] - M_dual = gem.ListTensor(inverse(M.T)) + M_dual = inverse(M.T) + if self.indices is not None: + M_dual = M_dual[numpy.ix_(self.indices, self.indices)] + M_dual = gem.ListTensor(M_dual) key = None - return MappedTabulation(M_dual, {key: Q}, indices=self.indices)[key] + return MappedTabulation(M_dual, {key: Q})[key] class DirectlyDefinedElement(NeedsCoordinateMappingElement): From b6e0652a4b588bdd26e18c39f5c02d346957dc6a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 16 Feb 2026 17:42:23 +0000 Subject: [PATCH 5/8] fix --- finat/physically_mapped.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index 9bdd1f503..bd66b2ff5 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -91,8 +91,9 @@ def dual_transformation(self, Q, coordinate_mapping=None): M = self.basis_transformation(coordinate_mapping) M = M.array - if M.shape[0] != M.shape[1]: - M = M[:, :self.space_dimension()] + if M.shape[1] > M.shape[0]: + M = M[:, :M.shape[0]] + M_dual = inverse(M.T) if self.indices is not None: M_dual = M_dual[numpy.ix_(self.indices, self.indices)] From e52200067811b6de053e3a8e7a465a533ac1cf62 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 17 Feb 2026 10:21:09 +0000 Subject: [PATCH 6/8] cleanup --- finat/physically_mapped.py | 9 ++++---- finat/restricted.py | 46 ++++++++++++++++++++++---------------- 2 files changed, 32 insertions(+), 23 deletions(-) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index bd66b2ff5..ea6931387 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -65,7 +65,7 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) cite("Kirby2018zany") cite("Kirby2019zany") - self.indices = None + self.restriction_indices = None @abstractmethod def basis_transformation(self, coordinate_mapping): @@ -77,7 +77,7 @@ def basis_transformation(self, coordinate_mapping): def map_tabulation(self, ref_tabulation, coordinate_mapping): assert coordinate_mapping is not None M = self.basis_transformation(coordinate_mapping) - return MappedTabulation(M, ref_tabulation, indices=self.indices) + return MappedTabulation(M, ref_tabulation, indices=self.restriction_indices) def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None): result = super().basis_evaluation(order, ps, entity=entity) @@ -95,8 +95,9 @@ def dual_transformation(self, Q, coordinate_mapping=None): M = M[:, :M.shape[0]] M_dual = inverse(M.T) - if self.indices is not None: - M_dual = M_dual[numpy.ix_(self.indices, self.indices)] + if self.restriction_indices is not None: + indices = self.restriction_indices + M_dual = M_dual[numpy.ix_(indices, indices)] M_dual = gem.ListTensor(M_dual) key = None diff --git a/finat/restricted.py b/finat/restricted.py index 4f6f8ebee..7bc9f62b9 100644 --- a/finat/restricted.py +++ b/finat/restricted.py @@ -1,6 +1,5 @@ from functools import singledispatch from itertools import chain -from copy import deepcopy import FIAT from FIAT.polynomial_set import mis @@ -15,21 +14,35 @@ class RestrictedPhysicallyMappedElement(PhysicallyMappedElement, FiatElement): - """A restricted PhysicallyMappedElement.""" - def __init__(self, element, indices, entity_dofs): + """A restricted PhysicallyMappedElement. + + :arg element: the finat.FiatElement to be restricted. + :arg indices: the indices of the degrees of freedom to be kept. + """ + def __init__(self, element, indices): super().__init__(element._element) - self.indices = indices - self._entity_dofs = entity_dofs - self.full_element = element + self.restriction_indices = indices + # Restrict the entity_dofs dict + edofs = element.entity_dofs() + rdofs = {d: {e: [indices.index(i) for i in edofs[d][e] if i in indices] + for e in edofs[d]} for d in edofs} + self.restriction_entity_dofs = rdofs + # Grab the basis transformation matrix from the parent element + if isinstance(element, PhysicallyMappedElement): + self.full_basis_transformation = element.basis_transformation + else: + self.full_basis_transformation = None def basis_transformation(self, coordinate_mapping): - return self.full_element.basis_transformation(coordinate_mapping) + if self.full_basis_transformation is None: + raise NotImplementedError("basis_transformation not implemented.") + return self.full_basis_transformation(coordinate_mapping) def space_dimension(self): - return len(self.indices) + return len(self.restriction_indices) def entity_dofs(self): - return self._entity_dofs + return self.restriction_entity_dofs @singledispatch @@ -68,16 +81,11 @@ def restrict_fiat(element, domain, take_closure): dual = element._element.get_dual_set() indices = dual.get_indices(domain, take_closure) if element.space_dimension() < element._element.space_dimension(): - # Throw away constraint dofs - indices = [i for i in indices if i < element.space_dimension()] - entity_dofs = deepcopy(element.entity_dofs()) - for dim in entity_dofs: - for entity in entity_dofs[dim]: - ids = entity_dofs[dim][entity] - entity_dofs[dim][entity] = [indices.index(i) for i in ids if i in indices] - else: - entity_dofs = re.entity_dofs() - return RestrictedPhysicallyMappedElement(element, indices, entity_dofs) + # Throw away constrained DOFs + space_dim = element.space_dimension() + indices = [i for i in indices if i < space_dim] + + return RestrictedPhysicallyMappedElement(element, indices) return FiatElement(re) From be2b0ef5681bd280918dc1b745401d4beeeb1a33 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 17 Feb 2026 10:38:03 +0000 Subject: [PATCH 7/8] Docs --- finat/physically_mapped.py | 6 +++++- finat/restricted.py | 4 ++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index c0733f9fe..df0d598ab 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -17,8 +17,12 @@ def dual_transformation(self, Q, coordinate_mapping=None): class MappedTabulation(Mapping): """A lazy tabulation dict that applies the basis transformation only - on the requested derivatives.""" + on the requested derivatives. + :arg M: a gem.ListTensor with the basis transformation matrix. + :arg ref_tabulation: a dict of tabulations on the reference cell. + :kwarg indices: an optional list of restriction indices on the basis functions. + """ def __init__(self, M, ref_tabulation, indices=None): self.M = M self.ref_tabulation = ref_tabulation diff --git a/finat/restricted.py b/finat/restricted.py index 7bc9f62b9..f07679c90 100644 --- a/finat/restricted.py +++ b/finat/restricted.py @@ -80,9 +80,9 @@ def restrict_fiat(element, domain, take_closure): if isinstance(element, PhysicallyMappedElement) and not (domain == "interior" and not take_closure): dual = element._element.get_dual_set() indices = dual.get_indices(domain, take_closure) - if element.space_dimension() < element._element.space_dimension(): + space_dim = element.space_dimension() + if space_dim < element._element.space_dimension(): # Throw away constrained DOFs - space_dim = element.space_dimension() indices = [i for i in indices if i < space_dim] return RestrictedPhysicallyMappedElement(element, indices) From ed6a0e3f8669e97a05f69e2cb4f422de86a91585 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 17 Feb 2026 14:09:20 +0000 Subject: [PATCH 8/8] cleanup --- finat/restricted.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/finat/restricted.py b/finat/restricted.py index f07679c90..b5144be9f 100644 --- a/finat/restricted.py +++ b/finat/restricted.py @@ -21,12 +21,19 @@ class RestrictedPhysicallyMappedElement(PhysicallyMappedElement, FiatElement): """ def __init__(self, element, indices): super().__init__(element._element) + # First sanitise the restriction indices + # Some finat elements, e.g. Bell, are the restriction of a FIAT element. + # Therefore, we need to compose the restrictions. + edofs = element.entity_dofs() + free_indices = set(chain.from_iterable(edofs[d][e] for d in edofs for e in edofs[d])) + indices = [i for i in indices if i in free_indices] self.restriction_indices = indices + # Restrict the entity_dofs dict - edofs = element.entity_dofs() rdofs = {d: {e: [indices.index(i) for i in edofs[d][e] if i in indices] for e in edofs[d]} for d in edofs} self.restriction_entity_dofs = rdofs + # Grab the basis transformation matrix from the parent element if isinstance(element, PhysicallyMappedElement): self.full_basis_transformation = element.basis_transformation @@ -62,7 +69,6 @@ def restrict(element, domain, take_closure): @restrict.register(FiatElement) -@restrict.register(PhysicallyMappedElement) def restrict_fiat(element, domain, take_closure): try: re = FIAT.RestrictedElement(element._element, @@ -78,14 +84,7 @@ def restrict_fiat(element, domain, take_closure): return element if isinstance(element, PhysicallyMappedElement) and not (domain == "interior" and not take_closure): - dual = element._element.get_dual_set() - indices = dual.get_indices(domain, take_closure) - space_dim = element.space_dimension() - if space_dim < element._element.space_dimension(): - # Throw away constrained DOFs - indices = [i for i in indices if i < space_dim] - - return RestrictedPhysicallyMappedElement(element, indices) + return RestrictedPhysicallyMappedElement(element, re._indices) return FiatElement(re)