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/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/FIAT/reference_element.py b/FIAT/reference_element.py index 606275656..edc902f0d 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 ( @@ -43,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. @@ -126,7 +129,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 +187,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 @@ -1131,6 +1137,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 @@ -1212,7 +1221,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 @@ -1255,6 +1264,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, TensorProductCell): @@ -1345,17 +1356,17 @@ def extrinsic_orientation_permutation_map(self): class Hypercube(Cell): - """ - For reference elements based on TensorProductCells""" - def __init__(self, product): - pt = product.get_topology() + """Abstract class for a reference hypercube""" + def __init__(self, dimension, product): + self.dimension = dimension + self.shape = hypercube_shapes[dimension] + + pt = product.get_topology() verts = product.get_vertices() topology = flatten_entities(pt) - # TODO this should be generalised ? - cube_types = {2: QUADRILATERAL, 3: HEXAHEDRON} - super().__init__(cube_types[sum(product.get_dimension())], verts, topology) + super().__init__(self.shape, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) @@ -1365,6 +1376,21 @@ def get_dimension(self): 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) + """ + sd = self.get_spatial_dimension() + if dimension > sd: + raise ValueError(f"Invalid dimension: {(dimension,)}") + elif dimension == sd: + return self + 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 `entity_i`-th subentity of dimension `dim` to the cell. @@ -1381,13 +1407,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 @@ -1405,14 +1432,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``).""" @@ -1431,7 +1459,33 @@ def __le__(self, other): return self.product <= other -class UFCQuadrilateral(Cell): +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().__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(f"Invalid dimension: {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). @@ -1476,207 +1530,18 @@ class UFCQuadrilateral(Cell): 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) - - 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 == 2: - return self - 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) - - def volume(self): - """Computes the volume in the appropriate dimensional measure.""" - return self.product.volume() - - def compute_reference_normal(self, facet_dim, facet_i): - """Returns the unit normal in infinity norm to facet_i.""" - assert facet_dim == 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', - 'taxicab' or rectilinear distance) to the cell. - - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. - - Returns - ------- - bool : True if the point is inside the cell, False otherwise. - - """ - 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 - 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] - - 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() - - def __gt__(self, other): - return self.product > other - - def __lt__(self, other): - return self.product < other - - def __ge__(self, other): - return self.product >= other - def __le__(self, other): - return self.product <= other + def __init__(self): + super().__init__(2) -class UFCHexahedron(Cell): +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): - product = TensorProductCell(UFCInterval(), UFCInterval(), UFCInterval()) - pt = product.get_topology() - - verts = product.get_vertices() - topology = flatten_entities(pt) - - super().__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) - - def volume(self): - """Computes the volume in the appropriate dimensional measure.""" - return self.product.volume() - - 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) - - 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. - - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. - - Returns - ------- - bool : True if the point is inside the cell, False otherwise. - - """ - 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 - 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, 48][dim] - - 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() - - def __gt__(self, other): - return self.product > other - - def __lt__(self, other): - return self.product < other - - def __ge__(self, other): - return self.product >= other - - def __le__(self, other): - return self.product <= other + super().__init__(3) def make_affine_mapping(xs, ys): @@ -1715,6 +1580,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.""" @@ -1727,7 +1607,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): @@ -1742,7 +1622,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): @@ -1782,7 +1662,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): @@ -1821,32 +1701,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): - res = False - if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)): - res = True + 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]) + return all(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 + return False def flatten_reference_cube(ref_el): """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" - from finat.element_factory import as_fiat_cell - flattened_cube = {2: as_fiat_cell("quadrilateral"), 3: as_fiat_cell("hexahedron")} - 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/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/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() 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__)) 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