From 9f2597ad3bc49bb3e55e608853e9be3ba0714839 Mon Sep 17 00:00:00 2001 From: ordinary-slim Date: Mon, 26 Aug 2024 09:53:22 +0000 Subject: [PATCH 1/9] Replacing np by npt in python function annotations --- python/dolfinx/fem/bcs.py | 5 +++-- python/dolfinx/fem/forms.py | 4 ++-- python/dolfinx/fem/function.py | 10 +++++----- python/dolfinx/graph.py | 4 +++- python/dolfinx/mesh.py | 10 +++++----- 5 files changed, 18 insertions(+), 15 deletions(-) diff --git a/python/dolfinx/fem/bcs.py b/python/dolfinx/fem/bcs.py index 1972740bfcb..09748d23e0b 100644 --- a/python/dolfinx/fem/bcs.py +++ b/python/dolfinx/fem/bcs.py @@ -17,6 +17,7 @@ from dolfinx.fem.function import Constant, Function import numpy as np +import numpy.typing as npt import dolfinx from dolfinx import cpp as _cpp @@ -57,7 +58,7 @@ def locate_dofs_geometrical( def locate_dofs_topological( V: typing.Union[dolfinx.fem.FunctionSpace, typing.Iterable[dolfinx.fem.FunctionSpace]], entity_dim: int, - entities: numpy.typing.NDArray[np.int32], + entities: npt.NDArray[np.int32], remote: bool = True, ) -> np.ndarray: """Locate degrees-of-freedom belonging to mesh entities topologically. @@ -135,7 +136,7 @@ def function_space(self) -> dolfinx.fem.FunctionSpace: def dirichletbc( value: typing.Union[Function, Constant, np.ndarray], - dofs: numpy.typing.NDArray[np.int32], + dofs: npt.NDArray[np.int32], V: typing.Optional[dolfinx.fem.FunctionSpace] = None, ) -> DirichletBC: """Create a representation of Dirichlet boundary condition which diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 6c4ed21311c..d732e382b85 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -102,7 +102,7 @@ def integral_types(self): def get_integration_domains( integral_type: IntegralType, - subdomain: typing.Optional[typing.Union[MeshTags, list[tuple[int, np.ndarray]]]], + subdomain: typing.Optional[typing.Union[MeshTags, list[tuple[int, npt.NDArray[np.int32]]]]], subdomain_ids: list[int], ) -> list[tuple[int, np.ndarray]]: """Get integration domains from subdomain data. @@ -190,7 +190,7 @@ def form( dtype: npt.DTypeLike = default_scalar_type, form_compiler_options: typing.Optional[dict] = None, jit_options: typing.Optional[dict] = None, - entity_maps: dict[Mesh, np.typing.NDArray[np.int32]] = {}, + entity_maps: dict[Mesh, npt.NDArray[np.int32]] = {}, ): """Create a Form or an array of Forms. diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index a04ce990c67..2e0f395d955 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -90,7 +90,7 @@ class Expression: def __init__( self, e: ufl.core.expr.Expr, - X: np.ndarray, + X: typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]], comm: typing.Optional[_MPI.Comm] = None, form_compiler_options: typing.Optional[dict] = None, jit_options: typing.Optional[dict] = None, @@ -196,7 +196,7 @@ def _create_expression(dtype): def eval( self, mesh: Mesh, - entities: np.ndarray, + entities: npt.NDArray[np.int32], values: typing.Optional[np.ndarray] = None, ) -> np.ndarray: """Evaluate Expression on entities. @@ -412,8 +412,8 @@ def interpolate_nonmatching( def interpolate( self, u0: typing.Union[typing.Callable, Expression, Function], - cells0: typing.Optional[np.ndarray] = None, - cells1: typing.Optional[np.ndarray] = None, + cells0: typing.Optional[npt.NDArray[np.int32]] = None, + cells1: typing.Optional[npt.NDArray[np.int32]] = None, ) -> None: """Interpolate an expression. @@ -582,7 +582,7 @@ def _create_dolfinx_element( comm: _MPI.Intracomm, cell_type: _cpp.mesh.CellType, ufl_e: ufl.FiniteElementBase, - dtype: np.dtype, + dtype: npt.DTypeLike, ) -> typing.Union[_cpp.fem.FiniteElement_float32, _cpp.fem.FiniteElement_float64]: """Create a DOLFINx element from a basix.ufl element.""" if np.issubdtype(dtype, np.float32): diff --git a/python/dolfinx/graph.py b/python/dolfinx/graph.py index dd16514a79a..fb8c7e83923 100644 --- a/python/dolfinx/graph.py +++ b/python/dolfinx/graph.py @@ -8,6 +8,8 @@ from __future__ import annotations import numpy as np +import typing +import numpy.typing as npt from dolfinx import cpp as _cpp from dolfinx.cpp.graph import partitioner @@ -31,7 +33,7 @@ __all__ = ["adjacencylist", "partitioner"] -def adjacencylist(data: np.ndarray, offsets=None): +def adjacencylist(data: typing.Union[npt.NDArray[np.int32], npt.NDArray[np.int64]], offsets: typing.Optional[npt.NDArray[np.int32]]=None): """Create an AdjacencyList for int32 or int64 datasets. Args: diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 5f1a0da94bb..a6c58240046 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -222,7 +222,7 @@ def compute_midpoints(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]): return _cpp.mesh.compute_midpoints(mesh._cpp_object, dim, entities) -def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: +def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> npt.NDArray[np.int32]: """Compute mesh entities satisfying a geometric marking function. Args: @@ -239,7 +239,7 @@ def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray return _cpp.mesh.locate_entities(mesh._cpp_object, dim, marker) -def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: +def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> npt.NDArray[np.int32]: """Compute mesh entities that are connected to an owned boundary facet and satisfy a geometric marking function. @@ -313,7 +313,7 @@ def transfer_meshtag( def refine( - mesh: Mesh, edges: typing.Optional[np.ndarray] = None, redistribute: bool = True + mesh: Mesh, edges: typing.Optional[npt.NDArray[np.int32]] = None, redistribute: bool = True ) -> Mesh: """Refine a mesh. @@ -336,7 +336,7 @@ def refine( def refine_interval( mesh: Mesh, - cells: typing.Optional[np.ndarray] = None, + cells: typing.Optional[npt.NDArray[np.int32]] = None, redistribute: bool = True, ghost_mode: GhostMode = GhostMode.shared_facet, ) -> tuple[Mesh, npt.NDArray[np.int32]]: @@ -367,7 +367,7 @@ def refine_interval( def refine_plaza( mesh: Mesh, - edges: typing.Optional[np.ndarray] = None, + edges: typing.Optional[npt.NDArray[np.int32]] = None, redistribute: bool = True, option: RefinementOption = RefinementOption.none, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int32]]: From 434aa5e269d8ca20f15520f465e7793eb09a5895 Mon Sep 17 00:00:00 2001 From: ordinary-slim Date: Mon, 26 Aug 2024 12:16:08 +0000 Subject: [PATCH 2/9] Improve python exposure determine_point_ownership * Python wrapper around nanobind exposed function * Extra optional arguments `cells` --- cpp/dolfinx/fem/interpolate.h | 7 +- cpp/dolfinx/geometry/utils.h | 20 +-- python/dolfinx/geometry.py | 43 +++++ python/dolfinx/wrappers/geometry.cpp | 9 +- .../unit/geometry/test_bounding_box_tree.py | 158 ++++++++++++++++++ 5 files changed, 222 insertions(+), 15 deletions(-) diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index 12956057cb7..55a4cee93bd 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -1070,7 +1070,12 @@ geometry::PointOwnershipData create_interpolation_data( x[3 * i + j] = coords[i + j * num_points]; // Determine ownership of each point - return geometry::determine_point_ownership(mesh1, x, padding); + const int tdim = mesh1.topology()->dim(); + auto cell_map1 = mesh1.topology()->index_map(tdim); + const std::int32_t num_cells1 = cell_map1->size_local(); + std::vector cells1(num_cells1, 0); + std::iota(cells1.begin(), cells1.end(), 0); + return geometry::determine_point_ownership(mesh1, x, cells1, padding); } /// @brief Interpolate a finite element Function defined on a mesh to a diff --git a/cpp/dolfinx/geometry/utils.h b/cpp/dolfinx/geometry/utils.h index 0f643e62643..bf8e13b000e 100644 --- a/cpp/dolfinx/geometry/utils.h +++ b/cpp/dolfinx/geometry/utils.h @@ -663,21 +663,21 @@ graph::AdjacencyList compute_colliding_cells( /// @param[in] mesh The mesh /// @param[in] points Points to check for collision (`shape=(num_points, /// 3)`). Storage is row-major. +/// @param[in] cells Cells to check for ownership /// @param[in] padding Amount of absolute padding of bounding boxes of the mesh. /// Each bounding box of the mesh is padded with this amount, to increase /// the number of candidates, avoiding rounding errors in determining the owner /// of a point if the point is on the surface of a cell in the mesh. -/// @return Tuple `(src_owner, dest_owner, dest_points, dest_cells)`, -/// where src_owner is a list of ranks corresponding to the input -/// points. dest_owner is a list of ranks corresponding to dest_points, -/// the points that this process owns. dest_cells contains the +/// @return Point ownership data with fields `src_owner`, `dest_owner`, `dest_points`, `dest_cells`, +/// where `src_owner` is a list of ranks corresponding to the input +/// points. `dest_owner` is a list of ranks corresponding to `dest_points`, +/// the points that this process owns. `dest_cells` contains the /// corresponding cell for each entry in dest_points. /// /// @note `dest_owner` is sorted -/// @note Returns -1 if no colliding process is found +/// @note `src_owner` is -1 if no colliding process is found /// @note dest_points is flattened row-major, shape `(dest_owner.size(), /// 3)` -/// @note Only looks through cells owned by the process /// @note A large padding value can increase the runtime of the function by /// orders of magnitude, because for non-colliding cells /// one has to determine the closest cell among all processes with an @@ -685,18 +685,14 @@ graph::AdjacencyList compute_colliding_cells( template PointOwnershipData determine_point_ownership(const mesh::Mesh& mesh, std::span points, - T padding) + std::span cells, + T padding = 0.0) { MPI_Comm comm = mesh.comm(); // Create a global bounding-box tree to find candidate processes with // cells that could collide with the points const int tdim = mesh.topology()->dim(); - auto cell_map = mesh.topology()->index_map(tdim); - const std::int32_t num_cells = cell_map->size_local(); - // NOTE: Should we send the cells in as input? - std::vector cells(num_cells, 0); - std::iota(cells.begin(), cells.end(), 0); BoundingBoxTree bb(mesh, tdim, cells, padding); BoundingBoxTree global_bbtree = bb.create_global_tree(comm); diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py index 03cd0870b24..1e58c020441 100644 --- a/python/dolfinx/geometry.py +++ b/python/dolfinx/geometry.py @@ -270,3 +270,46 @@ def compute_distance_gjk( """ return _cpp.geometry.compute_distance_gjk(p, q) + + +def determine_point_ownership( + mesh: Mesh, + points: npt.NDArray[np.floating], + cells: typing.Optional[npt.NDArray[np.int32]] = None, + padding: float = 0.0, +) -> PointOwnershipData: + """Build point ownership data for a mesh-points pair. + + First, potential collisions are found by computing intersections + between the bounding boxes of the cells and the set of points. + Then, actual containment pairs are determined using the GJK algorithm. + + Args: + mesh: The mesh + points: Points to check for collision (``shape=(num_points, gdim)``) + cells: Cells to check for ownership + If ``None`` then all cells are considered. + padding: Amount of absolute padding of bounding boxes of the mesh. + Each bounding box of the mesh is padded with this amount, to increase + the number of candidates, avoiding rounding errors in determining the owner + of a point if the point is on the surface of a cell in the mesh. + + Returns: + Point ownership data + + Note: + `dest_owner` is sorted + + `src_owner` is -1 if no colliding process is found + + A large padding value will increase the run-time of the code by orders + of magnitude. General advice is to use a padding on the scale of the + cell size. + + """ + if cells is None: + map = mesh.topology.index_map(mesh.topology.dim) + cells = np.arange(map.size_local, dtype=np.int32) + return PointOwnershipData( + _cpp.geometry.determine_point_ownership(mesh._cpp_object, points, cells, padding) + ) diff --git a/python/dolfinx/wrappers/geometry.cpp b/python/dolfinx/wrappers/geometry.cpp index d64e70cb584..ce0931134af 100644 --- a/python/dolfinx/wrappers/geometry.cpp +++ b/python/dolfinx/wrappers/geometry.cpp @@ -180,13 +180,18 @@ void declare_bbtree(nb::module_& m, std::string type) nb::arg("mesh"), nb::arg("dim"), nb::arg("indices"), nb::arg("points")); m.def("determine_point_ownership", [](const dolfinx::mesh::Mesh& mesh, - nb::ndarray points, const T padding) + nb::ndarray points, + nb::ndarray, nb::c_contig> cells, + const T padding) { const std::size_t p_s0 = points.ndim() == 1 ? 1 : points.shape(0); std::span _p(points.data(), 3 * p_s0); return dolfinx::geometry::determine_point_ownership(mesh, _p, + std::span(cells.data(), cells.size()), padding); - }); + }, + nb::arg("mesh"), nb::arg("points"), nb::arg("cells"), nb::arg("padding") = 0.0, + "Compute point ownership data for mesh-points pair."); std::string pod_pyclass_name = "PointOwnershipData_" + type; nb::class_>(m, diff --git a/python/test/unit/geometry/test_bounding_box_tree.py b/python/test/unit/geometry/test_bounding_box_tree.py index 2cc328300b6..f76d1ccd05f 100644 --- a/python/test/unit/geometry/test_bounding_box_tree.py +++ b/python/test/unit/geometry/test_bounding_box_tree.py @@ -12,6 +12,7 @@ from dolfinx import cpp as _cpp from dolfinx.geometry import ( + PointOwnershipData, bb_tree, compute_closest_entity, compute_colliding_cells, @@ -19,6 +20,7 @@ compute_collisions_trees, compute_distance_gjk, create_midpoint_tree, + determine_point_ownership, ) from dolfinx.mesh import ( CellType, @@ -496,3 +498,159 @@ def test_surface_bbtree_collision(dtype): collisions = compute_collisions_trees(bbtree1, bbtree2) assert len(collisions) == 1 + + +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +def test_determine_point_ownership(dim, dtype): + """Find point owners (ranks and cells) using bounding box trees + global communication + and compare to point ownership data results.""" + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + tdim = dim + num_cells_side = 4 + h = dtype(1.0 / num_cells_side) + if tdim == 2: + mesh = create_unit_square( + MPI.COMM_WORLD, num_cells_side, num_cells_side, CellType.quadrilateral, dtype=dtype + ) + else: + mesh = create_unit_cube( + MPI.COMM_WORLD, + num_cells_side, + num_cells_side, + num_cells_side, + CellType.hexahedron, + dtype=dtype, + ) + cell_map = mesh.topology.index_map(tdim) + + tree = bb_tree(mesh, mesh.topology.dim, np.arange(cell_map.size_local)) + num_global_cells = num_cells_side**dim + all_midpoints = np.zeros((num_global_cells, 3), dtype=dtype) + counter = 0 + Xs = [(i + 0.5) * h for i in range(num_cells_side)] + Ys = Xs + Zs = [0.0] if dim == 2 else Xs + for x in Xs: + for y in Ys: + for z in Zs: + all_midpoints[counter, :] = np.array([x, y, z], dtype=dtype) + counter += 1 + # Since ghost cells are left out and the points considered are midpoints + # of cells, they are only contained in a single process / cell + # Moreover, the bounding boxes of the elements correspond to the elements, + # hence the tree collisions are the exact collisions + tree_col = compute_collisions_points(tree, all_midpoints) + processwise_owners = np.zeros(2 * num_global_cells, dtype=np.int32) + owners = np.empty_like(processwise_owners) + for ipoint in range(num_global_cells): + cell = tree_col.links(ipoint) + if len(cell) > 0: + processwise_owners[2 * ipoint] = rank + processwise_owners[2 * ipoint + 1] = cell[0] + + # The value at a given index is null if it doesn't correspond + # to the current process. + # We can sum the processwise arrays to obtain a global array + comm.Allreduce(processwise_owners, owners, op=MPI.SUM) + owner_ranks = owners[[2 * i for i in range(num_global_cells)]] + owner_cells = owners[[2 * i + 1 for i in range(num_global_cells)]] + + # Reorganize ownership data (point, index, rank, cell) into dictionary + ownership_data = {} + for ipoint in range(num_global_cells): + ownership_data[tuple(all_midpoints[ipoint])] = ( + ipoint, + owner_ranks[ipoint], + owner_cells[ipoint], + ) + + def check_po(po: PointOwnershipData, src_points, ownership_data, global_dest_owners): + """ + Check point ownership data + + po: PointOwnershipData object to check + src_points: Points sent by process + ownership_data: {point:(global_index,rank,cell} + global_dest_owners: Rank who sent each point + """ + # Check src_owner: Check owner ranks of sent points + src_owner = po.src_owner() + for ipoint in range(src_points.shape[0]): + assert ownership_data[tuple(src_points[ipoint])][1] == src_owner[ipoint] + + dest_points = po.dest_points() + dest_owners = po.dest_owner() + dest_cells = po.dest_cells() + + # Check dest_points: All points that should have been found have been found + dest_points_indices = list(range(dest_points.shape[0])) + for point, data in ownership_data.items(): + (iglobal, processor, _) = data + if processor == rank: + found = False + point = np.array(point, dtype) + for jpoint in dest_points_indices: + found = np.allclose(point, dest_points[jpoint]) + if found: + break + assert found + dest_points_indices.remove(jpoint) + + # Check dest_owners and dest_cells + # dest_owners: Ranks that asked about the points we own + # dest_cells: Local index of cell that contains the points we own + for ipoint in range(dest_points.shape[0]): + iglobal = ownership_data[tuple(dest_points[ipoint])][0] + c = ownership_data[tuple(dest_points[ipoint])][2] + assert dest_owners[ipoint] == global_dest_owners[iglobal] + assert dest_cells[ipoint] == c + + def set_local_range(array): + N = array.shape[0] + n = N // comm.size + r = N % comm.size + # First r processes has one extra value + if rank < r: + (start, stop) = [rank * (n + 1), (rank + 1) * (n + 1)] + else: + (start, stop) = [rank * n + r, (rank + 1) * n + r] + return array[start:stop], start, stop + + def compute_global_owners(N, start, stop): + """Compute array of ranks who own each point""" + mask_points_owned = np.zeros(N, np.int32) + global_owners = np.empty_like(mask_points_owned) + mask_points_owned[start:stop] = rank + comm.Allreduce(mask_points_owned, global_owners, op=MPI.SUM) + return global_owners + + # All cells + points, start, stop = set_local_range(all_midpoints) + owners = compute_global_owners(np.int64(all_midpoints.shape[0]), start, stop) + all_cells = np.arange(cell_map.size_local, dtype=np.int32) + po = determine_point_ownership(mesh, points, all_cells) + + check_po(po, points, ownership_data, owners) + + # Left half + num_left_cells = np.rint((num_cells_side**dim) / 2).astype(np.int32) + left_midpoints = all_midpoints[np.arange(num_left_cells), :] + points, start, stop = set_local_range(left_midpoints) + owners = compute_global_owners(np.int64(all_midpoints.shape[0]), start, stop) + left_cells = locate_entities(mesh, tdim, lambda x: x[0] <= 0.5) + left_cells = np.array( + [cell for cell in left_cells if cell < cell_map.size_local], dtype=np.int32 + ) # Filter out ghost cells + lpo = determine_point_ownership(mesh, points, left_cells) + + left_ownership_data = {} + for ipoint in range(num_left_cells): + left_ownership_data[tuple(all_midpoints[ipoint])] = ( + ipoint, + owner_ranks[ipoint], + owner_cells[ipoint], + ) + check_po(lpo, points, left_ownership_data, owners) From 5b6726fbd2ce4df846f69b89c17c4e1477553d6f Mon Sep 17 00:00:00 2001 From: ordinary-slim Date: Thu, 29 Aug 2024 08:48:04 +0000 Subject: [PATCH 3/9] Ruff --- python/dolfinx/graph.py | 8 ++++++-- python/dolfinx/mesh.py | 4 +++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/graph.py b/python/dolfinx/graph.py index fb8c7e83923..df8d27dbfae 100644 --- a/python/dolfinx/graph.py +++ b/python/dolfinx/graph.py @@ -7,8 +7,9 @@ from __future__ import annotations -import numpy as np import typing + +import numpy as np import numpy.typing as npt from dolfinx import cpp as _cpp @@ -33,7 +34,10 @@ __all__ = ["adjacencylist", "partitioner"] -def adjacencylist(data: typing.Union[npt.NDArray[np.int32], npt.NDArray[np.int64]], offsets: typing.Optional[npt.NDArray[np.int32]]=None): +def adjacencylist( + data: typing.Union[npt.NDArray[np.int32], npt.NDArray[np.int64]], + offsets: typing.Optional[npt.NDArray[np.int32]] = None, +): """Create an AdjacencyList for int32 or int64 datasets. Args: diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 84a9b1a5648..d1edb0eb812 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -239,7 +239,9 @@ def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> npt.NDArra return _cpp.mesh.locate_entities(mesh._cpp_object, dim, marker) -def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> npt.NDArray[np.int32]: +def locate_entities_boundary( + mesh: Mesh, dim: int, marker: typing.Callable +) -> npt.NDArray[np.int32]: """Compute mesh entities that are connected to an owned boundary facet and satisfy a geometric marking function. From b06ab5a79dddab9efde79515daecbb0249dfb40b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 10 Sep 2024 11:52:40 +0200 Subject: [PATCH 4/9] Apply suggestions from code review --- python/dolfinx/wrappers/geometry.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/wrappers/geometry.cpp b/python/dolfinx/wrappers/geometry.cpp index ce0931134af..beeb6da681b 100644 --- a/python/dolfinx/wrappers/geometry.cpp +++ b/python/dolfinx/wrappers/geometry.cpp @@ -190,7 +190,7 @@ void declare_bbtree(nb::module_& m, std::string type) std::span(cells.data(), cells.size()), padding); }, - nb::arg("mesh"), nb::arg("points"), nb::arg("cells"), nb::arg("padding") = 0.0, + nb::arg("mesh"), nb::arg("points"), nb::arg("cells"), nb::arg("padding"), "Compute point ownership data for mesh-points pair."); std::string pod_pyclass_name = "PointOwnershipData_" + type; From 75a3bc52a499a8e6378dcd28d6d2038919dc7776 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Thu, 12 Sep 2024 12:49:53 +0200 Subject: [PATCH 5/9] Change to rst in docstring --- python/dolfinx/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py index 1e58c020441..8bcff028fea 100644 --- a/python/dolfinx/geometry.py +++ b/python/dolfinx/geometry.py @@ -298,9 +298,9 @@ def determine_point_ownership( Point ownership data Note: - `dest_owner` is sorted + ``dest_owner`` is sorted - `src_owner` is -1 if no colliding process is found + ``src_owner`` is -1 if no colliding process is found A large padding value will increase the run-time of the code by orders of magnitude. General advice is to use a padding on the scale of the From 942023804a2278f22b519d04d935dfab5ac3fde9 Mon Sep 17 00:00:00 2001 From: Mehdi Slimani Date: Tue, 1 Oct 2024 15:13:23 +0000 Subject: [PATCH 6/9] Wider `test_determine_point_ownership` coverage `test_determine_point_ownership` tests for triangle and tetra elements. --- cpp/dolfinx/fem/interpolate.h | 7 +- cpp/dolfinx/geometry/utils.h | 20 +-- python/dolfinx/geometry.py | 7 +- python/dolfinx/wrappers/geometry.cpp | 6 +- .../unit/geometry/test_bounding_box_tree.py | 121 +++++++++++++----- 5 files changed, 109 insertions(+), 52 deletions(-) diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index 2d3d245f99f..e25fdcbfa2b 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -1093,12 +1093,7 @@ geometry::PointOwnershipData create_interpolation_data( x[3 * i + j] = coords[i + j * num_points]; // Determine ownership of each point - const int tdim = mesh1.topology()->dim(); - auto cell_map1 = mesh1.topology()->index_map(tdim); - const std::int32_t num_cells1 = cell_map1->size_local(); - std::vector cells1(num_cells1, 0); - std::iota(cells1.begin(), cells1.end(), 0); - return geometry::determine_point_ownership(mesh1, x, cells1, padding); + return geometry::determine_point_ownership(mesh1, x, padding); } /// @brief Interpolate a finite element Function defined on a mesh to a diff --git a/cpp/dolfinx/geometry/utils.h b/cpp/dolfinx/geometry/utils.h index bf8e13b000e..0c772937f02 100644 --- a/cpp/dolfinx/geometry/utils.h +++ b/cpp/dolfinx/geometry/utils.h @@ -668,11 +668,7 @@ graph::AdjacencyList compute_colliding_cells( /// Each bounding box of the mesh is padded with this amount, to increase /// the number of candidates, avoiding rounding errors in determining the owner /// of a point if the point is on the surface of a cell in the mesh. -/// @return Point ownership data with fields `src_owner`, `dest_owner`, `dest_points`, `dest_cells`, -/// where `src_owner` is a list of ranks corresponding to the input -/// points. `dest_owner` is a list of ranks corresponding to `dest_points`, -/// the points that this process owns. `dest_cells` contains the -/// corresponding cell for each entry in dest_points. +/// @return Point ownership data. /// /// @note `dest_owner` is sorted /// @note `src_owner` is -1 if no colliding process is found @@ -685,14 +681,22 @@ graph::AdjacencyList compute_colliding_cells( template PointOwnershipData determine_point_ownership(const mesh::Mesh& mesh, std::span points, - std::span cells, - T padding = 0.0) + T padding, + std::span cells = {}) { MPI_Comm comm = mesh.comm(); + const int tdim = mesh.topology()->dim(); + + std::vector local_cells; + if (cells.empty()) { + auto cell_map = mesh.topology()->index_map(tdim); + local_cells.resize(cell_map->size_local()); + std::iota(local_cells.begin(), local_cells.end(), 0); + cells = std::span(local_cells.data(), local_cells.size()); + } // Create a global bounding-box tree to find candidate processes with // cells that could collide with the points - const int tdim = mesh.topology()->dim(); BoundingBoxTree bb(mesh, tdim, cells, padding); BoundingBoxTree global_bbtree = bb.create_global_tree(comm); diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py index 8bcff028fea..b356b346b4c 100644 --- a/python/dolfinx/geometry.py +++ b/python/dolfinx/geometry.py @@ -275,8 +275,8 @@ def compute_distance_gjk( def determine_point_ownership( mesh: Mesh, points: npt.NDArray[np.floating], + padding: float, cells: typing.Optional[npt.NDArray[np.int32]] = None, - padding: float = 0.0, ) -> PointOwnershipData: """Build point ownership data for a mesh-points pair. @@ -287,12 +287,12 @@ def determine_point_ownership( Args: mesh: The mesh points: Points to check for collision (``shape=(num_points, gdim)``) - cells: Cells to check for ownership - If ``None`` then all cells are considered. padding: Amount of absolute padding of bounding boxes of the mesh. Each bounding box of the mesh is padded with this amount, to increase the number of candidates, avoiding rounding errors in determining the owner of a point if the point is on the surface of a cell in the mesh. + cells: Cells to check for ownership + If ``None`` then all cells are considered. Returns: Point ownership data @@ -305,7 +305,6 @@ def determine_point_ownership( A large padding value will increase the run-time of the code by orders of magnitude. General advice is to use a padding on the scale of the cell size. - """ if cells is None: map = mesh.topology.index_map(mesh.topology.dim) diff --git a/python/dolfinx/wrappers/geometry.cpp b/python/dolfinx/wrappers/geometry.cpp index beeb6da681b..2419e7e511b 100644 --- a/python/dolfinx/wrappers/geometry.cpp +++ b/python/dolfinx/wrappers/geometry.cpp @@ -187,10 +187,10 @@ void declare_bbtree(nb::module_& m, std::string type) const std::size_t p_s0 = points.ndim() == 1 ? 1 : points.shape(0); std::span _p(points.data(), 3 * p_s0); return dolfinx::geometry::determine_point_ownership(mesh, _p, - std::span(cells.data(), cells.size()), - padding); + padding, + std::span(cells.data(), cells.size())); }, - nb::arg("mesh"), nb::arg("points"), nb::arg("cells"), nb::arg("padding"), + nb::arg("mesh"), nb::arg("points"), nb::arg("padding"), nb::arg("cells"), "Compute point ownership data for mesh-points pair."); std::string pod_pyclass_name = "PointOwnershipData_" + type; diff --git a/python/test/unit/geometry/test_bounding_box_tree.py b/python/test/unit/geometry/test_bounding_box_tree.py index f76d1ccd05f..adad166f0d6 100644 --- a/python/test/unit/geometry/test_bounding_box_tree.py +++ b/python/test/unit/geometry/test_bounding_box_tree.py @@ -24,6 +24,7 @@ ) from dolfinx.mesh import ( CellType, + compute_midpoints, create_box, create_unit_cube, create_unit_interval, @@ -501,56 +502,107 @@ def test_surface_bbtree_collision(dtype): @pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("affine", [True, False]) @pytest.mark.parametrize("dtype", [np.float32, np.float64]) -def test_determine_point_ownership(dim, dtype): +def test_determine_point_ownership(dim, affine, dtype): """Find point owners (ranks and cells) using bounding box trees + global communication and compare to point ownership data results.""" comm = MPI.COMM_WORLD rank = comm.Get_rank() + mpi_dtype = MPI.DOUBLE if dtype == np.float64 else MPI.FLOAT tdim = dim num_cells_side = 4 - h = dtype(1.0 / num_cells_side) if tdim == 2: - mesh = create_unit_square( - MPI.COMM_WORLD, num_cells_side, num_cells_side, CellType.quadrilateral, dtype=dtype - ) + ct = CellType.triangle if affine else CellType.quadrilateral + mesh = create_unit_square(MPI.COMM_WORLD, num_cells_side, num_cells_side, ct, dtype=dtype) else: + ct = CellType.tetrahedron if affine else CellType.hexahedron mesh = create_unit_cube( MPI.COMM_WORLD, num_cells_side, num_cells_side, num_cells_side, - CellType.hexahedron, + ct, dtype=dtype, ) cell_map = mesh.topology.index_map(tdim) tree = bb_tree(mesh, mesh.topology.dim, np.arange(cell_map.size_local)) - num_global_cells = num_cells_side**dim + num_global_cells = num_cells_side**tdim + if affine: + num_global_cells *= 2 * (3 ** (tdim - 2)) + local_midpoints = compute_midpoints( + mesh, tdim, np.arange(mesh.topology.index_map(tdim).size_local) + ) + midpoints_per_rank = np.zeros(comm.size, dtype=np.int32) + midpoints_offsets = np.zeros(comm.size, dtype=np.int32) + comm.Allgather(np.array([local_midpoints.shape[0]], dtype=np.int32), midpoints_per_rank) + midpoints_offsets[1:] = np.cumsum(midpoints_per_rank[:-1]) all_midpoints = np.zeros((num_global_cells, 3), dtype=dtype) - counter = 0 - Xs = [(i + 0.5) * h for i in range(num_cells_side)] - Ys = Xs - Zs = [0.0] if dim == 2 else Xs - for x in Xs: - for y in Ys: - for z in Zs: - all_midpoints[counter, :] = np.array([x, y, z], dtype=dtype) - counter += 1 - # Since ghost cells are left out and the points considered are midpoints - # of cells, they are only contained in a single process / cell - # Moreover, the bounding boxes of the elements correspond to the elements, - # hence the tree collisions are the exact collisions + comm.Allgatherv( + local_midpoints, [all_midpoints, midpoints_per_rank * 3, midpoints_offsets * 3, mpi_dtype] + ) + # Find potential owner cells tree_col = compute_collisions_points(tree, all_midpoints) + + mesh.topology.create_connectivity(tdim - 1, 0) + mesh.topology.create_connectivity(0, tdim) + cfc = mesh.topology.connectivity(tdim, tdim - 1) + fpc = mesh.topology.connectivity(tdim - 1, 0) + + # Narrow it down to a single owner cell + def is_inside(mesh, icell, point): + fdim = tdim - 1 + is_inside = True + cpoints = mesh.geometry.x[mesh.geometry.dofmap[icell, :]] # cell points + ccentroid = np.average(cpoints, axis=0) # cell centroid + for ifacet in cfc.links(icell): + fpoints_indices = _cpp.mesh.entities_to_geometry( + mesh._cpp_object, + 0, + fpc.links(ifacet), + False, + ) + fpoints_indices = fpoints_indices.reshape(fpoints_indices.size) + fpoints = mesh.geometry.x[fpoints_indices] + fcentroid = np.average(fpoints, axis=0) # facet centroid + # Compute facet normal pointing to outside of owner cell + normal = np.zeros(3, dtype=dtype) + facet_vector1 = fpoints[1, :] - fpoints[0, :] + if fdim == 1: + normal[0] = -facet_vector1[1] + normal[1] = +facet_vector1[0] + elif fdim == 2: + facet_vector2 = fpoints[2, :] - fpoints[0, :] + normal = np.cross(facet_vector1, facet_vector2) + else: + raise ValueError("Unexpected facet dimension.") + normal /= np.linalg.norm(normal) + # Re-align if pointing to inside the parent cell + normal = -normal if (np.dot((ccentroid - fcentroid), normal) > 0) else normal + # Test the point + signed_distance = np.dot((point - fcentroid), normal) + if signed_distance > 1e-9: + is_inside = False + break + return is_inside + processwise_owners = np.zeros(2 * num_global_cells, dtype=np.int32) owners = np.empty_like(processwise_owners) for ipoint in range(num_global_cells): - cell = tree_col.links(ipoint) - if len(cell) > 0: + potential_owners = tree_col.links(ipoint) + owner_cells = [] + for cell in potential_owners: + if is_inside(mesh, cell, all_midpoints[ipoint, :]): + owner_cells.append(cell) + if owner_cells: + assert len(owner_cells) == 1 processwise_owners[2 * ipoint] = rank - processwise_owners[2 * ipoint + 1] = cell[0] + processwise_owners[2 * ipoint + 1] = owner_cells[0] + # Since ghost cells are left out and the points considered are midpoints + # of cells, they are only contained in a single process / cell # The value at a given index is null if it doesn't correspond # to the current process. # We can sum the processwise arrays to obtain a global array @@ -591,7 +643,7 @@ def check_po(po: PointOwnershipData, src_points, ownership_data, global_dest_own (iglobal, processor, _) = data if processor == rank: found = False - point = np.array(point, dtype) + point = np.array(point, dtype=dtype) for jpoint in dest_points_indices: found = np.allclose(point, dest_points[jpoint]) if found: @@ -630,26 +682,33 @@ def compute_global_owners(N, start, stop): # All cells points, start, stop = set_local_range(all_midpoints) owners = compute_global_owners(np.int64(all_midpoints.shape[0]), start, stop) - all_cells = np.arange(cell_map.size_local, dtype=np.int32) - po = determine_point_ownership(mesh, points, all_cells) + all_cells = np.arange(cell_map.size_local, dtype=dtype) + po = determine_point_ownership(mesh, points, 0.0, all_cells) check_po(po, points, ownership_data, owners) # Left half - num_left_cells = np.rint((num_cells_side**dim) / 2).astype(np.int32) - left_midpoints = all_midpoints[np.arange(num_left_cells), :] + num_left_cells = np.rint(num_global_cells / 2).astype(np.int32) + left_midpoints = np.zeros((num_left_cells, 3), dtype=dtype) + counter = 0 + indices_left = [] + for ipoint in range(num_global_cells): + if all_midpoints[ipoint, 0] <= 0.5: + left_midpoints[counter] = all_midpoints[ipoint] + indices_left.append(ipoint) + counter += 1 points, start, stop = set_local_range(left_midpoints) owners = compute_global_owners(np.int64(all_midpoints.shape[0]), start, stop) left_cells = locate_entities(mesh, tdim, lambda x: x[0] <= 0.5) left_cells = np.array( [cell for cell in left_cells if cell < cell_map.size_local], dtype=np.int32 ) # Filter out ghost cells - lpo = determine_point_ownership(mesh, points, left_cells) + lpo = determine_point_ownership(mesh, points, 0.0, left_cells) left_ownership_data = {} - for ipoint in range(num_left_cells): + for idx, ipoint in enumerate(indices_left): left_ownership_data[tuple(all_midpoints[ipoint])] = ( - ipoint, + idx, owner_ranks[ipoint], owner_cells[ipoint], ) From 74b9a223fbd53120e8fe7863c4e2afd50c7535d5 Mon Sep 17 00:00:00 2001 From: Mehdi Slimani Date: Wed, 2 Oct 2024 09:07:00 +0000 Subject: [PATCH 7/9] Non-defaulted padding in BoundingBoxTree python constructor --- python/dolfinx/geometry.py | 2 +- python/test/unit/fem/test_function.py | 2 +- python/test/unit/fem/test_interpolation.py | 2 +- .../unit/geometry/test_bounding_box_tree.py | 40 +++++++++---------- python/test/unit/geometry/test_gjk.py | 2 +- .../unit/mesh/test_manifold_point_search.py | 2 +- 6 files changed, 25 insertions(+), 25 deletions(-) diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py index b356b346b4c..bd5282a9f39 100644 --- a/python/dolfinx/geometry.py +++ b/python/dolfinx/geometry.py @@ -103,8 +103,8 @@ def create_global_tree(self, comm) -> BoundingBoxTree: def bb_tree( mesh: Mesh, dim: int, + padding: float, entities: typing.Optional[npt.NDArray[np.int32]] = None, - padding: float = 0.0, ) -> BoundingBoxTree: """Create a bounding box tree for use in collision detection. diff --git a/python/test/unit/fem/test_function.py b/python/test/unit/fem/test_function.py index 550b0528e7a..c131aa47c4f 100644 --- a/python/test/unit/fem/test_function.py +++ b/python/test/unit/fem/test_function.py @@ -90,7 +90,7 @@ def e3(x): u3.interpolate(e3) x0 = (mesh.geometry.x[0] + mesh.geometry.x[1]) / 2.0 - tree = bb_tree(mesh, mesh.geometry.dim) + tree = bb_tree(mesh, mesh.geometry.dim, 0.0) cell_candidates = compute_collisions_points(tree, x0) cell = compute_colliding_cells(mesh, cell_candidates, x0).array assert len(cell) > 0 diff --git a/python/test/unit/fem/test_interpolation.py b/python/test/unit/fem/test_interpolation.py index 126e265fc88..994888b1a25 100644 --- a/python/test/unit/fem/test_interpolation.py +++ b/python/test/unit/fem/test_interpolation.py @@ -1025,7 +1025,7 @@ def f_test2(x): u1_exact.x.scatter_forward() # Find the single cell in mesh1 which is overlapped by mesh2 - tree1 = bb_tree(mesh1, mesh1.topology.dim) + tree1 = bb_tree(mesh1, mesh1.topology.dim, 0.0) cells_overlapped1 = compute_collisions_points( tree1, np.array([p0_mesh2, p0_mesh2, 0.0]) / 2 ).array diff --git a/python/test/unit/geometry/test_bounding_box_tree.py b/python/test/unit/geometry/test_bounding_box_tree.py index adad166f0d6..be485b60ff0 100644 --- a/python/test/unit/geometry/test_bounding_box_tree.py +++ b/python/test/unit/geometry/test_bounding_box_tree.py @@ -150,7 +150,7 @@ def rotation_matrix(axis, angle): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_empty_tree(dtype): mesh = create_unit_interval(MPI.COMM_WORLD, 16, dtype=dtype) - bbtree = bb_tree(mesh, mesh.topology.dim, np.array([], dtype=dtype)) + bbtree = bb_tree(mesh, mesh.topology.dim, 0.0, np.array([], dtype=dtype)) assert bbtree.num_bboxes == 0 @@ -167,7 +167,7 @@ def test_compute_collisions_point_1d(dtype): # Compute collision tdim = mesh.topology.dim - tree = bb_tree(mesh, tdim) + tree = bb_tree(mesh, tdim, 0.0) entities = compute_collisions_points(tree, p) assert len(entities.array) == 1 @@ -212,8 +212,8 @@ def locator_B(x): cells_B = np.sort(np.unique(np.hstack([v_to_c.links(vertex) for vertex in vertices_B]))) # Find colliding entities using bounding box trees - tree_A = bb_tree(mesh_A, mesh_A.topology.dim) - tree_B = bb_tree(mesh_B, mesh_B.topology.dim) + tree_A = bb_tree(mesh_A, mesh_A.topology.dim, 0.0) + tree_B = bb_tree(mesh_B, mesh_B.topology.dim, 0.0) entities = compute_collisions_trees(tree_A, tree_B) entities_A = np.sort(np.unique([q[0] for q in entities])) entities_B = np.sort(np.unique([q[1] for q in entities])) @@ -229,8 +229,8 @@ def test_compute_collisions_tree_2d(point, dtype): mesh_B = create_unit_square(MPI.COMM_WORLD, 5, 5, dtype=dtype) bgeom = mesh_B.geometry.x bgeom += point - tree_A = bb_tree(mesh_A, mesh_A.topology.dim) - tree_B = bb_tree(mesh_B, mesh_B.topology.dim) + tree_A = bb_tree(mesh_A, mesh_A.topology.dim, 0.0) + tree_B = bb_tree(mesh_B, mesh_B.topology.dim, 0.0) entities = compute_collisions_trees(tree_A, tree_B) entities_A = np.sort(np.unique([q[0] for q in entities])) @@ -251,8 +251,8 @@ def test_compute_collisions_tree_3d(point, dtype): bgeom = mesh_B.geometry.x bgeom += point - tree_A = bb_tree(mesh_A, mesh_A.topology.dim) - tree_B = bb_tree(mesh_B, mesh_B.topology.dim) + tree_A = bb_tree(mesh_A, mesh_A.topology.dim, 0.0) + tree_B = bb_tree(mesh_B, mesh_B.topology.dim, 0.0) entities = compute_collisions_trees(tree_A, tree_B) entities_A = np.sort(np.unique([q[0] for q in entities])) entities_B = np.sort(np.unique([q[1] for q in entities])) @@ -269,7 +269,7 @@ def test_compute_closest_entity_1d(dim, dtype): N = 16 points = np.array([[-ref_distance, 0, 0], [2 / N, 2 * ref_distance, 0]], dtype=dtype) mesh = create_unit_interval(MPI.COMM_WORLD, N, dtype=dtype) - tree = bb_tree(mesh, dim) + tree = bb_tree(mesh, dim, 0.0) num_entities_local = ( mesh.topology.index_map(dim).size_local + mesh.topology.index_map(dim).num_ghosts ) @@ -303,7 +303,7 @@ def test_compute_closest_entity_2d(dim, dtype): points = np.array([-1.0, -0.01, 0.0], dtype=dtype) mesh = create_unit_square(MPI.COMM_WORLD, 15, 15, dtype=dtype) mesh.topology.create_entities(dim) - tree = bb_tree(mesh, dim) + tree = bb_tree(mesh, dim, 0.0) num_entities_local = ( mesh.topology.index_map(dim).size_local + mesh.topology.index_map(dim).num_ghosts ) @@ -335,7 +335,7 @@ def test_compute_closest_entity_3d(dim, dtype): mesh = create_unit_cube(MPI.COMM_WORLD, 8, 8, 8, dtype=dtype) mesh.topology.create_entities(dim) - tree = bb_tree(mesh, dim) + tree = bb_tree(mesh, dim, 0.0) num_entities_local = ( mesh.topology.index_map(dim).size_local + mesh.topology.index_map(dim).num_ghosts ) @@ -368,7 +368,7 @@ def test_compute_closest_sub_entity(dim, dtype): mesh = create_unit_cube(MPI.COMM_WORLD, 8, 8, 8, dtype=dtype) mesh.topology.create_entities(dim) left_entities = locate_entities(mesh, dim, lambda x: x[0] <= xc) - tree = bb_tree(mesh, dim, left_entities) + tree = bb_tree(mesh, dim, 0.0, left_entities) midpoint_tree = create_midpoint_tree(mesh, dim, left_entities) closest_entities = compute_closest_entity(tree, midpoint_tree, mesh, points) @@ -396,7 +396,7 @@ def test_surface_bbtree(dtype): tdim = mesh.topology.dim f_to_c = mesh.topology.connectivity(tdim - 1, tdim) cells = np.array([f_to_c.links(f)[0] for f in sf], dtype=np.int32) - bbtree = bb_tree(mesh, tdim, cells) + bbtree = bb_tree(mesh, tdim, 0.0, cells) # test collision (should not collide with any) p = np.array([0.5, 0.5, 0.5]) @@ -413,7 +413,7 @@ def test_sub_bbtree_codim1(dtype): top_facets = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[2], 1)) f_to_c = mesh.topology.connectivity(tdim - 1, tdim) cells = np.array([f_to_c.links(f)[0] for f in top_facets], dtype=np.int32) - bbtree = bb_tree(mesh, tdim, cells) + bbtree = bb_tree(mesh, tdim, 0.0, cells) # Compute a BBtree for all processes process_bbtree = bbtree.create_global_tree(mesh.comm) @@ -441,7 +441,7 @@ def test_serial_global_bb_tree(dtype, comm): # entity tree with a serial mesh x = np.array([[2.0, 2.0, 3.0], [0.3, 0.2, 0.1]], dtype=dtype) - tree = bb_tree(mesh, mesh.topology.dim) + tree = bb_tree(mesh, mesh.topology.dim, 0.0) global_tree = tree.create_global_tree(mesh.comm) tree_col = compute_collisions_points(tree, x) @@ -465,12 +465,12 @@ def test_sub_bbtree_box(ct, N, dtype): facets = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[1], 1.0)) f_to_c = mesh.topology.connectivity(fdim, tdim) cells = np.int32(np.unique([f_to_c.links(f)[0] for f in facets])) - bbtree = bb_tree(mesh, tdim, cells) + bbtree = bb_tree(mesh, tdim, 0.0, cells) num_boxes = bbtree.num_bboxes if num_boxes > 0: bbox = bbtree.get_bbox(num_boxes - 1) assert np.isclose(bbox[0][1], (N - 1) / N) - tree = bb_tree(mesh, tdim) + tree = bb_tree(mesh, tdim, 0.0) assert num_boxes < tree.num_bboxes @@ -489,13 +489,13 @@ def test_surface_bbtree_collision(dtype): # Compute unique set of cells (some will be counted multiple times) cells = np.array(list(set([f_to_c.links(f)[0] for f in sf])), dtype=np.int32) - bbtree1 = bb_tree(mesh1, tdim, cells) + bbtree1 = bb_tree(mesh1, tdim, 0.0, cells) mesh2.topology.create_connectivity(mesh2.topology.dim - 1, mesh2.topology.dim) sf = exterior_facet_indices(mesh2.topology) f_to_c = mesh2.topology.connectivity(tdim - 1, tdim) cells = np.array(list(set([f_to_c.links(f)[0] for f in sf])), dtype=np.int32) - bbtree2 = bb_tree(mesh2, tdim, cells) + bbtree2 = bb_tree(mesh2, tdim, 0.0, cells) collisions = compute_collisions_trees(bbtree1, bbtree2) assert len(collisions) == 1 @@ -528,7 +528,7 @@ def test_determine_point_ownership(dim, affine, dtype): ) cell_map = mesh.topology.index_map(tdim) - tree = bb_tree(mesh, mesh.topology.dim, np.arange(cell_map.size_local)) + tree = bb_tree(mesh, mesh.topology.dim, 0.0, np.arange(cell_map.size_local)) num_global_cells = num_cells_side**tdim if affine: num_global_cells *= 2 * (3 ** (tdim - 2)) diff --git a/python/test/unit/geometry/test_gjk.py b/python/test/unit/geometry/test_gjk.py index 8b895a4633f..fb2dc706343 100644 --- a/python/test/unit/geometry/test_gjk.py +++ b/python/test/unit/geometry/test_gjk.py @@ -193,7 +193,7 @@ def test_collision_2nd_order_triangle(dtype): sample_points = np.array([[0.1, 0.3, 0.0], [0.2, 0.5, 0.0], [0.6, 0.6, 0.0]]) # Create boundingboxtree - tree = geometry.bb_tree(mesh, mesh.geometry.dim) + tree = geometry.bb_tree(mesh, mesh.geometry.dim, 0.0) cell_candidates = geometry.compute_collisions_points(tree, sample_points) colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, sample_points) # Check for collision diff --git a/python/test/unit/mesh/test_manifold_point_search.py b/python/test/unit/mesh/test_manifold_point_search.py index ec3428ca3f9..80704d2e0bc 100644 --- a/python/test/unit/mesh/test_manifold_point_search.py +++ b/python/test/unit/mesh/test_manifold_point_search.py @@ -18,7 +18,7 @@ def test_manifold_point_search(): cells = np.array([[0, 1, 2], [0, 1, 3]], dtype=np.int64) domain = ufl.Mesh(element("Lagrange", "triangle", 1, shape=(2,))) mesh = create_mesh(MPI.COMM_WORLD, cells, vertices, domain) - bb = bb_tree(mesh, mesh.topology.dim) + bb = bb_tree(mesh, mesh.topology.dim, 0.0) # Find cell colliding with point points = np.array([[0.5, 0.25, 0.75], [0.25, 0.5, 0.75]], dtype=default_real_type) From a96927116c09a3f94ac2255f7ea87872601387cd Mon Sep 17 00:00:00 2001 From: Mehdi Slimani Date: Wed, 2 Oct 2024 09:00:45 +0000 Subject: [PATCH 8/9] Non-defaulted padding in BoundingBoxTree c++ constructor - Changed argument order to match Python constructor --- cpp/dolfinx/geometry/BoundingBoxTree.h | 4 ++-- cpp/dolfinx/geometry/utils.h | 2 +- python/dolfinx/geometry.py | 4 ++-- python/dolfinx/wrappers/geometry.cpp | 11 +++++------ 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/cpp/dolfinx/geometry/BoundingBoxTree.h b/cpp/dolfinx/geometry/BoundingBoxTree.h index 64ede45a057..01e0db723f3 100644 --- a/cpp/dolfinx/geometry/BoundingBoxTree.h +++ b/cpp/dolfinx/geometry/BoundingBoxTree.h @@ -224,7 +224,7 @@ class BoundingBoxTree /// @param[in] padding Value to pad (extend) the the bounding box of /// each entity by. BoundingBoxTree(const mesh::Mesh& mesh, int tdim, - std::span entities, double padding = 0) + double padding, std::span entities) : _tdim(tdim) { if (tdim < 0 or tdim > mesh.topology()->dim()) @@ -266,7 +266,7 @@ class BoundingBoxTree /// build the bounding box tree for /// @param[in] padding Value to pad (extend) the the bounding box of /// each entity by. - BoundingBoxTree(const mesh::Mesh& mesh, int tdim, T padding = 0) + BoundingBoxTree(const mesh::Mesh& mesh, int tdim, T padding) : BoundingBoxTree::BoundingBoxTree( mesh, tdim, range(mesh.topology_mutable(), tdim), padding) { diff --git a/cpp/dolfinx/geometry/utils.h b/cpp/dolfinx/geometry/utils.h index 0c772937f02..c629c3ba598 100644 --- a/cpp/dolfinx/geometry/utils.h +++ b/cpp/dolfinx/geometry/utils.h @@ -697,7 +697,7 @@ PointOwnershipData determine_point_ownership(const mesh::Mesh& mesh, } // Create a global bounding-box tree to find candidate processes with // cells that could collide with the points - BoundingBoxTree bb(mesh, tdim, cells, padding); + BoundingBoxTree bb(mesh, tdim, padding, cells); BoundingBoxTree global_bbtree = bb.create_global_tree(comm); // Compute collisions: diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py index bd5282a9f39..b086dd55204 100644 --- a/python/dolfinx/geometry.py +++ b/python/dolfinx/geometry.py @@ -128,11 +128,11 @@ def bb_tree( dtype = mesh.geometry.x.dtype if np.issubdtype(dtype, np.float32): return BoundingBoxTree( - _cpp.geometry.BoundingBoxTree_float32(mesh._cpp_object, dim, entities, padding) + _cpp.geometry.BoundingBoxTree_float32(mesh._cpp_object, dim, padding, entities) ) elif np.issubdtype(dtype, np.float64): return BoundingBoxTree( - _cpp.geometry.BoundingBoxTree_float64(mesh._cpp_object, dim, entities, padding) + _cpp.geometry.BoundingBoxTree_float64(mesh._cpp_object, dim, padding, entities) ) else: raise NotImplementedError(f"Type {dtype} not supported.") diff --git a/python/dolfinx/wrappers/geometry.cpp b/python/dolfinx/wrappers/geometry.cpp index 2419e7e511b..cc81cd3d208 100644 --- a/python/dolfinx/wrappers/geometry.cpp +++ b/python/dolfinx/wrappers/geometry.cpp @@ -33,17 +33,16 @@ void declare_bbtree(nb::module_& m, std::string type) "__init__", [](dolfinx::geometry::BoundingBoxTree* bbt, const dolfinx::mesh::Mesh& mesh, int dim, + double padding, nb::ndarray, nb::c_contig> - entities, - double padding) + entities) { new (bbt) dolfinx::geometry::BoundingBoxTree( mesh, dim, - std::span(entities.data(), entities.size()), - padding); + padding, + std::span(entities.data(), entities.size())); }, - nb::arg("mesh"), nb::arg("dim"), nb::arg("entities"), - nb::arg("padding") = 0.0) + nb::arg("mesh"), nb::arg("dim"), nb::arg("padding"), nb::arg("entities")) .def_prop_ro("num_bboxes", &dolfinx::geometry::BoundingBoxTree::num_bboxes) .def( From 3e481f0e1c708bb0d2728be45719af315595ea48 Mon Sep 17 00:00:00 2001 From: Mehdi Slimani Date: Wed, 2 Oct 2024 09:33:59 +0000 Subject: [PATCH 9/9] added missing padding argument in python demo --- python/demo/demo_static-condensation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py index 93c7d29d690..d415c7aafed 100644 --- a/python/demo/demo_static-condensation.py +++ b/python/demo/demo_static-condensation.py @@ -183,7 +183,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL): A.assemble() # Create bounding box for function evaluation -bb_tree = geometry.bb_tree(msh, 2) +bb_tree = geometry.bb_tree(msh, 2, 0.0) # Check against standard table value p = np.array([[48.0, 52.0, 0.0]], dtype=np.float64)