Improve python exposure of determine_point_ownership
#3344
Conversation
I added a test comparing bounding box tree results to determine_point_ownership. Another possibility would be to look at post files and hardcode rank-cell data into a test. |
| def test_po_against_bb_tree(mesh,points,cells): | ||
| cells = np.array([cell for cell in cells if cell < cell_map.size_local]) | ||
| num_points = points.shape[0] | ||
| po = determine_point_ownership(mesh,points,cells) | ||
| tree = bb_tree(mesh, mesh.topology.dim, cells) |
There was a problem hiding this comment.
Run ruff format and ruff check --fix on the Python files to assure appropriate formatting.
|
Sorry, I made a mess with my branches in my fork and closed the pull request by mistake. Is it possible to reopen it pointing to https://github.com/ordinary-slim/dolfinx/tree/main ? |
I think you have to reopen it from your repository. See: https://gist.github.com/Farhan-Haseeb/fd4cbb476de4e2931c322d7bacb75e5f |
I can't see the PR from my repository. I deleted the original branch of the PR and it automatically closed it. The reopen button is greyed out: The link you provided says:
But I can't find this |
|
I cannot reopen. Please open a new PR:) |
|
Ok, I had to recover the deleted branch with the same name then I could reopen, as per this stackoverflow comment. Sorry for the inconvenience! |
| 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 |
There was a problem hiding this comment.
Can it be. parameterised across the mesh? (triangle, quad, tet, hex)?
There was a problem hiding this comment.
No, as the test is now it needs the elements to coincide with their bounding boxes. I'll commit a version that covers triangles and tet elements.
python/dolfinx/geometry.py
Outdated
| 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. | ||
|
|
cpp/dolfinx/geometry/utils.h
Outdated
| /// 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`, |
There was a problem hiding this comment.
I think this documentation should be moved out to the struct itself.
|
FYI: The geometry test |
`test_determine_point_ownership` tests for triangle and tetra elements.
| local_cells.resize(cell_map->size_local()); | ||
| std::iota(local_cells.begin(), local_cells.end(), 0); | ||
| cells = std::span<const std::int32_t>(local_cells.data(), local_cells.size()); | ||
| } |
There was a problem hiding this comment.
@jhale I made cells a default argument at the c++ level (not at the wrapper level).
Ideally this local_cells variable would be initialized inside of the empty conditional but I could not figure it out since the span that wraps it dangles if the underlying vector goes out of scope. I tried using the static keyword but didn't work.
- Changed argument order to match Python constructor
|
Let me know if I should move the last two commits (non-defaulting padding argument in BoundingBoxTree constructors) to a new PR. |
|
Bumping this @jhale |
| @@ -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<T>& mesh, int tdim, T padding = 0) | |||
There was a problem hiding this comment.
Could you switch this literal 0 to 0.0 so it's 100% clear T is a floating point type.
|
I'd like to get this merged - can you address the two final small comments? |
| /// @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 | ||
| /// intersecting bounding box, which is an expensive operation to perform. | ||
| template <std::floating_point T> | ||
| PointOwnershipData<T> determine_point_ownership(const mesh::Mesh<T>& mesh, | ||
| std::span<const T> points, | ||
| T padding) | ||
| T padding, | ||
| std::span<const std::int32_t> cells = {}) |
There was a problem hiding this comment.
Use std::optional rather than an empty span.
| between the bounding boxes of the cells and the set of points. | ||
| Then, actual containment pairs are determined using the GJK algorithm. | ||
|
|
||
| Args: |
There was a problem hiding this comment.
Formatting needs fixing - second and subsequent lines of docstring for an argument must be indented.
| of magnitude. General advice is to use a padding on the scale of the | ||
| cell size. | ||
| """ | ||
| if cells is None: |
There was a problem hiding this comment.
If you switch the C++ interface to std::optional, this if could be removed.
| @@ -182,13 +181,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<T>& mesh, | |||
| nb::ndarray<const T, nb::c_contig> points, const T padding) | |||
| nb::ndarray<const T, nb::c_contig> points, | |||
| nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig> cells, | |||
There was a problem hiding this comment.
Use std::optional for cells.
Python wrapper around nanobind exposure of
determine_point_ownershipand extra optional argumentscellsso that the search is performed only forcells.