Skip to content

Commit

Permalink
Vertex only mesh tolerance (#2449)
Browse files Browse the repository at this point in the history
Roundoff error can sometimes cause spurious loss of vertices when constructing a VertexOnlyMesh. This PR enables the user to work around this by setting the tolerance by which the local coordinates of a point will be allowed to lie outside a cell while still counting that point as within that cell.
  • Loading branch information
dham authored May 27, 2022
1 parent 10f92cb commit 6364be6
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 6 deletions.
18 changes: 12 additions & 6 deletions firedrake/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -2091,7 +2091,8 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', kern


@PETSc.Log.EventDecorator()
def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None):
def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None,
tolerance=None):
"""
Create a vertex only mesh, immersed in a given mesh, with vertices defined
by a list of coordinates.
Expand All @@ -2105,6 +2106,11 @@ def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None):
that they have the same list of vertices (else the test is not
possible): this operation scales with number of vertices and number of
ranks.
:kwarg tolerance: the amount by which the local coordinates of a point are
allowed to fall outside the cell while still having the point count as
in the cell. Increase the default (1.0e-14) somewhat if vertices
interior to the domain are being lost in the :class:`VertexOnlyMesh`
construction process.
.. note::
Expand Down Expand Up @@ -2175,7 +2181,7 @@ def VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour=None):
if pdim != gdim:
raise ValueError(f"Mesh geometric dimension {gdim} must match point list dimension {pdim}")

swarm = _pic_swarm_in_mesh(mesh, vertexcoords)
swarm = _pic_swarm_in_mesh(mesh, vertexcoords, tolerance=tolerance)

if missing_points_behaviour:

Expand Down Expand Up @@ -2249,7 +2255,7 @@ def compare_arrays(x, y, datatype):
return vmesh


def _pic_swarm_in_mesh(parent_mesh, coords, fields=None):
def _pic_swarm_in_mesh(parent_mesh, coords, fields=None, tolerance=None):
"""Create a Particle In Cell (PIC) DMSwarm immersed in a Mesh
This should only by used for meshes with straight edges. If not, the
Expand Down Expand Up @@ -2314,7 +2320,7 @@ def _pic_swarm_in_mesh(parent_mesh, coords, fields=None):
fields += [("parentcellnum", 1, IntType), ("refcoord", tdim, RealType)]

coords, reference_coords, parent_cell_nums = \
_parent_mesh_embedding(coords, parent_mesh)
_parent_mesh_embedding(coords, parent_mesh, tolerance)
# mesh.topology.cell_closure[:, -1] maps Firedrake cell numbers to plex numbers.
plex_parent_cell_nums = parent_mesh.topology.cell_closure[parent_cell_nums, -1]

Expand Down Expand Up @@ -2381,7 +2387,7 @@ def _pic_swarm_in_mesh(parent_mesh, coords, fields=None):
return swarm


def _parent_mesh_embedding(vertex_coords, parent_mesh):
def _parent_mesh_embedding(vertex_coords, parent_mesh, tolerance):
"""Find the parent cells and local coords for vertices on this rank.
Vertices not located in cells owned by this rank are discarded.
Expand All @@ -2405,7 +2411,7 @@ def _parent_mesh_embedding(vertex_coords, parent_mesh):
for i in range(max_num_vertices):
if i < num_vertices:
parent_cell_num, reference_coord = \
parent_mesh.locate_cell_and_reference_coordinate(vertex_coords[i])
parent_mesh.locate_cell_and_reference_coordinate(vertex_coords[i], tolerance)
# parent_cell_num >= parent_mesh.cell_set.size means the vertex is in the halo
# and is to be discarded.
if parent_cell_num is not None and parent_cell_num < parent_mesh.cell_set.size:
Expand Down
18 changes: 18 additions & 0 deletions tests/vertexonly/test_vertex_only_mesh_generation.py
Original file line number Diff line number Diff line change
Expand Up @@ -176,3 +176,21 @@ def test_extrude(parentmesh):
inputcoords, inputcoordslocal = cell_midpoints(parentmesh)
vm = VertexOnlyMesh(parentmesh, inputcoords)
ExtrudedMesh(vm, 1)


def test_point_tolerance():
"""Test the tolerance parameter to VertexOnlyMesh.
This test works by checking a point outside the domain. Tolerance does not
in fact promise to fix the problem of points outside the domain in the
general case. It is instead there to cope with losing points on internal
cell boundaries due to roundoff. The latter case is difficult to test in a
manner which is robust to roundoff behaviour in different environments."""
m = UnitSquareMesh(1, 1)
# Make the mesh non-axis-aligned.
m.coordinates.dat.data[1, :] = [1.1, 1]
coords = [[1.0501, 0.5]]
vm = VertexOnlyMesh(m, coords, tolerance=0.1)
assert vm.cell_set.size == 1
vm = VertexOnlyMesh(m, coords, tolerance=None)
assert vm.cell_set.size == 0

0 comments on commit 6364be6

Please sign in to comment.