diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 3f2fc20cd..85ed22916 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -507,11 +507,157 @@ def construct_subelement(self, dimension): def contains_point(self, point, epsilon=0): """Checks if reference cell contains given point - (with numerical tolerance).""" - result = (sum(point) - epsilon <= 1) - for c in point: - result &= (c + epsilon >= 0) - return result + (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.distance_to_point_l1(point) <= epsilon + + def distance_to_point_l1(self, point): + # noqa: D301 + """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear + distance) to a point with 0.0 if the point is inside the cell. + + Parameters + ---------- + point : numpy.ndarray or list + The coordinates of the point. + + Returns + ------- + float + The L1 distance, also known as taxicab, manhatten or rectilinear + distance, of the cell to the point. If 0.0 the point is inside the + cell. + + Notes + ----- + + This is done with the help of barycentric coordinates where the general + algorithm is to compute the most negative (i.e. minimum) barycentric + coordinate then return its negative. For implementation reasons we + return the sum of all the negative barycentric coordinates. In each of + the below examples the point coordinate is `X` with appropriate + dimensions. + + Consider, for example, a UFCInterval. We have two vertices which make + the interval, + `P0 = [0]` and + `P1 = [1]`. + Our point is + `X = [x]`. + Barycentric coordinates are defined as + `X = alpha * P0 + beta * P1` where + `alpha + beta = 1.0`. + The solution is + `alpha = 1 - X[0] = 1 - x` and + `beta = X[0] = x`. + If both `alpha` and `beta` are positive, the point is inside the + reference interval. + + `---regionA---P0=0------P1=1---regionB---` + + If we are in `regionA`, `alpha` is negative and + `-alpha = X[0] - 1.0` is the (positive) distance from `P0`. + If we are in `regionB`, `beta` is negative and `-beta = -X[0]` is + the exact (positive) distance from `P1`. Since we are in 1D the L1 + distance is the same as the L2 distance. If we are in the interval we + can just return 0.0. + + Things get more complicated when we consider higher dimensions. + Consider a UFCTriangle. We have three vertices which make the + reference triangle, + `P0 = (0, 0)`, + `P1 = (1, 0)` and + `P2 = (0, 1)`. + Our point is + `X = [x, y]`. + Below is a diagram of the cell (which may not render correctly in + sphinx): + + .. code-block:: text + ``` + y-axis + | + | + (0,1) P2 + | \\ + | \\ + | \\ + | \\ + | T \\ + | \\ + | \\ + | \\ + ---P0--------P1--- x-axis + (0,0) | (1,0) + ``` + + Barycentric coordinates are defined as + `X = alpha * P0 + beta * P1 + gamma * P2` where + `alpha + beta + gamma = 1.0`. + The solution is + `alpha = 1 - X[0] - X[1] = 1 - x - y`, + `beta = X[0] = x` and + `gamma = X[1] = y`. + If all three are positive, the point is inside the reference cell. + If any are negative, we are outside it. The absolute sum of any + negative barycentric coordinates usefully gives the L1 distance from + the cell to the point. For example the point (1,1) has L1 distance + 1 from the cell: on this case alpha = -1, beta = 1 and gamma = 1. + -alpha = 1 is the L1 distance. For comparison the L2 distance (the + length of the vector from the nearest point on the cell to the point) + is sqrt(0.5^2 + 0.5^2) = 0.707. Similarly the point (-1.0, -1.0) has + alpha = 3, beta = -1 and gamma = -1. The absolute sum of beta and gamma + 2 which is again the L1 distance. The L2 distance in this case is + sqrt(1^2 + 1^2) = 1.414. + + For a UFCTetrahedron we have four vertices + `P0 = (0,0,0)`, + `P1 = (1,0,0)`, + `P2 = (0,1,0)` and + `P3 = (0,0,1)`. + Our point is + `X = [x, y, z]`. + The barycentric coordinates are defined as + `X = alpha * P0 + beta * P1 + gamma * P2 + delta * P3` + where + `alpha + beta + gamma + delta = 1.0`. + The solution is + `alpha = 1 - X[0] - X[1] - X[2] = 1 - x - y - z`, + `beta = X[0] = x`, + `gamma = X[1] = y` and + `delta = X[2] = z`. + The rules are the same as for the triangle but with one extra + barycentric coordinate. Our approximate distance, the absolute sum of + the negative barycentric coordinates, is at worse around 4 times the + actual distance to the tetrahedron. + + """ + # bary = [alpha, beta, gamma, delta, ...] - see docstring + bary = [1.0 - sum(point)] + list(point) + # We avoid branching so that code can be generated from this (e.g. with + # sympy). bary-abs(bary) gets rid of the positive barycentric + # coordinates and doubles the negative distances. Summing, halving and + # taking the negative of these gives the L1 distance. So for example + # point = [-1, -1] + # bary = [3, -1, -1], + # bary-abs(bary) = [0, -2, -2], + # sum(bary-abs(bary)) = -4. + # - 0.5 * sum(bary-abs(bary)) = 2.0 + # which is the correct L1 distance from the cell to the point. + return - 0.5 * sum([b - abs(b) for b in bary]) class Point(Simplex): @@ -800,15 +946,43 @@ 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).""" - lengths = [c.get_spatial_dimension() for c in self.cells] - assert len(point) == sum(lengths) - slices = TensorProductCell._split_slices(lengths) + (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. + + """ + subcell_dimensions = [c.get_spatial_dimension() for c in self.cells] + assert len(point) == sum(subcell_dimensions) + point_slices = TensorProductCell._split_slices(subcell_dimensions) + subcell_points = [point[s] for s in point_slices] return reduce(operator.and_, - (c.contains_point(point[s], epsilon=epsilon) - for c, s in zip(self.cells, slices)), + (c.contains_point(p, epsilon=epsilon) + for c, p in zip(self.cells, subcell_points)), True) + 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. + + For more information see the docstring for the UFCSimplex method.""" + subcell_dimensions = [c.get_spatial_dimension() for c in self.cells] + assert len(point) == sum(subcell_dimensions) + point_slices = TensorProductCell._split_slices(subcell_dimensions) + subcell_points = [point[s] for s in point_slices] + subcell_distances = [c.distance_to_point_l1(p) + for c, p in zip(self.cells, subcell_points)] + return sum(subcell_distances) + def symmetry_group_size(self, dim): return tuple(c.symmetry_group_size(d) for d, c in zip(dim, self.cells)) @@ -912,9 +1086,30 @@ 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).""" + (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): + """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) + def symmetry_group_size(self, dim): return [1, 2, 8][dim] @@ -980,9 +1175,30 @@ 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).""" + (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): + """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) + def symmetry_group_size(self, dim): return [1, 2, 8, 48][dim] diff --git a/test/unit/test_reference_element.py b/test/unit/test_reference_element.py index f7b658c5d..2945508b2 100644 --- a/test/unit/test_reference_element.py +++ b/test/unit/test_reference_element.py @@ -18,6 +18,7 @@ import pytest import numpy as np import sys +from math import isclose from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron @@ -150,6 +151,252 @@ def test_reference_normal_vert(cell, normals): cell.compute_reference_normal(vert_dim, facet_number)) +@pytest.mark.parametrize(('cell', 'point', 'epsilon', 'expected'), + [(interval, [0.5], 0.0, True), + (interval, [0.0], 1e-14, True), + (interval, [1.0], 1e-14, True), + (interval, [-1e-12], 1e-11, True), + (interval, [1+1e-12], 1e-11, True), + (interval, [-1e-12], 1e-13, False), + (interval, [1+1e-12], 1e-13, False), + (triangle, [0.25, 0.25], 0.0, True), + (triangle, [0.0, 0.0], 1e-14, True), + (triangle, [1.0, 0.0], 1e-14, True), + (triangle, [0.0, 1.0], 1e-14, True), + (triangle, [0.5, 0.5], 1e-14, True), + (triangle, [-1e-12, 0.0], 1e-11, True), + (triangle, [1+1e-12, 0.0], 1e-11, True), + (triangle, [0.0, -1e-12], 1e-11, True), + (triangle, [0.0, 1+1e-12], 1e-11, True), + (triangle, [-1e-12, 0.0], 1e-13, False), + (triangle, [1+1e-12, 0.0], 1e-13, False), + (triangle, [0.0, -1e-12], 1e-13, False), + (triangle, [0.0, 1+1e-12], 1e-13, False), + (triangle, [0.5+1e-12, 0.5], 1e-13, False), + (triangle, [0.5, 0.5+1e-12], 1e-13, False), + (quadrilateral, [0.5, 0.5], 0.0, True), + (quadrilateral, [0.0, 0.0], 1e-14, True), + (quadrilateral, [1.0, 0.0], 1e-14, True), + (quadrilateral, [0.0, 1.0], 1e-14, True), + (quadrilateral, [1.0, 1.0], 1e-14, True), + (quadrilateral, [-1e-12, 0.5], 1e-11, True), + (quadrilateral, [1+1e-12, 0.5], 1e-11, True), + (quadrilateral, [0.5, -1e-12], 1e-11, True), + (quadrilateral, [0.5, 1+1e-12], 1e-11, True), + (quadrilateral, [-1e-12, 0.5], 1e-13, False), + (quadrilateral, [1+1e-12, 0.5], 1e-13, False), + (quadrilateral, [0.5, -1e-12], 1e-13, False), + (quadrilateral, [0.5, 1+1e-12], 1e-13, False), + (tetrahedron, [0.25, 0.25, 0.25], 0.0, True), + (tetrahedron, [1/3, 1/3, 1/3], 1e-14, True), + (tetrahedron, [0.0, 0.0, 0.0], 1e-14, True), + (tetrahedron, [1.0, 0.0, 0.0], 1e-14, True), + (tetrahedron, [0.0, 1.0, 0.0], 1e-14, True), + (tetrahedron, [0.0, 0.0, 1.0], 1e-14, True), + (tetrahedron, [0.0, 0.5, 0.5], 1e-14, True), + (tetrahedron, [0.5, 0.0, 0.5], 1e-14, True), + (tetrahedron, [0.5, 0.5, 0.0], 1e-14, True), + (tetrahedron, [-1e-12, 0.0, 0.0], 1e-11, True), + (tetrahedron, [1+1e-12, 0.0, 0.0], 1e-11, True), + (tetrahedron, [0.0, -1e-12, 0.0], 1e-11, True), + (tetrahedron, [0.0, 1+1e-12, 0.0], 1e-11, True), + (tetrahedron, [0.0, 0.0, -1e-12], 1e-11, True), + (tetrahedron, [0.0, 0.0, 1+1e-12], 1e-11, True), + (tetrahedron, [-1e-12, 0.0, 0.0], 1e-13, False), + (tetrahedron, [1+1e-12, 0.0, 0.0], 1e-13, False), + (tetrahedron, [0.0, -1e-12, 0.0], 1e-13, False), + (tetrahedron, [0.0, 1+1e-12, 0.0], 1e-13, False), + (tetrahedron, [0.0, 0.0, -1e-12], 1e-13, False), + (tetrahedron, [0.0, 0.0, 1+1e-12], 1e-13, False), + (tetrahedron, [0.5+1e-12, 0.5, 0.5], 1e-13, False), + (tetrahedron, [0.5, 0.5+1e-12, 0.5], 1e-13, False), + (tetrahedron, [0.5, 0.5, 0.5+1e-12], 1e-13, False), + (hexahedron, [0.5, 0.5, 0.5], 0.0, True), + (hexahedron, [0.0, 0.0, 0.0], 1e-14, True), + (hexahedron, [1.0, 0.0, 0.0], 1e-14, True), + (hexahedron, [0.0, 1.0, 0.0], 1e-14, True), + (hexahedron, [0.0, 0.0, 1.0], 1e-14, True), + (hexahedron, [1.0, 1.0, 0.0], 1e-14, True), + (hexahedron, [1.0, 0.0, 1.0], 1e-14, True), + (hexahedron, [0.0, 1.0, 1.0], 1e-14, True), + (hexahedron, [1.0, 1.0, 1.0], 1e-14, True), + (hexahedron, [-1e-12, 0.5, 0.5], 1e-11, True), + (hexahedron, [0.5, -1e-12, 0.5], 1e-11, True), + (hexahedron, [0.5, 0.5, -1e-12], 1e-11, True), + (hexahedron, [1+1e-12, 0.5, 0.5], 1e-11, True), + (hexahedron, [0.5, 1+1e-12, 0.5], 1e-11, True), + (hexahedron, [0.5, 0.5, 1+1e-12], 1e-11, True), + (hexahedron, [-1e-12, 0.5, 0.5], 1e-13, False), + (hexahedron, [0.5, -1e-12, 0.5], 1e-13, False), + (hexahedron, [0.5, 0.5, -1e-12], 1e-13, False), + (hexahedron, [1+1e-12, 0.5, 0.5], 1e-13, False), + (hexahedron, [0.5, 1+1e-12, 0.5], 1e-13, False), + (hexahedron, [0.5, 0.5, 1+1e-12], 1e-13, False), + (interval_x_interval, [0.5, 0.5], 0.0, True), + (interval_x_interval, [0.0, 0.0], 1e-14, True), + (interval_x_interval, [1.0, 0.0], 1e-14, True), + (interval_x_interval, [0.0, 1.0], 1e-14, True), + (interval_x_interval, [1.0, 1.0], 1e-14, True), + (interval_x_interval, [-1e-12, 0.5], 1e-11, True), + (interval_x_interval, [1+1e-12, 0.5], 1e-11, True), + (interval_x_interval, [0.5, -1e-12], 1e-11, True), + (interval_x_interval, [0.5, 1+1e-12], 1e-11, True), + (interval_x_interval, [-1e-12, 0.5], 1e-13, False), + (interval_x_interval, [1+1e-12, 0.5], 1e-13, False), + (interval_x_interval, [0.5, -1e-12], 1e-13, False), + (interval_x_interval, [0.5, 1+1e-12], 1e-13, False), + (triangle_x_interval, [0.25, 0.25, 0.5], 0.0, True), + (triangle_x_interval, [0.0, 0.0, 0.0], 1e-14, True), + (triangle_x_interval, [1.0, 0.0, 0.0], 1e-14, True), + (triangle_x_interval, [0.0, 1.0, 0.0], 1e-14, True), + (triangle_x_interval, [0.0, 0.0, 1.0], 1e-14, True), + (triangle_x_interval, [0.5, 0.5, 0.5], 1e-14, True), + (triangle_x_interval, [-1e-12, 0.0, 0.5], 1e-11, True), + (triangle_x_interval, [1+1e-12, 0.0, 0.5], 1e-11, True), + (triangle_x_interval, [0.0, -1e-12, 0.5], 1e-11, True), + (triangle_x_interval, [0.0, 1+1e-12, 0.5], 1e-11, True), + (triangle_x_interval, [0.0, 0.0, -1e-12], 1e-11, True), + (triangle_x_interval, [0.0, 0.0, 1+1e-12], 1e-11, True), + (triangle_x_interval, [-1e-12, 0.0, 0.5], 1e-13, False), + (triangle_x_interval, [1+1e-12, 0.0, 0.5], 1e-13, False), + (triangle_x_interval, [0.0, -1e-12, 0.5], 1e-13, False), + (triangle_x_interval, [0.0, 1+1e-12, 0.5], 1e-13, False), + (triangle_x_interval, [0.0, 0.0, -1e-12], 1e-13, False), + (triangle_x_interval, [0.0, 0.0, 1+1e-12], 1e-13, False), + (triangle_x_interval, [0.5+1e-12, 0.5, 0.5], 1e-13, False), + (triangle_x_interval, [0.5, 0.5+1e-12, 0.5], 1e-13, False), + (triangle_x_interval, [0.0, 0.0, -1e-12], 1e-13, False), + (triangle_x_interval, [0.0, 0.0, 1+1e-12], 1e-13, False), + (quadrilateral_x_interval, [0.5, 0.5, 0.5], 0.0, True), + (quadrilateral_x_interval, [0.0, 0.0, 0.0], 1e-14, True), + (quadrilateral_x_interval, [1.0, 0.0, 0.0], 1e-14, True), + (quadrilateral_x_interval, [0.0, 1.0, 0.0], 1e-14, True), + (quadrilateral_x_interval, [0.0, 0.0, 1.0], 1e-14, True), + (quadrilateral_x_interval, [-1e-12, 0.0, 0.0], 1e-11, True), + (quadrilateral_x_interval, [1+1e-12, 0.0, 0.0], 1e-11, True), + (quadrilateral_x_interval, [0.0, -1e-12, 0.0], 1e-11, True), + (quadrilateral_x_interval, [0.0, 1+1e-12, 0.0], 1e-11, True), + (quadrilateral_x_interval, [0.0, 0.0, -1e-12], 1e-11, True), + (quadrilateral_x_interval, [0.0, 0.0, 1+1e-12], 1e-11, True), + (quadrilateral_x_interval, [-1e-12, 0.0, 0.0], 1e-13, False), + (quadrilateral_x_interval, [1+1e-12, 0.0, 0.0], 1e-13, False), + (quadrilateral_x_interval, [0.0, -1e-12, 0.0], 1e-13, False), + (quadrilateral_x_interval, [0.0, 1+1e-12, 0.0], 1e-13, False), + (quadrilateral_x_interval, [0.0, 0.0, -1e-12], 1e-13, False), + (quadrilateral_x_interval, [0.0, 0.0, 1+1e-12], 1e-13, False)]) +def test_contains_point(cell, point, epsilon, expected): + assert cell.contains_point(point, epsilon) == expected + + +@pytest.mark.parametrize(('cell', 'point', 'expected'), + [(interval, [0.5], 0.0), + (interval, [0.0], 0.0), + (interval, [1.0], 0.0), + (interval, [-1e-12], 1e-12), + (interval, [1+1e-12], 1e-12), + (triangle, [0.25, 0.25], 0.0), + (triangle, [0.0, 0.0], 0.0), + (triangle, [1.0, 0.0], 0.0), + (triangle, [0.0, 1.0], 0.0), + (triangle, [0.5, 0.5], 0.0), + (triangle, [-1e-12, 0.0], 1e-12), + (triangle, [1+1e-12, 0.0], 1e-12), + (triangle, [0.0, -1e-12], 1e-12), + (triangle, [0.0, 1+1e-12], 1e-12), + (triangle, [0.5+1e-12, 0.5], 1e-12), + (triangle, [0.5, 0.5+1e-12], 1e-12), + (quadrilateral, [0.5, 0.5], 0.0), + (quadrilateral, [0.0, 0.0], 0.0), + (quadrilateral, [1.0, 0.0], 0.0), + (quadrilateral, [0.0, 1.0], 0.0), + (quadrilateral, [1.0, 1.0], 0.0), + (quadrilateral, [-1e-12, 0.5], 1e-12), + (quadrilateral, [1+1e-12, 0.5], 1e-12), + (quadrilateral, [0.5, -1e-12], 1e-12), + (quadrilateral, [0.5, 1+1e-12], 1e-12), + (quadrilateral, [-1e-12, 0.5], 1e-12), + (quadrilateral, [1+1e-12, 0.5], 1e-12), + (quadrilateral, [1+1e-12, 1+1e-12], 2e-12), + (tetrahedron, [0.25, 0.25, 0.25], 0.0), + (tetrahedron, [1/3, 1/3, 1/3], 0.0), + (tetrahedron, [0.0, 0.0, 0.0], 0.0), + (tetrahedron, [1.0, 0.0, 0.0], 0.0), + (tetrahedron, [0.0, 1.0, 0.0], 0.0), + (tetrahedron, [0.0, 0.0, 1.0], 0.0), + (tetrahedron, [0.0, 0.5, 0.5], 0.0), + (tetrahedron, [0.5, 0.0, 0.5], 0.0), + (tetrahedron, [0.5, 0.5, 0.0], 0.0), + (tetrahedron, [-1e-12, 0.0, 0.0], 1e-12), + (tetrahedron, [1+1e-12, 0.0, 0.0], 1e-12), + (tetrahedron, [0.0, -1e-12, 0.0], 1e-12), + (tetrahedron, [0.0, 1+1e-12, 0.0], 1e-12), + (tetrahedron, [0.0, 0.0, -1e-12], 1e-12), + (tetrahedron, [0.0, 0.0, 1+1e-12], 1e-12), + (tetrahedron, [1/3+1e-12, 1/3, 1/3], 1e-12), + (tetrahedron, [1/3, 1/3+1e-12, 1/3], 1e-12), + (tetrahedron, [1/3, 1/3, 1/3+1e-12], 1e-12), + (hexahedron, [0.5, 0.5, 0.5], 0.0), + (hexahedron, [0.0, 0.0, 0.0], 0.0), + (hexahedron, [1.0, 0.0, 0.0], 0.0), + (hexahedron, [0.0, 1.0, 0.0], 0.0), + (hexahedron, [0.0, 0.0, 1.0], 0.0), + (hexahedron, [1.0, 1.0, 0.0], 0.0), + (hexahedron, [1.0, 0.0, 1.0], 0.0), + (hexahedron, [0.0, 1.0, 1.0], 0.0), + (hexahedron, [1.0, 1.0, 1.0], 0.0), + (hexahedron, [-1e-12, 0.5, 0.5], 1e-12), + (hexahedron, [0.5, -1e-12, 0.5], 1e-12), + (hexahedron, [0.5, 0.5, -1e-12], 1e-12), + (hexahedron, [1+1e-12, 0.5, 0.5], 1e-12), + (hexahedron, [0.5, 1+1e-12, 0.5], 1e-12), + (hexahedron, [0.5, 0.5, 1+1e-12], 1e-12), + (hexahedron, [-1e-12, -1e-12, -1e-12], 3e-12), + (hexahedron, [1.0+1e-12, -1e-12, -1e-12], 3e-12), + (hexahedron, [-1e-12, 1.0+1e-12, -1e-12], 3e-12), + (hexahedron, [-1e-12, -1e-12, 1.0+1e-12], 3e-12), + (hexahedron, [1.0+1e-12, 1.0+1e-12, -1e-12], 3e-12), + (hexahedron, [1.0+1e-12, -1e-12, 1.0+1e-12], 3e-12), + (hexahedron, [-1e-12, 1.0+1e-12, 1.0+1e-12], 3e-12), + (hexahedron, [1.0+1e-12, 1.0+1e-12, 1.0+1e-12], 3e-12), + (interval_x_interval, [0.5, 0.5], 0.0), + (interval_x_interval, [0.0, 0.0], 0.0), + (interval_x_interval, [1.0, 0.0], 0.0), + (interval_x_interval, [0.0, 1.0], 0.0), + (interval_x_interval, [1.0, 1.0], 0.0), + (interval_x_interval, [-1e-12, 0.5], 1e-12), + (interval_x_interval, [1+1e-12, 0.5], 1e-12), + (interval_x_interval, [0.5, -1e-12], 1e-12), + (interval_x_interval, [0.5, 1+1e-12], 1e-12), + (interval_x_interval, [-1e-12, 0.5], 1e-12), + (interval_x_interval, [1+1e-12, 0.5], 1e-12), + (interval_x_interval, [1+1e-12, 1+1e-12], 2e-12), + (triangle_x_interval, [0.25, 0.25, 0.5], 0.0), + (triangle_x_interval, [0.0, 0.0, 0.0], 0.0), + (triangle_x_interval, [1.0, 0.0, 0.0], 0.0), + (triangle_x_interval, [0.0, 1.0, 0.0], 0.0), + (triangle_x_interval, [0.0, 0.0, 1.0], 0.0), + (triangle_x_interval, [0.5, 0.5, 0.5], 0.0), + (triangle_x_interval, [-1e-12, 0.0, 0.5], 1e-12), + (triangle_x_interval, [1+1e-12, 0.0, 0.5], 1e-12), + (triangle_x_interval, [0.0, -1e-12, 0.5], 1e-12), + (triangle_x_interval, [0.0, 1+1e-12, 0.5], 1e-12), + (triangle_x_interval, [0.0, 0.0, -1e-12], 1e-12), + (triangle_x_interval, [0.0, 0.0, 1+1e-12], 1e-12), + (quadrilateral_x_interval, [0.5, 0.5, 0.5], 0.0), + (quadrilateral_x_interval, [0.0, 0.0, 0.0], 0.0), + (quadrilateral_x_interval, [1.0, 0.0, 0.0], 0.0), + (quadrilateral_x_interval, [0.0, 1.0, 0.0], 0.0), + (quadrilateral_x_interval, [0.0, 0.0, 1.0], 0.0), + (quadrilateral_x_interval, [-1e-12, 0.0, 0.0], 1e-12), + (quadrilateral_x_interval, [1+1e-12, 0.0, 0.0], 1e-12), + (quadrilateral_x_interval, [0.0, -1e-12, 0.0], 1e-12), + (quadrilateral_x_interval, [0.0, 1+1e-12, 0.0], 1e-12), + (quadrilateral_x_interval, [0.0, 0.0, -1e-12], 1e-12), + (quadrilateral_x_interval, [0.0, 0.0, 1+1e-12], 1e-12)]) +def test_distance_to_point_l1(cell, point, expected): + assert isclose(cell.distance_to_point_l1(point), expected, rel_tol=1e-3) + + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__))