diff --git a/.gitignore b/.gitignore index b4109385..86451a5a 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,7 @@ /build/ /dist/ -/fenics_fiat.egg-info/ +/firedrake_fiat.egg-info/ /.cache/ /doc/sphinx/source/api-doc diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 87794b50..b00811e6 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -613,35 +613,56 @@ def get_dimension(self): spatial dimension.""" return self.get_spatial_dimension() - def compute_barycentric_coordinates(self, points, entity=None, rescale=False): + def compute_barycentric_coordinates(self, points, entity=None, rescale=False, facet_ordering=False): """Returns the barycentric coordinates of a list of points on the complex.""" - if len(points) == 0: - return points + if entity is None: entity = (self.get_spatial_dimension(), 0) + entity_dim, entity_id = entity top = self.get_topology() sd = self.get_spatial_dimension() - # 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) + indices = list(range(sd+1)) + global_indices = top[entity_dim][entity_id] # indices of all vertices that form the sub-entity + subcomplex = global_indices + if entity_dim != sd: + # get a sub-cell containing the sub-entity 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] + + # get the indices of all vertices that form the sub-cell subcomplex = top[sd][cell_id] - cell_verts = self.get_vertices_of_subcomplex(subcomplex) - ref_verts = numpy.eye(sd + 1) - A, b = make_affine_mapping(cell_verts, ref_verts) + if facet_ordering: + # Permute indices in facet order + # i.e., for each vertex in turn, get the position of the facet that excludes it + 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] + indices = [indices[p] for p in perm] + + 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 + A, b = make_affine_mapping(cell_verts, ref_verts) + A, b = A[indices], b[indices] if rescale: # rescale barycentric coordinates by the height wrt. to the facet h = 1 / numpy.linalg.norm(A, axis=1) b *= h A *= h[:, None] - out = numpy.dot(points, A.T) - return numpy.add(out, b, out=out) + + # 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 @@ -1406,6 +1427,58 @@ def extrinsic_orientation_permutation_map(self): def is_macrocell(self): return any(c.is_macrocell() for c in self.cells) + def compute_factor_barycentric_coordinates(self, points, entity=None, rescale=False): + """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. + + Returns + ------- + numpy.ndarray + 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 + + + axis_dims = [c.get_spatial_dimension() for c in self.cells] + point_slices = TensorProductCell._split_slices(axis_dims) + + # NOTE: entity is a tuple key of the tensor product cell's topology. It should be decomposed into per-factor sub-entities + # before being passed to factor.compute_barycentric_coordinates + if entity is not None: + raise NotImplementedError("Tensor-product sub-entity decomposition is not implemented yet.") + + 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) + + if isinstance(points, gem.Node): + # Flatten the array + # 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])]) + # returns a ListTensor wrapping the scalar GEM expr. of bary coords. + flat_result = gem.as_gem(flat_result) + else: + flat_result = numpy.concatenate(result) + + return flat_result + class Hypercube(Cell): """Abstract class for a reference hypercube""" @@ -1423,6 +1496,8 @@ def __init__(self, dimension, product): self.product = product self.unflattening_map = compute_unflattening_map(pt) + # 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 spatial dimension.""" @@ -1521,6 +1596,28 @@ 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. + + 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. + """ + if isinstance(points, numpy.ndarray) and len(points) == 0: + return points + + tp_bary_coords = self.product.compute_factor_barycentric_coordinates(points, entity, rescale) + + # return tp_bary_coords[self.facet_perm] + return tp_bary_coords + class UFCHypercube(Hypercube): """Reference UFC Hypercube @@ -1839,6 +1936,90 @@ def compute_unflattening_map(topology_dict): return unflattening_map +def compute_facet_permutation(unflattening_map, product): + """ + Returns a permutation mapping each facet of a Hypercube to the index of the + barycentric coordinate that vanishes on it. + + Let's take the example of a quad in 2D. Calling `compute_factor_barycentric_coordinates` returns + the barycentric coordinates on each of its 2 axes: + + axis 0 (x): lambda_x_1, lambda_x_2 + axis 1 (y): lambda_y_1, lambda_y_2 + + as a flat array bary_coords = [lambda_x_1, lambda_x_2, lambda_y_1, lambda_y_2] + + A quad has 4 facets which are numbered in UFC order as: + + facet 3 + ┌───────┐ + │ │ + facet 0 │ │ facet 1 + │ │ + │ │ + └───────┘ + facet 2 + + Since each axis is a UFCInterval (a simplex), with vertices at P1 = (0,) and P2 = (1,) its barycentric coordinates + are lambda_1 = 1 - t (vanishes at P2) and lambda_2 = t (vanishes at P1). This applies the rule for barycentric coordinates + on simplicies which is that lambda_i vanishes on the facet opposite vertex i. + + Therefore: + + - lambda_x_1 vanishes on facet 1 (x=1), lambda_x_2 vanishes on facet 0 (x=0) + - lambda_y_1 vanishes on facet 3 (y=1), lambda_y_2 vanishes on facet 2 (y=0) + + The permutation computed in this function reorders the array of barycentric coordinates such that: + + bary_coords[perm] = [lambda_x_2, lambda_x_1, lambda_y_2, lambda_y_1] + + where the i-th entry corresponds to the barycentric coordinate vanishing on facet i. + """ + # 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 integer permutation array + sd = len(product.cells) + num_facets = 2 * sd + 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 + dim_tuple, tp_entity = unflattening_map[(sd - 1, f)] + + # Determine the axis that's orthogonal to the facet + # 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 + # 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) + ) + tuple_ei = numpy.unravel_index(tp_entity, entity_shape) + local_facet = tuple_ei[axis] + + # 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] + bary_index + + return perm + + def max_complex(complexes): max_cell = max(complexes) if all(max_cell >= b for b in complexes): diff --git a/gem/gem.py b/gem/gem.py index 5c9231dc..9cb9f1f1 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 @@ -32,7 +36,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', @@ -81,12 +85,53 @@ def is_equal(self, other): self.children = other.children return result - def __getitem__(self, indices): - try: - indices = tuple(indices) - except TypeError: - indices = (indices, ) - return Indexed(self, indices) + def __getitem__( + self, + key: IndexT | tuple[IndexT, ...], + ) -> ComponentTensor | Indexed: + """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: + raise NotImplementedError("Multiple ellipses are not supported.") + ellipsis_pos = key.index(Ellipsis) + remaining_dims = len(self.shape) - (len(key) - 1) + if remaining_dims < 0: + raise IndexError("Too many indices provided.") + key = ( + key[:ellipsis_pos] + + (slice(None), ) * remaining_dims + + 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 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 + # Previously, built a ListTensor out of Indexed nodes, one for each element of the permutation + if has_array: + arr_pos = next(i for i, k in enumerate(key) if isinstance(k, (numpy.ndarray, list))) + list_index = ListIndex(key[arr_pos]) + new_key = key[:arr_pos] + (list_index,) + key[arr_pos+1:] + indexed = Indexed(self, new_key) # A[perm[i]] is a scalar + retval = ComponentTensor(indexed, (list_index.free_index,)) # CT(A[perm[i]], i) -> A[perm] + return retval + + # Point indexing + return Indexed(self, key) def __neg__(self): return componentwise(Product, minus, self) @@ -117,6 +162,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))) @@ -677,6 +723,52 @@ def __repr__(self): def __reduce__(self): return type(self), (self.expression,) +class ListIndex(IndexBase): + """ + Option 1: tensor[list_index] is tensor[index_array[list_index]] free index ranging over 0,1..., len(index_array) + + Option 2: tensor[list_index] is tensor[list_index.index_array[list_index.free_index]] + """ + + __slots__ = ('index_array', 'free_index', 'name') + + def __init__(self, index_array, name=None, free_index=None): + index_array = numpy.asarray(index_array) + assert numpy.issubdtype(index_array.dtype, numpy.integer) + # Wraps a free index together with the index array + self.index_array = index_array + if free_index is not None: + self.free_index = free_index + else: + self.free_index = Index(extent=len(index_array), name=name) + + def __str__(self): + if self.free_index.name: + return f"{self.index_array.tolist()}[i_{self.free_index.name}]" + return f"{self.index_array.tolist()}[i_{self.free_index.count}]" + + def __repr__(self): + if self.free_index.name: + return f"ListIndex({self.index_array.tolist()}, Index({self.free_index.name}))" + return f"ListIndex({self.index_array.tolist()}, Index({self.free_index.count}))" + + # 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 __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') @@ -716,9 +808,16 @@ def __new__(cls, aggregate, multiindex): kk = B.multiindex ff = C.free_indices if not any((j in ff) for j in jj): - # Only replace indices that are not present in C + # Only replace indices in kk that are in jj rep = dict(zip(jj, ii)) - ll = tuple(rep.get(k, k) for k in kk) + def _replace(k): + if isinstance(k, ListIndex) and k.free_index in rep: + replacement = rep[k.free_index] + if isinstance(replacement, int): + return int(k.index_array[replacement]) + return ListIndex(k.index_array, free_index=rep[k.free_index]) + return rep.get(k,k) + ll = tuple(_replace(k) for k in kk) aggregate = C multiindex = ll @@ -740,6 +839,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 @@ -752,6 +853,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) @@ -869,7 +972,7 @@ def __new__(cls, expression, multiindex): shape = tuple(index.extent for index in multiindex) assert all(s >= 0 for s in shape) - # Zero folding + # Zero foldingc if isinstance(expression, Zero): return Zero(shape, dtype=expression.dtype) @@ -1416,4 +1519,4 @@ def Piecewise(*args): expr = Literal(float("nan")) for v, c in reversed(pieces): expr = Conditional(c, v, expr) - return expr + return expr \ No newline at end of file diff --git a/gem/interpreter.py b/gem/interpreter.py index 13eeb44a..4def9bc8 100644 --- a/gem/interpreter.py +++ b/gem/interpreter.py @@ -262,28 +262,99 @@ def _evaluate_conditional(e, self): @_evaluate.register(gem.Indexed) 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)) + val = self(e.children[0]) # tensor being indexed + # fids = tuple(i for i in e.multiindex if isinstance(i, (gem.Index, gem.ListIndex))) # free indices in the multiindex + fids = tuple( + i.free_index if isinstance(i, gem.ListIndex) else i + 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)) + + # idx has one entry per axis but val.arr is a single array with shape (*val.fids extents, *tensor_shape) + # and the inner loop consumes the tensor_shape axes one by ones through indexing. + # When val.arr is indexed with idx as val[idx], NumPy broadcasts each entry of idx to operate + # on all axes simultaneously so we need to reshape each entry of idx appropriately. + + # First pick up all the free indices that's already part of the tensor + 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): + if 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)) + elif 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 + assert not result.tshape # int idx.append(result[()]) 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) +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) diff --git a/gem/optimise.py b/gem/optimise.py index caf254e0..528b01be 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) @@ -102,7 +102,14 @@ def _replace_indices_atomic(i, self, subst): return i if new_expr == i.expression else VariableIndex(new_expr) else: substitute = dict(subst) - return substitute.get(i, i) + def _replace(i): + if isinstance(i, ListIndex) and i.free_index in substitute: + replacement = substitute[i.free_index] + if isinstance(replacement, int): + return int(i.index_array[replacement]) + return ListIndex(i.index_array, free_index=substitute[i.free_index]) + return substitute.get(i,i) + return _replace(i) @replace_indices.register(Delta) @@ -120,6 +127,9 @@ def replace_indices_indexed(node, self, subst): multiindex = tuple(_replace_indices_atomic(i, self, subst) for i in node.multiindex) child, = node.children + # When node is Indexed(CT(A[j], (j,)), (i,)) with j and i as Index the below rule promotes the substitution j -> i + # so that A[j] becomes A[i] which eliminates the CT. This is wrong for when j a free index of ListIndex since A[j] is A[index_array[j]]; + # hence doing the subsitution effectively gets rid of index_array. if isinstance(child, ComponentTensor): # Indexing into ComponentTensor # Inline ComponentTensor and augment the substitution rules diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index d166ba6e..96fb361c 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -455,6 +455,173 @@ def test_flatten_maintains_ufc_status(cell): assert ufc_status == is_ufc(flat_cell) +@pytest.mark.parametrize(('cell', 'point'), + [(interval, [0.0]), # on facet 0 + (interval, [1.0]),]) # on facet 1 +def test_bary_coords_order_on_1d_simplex(cell, point, epsilon=1e-12): + facet_hits = [fid for fid, pts in cell.point_entity_ids([point])[0].items() if len(pts) > 0] + assert len(facet_hits) == 1 # lies on one facet + facet_id = facet_hits[0] + + # The default ordering of barycentric coords. is vertex ordering + # On a 1D interval: facets are vertices so they follow the vertex ordering + # But bary_coords[i] vanishes on the facet not containing i. In this case, this happens to be facet 1-i. + bary_default = cell.compute_barycentric_coordinates(point, facet_ordering=False) + assert np.isclose(bary_default[1 - facet_id], 0.0, atol=epsilon) + assert not np.isclose(bary_default[facet_id], 0.0, atol=epsilon) + + # Reverting the ordering of vertices inside `compute_barycentric_coordinates` gives us the facet ordering + # i.e, bary_coords[i] vanishes on the facet not containing i which happens to be facet i. + bary_ordered = cell.compute_barycentric_coordinates(point, facet_ordering=True) + assert np.isclose(bary_ordered[facet_id], 0.0, atol=epsilon) + + +@pytest.mark.parametrize(('cell', 'point'), + [(triangle, [0.5, 0.5]), # on facet 0 + (triangle, [0.0, 0.5]), # on facet 1 + (triangle, [0.5, 0.0]), # on facet 2 + (tetrahedron, [1/3, 1/3, 1/3]), # on facet 0 + (tetrahedron, [0.0, 0.25, 0.25]), # on facet 1 + (tetrahedron, [0.25, 0.0, 0.25]), # on facet 2 + (tetrahedron, [0.25, 0.25, 0.0]),]) # on facet 3 +def test_bary_coords_order_on_simplicies(cell, point, epsilon=1e-12): + sd = cell.get_spatial_dimension() + facet_hits = [fid for fid, pts in cell.point_entity_ids([point])[sd-1].items() if len(pts) > 0][0] + assert len(facet_hits) == 1 # lies on one facet + facet_id = facet_hits[0] + + # The default ordering of barycentric coords. is vertex ordering + # In triangles and tetrahedrons, facets are ordered by the vertex they include + # so: bary_coords[i] vanishes on facet not containing vertex i which is facet i + bary = cell.compute_barycentric_coordinates(point, facet_ordering=False) + assert np.isclose(bary[facet_id], 0.0, atol=epsilon) + + +@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_tp_bary_coords_are_in_factor_order(cell, point, epsilon=1e-12): + """Test that barycentric coordinates computed on tensor product cells are listed in factor order.""" + point = np.asarray(point) + axis_bary_coords = cell.compute_factor_barycentric_coordinates(point) + + assert type(axis_bary_coords) is np.ndarray + + offset = 0 + bary_offset = 0 + for factor in cell.cells: + sd = factor.get_spatial_dimension() + coords_k = point[offset: offset + sd] + offset += sd + if factor is interval: + expected = factor.compute_barycentric_coordinates(coords_k, facet_ordering=True) + else: + expected = factor.compute_barycentric_coordinates(coords_k) + 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'), + [(interval_x_interval, [0.0, 0.5]), + (interval_x_interval, [1.0, 0.5]), + (interval_x_interval, [0.5, 0.0]), + (interval_x_interval, [0.5, 1.0]), + (triangle_x_interval, [0.0, 0.5, 0.5]), + (triangle_x_interval, [0.5, 0.0, 0.5]), + (triangle_x_interval, [0.25, 0.25, 0.0]), + (triangle_x_interval, [0.25, 0.25, 1.0]), + (quadrilateral_x_interval, [0.0, 0.5, 0.5]), + (quadrilateral_x_interval, [1.0, 0.5, 0.5]), + (quadrilateral_x_interval, [0.5, 0.0, 0.5]), + (quadrilateral_x_interval, [0.5, 1.0, 0.5]), + (quadrilateral_x_interval, [0.5, 0.5, 0.0]), + (quadrilateral_x_interval, [0.5, 0.5, 1.0]),]) +def test_tp_bary_coords_are_in_facet_order(cell, point, epsilon=1e-12): + """Test that barycentric coordinates computed on tensor product cells are listed in facet order.""" + point = np.asarray(point) + bary_coords = cell.compute_factor_barycentric_coordinates(point) + + sd = cell.get_spatial_dimension() + top = cell.get_topology() + + # facets are entities where the sum of the dimension tuple is sd - 1 + ordered_facets = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim]) + if sum(dim) == sd - 1] + + point_entity_ids = cell.point_entity_ids([point]) + facet_hits = [i for i, (dim, entity) in enumerate(ordered_facets) + if len(point_entity_ids[dim][entity]) > 0] + + assert len(facet_hits) == 1 # lies on one facet + + facet_id = facet_hits[0] + 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) + + +@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_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) + + 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 # lies on one facet + + facet_id = facet_hits[0] + 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) + + +@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): + """Test the GEM expression for barycentric coordinates produces the expected numeric output.""" + 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) + + # If bary_gem is a numpy array convert to a GEM Node using before evaluating + # results, = evaluate((gem.as_gem(bary_gem),), bindings=bindings) + + results, = evaluate((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__))