From 28d5f1c8b3779128c68c84bd164bf90ad5675f37 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 1 May 2024 10:32:31 +0100 Subject: [PATCH 01/14] cleanup, add UFCHypercube class --- FIAT/reference_element.py | 231 +++++++++++++------------------------- 1 file changed, 81 insertions(+), 150 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 8ad5bd1a5..f04339fec 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -587,7 +587,7 @@ def construct_subelement(self, dimension): def contains_point(self, point, epsilon=0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -606,7 +606,7 @@ def contains_point(self, point, epsilon=0): def distance_to_point_l1(self, point): # noqa: D301 - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + """Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. Parameters @@ -617,7 +617,7 @@ def distance_to_point_l1(self, point): Returns ------- float - The L1 distance, also known as taxicab, manhatten or rectilinear + The L1 distance, also known as taxicab, manhattan or rectilinear distance, of the cell to the point. If 0.0 the point is inside the cell. @@ -1058,7 +1058,7 @@ def compute_reference_normal(self, facet_dim, facet_i): def contains_point(self, point, epsilon=0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -1083,7 +1083,7 @@ def contains_point(self, point, epsilon=0): True) def distance_to_point_l1(self, point): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + """Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" @@ -1103,6 +1103,8 @@ def cell_orientation_reflection_map(self): return make_cell_orientation_reflection_map_tensorproduct(self.cells) def compare(self, op, other): + """Parent-based comparison between simplicial complexes. + This is done dimension by dimension.""" if hasattr(other, "product"): other = other.product if isinstance(other, type(self)): @@ -1123,59 +1125,21 @@ def __le__(self, other): return self.compare(operator.le, other) -class UFCQuadrilateral(Cell): - r"""This is the reference quadrilateral with vertices - (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). - - Orientation of a physical cell is computed systematically - by comparing the canonical orderings of its facets and - the facets in the FIAT reference cell. - - As an example, we compute the orientation of a - quadrilateral cell: - - +---3---+ +--57---+ - | | | | - 0 1 43 55 - | | | | - +---2---+ +--42---+ - FIAT canonical Mapped example physical cell +class UFCHypercube(Cell): + """Abstract class for a reference unit hypercube [0, 1]^d with vertices in + lexicographical order.""" + dimension = None + shape = None - Suppose that the facets of the physical cell - are canonically ordered as: - - C = [55, 42, 43, 57] - - FIAT index to Physical index map must be such that - C[0] = 55 is mapped to a vertical facet; in this - example it is: - - M = [43, 55, 42, 57] - - C and M are decomposed into "vertical" and "horizontal" - parts, keeping the relative orders of numbers: - - C -> C0 = [55, 43], C1 = [42, 57] - M -> M0 = [43, 55], M1 = [42, 57] - - Then the orientation of the cell is computed as the - following: - - C0.index(M0[0]) = 1; C0.remove(M0[0]) - C0.index(M0[1]) = 0; C0.remove(M0[1]) - C1.index(M1[0]) = 0; C1.remove(M1[0]) - C1.index(M1[1]) = 0; C1.remove(M1[1]) - - o = 2 * 1 + 0 = 2 - """ def __init__(self): - product = TensorProductCell(UFCInterval(), UFCInterval()) + cells = [UFCInterval()] * self.dimension + product = TensorProductCell(*cells) pt = product.get_topology() verts = product.get_vertices() topology = flatten_entities(pt) - super(UFCQuadrilateral, self).__init__(QUADRILATERAL, verts, topology) + super(UFCHypercube, self).__init__(self.shape, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) @@ -1191,14 +1155,13 @@ def construct_subelement(self, dimension): :arg dimension: subentity dimension (integer) """ - if dimension == 2: + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError("Invalid dimension: %d" % (dimension,)) + elif dimension == sd: return self - elif dimension == 1: - return UFCInterval() - elif dimension == 0: - return Point() else: - raise ValueError("Invalid dimension: %d" % (dimension,)) + return ufc_hypercube(dimension) def get_entity_transform(self, dim, entity_i): """Returns a mapping of point coordinates from the @@ -1216,13 +1179,14 @@ def volume(self): def compute_reference_normal(self, facet_dim, facet_i): """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 1 + sd = self.get_spatial_dimension() + assert facet_dim == sd - 1 d, i = self.unflattening_map[(facet_dim, facet_i)] return self.product.compute_reference_normal(d, i) def contains_point(self, point, epsilon=0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -1240,14 +1204,14 @@ def contains_point(self, point, epsilon=0): return self.product.contains_point(point, epsilon=epsilon) def distance_to_point_l1(self, point): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + """Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" return self.product.distance_to_point_l1(point) def symmetry_group_size(self, dim): - return [1, 2, 8][dim] + return [1, 2, 8, 48][dim] def cell_orientation_reflection_map(self): """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" @@ -1266,109 +1230,61 @@ def __le__(self, other): return self.product <= other -class UFCHexahedron(Cell): - """This is the reference hexahedron with vertices - (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), - (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" - - def __init__(self): - product = TensorProductCell(UFCInterval(), UFCInterval(), UFCInterval()) - pt = product.get_topology() - - verts = product.get_vertices() - topology = flatten_entities(pt) - - super(UFCHexahedron, self).__init__(HEXAHEDRON, verts, topology) - - self.product = product - self.unflattening_map = compute_unflattening_map(pt) - - def get_dimension(self): - """Returns the subelement dimension of the cell. Same as the - spatial dimension.""" - return self.get_spatial_dimension() - - def construct_subelement(self, dimension): - """Constructs the reference element of a cell subentity - specified by subelement dimension. - - :arg dimension: subentity dimension (integer) - """ - if dimension == 3: - return self - elif dimension == 2: - return UFCQuadrilateral() - elif dimension == 1: - return UFCInterval() - elif dimension == 0: - return Point() - else: - raise ValueError("Invalid dimension: %d" % (dimension,)) - - def get_entity_transform(self, dim, entity_i): - """Returns a mapping of point coordinates from the - `entity_i`-th subentity of dimension `dim` to the cell. - - :arg dim: entity dimension (integer) - :arg entity_i: entity number (integer) - """ - d, e = self.unflattening_map[(dim, entity_i)] - return self.product.get_entity_transform(d, e) +class UFCQuadrilateral(UFCHypercube): + r"""This is the reference quadrilateral with vertices + (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). - def volume(self): - """Computes the volume in the appropriate dimensional measure.""" - return self.product.volume() + Orientation of a physical cell is computed systematically + by comparing the canonical orderings of its facets and + the facets in the FIAT reference cell. - def compute_reference_normal(self, facet_dim, facet_i): - """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 2 - d, i = self.unflattening_map[(facet_dim, facet_i)] - return self.product.compute_reference_normal(d, i) + As an example, we compute the orientation of a + quadrilateral cell: - def contains_point(self, point, epsilon=0): - """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', - 'taxicab' or rectilinear distance) to the cell. + +---3---+ +--57---+ + | | | | + 0 1 43 55 + | | | | + +---2---+ +--42---+ + FIAT canonical Mapped example physical cell - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. + Suppose that the facets of the physical cell + are canonically ordered as: - Returns - ------- - bool : True if the point is inside the cell, False otherwise. + C = [55, 42, 43, 57] - """ - return self.product.contains_point(point, epsilon=epsilon) + FIAT index to Physical index map must be such that + C[0] = 55 is mapped to a vertical facet; in this + example it is: - def distance_to_point_l1(self, point): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear - distance) to a point with 0.0 if the point is inside the cell. + M = [43, 55, 42, 57] - For more information see the docstring for the UFCSimplex method.""" - return self.product.distance_to_point_l1(point) + C and M are decomposed into "vertical" and "horizontal" + parts, keeping the relative orders of numbers: - def symmetry_group_size(self, dim): - return [1, 2, 8, 48][dim] + C -> C0 = [55, 43], C1 = [42, 57] + M -> M0 = [43, 55], M1 = [42, 57] - def cell_orientation_reflection_map(self): - """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" - return self.product.cell_orientation_reflection_map() + Then the orientation of the cell is computed as the + following: - def __gt__(self, other): - return self.product > other + C0.index(M0[0]) = 1; C0.remove(M0[0]) + C0.index(M0[1]) = 0; C0.remove(M0[1]) + C1.index(M1[0]) = 0; C1.remove(M1[0]) + C1.index(M1[1]) = 0; C1.remove(M1[1]) - def __lt__(self, other): - return self.product < other + o = 2 * 1 + 0 = 2 + """ + dimension = 2 + shape = QUADRILATERAL - def __ge__(self, other): - return self.product >= other - def __le__(self, other): - return self.product <= other +class UFCHexahedron(UFCHypercube): + """This is the reference hexahedron with vertices + (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), + (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" + dimension = 3 + shape = HEXAHEDRON def make_affine_mapping(xs, ys): @@ -1407,6 +1323,21 @@ def make_affine_mapping(xs, ys): return A, b +def ufc_hypercube(spatial_dim): + """Factory function that maps spatial dimension to an instance of + the UFC reference hypercube of that dimension.""" + if spatial_dim == 0: + return Point() + elif spatial_dim == 1: + return UFCInterval() + elif spatial_dim == 2: + return UFCQuadrilateral() + elif spatial_dim == 3: + return UFCHexahedron() + else: + raise RuntimeError("Can't create UFC hypercube of dimension %s." % str(spatial_dim)) + + def default_simplex(spatial_dim): """Factory function that maps spatial dimension to an instance of the default reference simplex of that dimension.""" From 8839f876732b1614e43b79fc7d7eb4b06987e391 Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Thu, 23 Jan 2025 22:47:18 +0000 Subject: [PATCH 02/14] Extend PolySet coefficients (#130) * Add extension of coefficients to allow union of polysets with different degrees * remove orthonormal as seems to be no longer necessary * Remove further orthonormals --- FIAT/nedelec.py | 4 +-- FIAT/polynomial_set.py | 40 +++++++++++++++++++++++++--- FIAT/raviart_thomas.py | 2 +- test/FIAT/unit/test_polynomial.py | 43 +++++++++++++++++++++++++++++++ 4 files changed, 83 insertions(+), 6 deletions(-) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c92fa06be..6cf122772 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -43,8 +43,8 @@ def NedelecSpace2D(ref_el, degree): CrossX = numpy.dot(numpy.array([[0.0, 1.0], [-1.0, 0.0]]), Qpts.T) PkHCrossX_at_Qpts = PkH_at_Qpts[:, None, :] * CrossX[None, :, :] - PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), Pkp1_at_Qpts.T) - + PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), + Pkp1_at_Qpts.T) PkHcrossX = polynomial_set.PolynomialSet(ref_el, k + 1, k + 1, diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 08e62b8d8..9686cf245 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -176,15 +176,49 @@ def polynomial_set_union_normalized(A, B): not contain any of the same members of the set, as we construct a span via SVD. """ - new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + assert A.get_reference_element() == B.get_reference_element() + new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B) + + deg = max(A.get_degree(), B.get_degree()) + em_deg = max(A.get_embedded_degree(), B.get_embedded_degree()) coeffs = spanning_basis(new_coeffs) return PolynomialSet(A.get_reference_element(), - A.get_degree(), - A.get_embedded_degree(), + deg, + em_deg, A.get_expansion_set(), coeffs) +def construct_new_coeffs(ref_el, A, B): + """ + Constructs new coefficients for the set union of A and B + If A and B are discontinuous and do not have the same degree the smaller one + is extended to match the larger. + + This does not handle the case that A and B have continuity but not the same degree. + """ + + if A.get_expansion_set().continuity != B.get_expansion_set().continuity: + raise ValueError("Continuity of expansion sets does not match.") + + if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None: + higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B + lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A + + diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] + + # pad only the 0th axis with the difference in size + padding = [(0, 0) for i in range(len(lower.coeffs.shape) - 1)] + [(0, diff)] + embedded_coeffs = numpy.pad(lower.coeffs, padding) + + new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0) + elif A.get_embedded_degree() == B.get_embedded_degree(): + new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + else: + raise NotImplementedError("Extending of coefficients is not implemented for PolynomialSets with continuity and different degrees") + return new_coeffs + + class ONSymTensorPolynomialSet(PolynomialSet): """Constructs an orthonormal basis for symmetric-tensor-valued polynomials on a reference element. diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index db0058d12..2c7d72c83 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -39,12 +39,12 @@ def RTSpace(ref_el, degree): # have to work on this through "tabulate" interface # first, tabulate PkH at quadrature points PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] x = Qpts.T PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) - PkHx = polynomial_set.PolynomialSet(ref_el, k, k + 1, diff --git a/test/FIAT/unit/test_polynomial.py b/test/FIAT/unit/test_polynomial.py index 2982cf11e..c211fb815 100644 --- a/test/FIAT/unit/test_polynomial.py +++ b/test/FIAT/unit/test_polynomial.py @@ -21,6 +21,7 @@ from FIAT import expansions, polynomial_set, reference_element from FIAT.quadrature_schemes import create_quadrature +from itertools import chain @pytest.fixture(params=(1, 2, 3)) @@ -107,3 +108,45 @@ def test_bubble_duality(cell, degree): results = scale * numpy.dot(numpy.multiply(phi_dual, qwts), phi.T) assert numpy.allclose(results, numpy.diag(numpy.diag(results))) assert numpy.allclose(numpy.diag(results), 1.0) + + +@pytest.mark.parametrize("degree", [10]) +def test_union_of_polysets(cell, degree): + """ demonstrates that polysets don't need to have the same degree for union + using RT space as an example""" + + sd = cell.get_spatial_dimension() + k = degree + vecPk = polynomial_set.ONPolynomialSet(cell, degree, (sd,)) + + vec_Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, (sd,), scale="orthonormal") + + dimPkp1 = expansions.polynomial_dimension(cell, k + 1) + dimPk = expansions.polynomial_dimension(cell, k) + dimPkm1 = expansions.polynomial_dimension(cell, k - 1) + + vec_Pk_indices = list(chain(*(range(i * dimPkp1, i * dimPkp1 + dimPk) + for i in range(sd)))) + vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) + + Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, scale="orthonormal") + PkH = Pkp1.take(list(range(dimPkm1, dimPk))) + + Q = create_quadrature(cell, 2 * (k + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() + + PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] + x = Qpts.T + PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] + PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) + PkHx = polynomial_set.PolynomialSet(cell, k, k + 1, vec_Pkp1.get_expansion_set(), PkHx_coeffs) + + same_deg = polynomial_set.polynomial_set_union_normalized(vec_Pk_from_Pkp1, PkHx) + different_deg = polynomial_set.polynomial_set_union_normalized(vecPk, PkHx) + + Q = create_quadrature(cell, 2*(degree)) + Qpts, _ = Q.get_points(), Q.get_weights() + same_vals = same_deg.tabulate(Qpts)[(0,) * sd] + diff_vals = different_deg.tabulate(Qpts)[(0,) * sd] + assert numpy.allclose(same_vals - diff_vals, 0) From 463c02d5c487491f938fb81f622fa8e51aa00a28 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Wed, 12 Feb 2025 13:53:21 +0000 Subject: [PATCH 03/14] Fix quadrature rule hash (#132) * Make hashing of quadrature rules safe * Also remove some uses of deprecated `@abstractproperty` and gives cells and point sets `repr()`s. * Add safe_repr function to handle floating point types --- FIAT/reference_element.py | 9 ++++- finat/point_set.py | 35 +++++++++++++--- finat/quadrature.py | 36 ++++++++++++++--- gem/utils.py | 55 ++++++++++++++++++++++++++ test/finat/test_create_fiat_element.py | 5 +++ test/finat/test_quadrature.py | 19 +++++++++ 6 files changed, 147 insertions(+), 12 deletions(-) create mode 100644 test/finat/test_quadrature.py diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 30fa94bb8..2cad80bf0 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -26,6 +26,7 @@ from math import factorial import numpy +from gem.utils import safe_repr from recursivenodes.nodes import _decode_family, _recursive from FIAT.orientation_utils import ( @@ -126,7 +127,7 @@ def linalg_subspace_intersection(A, B): return U[:, :rank_c] -class Cell(object): +class Cell: """Abstract class for a reference cell. Provides accessors for geometry (vertex coordinates) as well as topology (orderings of vertices that make up edges, faces, etc.""" @@ -184,6 +185,9 @@ def __init__(self, shape, vertices, topology): # Dictionary with derived cells self._split_cache = {} + def __repr__(self): + return f"{type(self).__name__}({self.shape!r}, {safe_repr(self.vertices)}, {self.topology!r})" + def _key(self): """Hashable object key data (excluding type).""" # Default: only type matters @@ -1130,6 +1134,9 @@ def __init__(self, *cells): super().__init__(TENSORPRODUCT, vertices, topology) self.cells = tuple(cells) + def __repr__(self): + return f"{type(self).__name__}({self.cells!r})" + def _key(self): return self.cells diff --git a/finat/point_set.py b/finat/point_set.py index 1497308c7..2a437c9c1 100644 --- a/finat/point_set.py +++ b/finat/point_set.py @@ -1,13 +1,15 @@ -from abc import ABCMeta, abstractproperty +import abc +import hashlib +from functools import cached_property from itertools import chain, product import numpy import gem -from gem.utils import cached_property +from gem.utils import safe_repr -class AbstractPointSet(metaclass=ABCMeta): +class AbstractPointSet(abc.ABC): """A way of specifying a known set of points, perhaps with some (tensor) structure. @@ -15,8 +17,15 @@ class AbstractPointSet(metaclass=ABCMeta): where point_set_shape is () for scalar, (N,) for N element vector, (N, M) for N x M matrix etc. """ + def __hash__(self): + return int.from_bytes(hashlib.md5(repr(self).encode()).digest(), byteorder="big") - @abstractproperty + @abc.abstractmethod + def __repr__(self): + pass + + @property + @abc.abstractmethod def points(self): """A flattened numpy array of points or ``UnknownPointsArray`` object with shape (# of points, point dimension).""" @@ -27,12 +36,14 @@ def dimension(self): _, dim = self.points.shape return dim - @abstractproperty + @property + @abc.abstractmethod def indices(self): """GEM indices with matching shape and extent to the structure of the point set.""" - @abstractproperty + @property + @abc.abstractmethod def expression(self): """GEM expression describing the points, with free indices ``self.indices`` and shape (point dimension,).""" @@ -53,6 +64,9 @@ def __init__(self, point): assert len(point.shape) == 1 self.point = point + def __repr__(self): + return f"{type(self).__name__}({safe_repr(self.point)})" + @cached_property def points(self): # Make sure we conform to the expected (# of points, point dimension) @@ -106,6 +120,9 @@ def __init__(self, points_expr): assert len(points_expr.shape) == 2 self._points_expr = points_expr + def __repr__(self): + return f"{type(self).__name__}({self._points_expr!r})" + @cached_property def points(self): return UnknownPointsArray(self._points_expr.shape) @@ -133,6 +150,9 @@ def __init__(self, points): assert len(points.shape) == 2 self.points = points + def __repr__(self): + return f"{type(self).__name__}({self.points!r})" + @cached_property def points(self): pass # set at initialisation @@ -177,6 +197,9 @@ class TensorPointSet(AbstractPointSet): def __init__(self, factors): self.factors = tuple(factors) + def __repr__(self): + return f"{type(self).__name__}({self.factors!r})" + @cached_property def points(self): return numpy.array([list(chain(*pt_tuple)) diff --git a/finat/quadrature.py b/finat/quadrature.py index ec6def127..10e4048b9 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -1,12 +1,13 @@ -from abc import ABCMeta, abstractproperty -from functools import reduce +import hashlib +from abc import ABCMeta, abstractmethod +from functools import cached_property, reduce import gem import numpy from FIAT.quadrature import GaussLegendreQuadratureLineRule from FIAT.quadrature_schemes import create_quadrature as fiat_scheme from FIAT.reference_element import LINE, QUADRILATERAL, TENSORPRODUCT -from gem.utils import cached_property +from gem.utils import safe_repr from finat.point_set import GaussLegendrePointSet, PointSet, TensorPointSet @@ -60,11 +61,23 @@ class AbstractQuadratureRule(metaclass=ABCMeta): """Abstract class representing a quadrature rule as point set and a corresponding set of weights.""" - @abstractproperty + def __hash__(self): + return int.from_bytes(hashlib.md5(repr(self).encode()).digest(), byteorder="big") + + def __eq__(self, other): + return type(other) is type(self) and repr(other) == repr(self) + + @abstractmethod + def __repr__(self): + pass + + @property + @abstractmethod def point_set(self): """Point set object representing the quadrature points.""" - @abstractproperty + @property + @abstractmethod def weight_expression(self): """GEM expression describing the weights, with the same free indices as the point set.""" @@ -110,6 +123,16 @@ def __init__(self, point_set, weights, ref_el=None, io_ornt_map_tuple=(None, )): self.weights = numpy.asarray(weights) self._intrinsic_orientation_permutation_map_tuple = io_ornt_map_tuple + def __repr__(self): + return ( + f"{type(self).__name__}(" + f"{self.point_set!r}, " + f"{safe_repr(self.weights)}, " + f"{self.ref_el!r}, " + f"{self._intrinsic_orientation_permutation_map_tuple!r}" + ")" + ) + @cached_property def point_set(self): pass # set at initialisation @@ -131,6 +154,9 @@ def __init__(self, factors, ref_el=None): for m in factor._intrinsic_orientation_permutation_map_tuple ) + def __repr__(self): + return f"{type(self).__name__}({self.factors!r}, {self.ref_el!r})" + @cached_property def point_set(self): return TensorPointSet(q.point_set for q in self.factors) diff --git a/gem/utils.py b/gem/utils.py index 7e855f768..a6bf97dfa 100644 --- a/gem/utils.py +++ b/gem/utils.py @@ -1,5 +1,10 @@ import collections +import functools +import numbers from functools import cached_property # noqa: F401 +from typing import Any + +import numpy as np def groupby(iterable, key=None): @@ -88,3 +93,53 @@ def __exit__(self, exc_type, exc_value, traceback): assert self.state is variable._head value, variable._head = variable._head self.state = None + + +@functools.singledispatch +def safe_repr(obj: Any) -> str: + """Return a 'safe' repr for an object, accounting for floating point error. + + Parameters + ---------- + obj : + The object to produce a repr for. + + Returns + ------- + str : + A repr for the object. + + """ + raise TypeError(f"Cannot provide a safe repr for {type(obj).__name__}") + + +@safe_repr.register(str) +def _(text: str) -> str: + return text + + +@safe_repr.register(numbers.Integral) +def _(num: numbers.Integral) -> str: + return repr(num) + + +@safe_repr.register(numbers.Real) +def _(num: numbers.Real) -> str: + # set roundoff to close-to-but-not-exactly machine epsilon + precision = np.finfo(num).precision - 2 + return "{:.{prec}}".format(num, prec=precision) + + +@safe_repr.register(np.ndarray) +def _(array: np.ndarray) -> str: + return f"{type(array).__name__}([{', '.join(map(safe_repr, array))}])" + + +@safe_repr.register(list) +def _(list_: list) -> str: + return f"[{', '.join(map(safe_repr, list_))}]" + + +@safe_repr.register(tuple) +def _(tuple_: tuple) -> str: + return f"({', '.join(map(safe_repr, tuple_))})" diff --git a/test/finat/test_create_fiat_element.py b/test/finat/test_create_fiat_element.py index f99053ab9..6530c8243 100644 --- a/test/finat/test_create_fiat_element.py +++ b/test/finat/test_create_fiat_element.py @@ -59,6 +59,11 @@ def tensor_name(request): ids=lambda x: x.cellname(), scope="module") def ufl_A(request, tensor_name): + if request.param == ufl.quadrilateral: + if tensor_name == "DG": + tensor_name = "DQ" + elif tensor_name == "DG L2": + tensor_name = "DQ L2" return finat.ufl.FiniteElement(tensor_name, request.param, 1) diff --git a/test/finat/test_quadrature.py b/test/finat/test_quadrature.py new file mode 100644 index 000000000..4d1acb9f7 --- /dev/null +++ b/test/finat/test_quadrature.py @@ -0,0 +1,19 @@ +import pytest + +from FIAT import ufc_cell +from finat.quadrature import make_quadrature + + +@pytest.mark.parametrize( + "cell_name", + ["interval", "triangle", "interval * interval", "triangle * interval"] +) +def test_quadrature_rules_are_hashable(cell_name): + ref_cell = ufc_cell(cell_name) + quadrature1 = make_quadrature(ref_cell, 3) + quadrature2 = make_quadrature(ref_cell, 3) + + assert quadrature1 is not quadrature2 + assert hash(quadrature1) == hash(quadrature2) + assert repr(quadrature1) == repr(quadrature2) + assert quadrature1 == quadrature2 From 637894a00d0257e00f9ce3254e4dfb4ec3e7aa42 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 13 Feb 2025 17:05:01 +0000 Subject: [PATCH 04/14] Fine-grained macro constraints (#133) * More general macro Constraints * Fine grained CkPolynomialSet * Test C2 on the double Alfeld split --- FIAT/dual_set.py | 3 ++ FIAT/expansions.py | 34 ++++++++++++++------- FIAT/finite_element.py | 3 ++ FIAT/macro.py | 59 +++++++++++++++++++++++------------- FIAT/polynomial_set.py | 3 ++ test/FIAT/unit/test_macro.py | 21 +++++++++++++ 6 files changed, 91 insertions(+), 32 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 74630363d..c09db974f 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -35,6 +35,9 @@ def __init__(self, nodes, ref_el, entity_ids, entity_permutations=None): def __iter__(self): return iter(self.nodes) + def __len__(self): + return len(self.nodes) + def get_nodes(self): return self.nodes diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 28eb561d4..cf0a67a3b 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -436,26 +436,38 @@ def tabulate_jumps(self, n, points, order=0): """ from FIAT.polynomial_set import mis + sd = self.ref_el.get_spatial_dimension() num_members = self.get_num_members(n) cell_node_map = self.get_cell_node_map(n) + cell_point_map = compute_cell_point_map(self.ref_el, points, unique=False) + + num_jumps = 0 + facet_point_map = {} + for facet in self.ref_el.get_interior_facets(sd-1): + try: + cells = self.ref_el.connectivity[(sd-1, sd)][facet] + ipts = list(set.intersection(*(set(cell_point_map[c]) for c in cells))) + if ipts != (): + facet_point_map[facet] = ipts + num_jumps += len(ipts) + except KeyError: + pass - top = self.ref_el.get_topology() - sd = self.ref_el.get_spatial_dimension() - interior_facets = self.ref_el.get_interior_facets(sd-1) derivs = {cell: self._tabulate_on_cell(n, points, order=order, cell=cell) - for cell in top[sd]} + for cell in cell_point_map} + jumps = {} for r in range(order+1): cur = 0 alphas = mis(sd, r) - jumps[r] = numpy.zeros((num_members, len(alphas)*len(interior_facets), len(points))) - for facet in interior_facets: - facet_verts = set(top[sd-1][facet]) - c0, c1 = tuple(k for k in top[sd] if facet_verts < set(top[sd][k])) + jumps[r] = numpy.zeros((num_members, len(alphas) * num_jumps)) + for facet, ipts in facet_point_map.items(): + c0, c1 = self.ref_el.connectivity[(sd-1, sd)][facet] for alpha in alphas: - jumps[r][cell_node_map[c1], cur] += derivs[c1][alpha] - jumps[r][cell_node_map[c0], cur] -= derivs[c0][alpha] - cur = cur + 1 + ijump = range(cur, cur + len(ipts)) + jumps[r][numpy.ix_(cell_node_map[c1], ijump)] += derivs[c1][alpha][:, ipts] + jumps[r][numpy.ix_(cell_node_map[c0], ijump)] -= derivs[c0][alpha][:, ipts] + cur += len(ipts) return jumps def get_dmats(self, degree, cell=0): diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 6f4ff24d5..96cf26286 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -134,6 +134,9 @@ def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref ref_complex = ref_complex or poly_set.get_reference_element() super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex) + if len(poly_set) != len(dual): + raise ValueError(f"Dimension of function space is {len(poly_set)}, but got {len(dual)} nodes.") + # build generalized Vandermonde matrix old_coeffs = poly_set.get_coeffs() dualmat = dual.to_riesz(poly_set) diff --git a/FIAT/macro.py b/FIAT/macro.py index 6662bdfcb..07a809acb 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -412,21 +412,35 @@ class CkPolynomialSet(polynomial_set.PolynomialSet): :arg ref_el: The simplicial complex. :arg degree: The polynomial degree. - :kwarg order: The order of continuity across subcells. + :kwarg order: The order of continuity across facets or a dict of dicts + mapping dimension and entity_id to the order of continuity at each entity. :kwarg vorder: The order of super-smoothness at interior vertices. :kwarg shape: The value shape. :kwarg variant: The variant for the underlying ExpansionSet. :kwarg scale: The scale for the underlying ExpansionSet. """ - def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): + def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): from FIAT.quadrature_schemes import create_quadrature + + if not isinstance(order, (int, dict)): + raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}") + + sd = ref_el.get_spatial_dimension() + if isinstance(order, int): + order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)} + if vorder is not None: + order[0] = dict.fromkeys(ref_el.get_interior_facets(0), vorder) + elif 0 not in order: + order[0] = {} + + if not all(k in {0, sd-1} for k in order): + raise NotImplementedError("Only face or vertex constraints have been implemented.") + expansion_set = expansions.ExpansionSet(ref_el, **kwargs) k = 1 if expansion_set.continuity == "C0" else 0 - num_members = expansion_set.get_num_members(degree) - sd = ref_el.get_spatial_dimension() + # Impose C^forder continuity across interior facets facet_el = ref_el.construct_subelement(sd-1) - phi_deg = 0 if sd == 1 else degree - k phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg) Q = create_quadrature(facet_el, 2 * phi_deg) @@ -435,30 +449,33 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): weights = numpy.multiply(phi_at_qpts, qwts) rows = [] - for facet in ref_el.get_interior_facets(sd-1): - jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) - for r in range(k, order+1): + for facet in order[sd-1]: + forder = order[sd-1][facet] + jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=forder) + for r in range(k, forder+1): num_wt = 1 if sd == 1 else expansions.polynomial_dimension(facet_el, degree-r) - rows.append(numpy.tensordot(weights[:num_wt], jumps[r], axes=(-1, -1)).reshape(-1, num_members)) + rows.append(numpy.tensordot(weights[:num_wt], jumps[r], axes=(-1, -1))) + + # Impose C^vorder super-smoothness at interior vertices + # C^forder automatically gives C^{forder+dim-1} at the interior vertex + verts = numpy.asarray(ref_el.get_vertices()) + for vorder in set(order[0].values()): + vids = [i for i in order[0] if order[0][i] == vorder] + facets = chain.from_iterable(ref_el.connectivity[(0, sd-1)][v] for v in vids) + forder = min(order[sd-1][f] for f in facets) + sorder = forder + sd - 1 + if vorder > sorder: + jumps = expansion_set.tabulate_jumps(degree, verts[vids], order=vorder) + rows.extend(numpy.vstack(jumps[r].T) for r in range(sorder+1, vorder+1)) if len(rows) > 0: + for row in rows: + row *= 1 / max(numpy.max(abs(row)), 1) dual_mat = numpy.vstack(rows) coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True) else: coeffs = numpy.eye(expansion_set.get_num_members(degree)) - # Impose C^vorder super-smoothness at interior vertices - # C^order automatically gives C^{order+dim-1} at the interior vertex - sorder = order + sd - 1 - if vorder > sorder: - verts = ref_el.get_vertices() - points = [verts[i] for i in ref_el.get_interior_facets(0)] - jumps = expansion_set.tabulate_jumps(degree, points, order=vorder) - for r in range(sorder+1, vorder+1): - dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T) - nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.dot(nsp, coeffs) - if shape != (): m, n = coeffs.shape ncomp = numpy.prod(shape) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 9686cf245..534c18642 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -104,6 +104,9 @@ def take(self, items): return PolynomialSet(self.ref_el, self.degree, self.embedded_degree, self.expansion_set, new_coeffs) + def __len__(self): + return self.num_members + class ONPolynomialSet(PolynomialSet): """Constructs an orthonormal basis out of expansion set by having an diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py index 1677d7779..8e367bf20 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -367,6 +367,27 @@ def test_Ck_basis(cell, order, degree, variant): assert numpy.allclose(local_phis, phis[:, ipts]) +def test_C2_double_alfeld(): + # Construct the quintic C2 spline on the double Alfeld split + # See Section 7.5 of Lai & Schumacher + K = ufc_simplex(2) + DCT = AlfeldSplit(AlfeldSplit(K)) + + degree = 5 + + # C3 on major split facets, C2 elsewhere + order = {} + order[1] = dict.fromkeys(DCT.get_interior_facets(1), 2) + order[1].update(dict.fromkeys(range(3, 6), 3)) + + # C4 at minor split barycenters, C3 at major split barycenter + order[0] = dict.fromkeys(DCT.get_interior_facets(0), 4) + order[0][3] = 3 + + P = CkPolynomialSet(DCT, degree, order=order, variant="bubble") + assert P.get_num_members() == 27 + + def test_distance_to_point_l1(cell): A = AlfeldSplit(cell) dim = A.get_spatial_dimension() From b90b09bac4fb63bda9743d89d7758fe74d018c70 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 14 Mar 2025 16:38:17 +0000 Subject: [PATCH 05/14] refactor to generalise beyond ufc --- FIAT/reference_element.py | 41 ++++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 64ad3c136..a1c8aaa60 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -44,6 +44,8 @@ HEXAHEDRON = 111 TENSORPRODUCT = 99 +hypercube_shapes = {2: QUADRILATERAL, 3: HEXAHEDRON} + def multiindex_equal(d, isum, imin=0): """A generator for d-tuple multi-indices whose sum is isum and minimum is imin. @@ -1352,21 +1354,26 @@ def extrinsic_orientation_permutation_map(self): return a -class UFCHypercube(Cell): - """Abstract class for a reference unit hypercube [0, 1]^d with vertices in - lexicographical order.""" - dimension = None - shape = None +class Hypercube(Cell): + """Abstract class for a reference hypercube + + if no tensor product is provided, it defaults to the + UFCHypercube: [0, 1]^d with vertices in + lexicographical order. """ - def __init__(self): - cells = [UFCInterval()] * self.dimension - product = TensorProductCell(*cells) - pt = product.get_topology() + def __init__(self, dimension, product=None): + self.dimension = dimension + self.shape = hypercube_shapes[dimension] + if product is None: + cells = [UFCInterval()] * self.dimension + product = TensorProductCell(*cells) + + pt = product.get_topology() verts = product.get_vertices() topology = flatten_entities(pt) - super(UFCHypercube, self).__init__(self.shape, verts, topology) + super(Hypercube, self).__init__(self.shape, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) @@ -1457,7 +1464,7 @@ def __le__(self, other): return self.product <= other -class UFCQuadrilateral(UFCHypercube): +class UFCQuadrilateral(Hypercube): r"""This is the reference quadrilateral with vertices (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). @@ -1502,16 +1509,18 @@ class UFCQuadrilateral(UFCHypercube): o = 2 * 1 + 0 = 2 """ - dimension = 2 - shape = QUADRILATERAL + + def __init__(self): + super(UFCQuadrilateral, self).__init__(2) -class UFCHexahedron(UFCHypercube): +class UFCHexahedron(Hypercube): """This is the reference hexahedron with vertices (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" - dimension = 3 - shape = HEXAHEDRON + + def __init__(self): + super(UFCHexahedron, self).__init__(3) def make_affine_mapping(xs, ys): From b75e2cc15e661027f79b7adeca8c6b30f0019c48 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 17 Mar 2025 09:09:37 +0000 Subject: [PATCH 06/14] Lint --- FIAT/reference_element.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index a1c8aaa60..95bfdb98a 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1356,7 +1356,7 @@ def extrinsic_orientation_permutation_map(self): class Hypercube(Cell): """Abstract class for a reference hypercube - + if no tensor product is provided, it defaults to the UFCHypercube: [0, 1]^d with vertices in lexicographical order. """ @@ -1368,7 +1368,7 @@ def __init__(self, dimension, product=None): if product is None: cells = [UFCInterval()] * self.dimension product = TensorProductCell(*cells) - + pt = product.get_topology() verts = product.get_vertices() topology = flatten_entities(pt) @@ -1518,7 +1518,7 @@ class UFCHexahedron(Hypercube): """This is the reference hexahedron with vertices (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" - + def __init__(self): super(UFCHexahedron, self).__init__(3) From 92be38f2e15aa4cfce0f1166eb505648361e5d25 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 17 Mar 2025 11:04:51 +0000 Subject: [PATCH 07/14] Modify symgroup size to be general formula --- FIAT/reference_element.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 95bfdb98a..0af7bd873 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1445,7 +1445,8 @@ def distance_to_point_l1(self, point, rescale=False): return self.product.distance_to_point_l1(point, rescale=rescale) def symmetry_group_size(self, dim): - return [1, 2, 8, 48][dim] + """ Size of hypercube symmetry group is d! * 2**d""" + return factorial(dim) * (2**dim) def cell_orientation_reflection_map(self): """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" From c1f1afd7f2d353dcbed7344b85b6edc9f8248ebf Mon Sep 17 00:00:00 2001 From: India Marsden Date: Mon, 17 Mar 2025 15:56:00 +0000 Subject: [PATCH 08/14] tidy comment and work around some places where UFC is assumed --- FIAT/reference_element.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 0af7bd873..2a434417d 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1364,8 +1364,10 @@ class Hypercube(Cell): def __init__(self, dimension, product=None): self.dimension = dimension self.shape = hypercube_shapes[dimension] + self.ufc = False if product is None: + self.ufc = True cells = [UFCInterval()] * self.dimension product = TensorProductCell(*cells) @@ -1394,8 +1396,11 @@ def construct_subelement(self, dimension): raise ValueError("Invalid dimension: %d" % (dimension,)) elif dimension == sd: return self - else: + elif self.ufc: return ufc_hypercube(dimension) + else: + sub_element = self.product.construct_subelement((dimension,) + (0,)*(len(self.product.cells) - 1)) + return flatten_reference_cube(sub_element) def get_entity_transform(self, dim, entity_i): """Returns a mapping of point coordinates from the @@ -1445,7 +1450,7 @@ def distance_to_point_l1(self, point, rescale=False): return self.product.distance_to_point_l1(point, rescale=rescale) def symmetry_group_size(self, dim): - """ Size of hypercube symmetry group is d! * 2**d""" + """Size of hypercube symmetry group is d! * 2**d""" return factorial(dim) * (2**dim) def cell_orientation_reflection_map(self): @@ -1682,7 +1687,7 @@ def tuple_sum(tree): def is_hypercube(cell): - if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)): + if isinstance(cell, (DefaultLine, UFCInterval, Hypercube)): return True elif isinstance(cell, TensorProductCell): return reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) @@ -1692,15 +1697,16 @@ def is_hypercube(cell): def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()} if numpy.sum(ref_el.get_dimension()) <= 1: # Just return point/interval cell arguments return ref_el else: # Handle cases where cell is a quad/cube constructed from a tensor product or # an already flattened element - if is_hypercube(ref_el): - return flattened_cube[numpy.sum(ref_el.get_dimension())] + if isinstance(ref_el, TensorProductCell): + return Hypercube(numpy.sum(ref_el.get_dimension()), ref_el) + elif is_hypercube(ref_el): + return ref_el else: raise TypeError('Invalid cell type') From 18e31ed03abe41259f75ac8517cc04022a75d91c Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Mon, 17 Mar 2025 16:12:57 +0000 Subject: [PATCH 09/14] Update FIAT/reference_element.py Co-authored-by: Pablo Brubeck --- FIAT/reference_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 2a434417d..3fa829a67 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1704,7 +1704,7 @@ def flatten_reference_cube(ref_el): # Handle cases where cell is a quad/cube constructed from a tensor product or # an already flattened element if isinstance(ref_el, TensorProductCell): - return Hypercube(numpy.sum(ref_el.get_dimension()), ref_el) + return Hypercube(ref_el.get_spatial_dimension(), ref_el) elif is_hypercube(ref_el): return ref_el else: From 9a0e0fd5360e46e10aefb97c8bc848db111f085b Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 18 Mar 2025 10:09:19 +0000 Subject: [PATCH 10/14] Refactor to make UFC hyper a separate class --- FIAT/reference_element.py | 48 +++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 3fa829a67..8c7cedc3d 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1355,21 +1355,11 @@ def extrinsic_orientation_permutation_map(self): class Hypercube(Cell): - """Abstract class for a reference hypercube + """Abstract class for a reference hypercube""" - if no tensor product is provided, it defaults to the - UFCHypercube: [0, 1]^d with vertices in - lexicographical order. """ - - def __init__(self, dimension, product=None): + def __init__(self, dimension, product): self.dimension = dimension self.shape = hypercube_shapes[dimension] - self.ufc = False - - if product is None: - self.ufc = True - cells = [UFCInterval()] * self.dimension - product = TensorProductCell(*cells) pt = product.get_topology() verts = product.get_vertices() @@ -1396,8 +1386,6 @@ def construct_subelement(self, dimension): raise ValueError("Invalid dimension: %d" % (dimension,)) elif dimension == sd: return self - elif self.ufc: - return ufc_hypercube(dimension) else: sub_element = self.product.construct_subelement((dimension,) + (0,)*(len(self.product.cells) - 1)) return flatten_reference_cube(sub_element) @@ -1470,7 +1458,33 @@ def __le__(self, other): return self.product <= other -class UFCQuadrilateral(Hypercube): +class UFCHypercube(Hypercube): + """Reference UFC Hypercube + + UFCHypercube: [0, 1]^d with vertices in + lexicographical order.""" + + def __init__(self, dim): + cells = [UFCInterval()] * dim + product = TensorProductCell(*cells) + super(UFCHypercube, self).__init__(dim, product) + + def construct_subelement(self, dimension): + """Constructs the reference element of a cell subentity + specified by subelement dimension. + + :arg dimension: subentity dimension (integer) + """ + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError("Invalid dimension: %d" % (dimension,)) + elif dimension == sd: + return self + else: + return ufc_hypercube(dimension) + + +class UFCQuadrilateral(UFCHypercube): r"""This is the reference quadrilateral with vertices (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). @@ -1520,7 +1534,7 @@ def __init__(self): super(UFCQuadrilateral, self).__init__(2) -class UFCHexahedron(Hypercube): +class UFCHexahedron(UFCHypercube): """This is the reference hexahedron with vertices (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" @@ -1697,7 +1711,7 @@ def is_hypercube(cell): def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - if numpy.sum(ref_el.get_dimension()) <= 1: + if ref_el.get_spatial_dimension() <= 1: # Just return point/interval cell arguments return ref_el else: From 94c93119cb17273115daffc84b574c0c47d87ab6 Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 18 Mar 2025 13:57:10 +0000 Subject: [PATCH 11/14] ensure if an element is ufc it is flattened to type ufc --- FIAT/reference_element.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 8c7cedc3d..a68e4d040 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1699,6 +1699,13 @@ def tuple_sum(tree): else: return tree +def is_ufc(cell): + if isinstance(cell, (Point, UFCInterval, UFCHypercube, UFCSimplex)): + return True + elif isinstance(cell, TensorProductCell): + return reduce(lambda a, b: a and b, [is_ufc(c) for c in cell.cells]) + else: + return False def is_hypercube(cell): if isinstance(cell, (DefaultLine, UFCInterval, Hypercube)): @@ -1718,6 +1725,8 @@ def flatten_reference_cube(ref_el): # Handle cases where cell is a quad/cube constructed from a tensor product or # an already flattened element if isinstance(ref_el, TensorProductCell): + if is_ufc(ref_el): + return ufc_hypercube(ref_el.get_spatial_dimension()) return Hypercube(ref_el.get_spatial_dimension(), ref_el) elif is_hypercube(ref_el): return ref_el From 8060e81ef850829aa0b921d65df2b3726aa79c4a Mon Sep 17 00:00:00 2001 From: India Marsden Date: Tue, 18 Mar 2025 16:01:26 +0000 Subject: [PATCH 12/14] lint and add tests --- FIAT/reference_element.py | 2 ++ test/FIAT/unit/test_reference_element.py | 41 ++++++++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index a68e4d040..831f6068c 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1699,6 +1699,7 @@ def tuple_sum(tree): else: return tree + def is_ufc(cell): if isinstance(cell, (Point, UFCInterval, UFCHypercube, UFCSimplex)): return True @@ -1707,6 +1708,7 @@ def is_ufc(cell): else: return False + def is_hypercube(cell): if isinstance(cell, (DefaultLine, UFCInterval, Hypercube)): return True diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index f551f4d69..d2153069a 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -21,6 +21,7 @@ from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron +from FIAT.reference_element import is_ufc, default_simplex, flatten_reference_cube, Hypercube point = Point() interval = UFCInterval() @@ -32,6 +33,12 @@ triangle_x_interval = TensorProductCell(triangle, interval) quadrilateral_x_interval = TensorProductCell(quadrilateral, interval) +default_interval = default_simplex(1) +default_triangle = default_simplex(2) +default_tetrahedron = default_simplex(3) +default_interval_x_interval = TensorProductCell(default_interval, default_interval) +default_hypercube = Hypercube(2, default_interval_x_interval) + ufc_tetrahedron_21connectivity = [(0, 1, 2), (0, 3, 4), (1, 3, 5), (2, 4, 5)] ufc_hexahedron_21connectivity = [(0, 1, 4, 5), (2, 3, 6, 7), (0, 2, 8, 9), (1, 3, 10, 11), (4, 6, 8, 10), (5, 7, 9, 11)] @@ -396,6 +403,40 @@ def test_distance_to_point_l1(cell, point, expected): assert np.isclose(cell.distance_to_point_l1(point), expected, rtol=1e-3) +@pytest.mark.parametrize(('cell', 'expected'), + [(interval, True), + (triangle, True), + (quadrilateral, True), + (tetrahedron, True), + (interval_x_interval, True), + (triangle_x_interval, True), + (quadrilateral_x_interval, True), + (hexahedron, True), + (default_interval, False), + (default_triangle, False), + (default_tetrahedron, False), + (default_interval_x_interval, False), + (default_hypercube, False),]) +def test_is_ufc(cell, expected): + assert is_ufc(cell) == expected + + +@pytest.mark.parametrize(('cell'), + [interval, + quadrilateral, + interval_x_interval, + triangle_x_interval, + quadrilateral_x_interval, + hexahedron, + default_interval, + default_interval_x_interval, + default_hypercube]) +def test_flatten_maintains_ufc_status(cell): + ufc_status = is_ufc(cell) + flat_cell = flatten_reference_cube(cell) + assert ufc_status == is_ufc(flat_cell) + + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__)) From 77ad0c0e747d4694896daba92ce71113d06549de Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Wed, 19 Mar 2025 16:33:42 +0000 Subject: [PATCH 13/14] General hypercube (#134) * cleanup, add UFCHypercube class * Modify symgroup size to be general formula * Refactor to make UFC hyper a separate class * ensure if an element is ufc it is flattened to type ufc --------- Co-authored-by: Pablo Brubeck Co-authored-by: Connor Ward --- FIAT/reference_element.py | 258 ++++++++++------------- test/FIAT/unit/test_reference_element.py | 59 ++++++ 2 files changed, 174 insertions(+), 143 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 2cad80bf0..547da2c04 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -44,6 +44,8 @@ HEXAHEDRON = 111 TENSORPRODUCT = 99 +hypercube_shapes = {2: QUADRILATERAL, 3: HEXAHEDRON} + def multiindex_equal(d, isum, imin=0): """A generator for d-tuple multi-indices whose sum is isum and minimum is imin. @@ -1218,7 +1220,7 @@ def compute_reference_normal(self, facet_dim, facet_i): def contains_point(self, point, epsilon=0.0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -1261,6 +1263,8 @@ def cell_orientation_reflection_map(self): return make_cell_orientation_reflection_map_tensorproduct(self.cells) def compare(self, op, other): + """Parent-based comparison between simplicial complexes. + This is done dimension by dimension.""" if hasattr(other, "product"): other = other.product if isinstance(other, type(self)): @@ -1350,59 +1354,18 @@ def extrinsic_orientation_permutation_map(self): return a -class UFCQuadrilateral(Cell): - r"""This is the reference quadrilateral with vertices - (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). - - Orientation of a physical cell is computed systematically - by comparing the canonical orderings of its facets and - the facets in the FIAT reference cell. - - As an example, we compute the orientation of a - quadrilateral cell: - - +---3---+ +--57---+ - | | | | - 0 1 43 55 - | | | | - +---2---+ +--42---+ - FIAT canonical Mapped example physical cell - - Suppose that the facets of the physical cell - are canonically ordered as: +class Hypercube(Cell): + """Abstract class for a reference hypercube""" - C = [55, 42, 43, 57] + def __init__(self, dimension, product): + self.dimension = dimension + self.shape = hypercube_shapes[dimension] - FIAT index to Physical index map must be such that - C[0] = 55 is mapped to a vertical facet; in this - example it is: - - M = [43, 55, 42, 57] - - C and M are decomposed into "vertical" and "horizontal" - parts, keeping the relative orders of numbers: - - C -> C0 = [55, 43], C1 = [42, 57] - M -> M0 = [43, 55], M1 = [42, 57] - - Then the orientation of the cell is computed as the - following: - - C0.index(M0[0]) = 1; C0.remove(M0[0]) - C0.index(M0[1]) = 0; C0.remove(M0[1]) - C1.index(M1[0]) = 0; C1.remove(M1[0]) - C1.index(M1[1]) = 0; C1.remove(M1[1]) - - o = 2 * 1 + 0 = 2 - """ - def __init__(self): - product = TensorProductCell(UFCInterval(), UFCInterval()) pt = product.get_topology() - verts = product.get_vertices() topology = flatten_entities(pt) - super().__init__(QUADRILATERAL, verts, topology) + super().__init__(self.shape, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) @@ -1418,14 +1381,14 @@ def construct_subelement(self, dimension): :arg dimension: subentity dimension (integer) """ - if dimension == 2: + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError(f"Invalid dimension: {(dimension,)}") + elif dimension == sd: return self - elif dimension == 1: - return UFCInterval() - elif dimension == 0: - return Point() else: - raise ValueError("Invalid dimension: %d" % (dimension,)) + sub_element = self.product.construct_subelement((dimension,) + (0,)*(len(self.product.cells) - 1)) + return flatten_reference_cube(sub_element) def get_entity_transform(self, dim, entity_i): """Returns a mapping of point coordinates from the @@ -1443,13 +1406,14 @@ def volume(self): def compute_reference_normal(self, facet_dim, facet_i): """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 1 + sd = self.get_spatial_dimension() + assert facet_dim == sd - 1 d, i = self.unflattening_map[(facet_dim, facet_i)] return self.product.compute_reference_normal(d, i) def contains_point(self, point, epsilon=0): """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', + (with numerical tolerance as given by the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to the cell. Parameters @@ -1467,14 +1431,15 @@ def contains_point(self, point, epsilon=0): return self.product.contains_point(point, epsilon=epsilon) def distance_to_point_l1(self, point, rescale=False): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + """Get the L1 distance (aka 'manhattan', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" return self.product.distance_to_point_l1(point, rescale=rescale) def symmetry_group_size(self, dim): - return [1, 2, 8][dim] + """Size of hypercube symmetry group is d! * 2**d""" + return factorial(dim) * (2**dim) def cell_orientation_reflection_map(self): """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" @@ -1493,27 +1458,16 @@ def __le__(self, other): return self.product <= other -class UFCHexahedron(Cell): - """This is the reference hexahedron with vertices - (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), - (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" - - def __init__(self): - product = TensorProductCell(UFCInterval(), UFCInterval(), UFCInterval()) - pt = product.get_topology() - - verts = product.get_vertices() - topology = flatten_entities(pt) +class UFCHypercube(Hypercube): + """Reference UFC Hypercube - super().__init__(HEXAHEDRON, verts, topology) + UFCHypercube: [0, 1]^d with vertices in + lexicographical order.""" - self.product = product - self.unflattening_map = compute_unflattening_map(pt) - - def get_dimension(self): - """Returns the subelement dimension of the cell. Same as the - spatial dimension.""" - return self.get_spatial_dimension() + def __init__(self, dim): + cells = [UFCInterval()] * dim + product = TensorProductCell(*cells) + super().__init__(dim, product) def construct_subelement(self, dimension): """Constructs the reference element of a cell subentity @@ -1521,81 +1475,72 @@ def construct_subelement(self, dimension): :arg dimension: subentity dimension (integer) """ - if dimension == 3: + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError(f"Invalid dimension: {dimension}") + elif dimension == sd: return self - elif dimension == 2: - return UFCQuadrilateral() - elif dimension == 1: - return UFCInterval() - elif dimension == 0: - return Point() else: - raise ValueError("Invalid dimension: %d" % (dimension,)) + return ufc_hypercube(dimension) - def get_entity_transform(self, dim, entity_i): - """Returns a mapping of point coordinates from the - `entity_i`-th subentity of dimension `dim` to the cell. - :arg dim: entity dimension (integer) - :arg entity_i: entity number (integer) - """ - d, e = self.unflattening_map[(dim, entity_i)] - return self.product.get_entity_transform(d, e) +class UFCQuadrilateral(UFCHypercube): + r"""This is the reference quadrilateral with vertices + (0.0, 0.0), (0.0, 1.0), (1.0, 0.0) and (1.0, 1.0). - def volume(self): - """Computes the volume in the appropriate dimensional measure.""" - return self.product.volume() + Orientation of a physical cell is computed systematically + by comparing the canonical orderings of its facets and + the facets in the FIAT reference cell. - def compute_reference_normal(self, facet_dim, facet_i): - """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 2 - d, i = self.unflattening_map[(facet_dim, facet_i)] - return self.product.compute_reference_normal(d, i) + As an example, we compute the orientation of a + quadrilateral cell: - def contains_point(self, point, epsilon=0): - """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', - 'taxicab' or rectilinear distance) to the cell. + +---3---+ +--57---+ + | | | | + 0 1 43 55 + | | | | + +---2---+ +--42---+ + FIAT canonical Mapped example physical cell - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. + Suppose that the facets of the physical cell + are canonically ordered as: - Returns - ------- - bool : True if the point is inside the cell, False otherwise. + C = [55, 42, 43, 57] - """ - return self.product.contains_point(point, epsilon=epsilon) + FIAT index to Physical index map must be such that + C[0] = 55 is mapped to a vertical facet; in this + example it is: - def distance_to_point_l1(self, point, rescale=False): - """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear - distance) to a point with 0.0 if the point is inside the cell. + M = [43, 55, 42, 57] - For more information see the docstring for the UFCSimplex method.""" - return self.product.distance_to_point_l1(point, rescale=rescale) + C and M are decomposed into "vertical" and "horizontal" + parts, keeping the relative orders of numbers: - def symmetry_group_size(self, dim): - return [1, 2, 8, 48][dim] + C -> C0 = [55, 43], C1 = [42, 57] + M -> M0 = [43, 55], M1 = [42, 57] - def cell_orientation_reflection_map(self): - """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" - return self.product.cell_orientation_reflection_map() + Then the orientation of the cell is computed as the + following: - def __gt__(self, other): - return self.product > other + C0.index(M0[0]) = 1; C0.remove(M0[0]) + C0.index(M0[1]) = 0; C0.remove(M0[1]) + C1.index(M1[0]) = 0; C1.remove(M1[0]) + C1.index(M1[1]) = 0; C1.remove(M1[1]) - def __lt__(self, other): - return self.product < other + o = 2 * 1 + 0 = 2 + """ - def __ge__(self, other): - return self.product >= other + def __init__(self): + super().__init__(2) - def __le__(self, other): - return self.product <= other + +class UFCHexahedron(UFCHypercube): + """This is the reference hexahedron with vertices + (0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), + (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).""" + + def __init__(self): + super().__init__(3) def make_affine_mapping(xs, ys): @@ -1634,6 +1579,21 @@ def make_affine_mapping(xs, ys): return A, b +def ufc_hypercube(spatial_dim): + """Factory function that maps spatial dimension to an instance of + the UFC reference hypercube of that dimension.""" + if spatial_dim == 0: + return Point() + elif spatial_dim == 1: + return UFCInterval() + elif spatial_dim == 2: + return UFCQuadrilateral() + elif spatial_dim == 3: + return UFCHexahedron() + else: + raise RuntimeError(f"Can't create UFC hypercube of dimension {spatial_dim}.") + + def default_simplex(spatial_dim): """Factory function that maps spatial dimension to an instance of the default reference simplex of that dimension.""" @@ -1646,7 +1606,7 @@ def default_simplex(spatial_dim): elif spatial_dim == 3: return DefaultTetrahedron() else: - raise RuntimeError("Can't create default simplex of dimension %s." % str(spatial_dim)) + raise RuntimeError(f"Can't create default simplex of dimension {spatial_dim}.") def ufc_simplex(spatial_dim): @@ -1661,7 +1621,7 @@ def ufc_simplex(spatial_dim): elif spatial_dim == 3: return UFCTetrahedron() else: - raise RuntimeError("Can't create UFC simplex of dimension %s." % str(spatial_dim)) + raise RuntimeError(f"Can't create UFC simplex of dimension {spatial_dim}.") def symmetric_simplex(spatial_dim): @@ -1701,7 +1661,7 @@ def ufc_cell(cell): elif celltype == "tetrahedron": return ufc_simplex(3) else: - raise RuntimeError("Don't know how to create UFC cell of type %s" % str(celltype)) + raise RuntimeError(f"Don't know how to create UFC cell of type {str(celltype)}") def volume(verts): @@ -1740,26 +1700,38 @@ def tuple_sum(tree): return tree +def is_ufc(cell): + if isinstance(cell, (Point, UFCInterval, UFCHypercube, UFCSimplex)): + return True + elif isinstance(cell, TensorProductCell): + return all(is_ufc(c) for c in cell.cells) + else: + return False + + def is_hypercube(cell): - if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)): + if isinstance(cell, (DefaultLine, UFCInterval, Hypercube)): return True elif isinstance(cell, TensorProductCell): - return reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) + return all(is_hypercube(c) for c in cell.cells) else: return False def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()} - if numpy.sum(ref_el.get_dimension()) <= 1: + if ref_el.get_spatial_dimension() <= 1: # Just return point/interval cell arguments return ref_el else: # Handle cases where cell is a quad/cube constructed from a tensor product or # an already flattened element - if is_hypercube(ref_el): - return flattened_cube[numpy.sum(ref_el.get_dimension())] + if isinstance(ref_el, TensorProductCell): + if is_ufc(ref_el): + return ufc_hypercube(ref_el.get_spatial_dimension()) + return Hypercube(ref_el.get_spatial_dimension(), ref_el) + elif is_hypercube(ref_el): + return ref_el else: raise TypeError('Invalid cell type') diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py index f551f4d69..d166ba6ef 100644 --- a/test/FIAT/unit/test_reference_element.py +++ b/test/FIAT/unit/test_reference_element.py @@ -21,6 +21,7 @@ from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron +from FIAT.reference_element import is_ufc, is_hypercube, default_simplex, flatten_reference_cube, Hypercube point = Point() interval = UFCInterval() @@ -32,6 +33,12 @@ triangle_x_interval = TensorProductCell(triangle, interval) quadrilateral_x_interval = TensorProductCell(quadrilateral, interval) +default_interval = default_simplex(1) +default_triangle = default_simplex(2) +default_tetrahedron = default_simplex(3) +default_interval_x_interval = TensorProductCell(default_interval, default_interval) +default_hypercube = Hypercube(2, default_interval_x_interval) + ufc_tetrahedron_21connectivity = [(0, 1, 2), (0, 3, 4), (1, 3, 5), (2, 4, 5)] ufc_hexahedron_21connectivity = [(0, 1, 4, 5), (2, 3, 6, 7), (0, 2, 8, 9), (1, 3, 10, 11), (4, 6, 8, 10), (5, 7, 9, 11)] @@ -396,6 +403,58 @@ def test_distance_to_point_l1(cell, point, expected): assert np.isclose(cell.distance_to_point_l1(point), expected, rtol=1e-3) +@pytest.mark.parametrize(('cell', 'expected'), + [(interval, True), + (triangle, True), + (quadrilateral, True), + (tetrahedron, True), + (interval_x_interval, True), + (triangle_x_interval, True), + (quadrilateral_x_interval, True), + (hexahedron, True), + (default_interval, False), + (default_triangle, False), + (default_tetrahedron, False), + (default_interval_x_interval, False), + (default_hypercube, False),]) +def test_is_ufc(cell, expected): + assert is_ufc(cell) == expected + + +@pytest.mark.parametrize(('cell', 'expected'), + [(interval, True), + (triangle, False), + (quadrilateral, True), + (tetrahedron, False), + (interval_x_interval, True), + (triangle_x_interval, False), + (quadrilateral_x_interval, True), + (hexahedron, True), + (default_interval, True), + (default_triangle, False), + (default_tetrahedron, False), + (default_interval_x_interval, True), + (default_hypercube, True),]) +def test_is_hypercube(cell, expected): + assert is_hypercube(cell) == expected + + +@pytest.mark.parametrize(('cell'), + [interval, + quadrilateral, + interval_x_interval, + triangle_x_interval, + quadrilateral_x_interval, + hexahedron, + default_interval, + default_interval_x_interval, + default_hypercube]) +def test_flatten_maintains_ufc_status(cell): + ufc_status = is_ufc(cell) + flat_cell = flatten_reference_cube(cell) + assert ufc_status == is_ufc(flat_cell) + + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__)) From 62206126e218af04c4fbc97bbca5d03b3adc59ab Mon Sep 17 00:00:00 2001 From: India Marsden Date: Fri, 21 Mar 2025 10:59:00 +0000 Subject: [PATCH 14/14] lint --- FIAT/discontinuous_pc.py | 10 ++-------- FIAT/reference_element.py | 13 ------------- 2 files changed, 2 insertions(+), 21 deletions(-) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index ba9894c91..8973a58f6 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -7,19 +7,13 @@ # Modified by David A. Ham (david.ham@imperial.ac.uk), 2018 from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.reference_element import (Point, - DefaultLine, - UFCInterval, - UFCQuadrilateral, - UFCHexahedron, - UFCTriangle, - UFCTetrahedron, - make_affine_mapping, +from FIAT.reference_element import (make_affine_mapping, flatten_reference_cube, cell_to_simplex) from FIAT.P0 import P0Dual import numpy as np + class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): flat_el = flatten_reference_cube(ref_el) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 028392807..edc902f0d 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1719,19 +1719,6 @@ def is_hypercube(cell): return False -def is_hypercube(cell): - if isinstance(cell, (DefaultLine, UFCInterval, Hypercube)): - return True - elif isinstance(cell, TensorProductCell): - res = reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) - else: - res = False - # breakpoint() - res2 = len(cell.vertices) == 2 ** cell.get_dimension() - # assert res == res2 - return res2 - - def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" if ref_el.get_spatial_dimension() <= 1: