Computing barycentric coordinates in Tensor Product Cells by achanbour · Pull Request #246 · firedrakeproject/fiat

achanbour · 2026-05-08T09:20:08Z

This is the first part of PR #245 concerning changes to FIAT only. In particular, I have extended FIAT to enable the computation of barycentric coordinates on tensor product cells.

The main distinction lies in how barycentric coordinates are represented in simplicies vs tensor product cells. For simplicies, there is one barycentric coordinate per vertex so the default ordering of these coordinates follows the vertex ordering. For tensor product cells, there is one barycentric coordinate per facet so the natural ordering is that of facets.

First, to enable facet ordering in simplicies, I have added a facet ordering permutation in SimplicialComplex.compute_barycentric_coordinates. Then, since barycentric coordinates in tensor product cells are computed recursively on each factor, having facet ordering on the simplicies that form each factor ensure that we obtain facet ordering on the tensor-product cell as a whole.

pbrubeck · 2026-05-08T11:45:17Z

Looks like you are including commits form #240. Are they needed?

…cells minor fixes modified compute_barycentric_coords to consistently handle input points of different sizes extended the tensor product barycentric coordinates computation function to handle nested tensor-product cells + added tests removed line enforcing points to be numpy arrays to ensure code generation works properly for symbolic points added conversion to numpy array only in axis_barycentric_coords generalised numpy operations in compute_barycentric_coordinates methods to work on different array shapes changed compute_barycentric_coordinates method to work on input GEM expressions representing point sets extended gem matmul to convert numpy arrays to Literals and deal separately with singleton contractions changes made to gem and FIAT to make barycentric coordinate computations work on GEM tensor expressions extended gem's slicing syntax and simplified the code in compute_axis_barycentric_coordinates in tp cells tidied up and added comments recent changes post merging compute barycentric coords symbolically in simplicies and hypercubes (to be tested) modified unit tests for computing bary coords on points passed as numpy arrays latest changes + gem tests implemented handler in gem.evaluate for FlexiblyIndexed nodes added a new ListIndex type to GEM for indexing GEM tensors using an index array final changes fixed formatting

changed ListIndex to be a subclass of Index

separated ListIndex from Index

…plicies and tensor product cells

… sub-entities of simplicies + added notes + modified test checking facet ordering of bary coords on tp cells modified the index substitution logic to correctly handle ListIndex changed index substitution optim isation

achanbour · 2026-05-12T17:53:39Z

Looks like you are including commits form #240. Are they needed?

No. I had accidentally merged some irrelevant changes into the old branch from which this PR was branched off. I've now reviewed the changes and removed unnecessary commits.

pbrubeck · 2026-05-13T10:10:24Z

+        subcomplex = global_indices
+
        if entity_dim != sd:


Suggested change

subcomplex = global_indices

if entity_dim != sd:

if entity_dim == sd:

subcomplex = global_indices

else:

pbrubeck · 2026-05-13T10:10:37Z

            cell_id = self.connectivity[(entity_dim, sd)][entity_id][0]
+
+            # restrict indices to the chosen sub-cell
            indices = [i for i, v in enumerate(top[sd][cell_id]) if v in subcomplex]


Suggested change

indices = [i for i, v in enumerate(top[sd][cell_id]) if v in subcomplex]

indices = [i for i, v in enumerate(top[sd][cell_id]) if v in global_indices]

pbrubeck · 2026-05-13T10:13:05Z

+            facet_ids = self.connectivity[(entity_dim, entity_dim-1)][entity_id]
+            facets = [top[entity_dim-1][f] for f in facet_ids] # facet as a tuple of vertex indices that form it
+
+            perm = [facets.index(tuple(sorted(set(global_indices) - {v}))) for v in global_indices]


I think we need to permute indices with the inverse of perm. facets.index returns a facet index, but indices refers to vertices.

pbrubeck · 2026-05-13T10:15:04Z

+        cell_verts = self.get_vertices_of_subcomplex(subcomplex) # get vertex coordinates of the sub-cell
+        ref_verts = numpy.eye(sd + 1) # get vertex coordinates of the standard reference cell


To keep a cleaner diff, it is better to not add extra comments on unmodified lines

Suggested change

cell_verts = self.get_vertices_of_subcomplex(subcomplex) # get vertex coordinates of the sub-cell

ref_verts = numpy.eye(sd + 1) # get vertex coordinates of the standard reference cell

# Get vertex coordinates of the physical and reference sub-cells

cell_verts = self.get_vertices_of_subcomplex(subcomplex)

ref_verts = numpy.eye(sd + 1)

pbrubeck · 2026-05-13T10:17:59Z

+        result = []
+        for factor, s in zip(self.cells, point_slices):
+            factor_sd = factor.get_spatial_dimension()
+            if factor_sd == 1:
+                factor_bary_coords = factor.compute_barycentric_coordinates(points[..., s], entity, rescale, facet_ordering=True)
+            else:
+                # For higher dimensional simplicies: vertex ordering is already guaranteed.
+                # NOTE: But what if we want to compute barycentric coordinates on a sub-entity and still want facet ordering?
+                # Then we need to pass facet_ordering=True even when factor_sd > 1
+                factor_bary_coords = factor.compute_barycentric_coordinates(points[..., s], entity, rescale)
+            result.append(factor_bary_coords)


Suggested change

result = []

for factor, s in zip(self.cells, point_slices):

factor_sd = factor.get_spatial_dimension()

if factor_sd == 1:

factor_bary_coords = factor.compute_barycentric_coordinates(points[..., s], entity, rescale, facet_ordering=True)

else:

# For higher dimensional simplicies: vertex ordering is already guaranteed.

# NOTE: But what if we want to compute barycentric coordinates on a sub-entity and still want facet ordering?

# Then we need to pass facet_ordering=True even when factor_sd > 1

factor_bary_coords = factor.compute_barycentric_coordinates(points[..., s], entity, rescale)

result.append(factor_bary_coords)

result = [factor.compute_barycentric_coordinates(points[..., s], entity, rescale, facet_ordering=True)

for factor, s in zip(self.cells, point_slices)]

pbrubeck · 2026-05-13T10:30:47Z

+            facets = [top[entity_dim-1][f] for f in facet_ids] # facet as a tuple of vertex indices that form it
+
+            perm = [facets.index(tuple(sorted(set(global_indices) - {v}))) for v in global_indices]
+            indices = [indices[p] for p in perm]


Suggested change

indices = [indices[p] for p in perm]

indices = [indices[p] for p in numpy.argsort(perm)]

pbrubeck · 2026-05-13T10:34:05Z

    return unflattening_map


+def compute_facet_permutation(unflattening_map, product):


We should get rid of this function if not used anymore

pbrubeck · 2026-05-13T10:34:41Z

    def __le__(self, other):
        return self.product <= other

+    def compute_barycentric_coordinates(self, points, entity=None, rescale=False):


Suggested change

def compute_barycentric_coordinates(self, points, entity=None, rescale=False):

def compute_barycentric_coordinates(self, points, entity=None, rescale=False, facet_ordering=True):

pbrubeck · 2026-05-13T10:36:28Z

+
+        Returns
+        -------
+        List of numpy.ndarray or GEM.ComponentTensor


Not quite true, so far we return a single numpy.ndarray

pbrubeck · 2026-05-13T10:41:30Z

-        # get a subcell containing the entity and the restriction indices of the entity
-        indices = slice(None)
-        subcomplex = top[entity_dim][entity_id]
+        # indices = slice(None)


Suggested change

# indices = slice(None)

# get a subcell containing the entity and the restriction indices of the entity

pbrubeck · 2026-05-13T10:43:38Z

+                          (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):


This test could pass at the moment, we just need to cast back and forth between gem and numpy

pbrubeck · 2026-05-13T10:44:29Z

+                          (hexahedron, [0.3, 0.4, 1.0]),])
+def test_hypercube_bary_coords_are_in_facet_order(cell, point, epsilon=1e-12):
+    """Test that the barycentric coordinates computed on hypercubes are listed in facet order."""
+    point = np.asarray(point)


we should support a single or multiple points as tuple or numpy.array

pbrubeck · 2026-05-13T10:45:33Z

+    coords = gem.Variable('X', (sd,))
+    bindings = {coords: point}
+
+    bary_gem = cell.compute_barycentric_coordinates(coords)


Suggested change

bary_gem = cell.compute_barycentric_coordinates(coords)

np_coords = tuple(coords[i] for i in range(sd))

bary_gem = gem.as_gem(cell.compute_barycentric_coordinates(np_coords))

pbrubeck · 2026-05-13T14:36:06Z

+                factor_bary_coords = factor.compute_barycentric_coordinates(points[..., s], entity, rescale)
+            result.append(factor_bary_coords)
+
+        flat_result = numpy.concatenate(result)


Suggested change

flat_result = numpy.concatenate(result)

flat_result = numpy.concatenate(result, axis=-1)

pbrubeck · 2026-05-13T14:37:10Z

+        if isinstance(points, numpy.ndarray) and len(points) == 0:
+            return points
+


Suggested change

if isinstance(points, numpy.ndarray) and len(points) == 0:

return points

pbrubeck · 2026-05-13T14:38:21Z

+        if isinstance(points, numpy.ndarray) and len(points) == 0:
+            return points
+
+        tp_bary_coords = self.product.compute_factor_barycentric_coordinates(points, entity, rescale)


Suggested change

tp_bary_coords = self.product.compute_factor_barycentric_coordinates(points, entity, rescale)

return self.product.compute_factor_barycentric_coordinates(points, entity, rescale)

achanbour requested a review from pbrubeck May 8, 2026 09:20

connorjward mentioned this pull request May 8, 2026

Extending the indexing of GEM nodes #247

Draft

achanbour force-pushed the achanbour/bary-coords-fiat branch from 448c540 to d71fa37 Compare May 12, 2026 17:38

achanbour added 9 commits May 12, 2026 18:45

added docstrings to tests

4564251

changed ListIndex to be a subclass of Index

renamed tp bary coords methods + minor fixes + updated tests

8b1add2

separated ListIndex from Index

added explanation to facet permutation to reorder barycentric coords

3a02fd6

enforced facet ordering in the barycentric coordinates across all sim…

ad9c358

…plicies and tensor product cells

undid changes to gem

d4fe50f

removed gem dependence from FIAT

7395468

restore gem to main

a3d51a1

achanbour force-pushed the achanbour/bary-coords-fiat branch from d71fa37 to a3d51a1 Compare May 12, 2026 17:49

pbrubeck reviewed May 13, 2026

View reviewed changes

	indices = [i for i, v in enumerate(top[sd][cell_id]) if v in subcomplex]
	indices = [i for i, v in enumerate(top[sd][cell_id]) if v in global_indices]

		cell_verts = self.get_vertices_of_subcomplex(subcomplex) # get vertex coordinates of the sub-cell
		ref_verts = numpy.eye(sd + 1) # get vertex coordinates of the standard reference cell

	indices = [indices[p] for p in perm]
	indices = [indices[p] for p in numpy.argsort(perm)]

		return unflattening_map


		def compute_facet_permutation(unflattening_map, product):

	def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
	def compute_barycentric_coordinates(self, points, entity=None, rescale=False, facet_ordering=True):

	# indices = slice(None)
	# get a subcell containing the entity and the restriction indices of the entity

	bary_gem = cell.compute_barycentric_coordinates(coords)
	np_coords = tuple(coords[i] for i in range(sd))
	bary_gem = gem.as_gem(cell.compute_barycentric_coordinates(np_coords))

	flat_result = numpy.concatenate(result)
	flat_result = numpy.concatenate(result, axis=-1)

		if isinstance(points, numpy.ndarray) and len(points) == 0:
		return points

	tp_bary_coords = self.product.compute_factor_barycentric_coordinates(points, entity, rescale)
	return self.product.compute_factor_barycentric_coordinates(points, entity, rescale)

Conversation

achanbour commented May 8, 2026

Uh oh!

pbrubeck commented May 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

achanbour commented May 12, 2026

Uh oh!

pbrubeck May 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pbrubeck May 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pbrubeck May 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pbrubeck May 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

pbrubeck commented May 8, 2026 •

edited

Loading

pbrubeck May 13, 2026 •

edited

Loading

pbrubeck May 13, 2026 •

edited

Loading

pbrubeck May 13, 2026 •

edited

Loading

pbrubeck May 13, 2026 •

edited

Loading