From 828691b2af6a86173ba646d03652c416900553ab Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 9 Feb 2026 12:12:08 +0000 Subject: [PATCH 01/21] added new methods for computing barycentric coords on tensor-product cells --- FIAT/reference_element.py | 103 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 103 insertions(+) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 87794b50..ff4682d2 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1405,6 +1405,25 @@ def extrinsic_orientation_permutation_map(self): def is_macrocell(self): return any(c.is_macrocell() for c in self.cells) + + def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=None): + """Compute axis-structured barycentric coordinates on a tensor-product cell. + + Returns a list of arrays, one per axis. The i-th entry has shape (npoints, nfacets_axis_i) and + contains the barycentric coordinates associated with facets normal to axis i. + """ + axis_dims = self.get_dimension() + point_slices = TensorProductCell._split_slices(axis_dims) + bary_coord_per_axis = [] + + # Compute barycentric coordinates independently on each axis + for cell, s in zip(self.cells, point_slices): + bary_axis = cell.compute_barycentric_coordinates( + points[..., s], entity=None, rescale=rescale + ) + bary_coord_per_axis.append(bary_axis) + + return bary_coord_per_axis class Hypercube(Cell): @@ -1521,6 +1540,22 @@ def __ge__(self, other): def __le__(self, other): return self.product <= other + def compute_barycentric_coordinates(self, points, entity=None, rescale=False): + """Returns the barycentric coordinates of a list of points on the hypercube.""" + if len(points) == 0: + return points + if entity is not None: + raise NotImplementedError( + "Sub-entity barycentric coordinates are not supported on tensor-product elements." + ) + + tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points) + tp_bary_coords = numpy.hstack(tp_bary_coords) # flatten barycentric coords. + + # Reorder the barycentric coords. in facet order + bary_coords = tp_bary_coords[:, self.facet_perm] + + return bary_coords class UFCHypercube(Hypercube): """Reference UFC Hypercube @@ -1838,6 +1873,65 @@ def compute_unflattening_map(topology_dict): return unflattening_map +def compute_facet_perm(unflattening_map, product): + """ + For a given tensor-product cell, this function returns a mapping between the cell's facets and its barycentric coords. indices. + + For each facet: + - determine which axis the facet is orthogonal to in the tensor-product representation of the cell (i.e., fixed at an endpoint) + - determine which endpoint on that axis the facet corresponds to + - map the endpoint to the correct barycentric coordinate index + """ + # First compute axis offsets into the flattened barycentric vector + axis_offsets = [] + offset = 0 + for axis_cell in product.cells: + axis_offsets.append(offset) + offset += axis_cell.get_dimension() + 1 + + # Initialise the permutation array + sd = len(product.cells) + num_facets = 2 * sd + perm = numpy.zeroes(num_facets, dtype=int) + + for f in range(num_facets): + # Recover the tensor-product representation of the facet + dim_tuple, tp_entity = unflattening_map[(sd - 1, f)] + + # Determine the axis that's orthogonal to the facet + # quad: dim_tuple = (0,1) -> facet is on the boundary of axis 0 (x fixed) -> x = 0 or x = 1 + # dim_tuple = (1,0) -> facet is on the boundary of axis 1 (y fixed) -> y = 0 or y = 1 + axis = next( + i for i, d in enumerate(dim_tuple) + if d == product.cells[i].get_dimension() - 1 + ) + + # Determine the index of the endpoint that produces the facet (local facet number in the axis space) + entity_shape = tuple( + len(c.get_topology()[d]) + for c, d in zip(product.cells, dim_tuple) + ) + tuple_ei = numpy.unravel_index(tp_entity, entity_shape) + local_facet = tuple_ei[axis] + + # Map the local_facet index to the corresponding barycentric coordinate index. + # + # For a 1D interval with endpoints 0 and 1, the barycentric coordinates are: + # lambda_0 = 1 - x + # lambda_1 = x + # + # Each endpoint facet is characterised by one barycentric coordinate vanishing: + # facet at x = 0 -> lambda_1 = 0 (opposite vertex 1) + # facet at x = 1 -> lambda_0 = 0 (opposite vertex 0) + # + # In general, FIAT's convention is: barycentric coordinate i vanishes on the + # facet opposite vertex i. + bary_index = get_bary_index_from_local_facet(axis_cell, local_facet) + + perm[f] = axis_offsets[axis] + bary_index + + return perm + def max_complex(complexes): max_cell = max(complexes) @@ -1845,3 +1939,12 @@ def max_complex(complexes): return max_cell else: raise ValueError("Cannot find the maximal complex") + +def get_bary_index_from_local_facet(simplicial_cell, local_facet): + """Return the barycentric coords. index that vanishes on a given facet. """ + + all_vertices = set(simplicial_cell.get_topology()[0].keys()) + facet_dim = simplicial_cell.get_dimension() - 1 + facet_vertices = set(simplicial_cell.get_topology()[facet_dim][local_facet]) + missing = all_vertices - facet_vertices + return next(iter(missing)) \ No newline at end of file From 85e17ef040c92ce884245c78eff59b25895ae59a Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Tue, 10 Feb 2026 09:37:39 +0000 Subject: [PATCH 02/21] minor fixes --- FIAT/reference_element.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index ff4682d2..13b7b2cf 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1442,6 +1442,8 @@ def __init__(self, dimension, product): self.product = product self.unflattening_map = compute_unflattening_map(pt) + self.facet_perm = compute_facet_perm(self.unflattening_map, self.product) + def get_dimension(self): """Returns the subelement dimension of the cell. Same as the spatial dimension.""" @@ -1892,7 +1894,7 @@ def compute_facet_perm(unflattening_map, product): # Initialise the permutation array sd = len(product.cells) num_facets = 2 * sd - perm = numpy.zeroes(num_facets, dtype=int) + perm = numpy.zeros(num_facets, dtype=int) for f in range(num_facets): # Recover the tensor-product representation of the facet From 0b38a3705094a90a4464c82efefb60b22f56475a Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Wed, 11 Feb 2026 11:21:18 +0000 Subject: [PATCH 03/21] modified compute_barycentric_coords to consistently handle input points of different sizes --- FIAT/reference_element.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 13b7b2cf..fb9d2e73 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -617,6 +617,10 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the complex.""" if len(points) == 0: return points + + if points.ndim == 1: + points = points[None, :] + if entity is None: entity = (self.get_spatial_dimension(), 0) entity_dim, entity_id = entity @@ -1546,6 +1550,12 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the hypercube.""" if len(points) == 0: return points + + # Accept a single point of shape (gdim, ) as well as a batch of points of shape (npoints, gdim) + points = numpy.asarray(points) + if points.ndim == 1: + points = points[None, :] + if entity is not None: raise NotImplementedError( "Sub-entity barycentric coordinates are not supported on tensor-product elements." From 78db95ed018b00330fd1419bff123c5421b3d8db Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 16 Feb 2026 11:44:55 +0000 Subject: [PATCH 04/21] extended the tensor product barycentric coordinates computation function to handle nested tensor-product cells + added tests --- FIAT/reference_element.py | 58 +++++++++++++++++++----- test/FIAT/unit/test_reference_element.py | 35 ++++++++++++++ 2 files changed, 81 insertions(+), 12 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index fb9d2e73..48d6cbfb 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -618,6 +618,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): if len(points) == 0: return points + points = numpy.asarray(points) if points.ndim == 1: points = points[None, :] @@ -1410,21 +1411,55 @@ def extrinsic_orientation_permutation_map(self): def is_macrocell(self): return any(c.is_macrocell() for c in self.cells) - def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=None): + @property + def simplex_cells(self): + """Return a cached flattened list of simplex axis factors. + + Example: + interval_x_interval -> (UFCInterval(), UFCInterval()) + quadrilateral_x_interval -> (UFCInterval(), UFCInterval(), UFCInterval()) + hexahedron -> (UFCInterval(), UFCInterval(), UFCInterval()) + """ + if not hasattr(self, "_simplex_cells"): + def _flatten(cell): + if cell.is_simplex(): + return [cell] + if hasattr(cell, "product"): + return _flatten(cell.product) + if isinstance(cell, TensorProductCell): + out = [] + for sub in cell.cells: + out.extend(_flatten(sub)) + return out + raise NotImplementedError(f"Cannot flatten factor {cell}") + out = [] + for c in self.cells: + out.extend(_flatten(c)) + self._simplex_cells = tuple(out) + return self._simplex_cells + + def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=False): """Compute axis-structured barycentric coordinates on a tensor-product cell. - Returns a list of arrays, one per axis. The i-th entry has shape (npoints, nfacets_axis_i) and - contains the barycentric coordinates associated with facets normal to axis i. + Returns a list of arrays, one per simplicial axis. + The i-th entry has shape (npoints, nfacets_axis_i) and contains the barycentric coordinates + associated with the tensor-product facets normal to axis i. """ - axis_dims = self.get_dimension() + if len(points) == 0: + return points + + points = numpy.asarray(points) + if points.ndim == 1: + points = numpy.asarray(points)[None, :] # convert to (1, dim) shape + + flat_factors = self.simplex_cells + axis_dims = [c.get_spatial_dimension() for c in flat_factors] point_slices = TensorProductCell._split_slices(axis_dims) bary_coord_per_axis = [] - # Compute barycentric coordinates independently on each axis - for cell, s in zip(self.cells, point_slices): - bary_axis = cell.compute_barycentric_coordinates( - points[..., s], entity=None, rescale=rescale - ) + # Compute barycentric coordinates axis-by-axis + for factor, s in zip(flat_factors, point_slices): + bary_axis = factor.compute_barycentric_coordinates(points[..., s], entity, rescale) bary_coord_per_axis.append(bary_axis) return bary_coord_per_axis @@ -1551,7 +1586,6 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): if len(points) == 0: return points - # Accept a single point of shape (gdim, ) as well as a batch of points of shape (npoints, gdim) points = numpy.asarray(points) if points.ndim == 1: points = points[None, :] @@ -1561,7 +1595,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): "Sub-entity barycentric coordinates are not supported on tensor-product elements." ) - tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points) + tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) tp_bary_coords = numpy.hstack(tp_bary_coords) # flatten barycentric coords. # Reorder the barycentric coords. in facet order @@ -1953,7 +1987,7 @@ def max_complex(complexes): raise ValueError("Cannot find the maximal complex") def get_bary_index_from_local_facet(simplicial_cell, local_facet): - """Return the barycentric coords. index that vanishes on a given facet. """ + """Return the barycentric coordinate index that vanishes on a given facet. """ all_vertices = set(simplicial_cell.get_topology()[0].keys()) facet_dim = simplicial_cell.get_dimension() - 1 diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index d166ba6e..736f84a6 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -455,6 +455,41 @@ def test_flatten_maintains_ufc_status(cell): assert ufc_status == is_ufc(flat_cell) +@pytest.mark.parametrize(('cell', 'point'), + [(interval_x_interval, [0.25, 0.6]), + (triangle_x_interval, [0.25, 0.25, 0.5]), + (quadrilateral_x_interval, [0.25, 0.25, 0.5])]) +def test_tensorproduct_barycentrics_match_factor_cells(cell, point): + bary_axes = cell.compute_axis_barycentric_coordinates(point) + + offset = 0 + for k, factor in enumerate(cell.simplex_cells): + sd = factor.get_spatial_dimension() + coords_k = point[offset: offset+sd] # coords on factor cell k + offset += sd + expected = factor.compute_barycentric_coordinates(coords_k) + assert np.allclose(bary_axes[k], expected) + +@pytest.mark.parametrize(('cell', 'point'), + [(quadrilateral, [0.0, 0.3]), + (quadrilateral, [1.0, 0.3]), + (quadrilateral, [0.3, 0.0]), + (quadrilateral, [0.3, 1.0]), + (hexahedron, [0.0, 0.3, 0.4]), + (hexahedron, [1.0, 0.3, 0.4]), + (hexahedron, [0.3, 0.0, 0.4]), + (hexahedron, [0.3, 1.0, 0.4]), + (hexahedron, [0.3, 0.4, 0.0]), + (hexahedron, [0.3, 0.4, 1.0]),]) +def test_hypercube_barycentrics_match_facet_order(cell, point, epsilon=1e-14): + facet_dim = cell.get_spatial_dimension() - 1 + point_entity_ids = cell.point_entity_ids([point]) + facet_hits = [fid for fid, pts in point_entity_ids[facet_dim].items() if len(pts) > 0] + assert len(facet_hits) == 1 + facet_id = facet_hits[0] + bary = cell.compute_barycentric_coordinates(point) + assert np.isclose(bary[0, facet_id], 0.0, atol=epsilon) + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__)) From 4f3d9cd89f436232e9ed5ae1f25e10c8d630dd00 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 16 Feb 2026 15:33:35 +0000 Subject: [PATCH 05/21] removed line enforcing points to be numpy arrays to ensure code generation works properly for symbolic points --- FIAT/reference_element.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 48d6cbfb..31c7cbba 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -617,10 +617,6 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the complex.""" if len(points) == 0: return points - - points = numpy.asarray(points) - if points.ndim == 1: - points = points[None, :] if entity is None: entity = (self.get_spatial_dimension(), 0) @@ -1448,10 +1444,6 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals if len(points) == 0: return points - points = numpy.asarray(points) - if points.ndim == 1: - points = numpy.asarray(points)[None, :] # convert to (1, dim) shape - flat_factors = self.simplex_cells axis_dims = [c.get_spatial_dimension() for c in flat_factors] point_slices = TensorProductCell._split_slices(axis_dims) @@ -1585,10 +1577,6 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the hypercube.""" if len(points) == 0: return points - - points = numpy.asarray(points) - if points.ndim == 1: - points = points[None, :] if entity is not None: raise NotImplementedError( From 4677313af96681ef4264157850978436698c7df7 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 16 Feb 2026 15:48:03 +0000 Subject: [PATCH 06/21] added conversion to numpy array only in axis_barycentric_coords --- FIAT/reference_element.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 31c7cbba..9027a387 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1444,6 +1444,10 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals if len(points) == 0: return points + points = numpy.asarray(points) + if points.ndim == 1: + points = points[None, :] + flat_factors = self.simplex_cells axis_dims = [c.get_spatial_dimension() for c in flat_factors] point_slices = TensorProductCell._split_slices(axis_dims) From c19b8bdce3e7900a8f444ed3c623c8d603aa6e17 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Thu, 19 Feb 2026 14:28:14 +0000 Subject: [PATCH 07/21] generalised numpy operations in compute_barycentric_coordinates methods to work on different array shapes --- FIAT/reference_element.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 9027a387..0d7e4d12 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -615,9 +615,12 @@ def get_dimension(self): def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the complex.""" - if len(points) == 0: + # if len(points) == 0: + # return points + # breakpoint() + if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: return points - + if entity is None: entity = (self.get_spatial_dimension(), 0) entity_dim, entity_id = entity @@ -642,7 +645,8 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): b *= h A *= h[:, None] out = numpy.dot(points, A.T) - return numpy.add(out, b, out=out) + # out = points @ A.T + return numpy.add(out, b) def compute_bubble(self, points, entity=None): """Returns the lowest-order bubble on an entity evaluated at the given @@ -1441,12 +1445,12 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals The i-th entry has shape (npoints, nfacets_axis_i) and contains the barycentric coordinates associated with the tensor-product facets normal to axis i. """ - if len(points) == 0: + if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: return points points = numpy.asarray(points) - if points.ndim == 1: - points = points[None, :] + # if points.ndim == 1: + # points = points[None, :] flat_factors = self.simplex_cells axis_dims = [c.get_spatial_dimension() for c in flat_factors] @@ -1591,7 +1595,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): tp_bary_coords = numpy.hstack(tp_bary_coords) # flatten barycentric coords. # Reorder the barycentric coords. in facet order - bary_coords = tp_bary_coords[:, self.facet_perm] + bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) return bary_coords From ac0d93e9c83471c7a497654548bbe525cd8ea164 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 23 Feb 2026 17:30:06 +0000 Subject: [PATCH 08/21] changed compute_barycentric_coordinates method to work on input GEM expressions representing point sets --- FIAT/reference_element.py | 81 ++++++++++++++++++++++++++------------- 1 file changed, 55 insertions(+), 26 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 0d7e4d12..4f1cbb84 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -617,7 +617,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the complex.""" # if len(points) == 0: # return points - # breakpoint() + if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: return points @@ -644,9 +644,11 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): h = 1 / numpy.linalg.norm(A, axis=1) b *= h A *= h[:, None] - out = numpy.dot(points, A.T) - # out = points @ A.T - return numpy.add(out, b) + # out = numpy.dot(points, A.T) + out = points @ A.T + + # return numpy.add(out, b) + return out + b def compute_bubble(self, points, entity=None): """Returns the lowest-order bubble on an entity evaluated at the given @@ -1413,7 +1415,7 @@ def is_macrocell(self): @property def simplex_cells(self): - """Return a cached flattened list of simplex axis factors. + """Return flattened list of simplex axis factors. Example: interval_x_interval -> (UFCInterval(), UFCInterval()) @@ -1443,23 +1445,27 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals Returns a list of arrays, one per simplicial axis. The i-th entry has shape (npoints, nfacets_axis_i) and contains the barycentric coordinates - associated with the tensor-product facets normal to axis i. + associated with the cell facets normal to axis i. """ + import gem + if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: return points - points = numpy.asarray(points) - # if points.ndim == 1: - # points = points[None, :] + # points = numpy.asarray(points) - flat_factors = self.simplex_cells + # NOTE: introduced a method that recursively finds all simplex factors of the tensor-product cells + flat_factors = self.simplex_cells + axis_dims = [c.get_spatial_dimension() for c in flat_factors] point_slices = TensorProductCell._split_slices(axis_dims) bary_coord_per_axis = [] # Compute barycentric coordinates axis-by-axis for factor, s in zip(flat_factors, point_slices): - bary_axis = factor.compute_barycentric_coordinates(points[..., s], entity, rescale) + # symbolic or numpy slicing + axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] + bary_axis = factor.compute_barycentric_coordinates(axis_points, entity, rescale) bary_coord_per_axis.append(bary_axis) return bary_coord_per_axis @@ -1583,7 +1589,9 @@ def __le__(self, other): def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the hypercube.""" - if len(points) == 0: + import gem + + if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: return points if entity is not None: @@ -1592,11 +1600,28 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): ) tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) - tp_bary_coords = numpy.hstack(tp_bary_coords) # flatten barycentric coords. - # Reorder the barycentric coords. in facet order - bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) + # NOTE: Previous operations (flattening + permuting) were only valid if tp_bary_coords returns a list of arrays + + # # Flatten barycentric coords. + # tp_bary_coords = numpy.hstack(tp_bary_coords) + # # Reorder the barycentric coords. in facet order + # bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) + # We now have self.facet_perm return a tuple (axis, axis_index) instead of an integer permutation + # in the flattened vector + + if isinstance(tp_bary_coords[0], gem.Node): + # breakpoint() + components = [ + gem.Indexed(tp_bary_coords[axis], (bary_index,)) + for axis, bary_index in self.facet_perm + ] + bary_coords = gem.ListTensor(components) # tensor built from scalar GEM expressions + else: + components = [tp_bary_coords[axis][..., bary_index] for axis, bary_index in self.facet_perm] + bary_coords = numpy.stack(components, axis=-1) + return bary_coords class UFCHypercube(Hypercube): @@ -1917,24 +1942,26 @@ def compute_unflattening_map(topology_dict): def compute_facet_perm(unflattening_map, product): """ - For a given tensor-product cell, this function returns a mapping between the cell's facets and its barycentric coords. indices. + For a given tensor-product cell, return a mapping between the cell's facets and the indices + in the axis-structured barycentric coordinate vector. For each facet: - - determine which axis the facet is orthogonal to in the tensor-product representation of the cell (i.e., fixed at an endpoint) - - determine which endpoint on that axis the facet corresponds to + - determine which axis the facet is orthogonal to + - determine which endpoint on that axis the facet corresponds to (facet is fixed at an endpoint) - map the endpoint to the correct barycentric coordinate index """ # First compute axis offsets into the flattened barycentric vector - axis_offsets = [] - offset = 0 - for axis_cell in product.cells: - axis_offsets.append(offset) - offset += axis_cell.get_dimension() + 1 + # axis_offsets = [] + # offset = 0 + # for axis_cell in product.cells: + # axis_offsets.append(offset) + # offset += axis_cell.get_dimension() + 1 # Initialise the permutation array sd = len(product.cells) num_facets = 2 * sd - perm = numpy.zeros(num_facets, dtype=int) + # perm = numpy.zeros(num_facets, dtype=int) + perm = [None] * num_facets for f in range(num_facets): # Recover the tensor-product representation of the facet @@ -1968,9 +1995,11 @@ def compute_facet_perm(unflattening_map, product): # # In general, FIAT's convention is: barycentric coordinate i vanishes on the # facet opposite vertex i. - bary_index = get_bary_index_from_local_facet(axis_cell, local_facet) + bary_index = get_bary_index_from_local_facet(product.cells[axis], local_facet) - perm[f] = axis_offsets[axis] + bary_index + # NOTE: this assumes that the barycentric coordinate vector is flattened + # perm[f] = axis_offsets[axis] + bary_index + perm[f] = (axis, bary_index) return perm From cf47fe9789466b1d03b85d9919acde9ce6aa5bae Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 23 Feb 2026 17:31:32 +0000 Subject: [PATCH 09/21] extended gem matmul to convert numpy arrays to Literals and deal separately with singleton contractions --- gem/gem.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/gem/gem.py b/gem/gem.py index 02aeaaf5..b70ebf0d 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -116,10 +116,19 @@ def __matmul__(self, other): raise ValueError("Both objects must have shape for matmul") elif self.shape[-1] != other.shape[0]: raise ValueError(f"Mismatching shapes {self.shape} and {other.shape} in matmul") - *i, k = indices(len(self.shape)) - _, *j = indices(len(other.shape)) - expr = Product(Indexed(self, (*i, k)), Indexed(other, (k, *j))) - return ComponentTensor(IndexSum(expr, (k, )), (*i, *j)) + + if self.shape[-1] == 1: + # Avoid generation of contraction loops over singleton dimensions (extent-1 indices) + # self: (*i, 1), other (1, *j) => result: (*i, *j) + i = indices(len(self.shape) - 1) + j = indices(len(other.shape) - 1) + expr = Product(Indexed(self, (*i, 0)), Indexed(other, (0, *j))) + return ComponentTensor(expr, (*i, *j)) + else: + *i, k = indices(len(self.shape)) + _, *j = indices(len(other.shape)) + expr = Product(Indexed(self, (*i, k)), Indexed(other, (k, *j))) + return ComponentTensor(IndexSum(expr, (k, )), (*i, *j)) def __rmatmul__(self, other): return as_gem(other).__matmul__(self) @@ -1324,7 +1333,7 @@ def as_gem(expr): """ if isinstance(expr, Node): return expr - elif isinstance(expr, Number): + elif isinstance(expr, Number) or isinstance(expr, numpy.ndarray): return Literal(expr) else: raise ValueError("Do not know how to convert %r to GEM" % expr) From c04b6fd361b28bfe495b956912466e616d39826a Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Wed, 25 Feb 2026 18:53:19 +0000 Subject: [PATCH 10/21] changes made to gem and FIAT to make barycentric coordinate computations work on GEM tensor expressions --- FIAT/reference_element.py | 9 +++++++-- gem/gem.py | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 4f1cbb84..b3086701 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1464,10 +1464,13 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals # Compute barycentric coordinates axis-by-axis for factor, s in zip(flat_factors, point_slices): # symbolic or numpy slicing + # TODO: implement slicing in GEM + # implement __getitem__ in GEM + # does gem.view() slice on the last dimension by default? what if we want to slice on multiple dimensions? axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] bary_axis = factor.compute_barycentric_coordinates(axis_points, entity, rescale) bary_coord_per_axis.append(bary_axis) - + return bary_coord_per_axis @@ -1605,12 +1608,14 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): # # Flatten barycentric coords. # tp_bary_coords = numpy.hstack(tp_bary_coords) + # # Reorder the barycentric coords. in facet order + # bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) # We now have self.facet_perm return a tuple (axis, axis_index) instead of an integer permutation # in the flattened vector - + if isinstance(tp_bary_coords[0], gem.Node): # breakpoint() components = [ diff --git a/gem/gem.py b/gem/gem.py index b70ebf0d..5e5655ac 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -124,6 +124,7 @@ def __matmul__(self, other): j = indices(len(other.shape) - 1) expr = Product(Indexed(self, (*i, 0)), Indexed(other, (0, *j))) return ComponentTensor(expr, (*i, *j)) + else: *i, k = indices(len(self.shape)) _, *j = indices(len(other.shape)) From f339b0195a426789a797f395877af0873181eace Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Thu, 26 Feb 2026 14:50:23 +0000 Subject: [PATCH 11/21] extended gem's slicing syntax and simplified the code in compute_axis_barycentric_coordinates in tp cells --- FIAT/reference_element.py | 62 +++++++++++++++++++++------------------ gem/gem.py | 57 ++++++++++++++++++++++++++++------- 2 files changed, 80 insertions(+), 39 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index b3086701..2f552cfe 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1443,31 +1443,37 @@ def _flatten(cell): def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=False): """Compute axis-structured barycentric coordinates on a tensor-product cell. - Returns a list of arrays, one per simplicial axis. - The i-th entry has shape (npoints, nfacets_axis_i) and contains the barycentric coordinates - associated with the cell facets normal to axis i. - """ - import gem + Parameters + ---------- + points: numpy.ndarray or GEM.Node + The reference coordinates of the points. - if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: + NOTE: + Should we support computing barycentric coordinates on cell entities? + Should we support rescaling the barycentric coordinates? + + Returns + ------- + List of numpy.ndarray or GEM.ComponentTensor + Returns a list of barycentric coordinates computed on each simplicial axis of the tensor-product cell. + The i-th array has shape (npoints, nfacets_axis_i) and contains the barycentric coordinates + associated with the cell facets normal to axis i. + """ + if isinstance(points, numpy.ndarray) and len(points) == 0: return points - # points = numpy.asarray(points) - - # NOTE: introduced a method that recursively finds all simplex factors of the tensor-product cells + # Find all simplex factors of the tensor-product cells flat_factors = self.simplex_cells axis_dims = [c.get_spatial_dimension() for c in flat_factors] point_slices = TensorProductCell._split_slices(axis_dims) bary_coord_per_axis = [] - # Compute barycentric coordinates axis-by-axis + # Compute barycentric coordinates on each axis for factor, s in zip(flat_factors, point_slices): # symbolic or numpy slicing - # TODO: implement slicing in GEM - # implement __getitem__ in GEM - # does gem.view() slice on the last dimension by default? what if we want to slice on multiple dimensions? - axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] + # axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] + axis_points = points[..., s] bary_axis = factor.compute_barycentric_coordinates(axis_points, entity, rescale) bary_coord_per_axis.append(bary_axis) @@ -1599,30 +1605,30 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): if entity is not None: raise NotImplementedError( - "Sub-entity barycentric coordinates are not supported on tensor-product elements." + "Sub-entity barycentric coordinates are not supported on hypercubes." ) tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) # NOTE: Previous operations (flattening + permuting) were only valid if tp_bary_coords returns a list of arrays + """ + # Flatten barycentric coords. + tp_bary_coords = numpy.hstack(tp_bary_coords) - # # Flatten barycentric coords. - # tp_bary_coords = numpy.hstack(tp_bary_coords) - - # # Reorder the barycentric coords. in facet order - - # bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) - - # We now have self.facet_perm return a tuple (axis, axis_index) instead of an integer permutation - # in the flattened vector - + # Reorder the barycentric coords. in facet order + bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) + """ + # We now have `self.facet_perm` return a list of tuple (axis, axis_index) instead of an integer permutation + # describing for each local facet `lf` which axis of the TP cell it corresponds to and which barycentric coordinate + # on that axis vanishes on the said facet. + if isinstance(tp_bary_coords[0], gem.Node): # breakpoint() components = [ gem.Indexed(tp_bary_coords[axis], (bary_index,)) for axis, bary_index in self.facet_perm ] - bary_coords = gem.ListTensor(components) # tensor built from scalar GEM expressions + bary_coords = gem.ListTensor(components) else: components = [tp_bary_coords[axis][..., bary_index] for axis, bary_index in self.facet_perm] bary_coords = numpy.stack(components, axis=-1) @@ -1947,7 +1953,7 @@ def compute_unflattening_map(topology_dict): def compute_facet_perm(unflattening_map, product): """ - For a given tensor-product cell, return a mapping between the cell's facets and the indices + For a given hypercube, return a mapping between the cell's facets and the indices in the axis-structured barycentric coordinate vector. For each facet: @@ -2017,7 +2023,7 @@ def max_complex(complexes): raise ValueError("Cannot find the maximal complex") def get_bary_index_from_local_facet(simplicial_cell, local_facet): - """Return the barycentric coordinate index that vanishes on a given facet. """ + """Return the barycentric coordinate index that vanishes on a given facet of a simplicial cell. """ all_vertices = set(simplicial_cell.get_topology()[0].keys()) facet_dim = simplicial_cell.get_dimension() - 1 diff --git a/gem/gem.py b/gem/gem.py index 5e5655ac..0ded808f 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -80,13 +80,49 @@ def is_equal(self, other): if result: self.children = other.children return result + + def __getitem__(self, key): + """ + Generalised indexing interface. - def __getitem__(self, indices): - try: - indices = tuple(indices) - except TypeError: - indices = (indices, ) - return Indexed(self, indices) + Parameters + ---------- + key: int, Index, VariableIndex, slice, Ellipsis or tuple + """ + # Normalize to tuple so we can inspect the key + if not isinstance(key, tuple): + key = (key,) + + # Expand ellipsis -> fill remaining dimensions with slice(None) + if Ellipsis in key: + if key.count(Ellipsis) > 1: + raise NotImplementedError("Multiple ellipses are not supported when indexing a gem.Node tensor") + ellipsis_pos = key.index(Ellipsis) + num_missing = len(self.shape) - (len(key) - 1) + if num_missing < 0: + raise IndexError("Too many indices provided when indexing the gem.Node tensor") + key = ( + key[:ellipsis_pos] + + (slice(None), ) * num_missing + + key[ellipsis_pos + 1:] + ) + # Slice indexing -> delegate to view() + if any(isinstance(k, slice) for k in key): + # view expects one slice for each axis/dim of the tensor + if len(key) != len(self.shape): + raise IndexError("Slice indexing expects the number of slices to match the gem.Node tensor rank") + return view(self, *key) + # Point indexing -> delegate to Index + else: + return Indexed(self, key) + + # OLD + # def __getitem__(self, indices): + # try: + # indices = tuple(indices) + # except TypeError: + # indices = (indices, ) + # return Indexed(self, indices) def __add__(self, other): return componentwise(Sum, self, as_gem(other)) @@ -119,11 +155,10 @@ def __matmul__(self, other): if self.shape[-1] == 1: # Avoid generation of contraction loops over singleton dimensions (extent-1 indices) - # self: (*i, 1), other (1, *j) => result: (*i, *j) - i = indices(len(self.shape) - 1) - j = indices(len(other.shape) - 1) - expr = Product(Indexed(self, (*i, 0)), Indexed(other, (0, *j))) - return ComponentTensor(expr, (*i, *j)) + i = indices(len(self.shape) - 1) # non-contracted axes of self + j = indices(len(other.shape) - 1) # non-contracted axes of other + expr = Product(Indexed(self, (*i, 0)), Indexed(other, (0, *j))) # scalar expr with free indices i and j + return ComponentTensor(expr, (*i, *j)) else: *i, k = indices(len(self.shape)) From ddf6eb411f13d5fc3ecaba53006c3f43353255e1 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Thu, 26 Feb 2026 14:58:32 +0000 Subject: [PATCH 12/21] tidied up and added comments --- FIAT/reference_element.py | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 2f552cfe..cc5841d8 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1415,7 +1415,7 @@ def is_macrocell(self): @property def simplex_cells(self): - """Return flattened list of simplex axis factors. + """Return a flattened list of simplex axis factors. Example: interval_x_interval -> (UFCInterval(), UFCInterval()) @@ -1597,10 +1597,22 @@ def __le__(self, other): return self.product <= other def compute_barycentric_coordinates(self, points, entity=None, rescale=False): - """Returns the barycentric coordinates of a list of points on the hypercube.""" + """Returns the barycentric coordinates of a list of points on the hypercube. + + Parameters + ---------- + points: numpy.ndarray or GEM.Node + The reference coordinates of the points. + + Returns + ------- + List of numpy.ndarray or GEM.ComponentTensor + Returns a list of barycentric coordinates in local facet order such that for any point + lying on local facet `lf` of the cell, the barycentric coordinate at index `lf` vanishes. + """ import gem - if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: + if isinstance(points, numpy.ndarray) and len(points) == 0: return points if entity is not None: @@ -1610,18 +1622,18 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) - # NOTE: Previous operations (flattening + permuting) were only valid if tp_bary_coords returns a list of arrays - """ + # Previous operations involved numpy and are only valid if tp_bary_coords returns a list of arrays # Flatten barycentric coords. - tp_bary_coords = numpy.hstack(tp_bary_coords) + #tp_bary_coords = numpy.hstack(tp_bary_coords) # Reorder the barycentric coords. in facet order - bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) - """ + #bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) + # We now have `self.facet_perm` return a list of tuple (axis, axis_index) instead of an integer permutation # describing for each local facet `lf` which axis of the TP cell it corresponds to and which barycentric coordinate # on that axis vanishes on the said facet. - + + # We use `self.facet_perm` to reorder the barycentric coordinates in facet order if isinstance(tp_bary_coords[0], gem.Node): # breakpoint() components = [ From 7b35ca3dabf271cdde2dc855913ba527c0b7b8bb Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 20 Apr 2026 12:28:18 +0100 Subject: [PATCH 13/21] recent changes post merging --- FIAT/reference_element.py | 71 +++---- firedrake_fiat.egg-info/PKG-INFO | 75 +++++++ firedrake_fiat.egg-info/SOURCES.txt | 212 +++++++++++++++++++ firedrake_fiat.egg-info/dependency_links.txt | 1 + firedrake_fiat.egg-info/requires.txt | 13 ++ firedrake_fiat.egg-info/top_level.txt | 3 + gem/gem.py | 16 +- 7 files changed, 340 insertions(+), 51 deletions(-) create mode 100644 firedrake_fiat.egg-info/PKG-INFO create mode 100644 firedrake_fiat.egg-info/SOURCES.txt create mode 100644 firedrake_fiat.egg-info/dependency_links.txt create mode 100644 firedrake_fiat.egg-info/requires.txt create mode 100644 firedrake_fiat.egg-info/top_level.txt diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index cc5841d8..65fbff95 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1413,45 +1413,41 @@ def extrinsic_orientation_permutation_map(self): def is_macrocell(self): return any(c.is_macrocell() for c in self.cells) - @property - def simplex_cells(self): - """Return a flattened list of simplex axis factors. + # @property + # def simplex_cells(self): + # """Return a flattened list of simplex axis factors. - Example: - interval_x_interval -> (UFCInterval(), UFCInterval()) - quadrilateral_x_interval -> (UFCInterval(), UFCInterval(), UFCInterval()) - hexahedron -> (UFCInterval(), UFCInterval(), UFCInterval()) - """ - if not hasattr(self, "_simplex_cells"): - def _flatten(cell): - if cell.is_simplex(): - return [cell] - if hasattr(cell, "product"): - return _flatten(cell.product) - if isinstance(cell, TensorProductCell): - out = [] - for sub in cell.cells: - out.extend(_flatten(sub)) - return out - raise NotImplementedError(f"Cannot flatten factor {cell}") - out = [] - for c in self.cells: - out.extend(_flatten(c)) - self._simplex_cells = tuple(out) - return self._simplex_cells + # Example: + # interval_x_interval -> (UFCInterval(), UFCInterval()) + # quadrilateral_x_interval -> (UFCInterval(), UFCInterval(), UFCInterval()) + # hexahedron -> (UFCInterval(), UFCInterval(), UFCInterval()) + # """ + # if not hasattr(self, "_simplex_cells"): + # def _flatten(cell): + # if cell.is_simplex(): + # return [cell] + # if hasattr(cell, "product"): + # return _flatten(cell.product) + # if isinstance(cell, TensorProductCell): + # out = [] + # for sub in cell.cells: + # out.extend(_flatten(sub)) + # return out + # raise NotImplementedError(f"Cannot flatten factor {cell}") + # out = [] + # for c in self.cells: + # out.extend(_flatten(c)) + # self._simplex_cells = tuple(out) + # return self._simplex_cells def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=False): - """Compute axis-structured barycentric coordinates on a tensor-product cell. + """Compute barycentric coordinates on each axis (factor) of a tensor-product cell. Parameters ---------- points: numpy.ndarray or GEM.Node The reference coordinates of the points. - NOTE: - Should we support computing barycentric coordinates on cell entities? - Should we support rescaling the barycentric coordinates? - Returns ------- List of numpy.ndarray or GEM.ComponentTensor @@ -1462,17 +1458,14 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals if isinstance(points, numpy.ndarray) and len(points) == 0: return points - # Find all simplex factors of the tensor-product cells - flat_factors = self.simplex_cells - - axis_dims = [c.get_spatial_dimension() for c in flat_factors] + axis_dims = [c.get_spatial_dimension() for c in self.cells] point_slices = TensorProductCell._split_slices(axis_dims) bary_coord_per_axis = [] - # Compute barycentric coordinates on each axis - for factor, s in zip(flat_factors, point_slices): - # symbolic or numpy slicing - # axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] + # Compute barycentric coordinates on each factor + for factor, s in zip(self.cells, point_slices): + # support symbolic or numpy slicing in the same way + # previously did: axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] axis_points = points[..., s] bary_axis = factor.compute_barycentric_coordinates(axis_points, entity, rescale) bary_coord_per_axis.append(bary_axis) @@ -1633,7 +1626,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): # describing for each local facet `lf` which axis of the TP cell it corresponds to and which barycentric coordinate # on that axis vanishes on the said facet. - # We use `self.facet_perm` to reorder the barycentric coordinates in facet order + # We use `self.facet_perm` to reorder the barycentric coordinates in facet order if isinstance(tp_bary_coords[0], gem.Node): # breakpoint() components = [ diff --git a/firedrake_fiat.egg-info/PKG-INFO b/firedrake_fiat.egg-info/PKG-INFO new file mode 100644 index 00000000..1eaca9f4 --- /dev/null +++ b/firedrake_fiat.egg-info/PKG-INFO @@ -0,0 +1,75 @@ +Metadata-Version: 2.4 +Name: firedrake-fiat +Version: 2025.11.0.dev0 +Summary: FInite element Automatic Tabulator +Author-email: "Robert C. Kirby et al." , Imperial College London and others +Project-URL: Repository, https://github.com/firedrakeproject/fiat.git +Classifier: Programming Language :: Python +Requires-Python: >=3.10 +Description-Content-Type: text/x-rst +License-File: COPYING +License-File: COPYING.LESSER +License-File: AUTHORS +Requires-Dist: fenics-ufl @ git+https://github.com/FEniCS/ufl.git@main +Requires-Dist: numpy>=1.16 +Requires-Dist: recursivenodes +Requires-Dist: scipy +Requires-Dist: symengine +Requires-Dist: sympy +Provides-Extra: doc +Requires-Dist: setuptools; extra == "doc" +Requires-Dist: sphinx; extra == "doc" +Provides-Extra: test +Requires-Dist: pytest; extra == "test" +Dynamic: license-file + +======================================== +FIAT: FInite element Automatic Tabulator +======================================== + +The FInite element Automatic Tabulator FIAT supports generation of +arbitrary order instances of the Lagrange elements on lines, +triangles, and tetrahedra. It is also capable of generating arbitrary +order instances of Jacobi-type quadrature rules on the same element +shapes. Further, H(div) and H(curl) conforming finite element spaces +such as the families of Raviart-Thomas, Brezzi-Douglas-Marini and +Nedelec are supported on triangles and tetrahedra. Upcoming versions +will also support Hermite and nonconforming elements. + +FIAT is part of the FEniCS Project. + +For more information, visit http://www.fenicsproject.org + + +.. image:: https://github.com/FEniCS/fiat/workflows/FIAT%20CI/badge.svg + :target: https://github.com/FEniCS/fiat/actions?query=workflow%3A%22FIAT+CI%22 + :alt: Build Status + +.. image:: https://coveralls.io/repos/github/FEniCS/fiat/badge.svg?branch=main + :target: https://coveralls.io/github/FEniCS/fiat?branch=main + :alt: Coverage Status + +.. image:: https://readthedocs.org/projects/fenics-fiat/badge/?version=latest + :target: http://fenics.readthedocs.io/projects/fiat/en/latest/?badge=latest + :alt: Documentation Status + +Documentation +============= + +Documentation can be viewed at http://fenics-fiat.readthedocs.org/. + +License +======= + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program. If not, see . diff --git a/firedrake_fiat.egg-info/SOURCES.txt b/firedrake_fiat.egg-info/SOURCES.txt new file mode 100644 index 00000000..e1867ace --- /dev/null +++ b/firedrake_fiat.egg-info/SOURCES.txt @@ -0,0 +1,212 @@ +AUTHORS +COPYING +COPYING.LESSER +MANIFEST.in +README.rst +pyproject.toml +setup.cfg +FIAT/P0.py +FIAT/Sminus.py +FIAT/SminusCurl.py +FIAT/SminusDiv.py +FIAT/__init__.py +FIAT/alfeld_sorokina.py +FIAT/argyris.py +FIAT/arnold_qin.py +FIAT/arnold_winther.py +FIAT/barycentric_interpolation.py +FIAT/bell.py +FIAT/bernardi_raugel.py +FIAT/bernstein.py +FIAT/brezzi_douglas_fortin_marini.py +FIAT/brezzi_douglas_marini.py +FIAT/brezzi_douglas_marini_cube.py +FIAT/bubble.py +FIAT/c2_elements.py +FIAT/check_format_variant.py +FIAT/christiansen_hu.py +FIAT/crouzeix_raviart.py +FIAT/discontinuous.py +FIAT/discontinuous_lagrange.py +FIAT/discontinuous_pc.py +FIAT/discontinuous_raviart_thomas.py +FIAT/discontinuous_taylor.py +FIAT/dual_set.py +FIAT/enriched.py +FIAT/expansions.py +FIAT/fdm_element.py +FIAT/finite_element.py +FIAT/functional.py +FIAT/gauss_legendre.py +FIAT/gauss_lobatto_legendre.py +FIAT/gauss_radau.py +FIAT/gopalakrishnan_lederer_schoberl.py +FIAT/guzman_neilan.py +FIAT/hct.py +FIAT/hdiv_trace.py +FIAT/hdivcurl.py +FIAT/hellan_herrmann_johnson.py +FIAT/hermite.py +FIAT/hierarchical.py +FIAT/histopolation.py +FIAT/hu_zhang.py +FIAT/jacobi.py +FIAT/johnson_mercier.py +FIAT/kong_mulder_veldhuizen.py +FIAT/lagrange.py +FIAT/macro.py +FIAT/mardal_tai_winther.py +FIAT/mixed.py +FIAT/morley.py +FIAT/nedelec.py +FIAT/nedelec_second_kind.py +FIAT/nodal_enriched.py +FIAT/orientation_utils.py +FIAT/orthopoly.py +FIAT/pointwise_dual.py +FIAT/polynomial_set.py +FIAT/powell_sabin.py +FIAT/quadrature.py +FIAT/quadrature_element.py +FIAT/quadrature_schemes.py +FIAT/raviart_thomas.py +FIAT/reference_element.py +FIAT/regge.py +FIAT/restricted.py +FIAT/serendipity.py +FIAT/tensor_product.py +FIAT/walkington.py +FIAT/wuxu.py +FIAT/xg_quad_data.py +finat/__init__.py +finat/alfeld_sorokina.py +finat/argyris.py +finat/arnold_qin.py +finat/aw.py +finat/bell.py +finat/bernardi_raugel.py +finat/c2_elements.py +finat/cell_tools.py +finat/christiansen_hu.py +finat/citations.py +finat/cube.py +finat/direct_serendipity.py +finat/discontinuous.py +finat/element_factory.py +finat/enriched.py +finat/fiat_elements.py +finat/finiteelementbase.py +finat/guzman_neilan.py +finat/hct.py +finat/hdivcurl.py +finat/hermite.py +finat/hz.py +finat/johnson_mercier.py +finat/mixed.py +finat/morley.py +finat/mtw.py +finat/nodal_enriched.py +finat/physically_mapped.py +finat/piola_mapped.py +finat/point_set.py +finat/powell_sabin.py +finat/quadrature.py +finat/quadrature_element.py +finat/restricted.py +finat/runtime_tabulated.py +finat/spectral.py +finat/sympy2gem.py +finat/tensor_product.py +finat/tensorfiniteelement.py +finat/walkington.py +finat/wuxu.py +finat/ufl/__init__.py +finat/ufl/brokenelement.py +finat/ufl/elementlist.py +finat/ufl/enrichedelement.py +finat/ufl/finiteelement.py +finat/ufl/finiteelementbase.py +finat/ufl/hdivcurl.py +finat/ufl/mixedelement.py +finat/ufl/restrictedelement.py +finat/ufl/tensorproductelement.py +firedrake_fiat.egg-info/PKG-INFO +firedrake_fiat.egg-info/SOURCES.txt +firedrake_fiat.egg-info/dependency_links.txt +firedrake_fiat.egg-info/requires.txt +firedrake_fiat.egg-info/top_level.txt +gem/__init__.py +gem/coffee.py +gem/flop_count.py +gem/gem.py +gem/impero.py +gem/impero_utils.py +gem/interpreter.py +gem/node.py +gem/optimise.py +gem/pprint.py +gem/refactorise.py +gem/scheduling.py +gem/unconcatenate.py +gem/utils.py +test/conftest.py +test/FIAT/regression/.gitignore +test/FIAT/regression/README.rst +test/FIAT/regression/conftest.py +test/FIAT/regression/fiat-reference-data-id +test/FIAT/regression/test_regression.py +test/FIAT/regression/scripts/download +test/FIAT/regression/scripts/getdata +test/FIAT/regression/scripts/getreferencerepo +test/FIAT/regression/scripts/parameters +test/FIAT/regression/scripts/upload +test/FIAT/unit/test_argyris.py +test/FIAT/unit/test_awc.py +test/FIAT/unit/test_awnc.py +test/FIAT/unit/test_bdfm.py +test/FIAT/unit/test_bernstein.py +test/FIAT/unit/test_discontinuous_pc.py +test/FIAT/unit/test_discontinuous_taylor.py +test/FIAT/unit/test_facet_support_dofs.py +test/FIAT/unit/test_fdm.py +test/FIAT/unit/test_fiat.py +test/FIAT/unit/test_gauss_legendre.py +test/FIAT/unit/test_gauss_lobatto_legendre.py +test/FIAT/unit/test_gauss_radau.py +test/FIAT/unit/test_gopalakrishnan_lederer_schoberl.py +test/FIAT/unit/test_hct.py +test/FIAT/unit/test_hdivtrace.py +test/FIAT/unit/test_hierarchical.py +test/FIAT/unit/test_histopolation.py +test/FIAT/unit/test_johnson_mercier.py +test/FIAT/unit/test_kong_mulder_veldhuizen.py +test/FIAT/unit/test_macro.py +test/FIAT/unit/test_mtw.py +test/FIAT/unit/test_orientation.py +test/FIAT/unit/test_pointwise_dual.py +test/FIAT/unit/test_polynomial.py +test/FIAT/unit/test_powell_sabin.py +test/FIAT/unit/test_quadrature.py +test/FIAT/unit/test_quadrature_element.py +test/FIAT/unit/test_reference_element.py +test/FIAT/unit/test_regge_hhj.py +test/FIAT/unit/test_serendipity.py +test/FIAT/unit/test_stokes_complex.py +test/FIAT/unit/test_tensor_product.py +test/FIAT/unit/test_walkington.py +test/finat/conftest.py +test/finat/test_create_broken_element.py +test/finat/test_create_fiat_element.py +test/finat/test_create_finat_element.py +test/finat/test_create_restricted_element.py +test/finat/test_direct_serendipity.py +test/finat/test_dual_basis.py +test/finat/test_hash.py +test/finat/test_mass_conditioning.py +test/finat/test_point_evaluation.py +test/finat/test_quadrature.py +test/finat/test_quadrature_element.py +test/finat/test_restriction.py +test/finat/test_ufl_elements.py +test/finat/test_zany_mapping.py +test/gem/test_simplify.py \ No newline at end of file diff --git a/firedrake_fiat.egg-info/dependency_links.txt b/firedrake_fiat.egg-info/dependency_links.txt new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/firedrake_fiat.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/firedrake_fiat.egg-info/requires.txt b/firedrake_fiat.egg-info/requires.txt new file mode 100644 index 00000000..08d00a07 --- /dev/null +++ b/firedrake_fiat.egg-info/requires.txt @@ -0,0 +1,13 @@ +fenics-ufl @ git+https://github.com/FEniCS/ufl.git@main +numpy>=1.16 +recursivenodes +scipy +symengine +sympy + +[doc] +setuptools +sphinx + +[test] +pytest diff --git a/firedrake_fiat.egg-info/top_level.txt b/firedrake_fiat.egg-info/top_level.txt new file mode 100644 index 00000000..6bc74dc3 --- /dev/null +++ b/firedrake_fiat.egg-info/top_level.txt @@ -0,0 +1,3 @@ +FIAT +finat +gem diff --git a/gem/gem.py b/gem/gem.py index fc92d78a..070498fb 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -154,18 +154,10 @@ def __matmul__(self, other): elif self.shape[-1] != other.shape[0]: raise ValueError(f"Mismatching shapes {self.shape} and {other.shape} in matmul") - if self.shape[-1] == 1: - # Avoid generation of contraction loops over singleton dimensions (extent-1 indices) - i = indices(len(self.shape) - 1) # non-contracted axes of self - j = indices(len(other.shape) - 1) # non-contracted axes of other - expr = Product(Indexed(self, (*i, 0)), Indexed(other, (0, *j))) # scalar expr with free indices i and j - return ComponentTensor(expr, (*i, *j)) - - else: - *i, k = indices(len(self.shape)) - _, *j = indices(len(other.shape)) - expr = Product(Indexed(self, (*i, k)), Indexed(other, (k, *j))) - return ComponentTensor(IndexSum(expr, (k, )), (*i, *j)) + *i, k = indices(len(self.shape)) + _, *j = indices(len(other.shape)) + expr = Product(Indexed(self, (*i, k)), Indexed(other, (k, *j))) + return ComponentTensor(IndexSum(expr, (k, )), (*i, *j)) def __rmatmul__(self, other): return as_gem(other).__matmul__(self) From 5919d6148133d8d0eb79b7ad0e80b50a7c2d0625 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 20 Apr 2026 15:16:15 +0100 Subject: [PATCH 14/21] compute barycentric coords symbolically in simplicies and hypercubes (to be tested) --- FIAT/reference_element.py | 156 +++++++++----------------------------- gem/gem.py | 35 ++++----- 2 files changed, 50 insertions(+), 141 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 65fbff95..5028be94 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1412,33 +1412,6 @@ def extrinsic_orientation_permutation_map(self): def is_macrocell(self): return any(c.is_macrocell() for c in self.cells) - - # @property - # def simplex_cells(self): - # """Return a flattened list of simplex axis factors. - - # Example: - # interval_x_interval -> (UFCInterval(), UFCInterval()) - # quadrilateral_x_interval -> (UFCInterval(), UFCInterval(), UFCInterval()) - # hexahedron -> (UFCInterval(), UFCInterval(), UFCInterval()) - # """ - # if not hasattr(self, "_simplex_cells"): - # def _flatten(cell): - # if cell.is_simplex(): - # return [cell] - # if hasattr(cell, "product"): - # return _flatten(cell.product) - # if isinstance(cell, TensorProductCell): - # out = [] - # for sub in cell.cells: - # out.extend(_flatten(sub)) - # return out - # raise NotImplementedError(f"Cannot flatten factor {cell}") - # out = [] - # for c in self.cells: - # out.extend(_flatten(c)) - # self._simplex_cells = tuple(out) - # return self._simplex_cells def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=False): """Compute barycentric coordinates on each axis (factor) of a tensor-product cell. @@ -1450,28 +1423,23 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals Returns ------- - List of numpy.ndarray or GEM.ComponentTensor - Returns a list of barycentric coordinates computed on each simplicial axis of the tensor-product cell. - The i-th array has shape (npoints, nfacets_axis_i) and contains the barycentric coordinates - associated with the cell facets normal to axis i. + numpy.ndarray + An array of shape ``(nfactors,)`` and dtype object if points are GEM nodes, + otherwise dtype numeric. The i-th entry contains the barycentric coordinates + on the i-th factor cell. If factor i is a simplex of dimension d, this will + have shape ``(npoints, d+1)``. If factor i is a hypercube of dimension d, + this will have shape ``(npoints, 2*d)``. """ if isinstance(points, numpy.ndarray) and len(points) == 0: return points axis_dims = [c.get_spatial_dimension() for c in self.cells] point_slices = TensorProductCell._split_slices(axis_dims) - bary_coord_per_axis = [] - - # Compute barycentric coordinates on each factor - for factor, s in zip(self.cells, point_slices): - # support symbolic or numpy slicing in the same way - # previously did: axis_points = gem.view(points, s) if isinstance(points, gem.Node) else numpy.asarray(points)[...,s] - axis_points = points[..., s] - bary_axis = factor.compute_barycentric_coordinates(axis_points, entity, rescale) - bary_coord_per_axis.append(bary_axis) - - return bary_coord_per_axis + return numpy.array([ + factor.compute_barycentric_coordinates(points[..., s], entity, rescale) + for factor, s in zip(self.cells, point_slices) + ]) class Hypercube(Cell): """Abstract class for a reference hypercube""" @@ -1489,7 +1457,7 @@ def __init__(self, dimension, product): self.product = product self.unflattening_map = compute_unflattening_map(pt) - self.facet_perm = compute_facet_perm(self.unflattening_map, self.product) + self.facet_perm = compute_facet_permutation(self.unflattening_map, self.product) def get_dimension(self): """Returns the subelement dimension of the cell. Same as the @@ -1603,42 +1571,12 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): Returns a list of barycentric coordinates in local facet order such that for any point lying on local facet `lf` of the cell, the barycentric coordinate at index `lf` vanishes. """ - import gem - if isinstance(points, numpy.ndarray) and len(points) == 0: return points - if entity is not None: - raise NotImplementedError( - "Sub-entity barycentric coordinates are not supported on hypercubes." - ) - tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) - - # Previous operations involved numpy and are only valid if tp_bary_coords returns a list of arrays - # Flatten barycentric coords. - #tp_bary_coords = numpy.hstack(tp_bary_coords) - - # Reorder the barycentric coords. in facet order - #bary_coords = numpy.take(tp_bary_coords, self.facet_perm, axis=-1) - - # We now have `self.facet_perm` return a list of tuple (axis, axis_index) instead of an integer permutation - # describing for each local facet `lf` which axis of the TP cell it corresponds to and which barycentric coordinate - # on that axis vanishes on the said facet. - - # We use `self.facet_perm` to reorder the barycentric coordinates in facet order - if isinstance(tp_bary_coords[0], gem.Node): - # breakpoint() - components = [ - gem.Indexed(tp_bary_coords[axis], (bary_index,)) - for axis, bary_index in self.facet_perm - ] - bary_coords = gem.ListTensor(components) - else: - components = [tp_bary_coords[axis][..., bary_index] for axis, bary_index in self.facet_perm] - bary_coords = numpy.stack(components, axis=-1) - return bary_coords + return tp_bary_coords[..., self.facet_perm] class UFCHypercube(Hypercube): """Reference UFC Hypercube @@ -1956,42 +1894,38 @@ def compute_unflattening_map(topology_dict): return unflattening_map -def compute_facet_perm(unflattening_map, product): +def compute_facet_permutation(unflattening_map, product): """ - For a given hypercube, return a mapping between the cell's facets and the indices - in the axis-structured barycentric coordinate vector. - - For each facet: - - determine which axis the facet is orthogonal to - - determine which endpoint on that axis the facet corresponds to (facet is fixed at an endpoint) - - map the endpoint to the correct barycentric coordinate index + Return a permutation mapping each hypercube facet to the index of its + vanishing barycentric coordinate in the axis-structured barycentric array. """ - # First compute axis offsets into the flattened barycentric vector - # axis_offsets = [] - # offset = 0 - # for axis_cell in product.cells: - # axis_offsets.append(offset) - # offset += axis_cell.get_dimension() + 1 + # First compute axis offsets into the flattened barycentric coordinate array. + axis_offsets = [] + offset = 0 + for axis_cell in product.cells: + axis_offsets.append(offset) + offset += axis_cell.get_dimension() + 1 - # Initialise the permutation array + # Initialise the integer permutation array sd = len(product.cells) num_facets = 2 * sd - # perm = numpy.zeros(num_facets, dtype=int) - perm = [None] * num_facets + perm = numpy.zeroes(num_facets, dtype=int) for f in range(num_facets): - # Recover the tensor-product representation of the facet + # Recover the tensor-product representation of the facet as given by the unflattening map dim_tuple, tp_entity = unflattening_map[(sd - 1, f)] # Determine the axis that's orthogonal to the facet - # quad: dim_tuple = (0,1) -> facet is on the boundary of axis 0 (x fixed) -> x = 0 or x = 1 - # dim_tuple = (1,0) -> facet is on the boundary of axis 1 (y fixed) -> y = 0 or y = 1 + # E.g., in a quad: + # if dim_tuple = (0,1) -> facet has dimension 0 on the first component -> fixed at x = 0 or x = 1 + # if dim_tuple = (1,0) -> facet has dimension 0 on the second component -> fixed at y = 0 or y = 1 axis = next( i for i, d in enumerate(dim_tuple) if d == product.cells[i].get_dimension() - 1 ) - # Determine the index of the endpoint that produces the facet (local facet number in the axis space) + # Determine the index of the endpoint that produces the facet + # which gives the local facet number in the axis space entity_shape = tuple( len(c.get_topology()[d]) for c, d in zip(product.cells, dim_tuple) @@ -2000,22 +1934,11 @@ def compute_facet_perm(unflattening_map, product): local_facet = tuple_ei[axis] # Map the local_facet index to the corresponding barycentric coordinate index. - # - # For a 1D interval with endpoints 0 and 1, the barycentric coordinates are: - # lambda_0 = 1 - x - # lambda_1 = x - # - # Each endpoint facet is characterised by one barycentric coordinate vanishing: - # facet at x = 0 -> lambda_1 = 0 (opposite vertex 1) - # facet at x = 1 -> lambda_0 = 0 (opposite vertex 0) - # - # In general, FIAT's convention is: barycentric coordinate i vanishes on the - # facet opposite vertex i. - bary_index = get_bary_index_from_local_facet(product.cells[axis], local_facet) - - # NOTE: this assumes that the barycentric coordinate vector is flattened - # perm[f] = axis_offsets[axis] + bary_index - perm[f] = (axis, bary_index) + # For a UFCInterval, the local facet index equals the index of the opposite vertex, + # (e.g., local facet 1 is the facet that does not contain vertex 1) which also corresponds + # to the index of the barycentric coordinate vanishing on that facet. + + perm[f] = axis_offsets[axis] + local_facet return perm @@ -2025,13 +1948,4 @@ def max_complex(complexes): if all(max_cell >= b for b in complexes): return max_cell else: - raise ValueError("Cannot find the maximal complex") - -def get_bary_index_from_local_facet(simplicial_cell, local_facet): - """Return the barycentric coordinate index that vanishes on a given facet of a simplicial cell. """ - - all_vertices = set(simplicial_cell.get_topology()[0].keys()) - facet_dim = simplicial_cell.get_dimension() - 1 - facet_vertices = set(simplicial_cell.get_topology()[facet_dim][local_facet]) - missing = all_vertices - facet_vertices - return next(iter(missing)) \ No newline at end of file + raise ValueError("Cannot find the maximal complex") \ No newline at end of file diff --git a/gem/gem.py b/gem/gem.py index 070498fb..1d2a49b2 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -1,3 +1,5 @@ +from __future__ import annotations + """GEM is the intermediate language of TSFC for describing tensor-valued mathematical expressions and tensor operations. It is similar to Einstein's notation. @@ -20,6 +22,8 @@ from operator import attrgetter from numbers import Integral, Number +from types import EllipsisType + import numpy from numpy import asarray @@ -80,20 +84,18 @@ def is_equal(self, other): if result: self.children = other.children return result - - def __getitem__(self, key): - """ - Generalised indexing interface. - Parameters - ---------- - key: int, Index, VariableIndex, slice, Ellipsis or tuple - """ - # Normalize to tuple so we can inspect the key + def __getitem__( + self, + key: int | Index | VariableIndex | slice | EllipsisType | + tuple[int | Index | VariableIndex | slice | EllipsisType, ...] + ) -> ComponentTensor | Indexed: + """A generalised interface for indexing a Gem.Node""" + if not isinstance(key, tuple): key = (key,) - # Expand ellipsis -> fill remaining dimensions with slice(None) + # Expand ellipsis -> fill in remaining dimensions with slice(None) if Ellipsis in key: if key.count(Ellipsis) > 1: raise NotImplementedError("Multiple ellipses are not supported when indexing a gem.Node tensor") @@ -106,6 +108,7 @@ def __getitem__(self, key): + (slice(None), ) * num_missing + key[ellipsis_pos + 1:] ) + # Slice indexing -> delegate to view() if any(isinstance(k, slice) for k in key): # view expects one slice for each axis/dim of the tensor @@ -116,14 +119,6 @@ def __getitem__(self, key): else: return Indexed(self, key) - # OLD - # def __getitem__(self, indices): - # try: - # indices = tuple(indices) - # except TypeError: - # indices = (indices, ) - # return Indexed(self, indices) - def __neg__(self): return componentwise(Product, minus, self) @@ -153,7 +148,7 @@ def __matmul__(self, other): raise ValueError("Both objects must have shape for matmul") elif self.shape[-1] != other.shape[0]: raise ValueError(f"Mismatching shapes {self.shape} and {other.shape} in matmul") - + *i, k = indices(len(self.shape)) _, *j = indices(len(other.shape)) expr = Product(Indexed(self, (*i, k)), Indexed(other, (k, *j))) @@ -1382,7 +1377,7 @@ def as_gem(expr): """ if isinstance(expr, Node): return expr - elif isinstance(expr, Number) or isinstance(expr, numpy.ndarray): + elif isinstance(expr, Number): return Literal(expr) elif isinstance(expr, (bool, numpy.bool)): return Literal(bool(expr)) From 3d87fc1118a0dddc431b53d5306bd6b44e313da6 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Mon, 20 Apr 2026 16:27:24 +0100 Subject: [PATCH 15/21] modified unit tests for computing bary coords on points passed as numpy arrays --- FIAT/reference_element.py | 48 ++++++++++++------------ test/FIAT/unit/test_reference_element.py | 33 +++++++++++----- 2 files changed, 49 insertions(+), 32 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 5028be94..8ab8c93a 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -615,12 +615,10 @@ def get_dimension(self): def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the complex.""" - # if len(points) == 0: - # return points - if isinstance(points, (list, tuple, numpy.ndarray)) and len(points) == 0: + if isinstance(points, numpy.ndarray) and len(points) == 0: return points - + if entity is None: entity = (self.get_spatial_dimension(), 0) entity_dim, entity_id = entity @@ -1412,7 +1410,7 @@ def extrinsic_orientation_permutation_map(self): def is_macrocell(self): return any(c.is_macrocell() for c in self.cells) - + def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=False): """Compute barycentric coordinates on each axis (factor) of a tensor-product cell. @@ -1424,7 +1422,7 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals Returns ------- numpy.ndarray - An array of shape ``(nfactors,)`` and dtype object if points are GEM nodes, + An array of shape ``(nfactors, total_bary_coords)`` and dtype object if points are GEM nodes, otherwise dtype numeric. The i-th entry contains the barycentric coordinates on the i-th factor cell. If factor i is a simplex of dimension d, this will have shape ``(npoints, d+1)``. If factor i is a hypercube of dimension d, @@ -1432,14 +1430,15 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals """ if isinstance(points, numpy.ndarray) and len(points) == 0: return points - + axis_dims = [c.get_spatial_dimension() for c in self.cells] point_slices = TensorProductCell._split_slices(axis_dims) - return numpy.array([ - factor.compute_barycentric_coordinates(points[..., s], entity, rescale) - for factor, s in zip(self.cells, point_slices) - ]) + result = numpy.empty(len(self.cells), dtype=object) + for k, (factor, s) in enumerate(zip(self.cells, point_slices)): + result[k] = factor.compute_barycentric_coordinates(points[..., s], entity, rescale) + return numpy.concatenate(result) + class Hypercube(Cell): """Abstract class for a reference hypercube""" @@ -1559,7 +1558,7 @@ def __le__(self, other): def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on the hypercube. - + Parameters ---------- points: numpy.ndarray or GEM.Node @@ -1575,9 +1574,10 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): return points tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) - + return tp_bary_coords[..., self.facet_perm] + class UFCHypercube(Hypercube): """Reference UFC Hypercube @@ -1894,22 +1894,23 @@ def compute_unflattening_map(topology_dict): return unflattening_map + def compute_facet_permutation(unflattening_map, product): """ Return a permutation mapping each hypercube facet to the index of its vanishing barycentric coordinate in the axis-structured barycentric array. """ - # First compute axis offsets into the flattened barycentric coordinate array. + # First compute axis offsets into the flattened barycentric coordinate array. axis_offsets = [] offset = 0 for axis_cell in product.cells: axis_offsets.append(offset) - offset += axis_cell.get_dimension() + 1 - + offset += axis_cell.get_dimension() + 1 + # Initialise the integer permutation array sd = len(product.cells) num_facets = 2 * sd - perm = numpy.zeroes(num_facets, dtype=int) + perm = numpy.zeros(num_facets, dtype=int) for f in range(num_facets): # Recover the tensor-product representation of the facet as given by the unflattening map @@ -1933,12 +1934,13 @@ def compute_facet_permutation(unflattening_map, product): tuple_ei = numpy.unravel_index(tp_entity, entity_shape) local_facet = tuple_ei[axis] - # Map the local_facet index to the corresponding barycentric coordinate index. - # For a UFCInterval, the local facet index equals the index of the opposite vertex, - # (e.g., local facet 1 is the facet that does not contain vertex 1) which also corresponds - # to the index of the barycentric coordinate vanishing on that facet. + # For a simplex (UFCInterval, UFCTriangle), the barycentric coordinate that vanishes on local facet i + # corresponds to the ID of the vertex that doesn't belong to that facet + all_vertices = set(product.cells[axis].get_topology()[0].keys()) + facet_vertices = set(product.cells[axis].get_topology()[0][local_facet]) + bary_index = next(iter(all_vertices - facet_vertices)) - perm[f] = axis_offsets[axis] + local_facet + perm[f] = axis_offsets[axis] + bary_index return perm @@ -1948,4 +1950,4 @@ def max_complex(complexes): if all(max_cell >= b for b in complexes): return max_cell else: - raise ValueError("Cannot find the maximal complex") \ No newline at end of file + raise ValueError("Cannot find the maximal complex") diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index 736f84a6..7c8304d5 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -459,16 +459,23 @@ def test_flatten_maintains_ufc_status(cell): [(interval_x_interval, [0.25, 0.6]), (triangle_x_interval, [0.25, 0.25, 0.5]), (quadrilateral_x_interval, [0.25, 0.25, 0.5])]) -def test_tensorproduct_barycentrics_match_factor_cells(cell, point): - bary_axes = cell.compute_axis_barycentric_coordinates(point) - +def test_tp_axis_bary_coords(cell, point, epsilon=1e-12): + point = np.asarray(point) + axis_bary_coords = cell.compute_axis_barycentric_coordinates(point) + + assert type(axis_bary_coords) is np.ndarray + offset = 0 - for k, factor in enumerate(cell.simplex_cells): + bary_offset = 0 + for factor in cell.cells: sd = factor.get_spatial_dimension() - coords_k = point[offset: offset+sd] # coords on factor cell k + coords_k = point[offset: offset + sd] offset += sd expected = factor.compute_barycentric_coordinates(coords_k) - assert np.allclose(bary_axes[k], expected) + n_bary = len(expected) + assert np.allclose(axis_bary_coords[bary_offset: bary_offset + n_bary], expected, atol=epsilon) + bary_offset += n_bary + @pytest.mark.parametrize(('cell', 'point'), [(quadrilateral, [0.0, 0.3]), @@ -481,14 +488,22 @@ def test_tensorproduct_barycentrics_match_factor_cells(cell, point): (hexahedron, [0.3, 1.0, 0.4]), (hexahedron, [0.3, 0.4, 0.0]), (hexahedron, [0.3, 0.4, 1.0]),]) -def test_hypercube_barycentrics_match_facet_order(cell, point, epsilon=1e-14): +def test_hypercube_bary_coords_are_in_facet_order(cell, point, epsilon=1e-12): + point = np.asarray(point) + facet_dim = cell.get_spatial_dimension() - 1 point_entity_ids = cell.point_entity_ids([point]) facet_hits = [fid for fid, pts in point_entity_ids[facet_dim].items() if len(pts) > 0] assert len(facet_hits) == 1 + facet_id = facet_hits[0] - bary = cell.compute_barycentric_coordinates(point) - assert np.isclose(bary[0, facet_id], 0.0, atol=epsilon) + bary_coords = cell.compute_barycentric_coordinates(point) + assert np.isclose(bary_coords[facet_id], 0.0, atol=epsilon) + + mask = np.ones(len(bary_coords), dtype=bool) + mask[facet_id] = False + assert np.all(bary_coords[mask] > epsilon) + if __name__ == '__main__': import os From 76e63f73befbf2f0daeddd99fb7f16161d00e144 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Tue, 21 Apr 2026 10:46:22 +0100 Subject: [PATCH 16/21] latest changes + gem tests --- FIAT/reference_element.py | 19 ++++++-- gem/gem.py | 59 ++++++++++++++++++------ test/FIAT/unit/test_reference_element.py | 24 ++++++++++ 3 files changed, 84 insertions(+), 18 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 8ab8c93a..c20ece9c 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1422,12 +1422,14 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals Returns ------- numpy.ndarray - An array of shape ``(nfactors, total_bary_coords)`` and dtype object if points are GEM nodes, + A flattened array of shape ``(total_bary_coords, )`` and dtype object if points are GEM nodes, otherwise dtype numeric. The i-th entry contains the barycentric coordinates on the i-th factor cell. If factor i is a simplex of dimension d, this will have shape ``(npoints, d+1)``. If factor i is a hypercube of dimension d, this will have shape ``(npoints, 2*d)``. """ + import gem + if isinstance(points, numpy.ndarray) and len(points) == 0: return points @@ -1437,7 +1439,18 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals result = numpy.empty(len(self.cells), dtype=object) for k, (factor, s) in enumerate(zip(self.cells, point_slices)): result[k] = factor.compute_barycentric_coordinates(points[..., s], entity, rescale) - return numpy.concatenate(result) + + # Flatten the array + # NOTE: cannot construct the flat array directly since we may not know upfront the total number + # of barycentric coordinates (e.g., in a simplex: d+1, in a hypercube: 2*d) + flat_result = numpy.array([bary[j] for bary in result for j in range(bary.shape[0])]) + + if isinstance(points, numpy.ndarray): + return flat_result + elif isinstance(points, gem.Node): + return gem.as_gem(flat_result) + else: + raise ValueError(f"Expected a numpy.ndarray or gem.Node, got {type(points)}") class Hypercube(Cell): @@ -1575,7 +1588,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): tp_bary_coords = self.product.compute_axis_barycentric_coordinates(points, entity, rescale) - return tp_bary_coords[..., self.facet_perm] + return tp_bary_coords[self.facet_perm] class UFCHypercube(Hypercube): diff --git a/gem/gem.py b/gem/gem.py index 1d2a49b2..b9341c0a 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -23,6 +23,7 @@ from numbers import Integral, Number from types import EllipsisType +import itertools import numpy from numpy import asarray @@ -87,37 +88,65 @@ def is_equal(self, other): def __getitem__( self, - key: int | Index | VariableIndex | slice | EllipsisType | - tuple[int | Index | VariableIndex | slice | EllipsisType, ...] - ) -> ComponentTensor | Indexed: + key: int | Index | VariableIndex | slice | EllipsisType | list | numpy.ndarray | + tuple[int | Index | VariableIndex | slice | EllipsisType, ...], + ) -> ComponentTensor | Indexed | ListTensor: """A generalised interface for indexing a Gem.Node""" - + # Normalize to tuple for inspection if not isinstance(key, tuple): key = (key,) + # Support a list or array of integer indices + if any(isinstance(k, (list, numpy.ndarray)) for k in key): + pos = next(i for i, k in enumerate(key) if isinstance(k, (list, numpy.ndarray))) + components = numpy.array( + [Indexed(self, key[:pos] + (int(i),) + key[pos + 1:]) for i in key[pos]], + dtype=object + ) + return as_gem(components) + # Expand ellipsis -> fill in remaining dimensions with slice(None) if Ellipsis in key: if key.count(Ellipsis) > 1: - raise NotImplementedError("Multiple ellipses are not supported when indexing a gem.Node tensor") + raise NotImplementedError("Multiple ellipses are not supported.") ellipsis_pos = key.index(Ellipsis) - num_missing = len(self.shape) - (len(key) - 1) - if num_missing < 0: - raise IndexError("Too many indices provided when indexing the gem.Node tensor") + remaining_dims = len(self.shape) - (len(key) - 1) + if remaining_dims < 0: + raise IndexError("Too many indices provided.") key = ( key[:ellipsis_pos] - + (slice(None), ) * num_missing + + (slice(None), ) * remaining_dims + key[ellipsis_pos + 1:] ) # Slice indexing -> delegate to view() + # NOTE: view() produces a ComponentTensor with a FlexiblyIndexed Node that gem's evaluator cannot evaluate + # if any(isinstance(k, slice) for k in key): + # # view expects one slice for each axis/dim of the tensor + # if len(key) != len(self.shape): + # raise IndexError("Expects the number of slices to match the gem.Node tensor rank") + # return view(self, *key) + + # Slice indexing -> build ListTensor of Indexed nodes if any(isinstance(k, slice) for k in key): - # view expects one slice for each axis/dim of the tensor if len(key) != len(self.shape): - raise IndexError("Slice indexing expects the number of slices to match the gem.Node tensor rank") - return view(self, *key) - # Point indexing -> delegate to Index - else: - return Indexed(self, key) + raise IndexError("Expects one key per dimension.") + # Unpack slices into ranges + ranges = [ + range( + k.start if k.start is not None else 0, + k.stop if k.stop is not None else s, + k.step if k.step is not None else 1 + ) for k, s in zip(key, self.shape) + ] + # Extract tensor indices from the cartesian product of ranges + components = numpy.array( + [Indexed(self, idx) for idx in itertools.product(*ranges)], + dtype=object).reshape([len(r) for r in ranges]) + return as_gem(components) + + # Point indexing + return Indexed(self, key) def __neg__(self): return componentwise(Product, minus, self) diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index 7c8304d5..e61b6a55 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -505,6 +505,30 @@ def test_hypercube_bary_coords_are_in_facet_order(cell, point, epsilon=1e-12): assert np.all(bary_coords[mask] > epsilon) +@pytest.mark.parametrize(('cell', 'point'), + [(interval, [0.5]), + (triangle, [0.25, 0.25]), + (tetrahedron, [0.25, 0.25, 0.25]), + (quadrilateral, [0.25, 0.5]), + (hexahedron, [0.25, 0.5, 0.25]),]) +def test_bary_coords_gem(cell, point): + import gem + from gem.interpreter import evaluate + + point = np.asarray(point) + sd = cell.get_spatial_dimension() + + coords = gem.Variable('X', (sd,)) + bindings = {coords: point} + + bary_gem = cell.compute_barycentric_coordinates(coords) + results, = evaluate((gem.as_gem(bary_gem),), bindings=bindings) + results = results.arr + + bary_numpy = cell.compute_barycentric_coordinates(point) + assert np.allclose(results, bary_numpy) + + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__)) From e2403c4f2b615a780fd34e2a36d75f92c93a8da0 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Tue, 21 Apr 2026 10:53:34 +0100 Subject: [PATCH 17/21] removed egg-info from tracking --- .gitignore | 2 + firedrake_fiat.egg-info/PKG-INFO | 75 ------- firedrake_fiat.egg-info/SOURCES.txt | 212 ------------------- firedrake_fiat.egg-info/dependency_links.txt | 1 - firedrake_fiat.egg-info/requires.txt | 13 -- firedrake_fiat.egg-info/top_level.txt | 3 - 6 files changed, 2 insertions(+), 304 deletions(-) delete mode 100644 firedrake_fiat.egg-info/PKG-INFO delete mode 100644 firedrake_fiat.egg-info/SOURCES.txt delete mode 100644 firedrake_fiat.egg-info/dependency_links.txt delete mode 100644 firedrake_fiat.egg-info/requires.txt delete mode 100644 firedrake_fiat.egg-info/top_level.txt diff --git a/.gitignore b/.gitignore index b4109385..7557ed0b 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,5 @@ release/ /docs/source/finat.rst /docs/source/finat.ufl.rst /docs/source/gem.rst + +firedrake_fiat.egg-info/ diff --git a/firedrake_fiat.egg-info/PKG-INFO b/firedrake_fiat.egg-info/PKG-INFO deleted file mode 100644 index 1eaca9f4..00000000 --- a/firedrake_fiat.egg-info/PKG-INFO +++ /dev/null @@ -1,75 +0,0 @@ -Metadata-Version: 2.4 -Name: firedrake-fiat -Version: 2025.11.0.dev0 -Summary: FInite element Automatic Tabulator -Author-email: "Robert C. Kirby et al." , Imperial College London and others -Project-URL: Repository, https://github.com/firedrakeproject/fiat.git -Classifier: Programming Language :: Python -Requires-Python: >=3.10 -Description-Content-Type: text/x-rst -License-File: COPYING -License-File: COPYING.LESSER -License-File: AUTHORS -Requires-Dist: fenics-ufl @ git+https://github.com/FEniCS/ufl.git@main -Requires-Dist: numpy>=1.16 -Requires-Dist: recursivenodes -Requires-Dist: scipy -Requires-Dist: symengine -Requires-Dist: sympy -Provides-Extra: doc -Requires-Dist: setuptools; extra == "doc" -Requires-Dist: sphinx; extra == "doc" -Provides-Extra: test -Requires-Dist: pytest; extra == "test" -Dynamic: license-file - -======================================== -FIAT: FInite element Automatic Tabulator -======================================== - -The FInite element Automatic Tabulator FIAT supports generation of -arbitrary order instances of the Lagrange elements on lines, -triangles, and tetrahedra. It is also capable of generating arbitrary -order instances of Jacobi-type quadrature rules on the same element -shapes. Further, H(div) and H(curl) conforming finite element spaces -such as the families of Raviart-Thomas, Brezzi-Douglas-Marini and -Nedelec are supported on triangles and tetrahedra. Upcoming versions -will also support Hermite and nonconforming elements. - -FIAT is part of the FEniCS Project. - -For more information, visit http://www.fenicsproject.org - - -.. image:: https://github.com/FEniCS/fiat/workflows/FIAT%20CI/badge.svg - :target: https://github.com/FEniCS/fiat/actions?query=workflow%3A%22FIAT+CI%22 - :alt: Build Status - -.. image:: https://coveralls.io/repos/github/FEniCS/fiat/badge.svg?branch=main - :target: https://coveralls.io/github/FEniCS/fiat?branch=main - :alt: Coverage Status - -.. image:: https://readthedocs.org/projects/fenics-fiat/badge/?version=latest - :target: http://fenics.readthedocs.io/projects/fiat/en/latest/?badge=latest - :alt: Documentation Status - -Documentation -============= - -Documentation can be viewed at http://fenics-fiat.readthedocs.org/. - -License -======= - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with this program. If not, see . diff --git a/firedrake_fiat.egg-info/SOURCES.txt b/firedrake_fiat.egg-info/SOURCES.txt deleted file mode 100644 index e1867ace..00000000 --- a/firedrake_fiat.egg-info/SOURCES.txt +++ /dev/null @@ -1,212 +0,0 @@ -AUTHORS -COPYING -COPYING.LESSER -MANIFEST.in -README.rst -pyproject.toml -setup.cfg -FIAT/P0.py -FIAT/Sminus.py -FIAT/SminusCurl.py -FIAT/SminusDiv.py -FIAT/__init__.py -FIAT/alfeld_sorokina.py -FIAT/argyris.py -FIAT/arnold_qin.py -FIAT/arnold_winther.py -FIAT/barycentric_interpolation.py -FIAT/bell.py -FIAT/bernardi_raugel.py -FIAT/bernstein.py -FIAT/brezzi_douglas_fortin_marini.py -FIAT/brezzi_douglas_marini.py -FIAT/brezzi_douglas_marini_cube.py -FIAT/bubble.py -FIAT/c2_elements.py -FIAT/check_format_variant.py -FIAT/christiansen_hu.py -FIAT/crouzeix_raviart.py -FIAT/discontinuous.py -FIAT/discontinuous_lagrange.py -FIAT/discontinuous_pc.py -FIAT/discontinuous_raviart_thomas.py -FIAT/discontinuous_taylor.py -FIAT/dual_set.py -FIAT/enriched.py -FIAT/expansions.py -FIAT/fdm_element.py -FIAT/finite_element.py -FIAT/functional.py -FIAT/gauss_legendre.py -FIAT/gauss_lobatto_legendre.py -FIAT/gauss_radau.py -FIAT/gopalakrishnan_lederer_schoberl.py -FIAT/guzman_neilan.py -FIAT/hct.py -FIAT/hdiv_trace.py -FIAT/hdivcurl.py -FIAT/hellan_herrmann_johnson.py -FIAT/hermite.py -FIAT/hierarchical.py -FIAT/histopolation.py -FIAT/hu_zhang.py -FIAT/jacobi.py -FIAT/johnson_mercier.py -FIAT/kong_mulder_veldhuizen.py -FIAT/lagrange.py -FIAT/macro.py -FIAT/mardal_tai_winther.py -FIAT/mixed.py -FIAT/morley.py -FIAT/nedelec.py -FIAT/nedelec_second_kind.py -FIAT/nodal_enriched.py -FIAT/orientation_utils.py -FIAT/orthopoly.py -FIAT/pointwise_dual.py -FIAT/polynomial_set.py -FIAT/powell_sabin.py -FIAT/quadrature.py -FIAT/quadrature_element.py -FIAT/quadrature_schemes.py -FIAT/raviart_thomas.py -FIAT/reference_element.py -FIAT/regge.py -FIAT/restricted.py -FIAT/serendipity.py -FIAT/tensor_product.py -FIAT/walkington.py -FIAT/wuxu.py -FIAT/xg_quad_data.py -finat/__init__.py -finat/alfeld_sorokina.py -finat/argyris.py -finat/arnold_qin.py -finat/aw.py -finat/bell.py -finat/bernardi_raugel.py -finat/c2_elements.py -finat/cell_tools.py -finat/christiansen_hu.py -finat/citations.py -finat/cube.py -finat/direct_serendipity.py -finat/discontinuous.py -finat/element_factory.py -finat/enriched.py -finat/fiat_elements.py -finat/finiteelementbase.py -finat/guzman_neilan.py -finat/hct.py -finat/hdivcurl.py -finat/hermite.py -finat/hz.py -finat/johnson_mercier.py -finat/mixed.py -finat/morley.py -finat/mtw.py -finat/nodal_enriched.py -finat/physically_mapped.py -finat/piola_mapped.py -finat/point_set.py -finat/powell_sabin.py -finat/quadrature.py -finat/quadrature_element.py -finat/restricted.py -finat/runtime_tabulated.py -finat/spectral.py -finat/sympy2gem.py -finat/tensor_product.py -finat/tensorfiniteelement.py -finat/walkington.py -finat/wuxu.py -finat/ufl/__init__.py -finat/ufl/brokenelement.py -finat/ufl/elementlist.py -finat/ufl/enrichedelement.py -finat/ufl/finiteelement.py -finat/ufl/finiteelementbase.py -finat/ufl/hdivcurl.py -finat/ufl/mixedelement.py -finat/ufl/restrictedelement.py -finat/ufl/tensorproductelement.py -firedrake_fiat.egg-info/PKG-INFO -firedrake_fiat.egg-info/SOURCES.txt -firedrake_fiat.egg-info/dependency_links.txt -firedrake_fiat.egg-info/requires.txt -firedrake_fiat.egg-info/top_level.txt -gem/__init__.py -gem/coffee.py -gem/flop_count.py -gem/gem.py -gem/impero.py -gem/impero_utils.py -gem/interpreter.py -gem/node.py -gem/optimise.py -gem/pprint.py -gem/refactorise.py -gem/scheduling.py -gem/unconcatenate.py -gem/utils.py -test/conftest.py -test/FIAT/regression/.gitignore -test/FIAT/regression/README.rst -test/FIAT/regression/conftest.py -test/FIAT/regression/fiat-reference-data-id -test/FIAT/regression/test_regression.py -test/FIAT/regression/scripts/download -test/FIAT/regression/scripts/getdata -test/FIAT/regression/scripts/getreferencerepo -test/FIAT/regression/scripts/parameters -test/FIAT/regression/scripts/upload -test/FIAT/unit/test_argyris.py -test/FIAT/unit/test_awc.py -test/FIAT/unit/test_awnc.py -test/FIAT/unit/test_bdfm.py -test/FIAT/unit/test_bernstein.py -test/FIAT/unit/test_discontinuous_pc.py -test/FIAT/unit/test_discontinuous_taylor.py -test/FIAT/unit/test_facet_support_dofs.py -test/FIAT/unit/test_fdm.py -test/FIAT/unit/test_fiat.py -test/FIAT/unit/test_gauss_legendre.py -test/FIAT/unit/test_gauss_lobatto_legendre.py -test/FIAT/unit/test_gauss_radau.py -test/FIAT/unit/test_gopalakrishnan_lederer_schoberl.py -test/FIAT/unit/test_hct.py -test/FIAT/unit/test_hdivtrace.py -test/FIAT/unit/test_hierarchical.py -test/FIAT/unit/test_histopolation.py -test/FIAT/unit/test_johnson_mercier.py -test/FIAT/unit/test_kong_mulder_veldhuizen.py -test/FIAT/unit/test_macro.py -test/FIAT/unit/test_mtw.py -test/FIAT/unit/test_orientation.py -test/FIAT/unit/test_pointwise_dual.py -test/FIAT/unit/test_polynomial.py -test/FIAT/unit/test_powell_sabin.py -test/FIAT/unit/test_quadrature.py -test/FIAT/unit/test_quadrature_element.py -test/FIAT/unit/test_reference_element.py -test/FIAT/unit/test_regge_hhj.py -test/FIAT/unit/test_serendipity.py -test/FIAT/unit/test_stokes_complex.py -test/FIAT/unit/test_tensor_product.py -test/FIAT/unit/test_walkington.py -test/finat/conftest.py -test/finat/test_create_broken_element.py -test/finat/test_create_fiat_element.py -test/finat/test_create_finat_element.py -test/finat/test_create_restricted_element.py -test/finat/test_direct_serendipity.py -test/finat/test_dual_basis.py -test/finat/test_hash.py -test/finat/test_mass_conditioning.py -test/finat/test_point_evaluation.py -test/finat/test_quadrature.py -test/finat/test_quadrature_element.py -test/finat/test_restriction.py -test/finat/test_ufl_elements.py -test/finat/test_zany_mapping.py -test/gem/test_simplify.py \ No newline at end of file diff --git a/firedrake_fiat.egg-info/dependency_links.txt b/firedrake_fiat.egg-info/dependency_links.txt deleted file mode 100644 index 8b137891..00000000 --- a/firedrake_fiat.egg-info/dependency_links.txt +++ /dev/null @@ -1 +0,0 @@ - diff --git a/firedrake_fiat.egg-info/requires.txt b/firedrake_fiat.egg-info/requires.txt deleted file mode 100644 index 08d00a07..00000000 --- a/firedrake_fiat.egg-info/requires.txt +++ /dev/null @@ -1,13 +0,0 @@ -fenics-ufl @ git+https://github.com/FEniCS/ufl.git@main -numpy>=1.16 -recursivenodes -scipy -symengine -sympy - -[doc] -setuptools -sphinx - -[test] -pytest diff --git a/firedrake_fiat.egg-info/top_level.txt b/firedrake_fiat.egg-info/top_level.txt deleted file mode 100644 index 6bc74dc3..00000000 --- a/firedrake_fiat.egg-info/top_level.txt +++ /dev/null @@ -1,3 +0,0 @@ -FIAT -finat -gem From 26e9f39cba6807e9318d4d9c51b530f39b832818 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Thu, 23 Apr 2026 14:42:06 +0100 Subject: [PATCH 18/21] implemented handler in gem.evaluate for FlexiblyIndexed nodes --- gem/gem.py | 28 ++++----------------------- gem/interpreter.py | 48 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 24 deletions(-) diff --git a/gem/gem.py b/gem/gem.py index b9341c0a..8810fb26 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -91,8 +91,7 @@ def __getitem__( key: int | Index | VariableIndex | slice | EllipsisType | list | numpy.ndarray | tuple[int | Index | VariableIndex | slice | EllipsisType, ...], ) -> ComponentTensor | Indexed | ListTensor: - """A generalised interface for indexing a Gem.Node""" - # Normalize to tuple for inspection + """A generalised interface for indexing GEM tensors""" if not isinstance(key, tuple): key = (key,) @@ -120,30 +119,11 @@ def __getitem__( ) # Slice indexing -> delegate to view() - # NOTE: view() produces a ComponentTensor with a FlexiblyIndexed Node that gem's evaluator cannot evaluate - # if any(isinstance(k, slice) for k in key): - # # view expects one slice for each axis/dim of the tensor - # if len(key) != len(self.shape): - # raise IndexError("Expects the number of slices to match the gem.Node tensor rank") - # return view(self, *key) - - # Slice indexing -> build ListTensor of Indexed nodes if any(isinstance(k, slice) for k in key): + # view expects one slice for each axis/dim of the tensor if len(key) != len(self.shape): - raise IndexError("Expects one key per dimension.") - # Unpack slices into ranges - ranges = [ - range( - k.start if k.start is not None else 0, - k.stop if k.stop is not None else s, - k.step if k.step is not None else 1 - ) for k, s in zip(key, self.shape) - ] - # Extract tensor indices from the cartesian product of ranges - components = numpy.array( - [Indexed(self, idx) for idx in itertools.product(*ranges)], - dtype=object).reshape([len(r) for r in ranges]) - return as_gem(components) + raise IndexError("Expects the number of slices to match the gem.Node tensor rank") + return view(self, *key) # Point indexing return Indexed(self, key) diff --git a/gem/interpreter.py b/gem/interpreter.py index 13eeb44a..70094c10 100644 --- a/gem/interpreter.py +++ b/gem/interpreter.py @@ -286,6 +286,54 @@ def _evaluate_indexed(e, self): return Result(val[idx], val.fids + fids) +@_evaluate.register(gem.FlexiblyIndexed) +def _evaluate_flexiblyindexed(e, self): + val = self(e.children[0]) + + all_fids = val.fids + e.free_indices + + idx = [] + + # Collect existing free indices + for fid in val.fids: + shape = tuple(fid.extent if f is fid else 1 for f in all_fids) + fidx_array = numpy.arange(fid.extent).reshape(shape) + idx.append(fidx_array) + + # Compute a flat index for each dim from dim2idxs + for dim, (offset, idxs) in zip(e.children[0].shape, e.dim2idxs): + if isinstance(offset, gem.Node): + offset_result = self(offset) + assert not offset_result.tshape + offset_val = offset_result[()] + else: + offset_val = offset + + dim_flat_idx_components = [] + for i, s in idxs: + # Index i may be one of: Index, VariableIndex, int + if isinstance(i, gem.Index): + # Free index contributes an array index given by the extent + shape = tuple(i.extent if f is i else 1 for f in all_fids) + i_vals = numpy.arange(i.extent).reshape(shape) * s + dim_flat_idx_components.append(i_vals) + elif isinstance(i, gem.VariableIndex): + # Variable index gives a single integer index + # obtained by evaluating its inner expression + result, = self(i.expression) + assert not result.tshape + dim_flat_idx_components.append(result[()]*s) + else: + # Fixed index is just an integer so use that + dim_flat_idx_components.append(i*s) + + # From dim_flat_idx_components compute the flat index into dim + dim_idx_array = offset_val + sum(flat_index_component for flat_index_component in dim_flat_idx_components) + idx.append(dim_idx_array) + + return Result(val[idx], all_fids) + + @_evaluate.register(gem.ComponentTensor) def _evaluate_componenttensor(e, self): """Component tensors map free indices to shape.""" From d9a5e4b2748e2bda709b6429303f5c115ad54035 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Thu, 23 Apr 2026 16:54:12 +0100 Subject: [PATCH 19/21] added a new ListIndex type to GEM for indexing GEM tensors using an index array --- FIAT/reference_element.py | 13 ++----- gem/gem.py | 78 ++++++++++++++++++++++++++++++--------- gem/interpreter.py | 24 +++++++++--- 3 files changed, 83 insertions(+), 32 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index c20ece9c..51f9581b 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1428,7 +1428,6 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals have shape ``(npoints, d+1)``. If factor i is a hypercube of dimension d, this will have shape ``(npoints, 2*d)``. """ - import gem if isinstance(points, numpy.ndarray) and len(points) == 0: return points @@ -1439,19 +1438,13 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals result = numpy.empty(len(self.cells), dtype=object) for k, (factor, s) in enumerate(zip(self.cells, point_slices)): result[k] = factor.compute_barycentric_coordinates(points[..., s], entity, rescale) - + # Flatten the array # NOTE: cannot construct the flat array directly since we may not know upfront the total number - # of barycentric coordinates (e.g., in a simplex: d+1, in a hypercube: 2*d) + # of barycentric coordinates (e.g., in a simplex it is d+1, in a hypercube it is 2*d) flat_result = numpy.array([bary[j] for bary in result for j in range(bary.shape[0])]) - if isinstance(points, numpy.ndarray): - return flat_result - elif isinstance(points, gem.Node): - return gem.as_gem(flat_result) - else: - raise ValueError(f"Expected a numpy.ndarray or gem.Node, got {type(points)}") - + return flat_result class Hypercube(Cell): """Abstract class for a reference hypercube""" diff --git a/gem/gem.py b/gem/gem.py index 8810fb26..dc921438 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -37,7 +37,7 @@ 'Variable', 'Sum', 'Product', 'Division', 'FloorDiv', 'Remainder', 'Power', 'MathFunction', 'MinValue', 'MaxValue', 'Comparison', 'LogicalNot', 'LogicalAnd', 'LogicalOr', 'Conditional', - 'Index', 'VariableIndex', 'Indexed', 'ComponentTensor', + 'Index', 'VariableIndex', 'ListIndex', 'Indexed', 'ComponentTensor', 'IndexSum', 'ListTensor', 'Concatenate', 'Delta', 'OrientationVariableIndex', 'index_sum', 'partial_indexed', 'reshape', 'view', 'indices', 'as_gem', 'FlexiblyIndexed', @@ -46,7 +46,6 @@ uint_type = numpy.dtype(numpy.uintc) - class NodeMeta(type): """Metaclass of GEM nodes. @@ -88,24 +87,14 @@ def is_equal(self, other): def __getitem__( self, - key: int | Index | VariableIndex | slice | EllipsisType | list | numpy.ndarray | - tuple[int | Index | VariableIndex | slice | EllipsisType, ...], - ) -> ComponentTensor | Indexed | ListTensor: + key: IndexT | tuple[IndexT, ...], + ) -> ComponentTensor | Indexed: """A generalised interface for indexing GEM tensors""" if not isinstance(key, tuple): key = (key,) - - # Support a list or array of integer indices - if any(isinstance(k, (list, numpy.ndarray)) for k in key): - pos = next(i for i, k in enumerate(key) if isinstance(k, (list, numpy.ndarray))) - components = numpy.array( - [Indexed(self, key[:pos] + (int(i),) + key[pos + 1:]) for i in key[pos]], - dtype=object - ) - return as_gem(components) - + # Expand ellipsis -> fill in remaining dimensions with slice(None) - if Ellipsis in key: + if any(k is Ellipsis for k in key): if key.count(Ellipsis) > 1: raise NotImplementedError("Multiple ellipses are not supported.") ellipsis_pos = key.index(Ellipsis) @@ -118,13 +107,29 @@ def __getitem__( + key[ellipsis_pos + 1:] ) + has_slice = any(isinstance(k, slice) for k in key) + has_array = any(isinstance(k, (numpy.ndarray, list)) for k in key) + + if has_slice and has_array: + raise NotImplementedError("Mixed slice and array indexing is not supported.") + # Slice indexing -> delegate to view() - if any(isinstance(k, slice) for k in key): + if has_slice: # view expects one slice for each axis/dim of the tensor if len(key) != len(self.shape): raise IndexError("Expects the number of slices to match the gem.Node tensor rank") return view(self, *key) + # Support a list or array of integer indices + # Old approach: build a ListTensor out of Indexed nodes, one for each element of the permutation + if has_array: + pos = next(i for i, k in enumerate(key) if isinstance(k, (numpy.ndarray, list))) + arr = numpy.asarray(key[pos]) + list_index = ListIndex(arr) + new_key = key[:pos] + (list_index,) + key[pos+1:] + indexed = Indexed(self, new_key) + return ComponentTensor(indexed, (list_index.free_index,)) # convert free index back to shape + # Point indexing return Indexed(self, key) @@ -718,6 +723,38 @@ def __repr__(self): def __reduce__(self): return type(self), (self.expression,) +class ListIndex(IndexBase): + """A lookup index in the form of an index array""" + + __slots__ = ('index_array', 'free_index',) + + def __init__(self, index_array): + assert isinstance(index_array, numpy.ndarray) + assert numpy.issubdtype(index_array.dtype, numpy.integer) + self.index_array = index_array + self.free_index = Index(extent=len(self.index_array)) + + def __eq__(self, other): + if type(self) is not type(other): + return False + return numpy.array_equal(self.index_array, other.index_array) + + def __ne__(self, other): + return not self.__eq__(other) + + def __hash__(self): + return hash((type(self), self.index_array.tobytes())) + + def __str__(self): + return str(self.index_array) + + def __repr__(self): + return "%r(%s)" % (type(self), self.index_array) + + def __reduce__(self): + return type(self), (self.index_array, ) + +IndexT = int | Index | VariableIndex | ListIndex | slice | EllipsisType | list | numpy.ndarray class Indexed(Scalar): __slots__ = ('children', 'multiindex', 'indirect_children') @@ -734,6 +771,9 @@ def __new__(cls, aggregate, multiindex): assert isinstance(index, IndexBase) if isinstance(index, Index): index.set_extent(extent) + elif isinstance(index, ListIndex): + if numpy.any(index.index_array < 0) or numpy.any(index.index_array >= extent): + raise IndexError("Invalid index in ListIndex") elif isinstance(index, int) and not (0 <= index < extent): raise IndexError("Invalid literal index") @@ -781,6 +821,8 @@ def __new__(cls, aggregate, multiindex): new_indices.append(i) elif isinstance(i, VariableIndex): new_indices.extend(i.expression.free_indices) + elif isinstance(i, ListIndex): + new_indices.append(i.free_index) self.free_indices = unique(aggregate.free_indices + tuple(new_indices)) return self @@ -793,6 +835,8 @@ def index_ordering(self): free_indices.append(i) elif isinstance(i, VariableIndex): free_indices.extend(i.expression.free_indices) + elif isinstance(i, ListIndex): + free_indices.append(i.free_index) return tuple(free_indices) diff --git a/gem/interpreter.py b/gem/interpreter.py index 70094c10..a55db78a 100644 --- a/gem/interpreter.py +++ b/gem/interpreter.py @@ -263,27 +263,41 @@ def _evaluate_conditional(e, self): def _evaluate_indexed(e, self): """Indexing maps shape to free indices""" val = self(e.children[0]) - fids = tuple(i for i in e.multiindex if isinstance(i, gem.Index)) + fids = tuple( + i if isinstance(i, gem.Index) else i.free_index + for i in e.multiindex + if isinstance(i, (gem.Index, gem.ListIndex)) + ) + all_fids = val.fids + fids idx = [] # First pick up all the existing free indices - for _ in val.fids: - idx.append(slice(None)) + for fid in val.fids: + # idx.append(slice(None)) + shape = tuple(fid.extent if f is fid else 1 for f in all_fids) + idx.append(numpy.arange(fid.extent).reshape(shape)) + # Now grab the shape axes for i in e.multiindex: if isinstance(i, gem.Index): # Free index, want entire extent - idx.append(slice(None)) + # idx.append(slice(None)) + shape = tuple(i.extent if f is i else 1 for f in all_fids) + idx.append(numpy.arange(i.extent).reshape(shape)) elif isinstance(i, gem.VariableIndex): # Variable index, evaluate inner expression result, = self(i.expression) assert not result.tshape idx.append(result[()]) + elif isinstance(i, gem.ListIndex): + # ListIndex, use the index array + shape = tuple(i.free_index.extent if f is i.free_index else 1 for f in all_fids) + idx.append(i.index_array.reshape(shape)) else: # Fixed index, just pick that value idx.append(i) assert len(idx) == len(val.tshape) - return Result(val[idx], val.fids + fids) + return Result(val[idx], all_fids) @_evaluate.register(gem.FlexiblyIndexed) From 05ea454315c34d9fd8c99c0cec57d3b71f83e99c Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Fri, 24 Apr 2026 15:43:31 +0100 Subject: [PATCH 20/21] final changes --- FIAT/reference_element.py | 13 +++++++++---- gem/optimise.py | 9 ++++++++- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 51f9581b..3e41ee3f 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1438,14 +1438,15 @@ def compute_axis_barycentric_coordinates(self, points, entity=None, rescale=Fals result = numpy.empty(len(self.cells), dtype=object) for k, (factor, s) in enumerate(zip(self.cells, point_slices)): result[k] = factor.compute_barycentric_coordinates(points[..., s], entity, rescale) - + # Flatten the array - # NOTE: cannot construct the flat array directly since we may not know upfront the total number + # We cannot construct the flat array directly since we may not know upfront the total number # of barycentric coordinates (e.g., in a simplex it is d+1, in a hypercube it is 2*d) flat_result = numpy.array([bary[j] for bary in result for j in range(bary.shape[0])]) return flat_result + class Hypercube(Cell): """Abstract class for a reference hypercube""" @@ -1903,8 +1904,12 @@ def compute_unflattening_map(topology_dict): def compute_facet_permutation(unflattening_map, product): """ - Return a permutation mapping each hypercube facet to the index of its - vanishing barycentric coordinate in the axis-structured barycentric array. + Returns a permutation mapping each facet of a `~.Hypercube` to the index of the + barycentric coordinate that vanishes on it. + + The order of barycentric coordinates returned by `compute_axis_barycentric_coordinates` + is determined by axis structure, not by facet numbering. Reordering them by this permutation + yields the invariant: the i-th barycentric coordinate vanishes on the i-th facet. """ # First compute axis offsets into the flattened barycentric coordinate array. axis_offsets = [] diff --git a/gem/optimise.py b/gem/optimise.py index caf254e0..d57f3ce1 100644 --- a/gem/optimise.py +++ b/gem/optimise.py @@ -13,7 +13,7 @@ reuse_if_untouched_arg, traversal) from gem.gem import (Node, Failure, Identity, Constant, Literal, Zero, Product, Sum, Comparison, Conditional, Division, - Index, VariableIndex, Indexed, FlexiblyIndexed, + Index, VariableIndex, ListIndex, Indexed, FlexiblyIndexed, IndexSum, ComponentTensor, ListTensor, Delta, partial_indexed, one) @@ -145,6 +145,13 @@ def replace_indices_indexed(node, self, subst): child = Literal(sub, dtype=child.dtype) if isinstance(child, Constant) else ListTensor(sub) multiindex = tuple(i for i in multiindex if not isinstance(i, Integral)) + elif isinstance(child, ListTensor) and len(multiindex) == 1 and isinstance(multiindex[0], ListIndex): + # ListIndex into a ListTensor: apply the permutation at compile time + # by reordering the ListTensor elements, eliminating the free index. + list_index = multiindex[0] + child = ListTensor(child.array[list_index.index_array]) + multiindex = (list_index.free_index,) + if multiindex == node.multiindex and child == node.children[0]: return node else: From 0cb039843cf8a66ffca5f123a15ff4a2b6165d21 Mon Sep 17 00:00:00 2001 From: Anastasia Chanbour Date: Fri, 24 Apr 2026 15:47:40 +0100 Subject: [PATCH 21/21] fixed formatting --- gem/gem.py | 23 +++++++++++++---------- gem/interpreter.py | 4 ++-- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/gem/gem.py b/gem/gem.py index dc921438..54fb060c 100644 --- a/gem/gem.py +++ b/gem/gem.py @@ -23,7 +23,6 @@ from numbers import Integral, Number from types import EllipsisType -import itertools import numpy from numpy import asarray @@ -46,6 +45,7 @@ uint_type = numpy.dtype(numpy.uintc) + class NodeMeta(type): """Metaclass of GEM nodes. @@ -92,7 +92,7 @@ def __getitem__( """A generalised interface for indexing GEM tensors""" if not isinstance(key, tuple): key = (key,) - + # Expand ellipsis -> fill in remaining dimensions with slice(None) if any(k is Ellipsis for k in key): if key.count(Ellipsis) > 1: @@ -128,7 +128,7 @@ def __getitem__( list_index = ListIndex(arr) new_key = key[:pos] + (list_index,) + key[pos+1:] indexed = Indexed(self, new_key) - return ComponentTensor(indexed, (list_index.free_index,)) # convert free index back to shape + return ComponentTensor(indexed, (list_index.free_index,)) # convert free index back to shape # Point indexing return Indexed(self, key) @@ -723,9 +723,10 @@ def __repr__(self): def __reduce__(self): return type(self), (self.expression,) + class ListIndex(IndexBase): """A lookup index in the form of an index array""" - + __slots__ = ('index_array', 'free_index',) def __init__(self, index_array): @@ -733,28 +734,30 @@ def __init__(self, index_array): assert numpy.issubdtype(index_array.dtype, numpy.integer) self.index_array = index_array self.free_index = Index(extent=len(self.index_array)) - + def __eq__(self, other): if type(self) is not type(other): return False return numpy.array_equal(self.index_array, other.index_array) - + def __ne__(self, other): return not self.__eq__(other) - + def __hash__(self): return hash((type(self), self.index_array.tobytes())) - + def __str__(self): return str(self.index_array) def __repr__(self): return "%r(%s)" % (type(self), self.index_array) - + def __reduce__(self): return type(self), (self.index_array, ) -IndexT = int | Index | VariableIndex | ListIndex | slice | EllipsisType | list | numpy.ndarray + +IndexT = int | Index | VariableIndex | ListIndex | slice | EllipsisType | list | numpy.ndarray + class Indexed(Scalar): __slots__ = ('children', 'multiindex', 'indirect_children') diff --git a/gem/interpreter.py b/gem/interpreter.py index a55db78a..d3f71f75 100644 --- a/gem/interpreter.py +++ b/gem/interpreter.py @@ -322,7 +322,7 @@ def _evaluate_flexiblyindexed(e, self): offset_val = offset_result[()] else: offset_val = offset - + dim_flat_idx_components = [] for i, s in idxs: # Index i may be one of: Index, VariableIndex, int @@ -332,7 +332,7 @@ def _evaluate_flexiblyindexed(e, self): i_vals = numpy.arange(i.extent).reshape(shape) * s dim_flat_idx_components.append(i_vals) elif isinstance(i, gem.VariableIndex): - # Variable index gives a single integer index + # Variable index gives a single integer index # obtained by evaluating its inner expression result, = self(i.expression) assert not result.tshape