Geometry (dolfinx::geometry)#

namespace geometry#

Geometry data structures and algorithms.

Tools for geometric data structures and operations, e.g. searching.

Functions

template<std::floating_point T, typename U = boost::multiprecision::cpp_bin_float_double_extended>
std::array<T, 3> compute_distance_gjk(std::span<const T> p0, std::span<const T> q0)#

Compute the distance between two convex bodies p0 and q0, each defined by a set of points.

Uses the Gilbert–Johnson–Keerthi (GJK) distance algorithm.

Parameters:

  • p0 – [in] Body 1 list of points, shape=(num_points, 3). Row-major storage.

  • q0 – [in] Body 2 list of points, shape=(num_points, 3). Row-major storage.

Template Parameters:

  • T – Floating point type

  • U – Floating point type used for geometry computations internally, which should be higher precision than T, to maintain accuracy.

Returns:

shortest vector between bodies

template<std::floating_point T, typename U = boost::multiprecision::cpp_bin_float_double_extended>
std::vector<T> compute_distances_gjk(const std::vector<std::span<const T>> &bodies, std::span<const T> q, int num_threads)#

Compute the distance between a sequence of convex bodies p0, ..., pN and q, each defined by a set of points.

Uses the Gilbert–Johnson–Keerthi (GJK) distance algorithm.

Parameters:

  • bodies – [in] List of the list of points that make up each of N bodies considered as body 1. shape=(num_bodies, (num_points_body_j, 3). Row-major storage.

  • q – [in] Body 2 list of points, shape=(num_points, 3). Row-major storage.

  • num_threads – [in] Number of threads to use. Must be >= 1.

Template Parameters:

  • T – Floating point type

  • U – Floating point type used for geometry computations internally, which should be higher precision than T, to maintain accuracy.

Returns:

For each body in p_j, return the shortest distance vector to body 2. Shape (num_points, 3).

template<std::floating_point T>
std::vector<T> shortest_vector(const mesh::Mesh<T> &mesh, int dim, std::span<const std::int32_t> entities, std::span<const T> points)#

Compute the shortest vector from a mesh entity to a point.

Parameters:

  • mesh – [in] The mesh

  • dim – [in] Topological dimension of the mesh entity

  • entities – [in] List of entities

  • points – [in] Set of points (shape=(num_points, 3)), using row-major storage.

Returns:

An array of vectors (shape=(num_points, 3)) where the ith row is the shortest vector between the ith entity and the ith point. Storage is row-major.

template<std::floating_point T>
T compute_squared_distance_bbox(std::span<const T, 6> b, std::span<const T, 3> x)#

Compute squared distance between point and bounding box.

Parameters:

  • b – [in] Bounding box coordinates

  • x – [in] A point

Returns:

The shortest distance between the bounding box b and the point x. Returns zero if x is inside box.

template<std::floating_point T>
std::vector<T> squared_distance(const mesh::Mesh<T> &mesh, int dim, std::span<const std::int32_t> entities, std::span<const T> points)#

Compute the squared distance between a point and a mesh entity.

The distance is computed between the ith input points and the ith input entity.

Note

Uses the GJK algorithm, see geometry::compute_distance_gjk for details.

Note

Uses a convex hull approximation of linearized geometry

Parameters:

  • mesh – [in] Mesh containing the entities

  • dim – [in] The topological dimension of the mesh entities

  • entities – [in] The indices of the mesh entities (local to process)

  • points – [in] The set points from which to computed the shortest (shape=(num_points, 3)). Storage is row-major.

Returns:

Squared shortest distance from points[i] to entities[i]

template<std::floating_point T>
BoundingBoxTree<T> create_midpoint_tree(const mesh::Mesh<T> &mesh, int tdim, std::span<const std::int32_t> entities)#

Create a bounding box tree for the midpoints of a subset of entities.

Parameters:

  • mesh – [in] The mesh

  • tdim – [in] The topological dimension of the entity

  • entities – [in] List of local entity indices

Returns:

Bounding box tree for midpoints of entities

template<std::floating_point T>
std::vector<std::int32_t> compute_collisions(const BoundingBoxTree<T> &tree0, const BoundingBoxTree<T> &tree1)#

Compute all collisions between two bounding box trees.

Parameters:
Returns:

List of pairs of intersecting box indices from each tree, flattened as a vector of size num_intersections*2

template<std::floating_point T>
graph::AdjacencyList<std::int32_t> compute_collisions(const BoundingBoxTree<T> &tree, std::span<const T> points)#

Compute collisions between points and leaf bounding boxes.

Bounding boxes can overlap, therefore points can collide with more than one box.

Parameters:

  • tree – [in] The bounding box tree

  • points – [in] The points (shape=(num_points, 3)). Storage is row-major.

Returns:

For each point, the bounding box leaves that collide with the point.

template<std::floating_point T>
std::int32_t compute_first_colliding_cell(const mesh::Mesh<T> &mesh, std::span<const std::int32_t> cells, std::array<T, 3> point, T tol, std::span<T> coordinate_dofs)#

Given a set of cells, find the first one that collides with a point.

A point can collide with more than one cell. The first cell detected to collide with the point is returned. If no collision is detected, -1 is returned.

Note

cells can for instance be found by using geometry::compute_collisions between a bounding box tree for the cells of the mesh and the point.

Parameters:

  • mesh – [in] The mesh.

  • cells – [in] Candidate cells.

  • point – [in] The point (shape=(3,)).

  • tol – [in] Tolerance for accepting a collision (in the squared distance).

  • coordinate_dofs – [inout] Scratch buffer, sized at least (number of nodes per cell) * 3, i.e. mesh.geometry().dofmaps().front().extent(1) * 3. A larger buffer is permitted (e.g. sized for the largest cell type in a mixed-topology mesh); only the first (number of nodes per cell) * 3 entries are used. Callers invoking this function once per point in a loop should hoist a single buffer outside the loop and pass it here to avoid a per-call allocation. Its contents on return are unspecified.

Returns:

Local cell index, -1 if not found.

template<std::floating_point T>
std::vector<std::int32_t> compute_closest_entity(const BoundingBoxTree<T> &tree, const BoundingBoxTree<T> &midpoint_tree, const mesh::Mesh<T> &mesh, std::span<const T> points)#

Compute closest mesh entity to a point.

Note

Returns a vector filled with index -1 if the bounding box tree is empty.

Parameters:

  • tree – [in] The bounding box tree for the entities

  • midpoint_tree – [in] A bounding box tree with the midpoints of all the mesh entities. This is used to accelerate the search.

  • mesh – [in] The mesh

  • points – [in] The set of points (shape=(num_points, 3)). Storage is row-major.

Returns:

For each point, the index of the closest mesh entity.

template<std::floating_point T>
graph::AdjacencyList<std::int32_t> compute_colliding_cells(const mesh::Mesh<T> &mesh, const graph::AdjacencyList<std::int32_t> &candidate_cells, std::span<const T> points)#

Compute which cells collide with a point.

Note

Uses the GJK algorithm, see geometry::compute_distance_gjk for details.

Note

candidate_cells can for instance be found by using geometry::compute_collisions between a bounding box tree and the set of points.

Parameters:

  • mesh – [in] The mesh

  • candidate_cells – [in] List of candidate colliding cells for the ith point in points

  • points – [in] Points to check for collision (shape=(num_points, 3)). Storage is row-major.

Returns:

For each point, the cells that collide with the point.

template<std::floating_point T>
PointOwnershipData<T> determine_point_ownership(const mesh::Mesh<T> &mesh, std::span<const T> points, T padding, std::optional<std::span<const std::int32_t>> cells, bool find_closest_cell = true)#

Determine, for a set of points, the owning process of the cell (if any) that contains each point.

A cell is a candidate for a point if the cell’s bounding box, padded by padding, contains the point. Each candidate is then tested for actual containment of the point with the GJK algorithm. If no candidate actually contains a point, the point is either left unowned or, if find_closest_cell is true, assigned to the candidate cell closest to it (by GJK distance).

Note

dest_owner is sorted

Note

An entry of src_owner is -1 if the corresponding point was not contained in any candidate cell and, if find_closest_cell is true, had no candidate cell to fall back on either (e.g. because padding was too small).

Note

dest_points is flattened row-major, shape (dest_owner.size(), 3)

Note

With find_closest_cell set to true, a large padding value can increase the runtime of the function by orders of magnitude, since a point not contained in any candidate then requires a GJK distance computation against every candidate cell across all processes with an intersecting bounding box.

Note

find_closest_cell must be the same on every rank of mesh.comm(): it gates collective MPI calls, so ranks disagreeing on its value will deadlock.

Parameters:

  • mesh – [in] The mesh

  • points – [in] Points to check for collision (shape=(num_points, 3)). Storage is row-major.

  • padding – [in] Amount of absolute padding applied to each cell’s bounding box before searching for candidate cells/processes. Increasing padding increases the number of cells considered as candidates for a point; it does not by itself decide whether a point with no actually-containing cell is assigned an owner, which is controlled by find_closest_cell.

  • cells – [in] Cells to check for ownership

  • find_closest_cell – [in] If true (default), a point not actually contained in any candidate cell is instead assigned to the process owning the candidate cell closest to it. If false, such a point is left unowned.

Returns:

Point ownership data.

template<std::floating_point T>
class BoundingBoxTree#
#include <BoundingBoxTree.h>

Axis-Aligned bounding box binary tree. It is used to find entities in a collection (often a mesh::Mesh).

template<std::floating_point T>
struct PointOwnershipData#
#include <utils.h>

Information on the ownership of points distributed across processes.

Template Parameters:

T – Mesh geometry floating type.

namespace impl_bb#

Functions

template<std::floating_point T>
std::array<T, 6> compute_bbox_of_entity(const mesh::Mesh<T> &mesh, int dim, std::int32_t index)#
template<std::floating_point T>
std::array<T, 6> compute_bbox_of_bboxes(std::span<const std::pair<std::array<T, 6>, std::int32_t>> leaf_bboxes)#
template<std::floating_point T>
std::pair<std::vector<std::int32_t>, std::vector<T>> build_from_leaf(std::vector<std::pair<std::array<T, 6>, std::int32_t>> &leaf_bboxes)#
namespace impl_gjk#

Functions

template<typename T>
inline std::array<T, 4> det4(const std::array<T, 12> &s)#

Compute the four determinants, of 3x3 matrices given by all the combinations of four vectors (0, 1, 2), (0, 1, 3), (0, 2, 3) and (1, 2, 3). Equivalent to the determinant cofactors of an augmented 4x4 matrix of the same vertices.

Parameters:

s – Flattened array of four 3D vertices (row-major) 4x3.

Returns:

Cofactors/determinants for each set of three vertices

template<typename Vec>
inline Vec::value_type dot3(const Vec &a, const Vec &b)#

Dot product of vectors a and b, both size 3.

Parameters:

  • a – Vector of size 3

  • b – Vector of size 3

Returns:

a.b

template<typename T, std::size_t simplex_size>
void nearest_simplex(const std::array<T, 12> &s, std::array<T, 4> &coordinates)#

Find the barycentric coordinates in the simplex s, of the point in s which is closest to the origin.

Note

s may be an interval, a triangle or a tetrahedron.

Template Parameters:

  • T – The scalar type of the coordinates.

  • simplex_size – The number of points in the simplex (2,3 or 4)

Parameters:

  • s – Simplex described by a set of points in 3D, row-major, flattened.

  • coordinates – [inout] Barycentric coordinates of the point in s closest to the origin.

template<typename T>
inline int support(std::span<const T> bd, const std::array<T, 3> &v)#

‘support’ function, finds point p in bd which maximises p.v

Parameters:

  • bd – Body described by set of 3D points, flattened

  • v – A point in 3D

Returns:

Index of point p in bd which maximises p.v

namespace impl#

Functions

bool is_leaf(std::array<int, 2> bbox)#

Check whether bounding box is a leaf node.

template<std::floating_point T>
bool point_in_bbox(std::span<const T, 6> b, std::span<const T, 3> x)#

A point x is inside a bounding box b if each component of its coordinates lies within the range [b(0,i), b(1,i)] that defines the bounds of the bounding box, b(0,i) <= x[i] <= b(1,i) for i = 0, 1, 2

template<std::floating_point T>
bool bbox_in_bbox(std::span<const T, 6> a, std::span<const T, 6> b)#

A bounding box “a” is contained inside another bounding box “b”, if each of its intervals [a(0,i), a(1,i)] is contained in [b(0,i), b(1,i)], a(0,i) <= b(1, i) and a(1,i) >= b(0, i)