Skip to content

Commit

Permalink
Merge pull request #1212 from LLNL/bugfix/whitlock/kdtree_pointquery_fix
Browse files Browse the repository at this point in the history 
Fix kdtree extents to better support queries of planar arrangements of points in 3D.
  • Loading branch information
@BradWhitlock
BradWhitlock committed Dec 13, 2023
2 parents 11bcdf6 + 5a01679 commit c8ea7af
Show file tree
Hide file tree
Showing 3 changed files with 96 additions and 22 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ and this project aspires to adhere to [Semantic Versioning](https://semver.org/s
- The `conduit::blueprint::mesh::partition()` function no longer issues an error when it receives a "maxshare" adjset.
- The partitioner is better about outputting a "material_map" node for matsets. The "material_map" node is optional for some varieties of matset but they can also help the `conduit::blueprint::mesh::matset::to_silo()` function generate the right material numbers when a domain does not contain all materials.
- The `conduit::Node::swap()` and `conduit::Node::move()` functions no longer cause node names to disappear.
- The `conduit::blueprint::mesh::utils::kdtree` could erroneously return that points were not found when one of the coordset dimensions had a very narrow range of values. This could happen with planar 2D geometries embedded in 3D, such as inside a `MatchQuery` during adjacency set creation.

## [0.8.8] - Released 2023-05-18

Expand Down
12 changes: 9 additions & 3 deletions src/libs/blueprint/conduit_blueprint_mesh_kdtree.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ class kdtree
public:
static const conduit::index_t NoChild;
static const conduit::index_t NotFound;
static const CoordinateType DEFAULT_POINT_TOLERANCE;

using IndexableType = Indexable;
using BoxType = CoordinateType[NDIMS][2];
Expand Down Expand Up @@ -181,6 +182,9 @@ const conduit::index_t kdtree<Indexable, CoordinateType, NDIMS>::NoChild = -1;
template <typename Indexable, typename CoordinateType, int NDIMS>
const conduit::index_t kdtree<Indexable, CoordinateType, NDIMS>::NotFound = -1;

template <typename Indexable, typename CoordinateType, int NDIMS>
const CoordinateType kdtree<Indexable, CoordinateType, NDIMS>::DEFAULT_POINT_TOLERANCE = static_cast<CoordinateType>(1.e-9);

//---------------------------------------------------------------------------
template <typename Indexable, typename CoordinateType, int NDIMS>
kdtree<Indexable, CoordinateType, NDIMS>::kdtree() : boxes(), index()
Expand All @@ -192,7 +196,6 @@ kdtree<Indexable, CoordinateType, NDIMS>::kdtree() : boxes(), index()
coords[i] = Indexable{};
}
coordlen = 0;
constexpr auto DEFAULT_POINT_TOLERANCE = static_cast<CoordinateType>(1.e-9);
setPointTolerance(DEFAULT_POINT_TOLERANCE);
}

Expand Down Expand Up @@ -458,10 +461,13 @@ void kdtree<Indexable, CoordinateType, NDIMS>::calculateExtents()
}
}

// Expand the box a little
// Expand the box a little. Note that we have a minimum expansion to account
// for when a dimension has all the same values.
const CoordinateType minExpansion = DEFAULT_POINT_TOLERANCE * 200.;
for(int i = 0; i < dims(); i++)
{
CoordinateType d = (box[i][1] - box[i][0]) / static_cast<CoordinateType>(200.);
CoordinateType side = std::max(box[i][1] - box[i][0], minExpansion);
CoordinateType d = side / static_cast<CoordinateType>(200.);
box[i][0] -= d;
box[i][1] += d;
}
Expand Down
105 changes: 86 additions & 19 deletions src/tests/blueprint/t_blueprint_mesh_query.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "conduit_log.hpp"

#include <cmath>
#include <cstdlib>
#include <set>
#include <vector>
#include <string>
Expand Down Expand Up @@ -279,26 +280,55 @@ TEST(conduit_blueprint_mesh_query, point_query)
}
}
}
//-----------------------------------------------------------------------------
template <typename Precision>
Precision rand_float(Precision scale = Precision{1})
{
Precision num01 = static_cast<Precision>(rand()) / static_cast<Precision>(RAND_MAX);
Precision num02 = Precision{2} * num01;
Precision num_neg1_pos1 = num02 - Precision{1};
Precision retval = scale * num_neg1_pos1;
return retval;
}

//---------------------------------------------------------------------------
template <typename T>
void make_coords_3d(T **coords, int dims[3])
void make_coords_3d(T **coords, int dims[3], int fuzz = -1)
{
// Initialization.
int npts = dims[0] * dims[1] * dims[2];
for(int dim = 0; dim < 3; dim++)
coords[dim] = new T[npts];

// Make coordinates.
int idx = 0;
for(int k = 0; k < dims[2]; k++)
for(int j = 0; j < dims[1]; j++)
for(int i = 0; i < dims[0]; i++, idx++)
if(fuzz < 0)
{
coords[0][idx] = static_cast<T>(i);
coords[1][idx] = static_cast<T>(j);
coords[2][idx] = static_cast<T>(k);
}
int idx = 0;
for(int k = 0; k < dims[2]; k++)
for(int j = 0; j < dims[1]; j++)
for(int i = 0; i < dims[0]; i++, idx++)
{
coords[0][idx] = static_cast<T>(i);
coords[1][idx] = static_cast<T>(j);
coords[2][idx] = static_cast<T>(k);
}
}
else
{
// Fuzz the coordinate in one of the dimensions.
const T maxScale = T{1.e-12};
int idx = 0;
for(int k = 0; k < dims[2]; k++)
for(int j = 0; j < dims[1]; j++)
for(int i = 0; i < dims[0]; i++, idx++)
{
const T f = rand_float<T>(maxScale);
const T z = 0;
coords[0][idx] = static_cast<T>(i) + ((fuzz == 0) ? f : z);
coords[1][idx] = static_cast<T>(j) + ((fuzz == 1) ? f : z);
coords[2][idx] = static_cast<T>(k) + ((fuzz == 2) ? f : z);
}
}
}

//---------------------------------------------------------------------------
Expand Down Expand Up @@ -339,11 +369,12 @@ void free_coords_2d(T **coords)

//---------------------------------------------------------------------------
template <typename T>
void kdtree_3d(int dims[3])
void kdtree_3d(int dims[3], int fuzz = -1)
{
// Initialization.
T *coords[3];
T *coords[3], *query_coords[3];
make_coords_3d(coords, dims);
make_coords_3d(query_coords, dims, fuzz);
int npts = dims[0] * dims[1] * dims[2];

// Make sure we can identify all of the coordinates.
Expand All @@ -353,9 +384,9 @@ void kdtree_3d(int dims[3])
for(int i = 0; i < npts; i++)
{
T pt[3];
pt[0] = coords[0][i];
pt[1] = coords[1][i];
pt[2] = coords[2][i];
pt[0] = query_coords[0][i];
pt[1] = query_coords[1][i];
pt[2] = query_coords[2][i];

found = search.findPoint(pt);
EXPECT_EQ(found, i);
Expand All @@ -370,6 +401,7 @@ void kdtree_3d(int dims[3])
EXPECT_EQ(found, search.NotFound);

free_coords_3d(coords);
free_coords_3d(query_coords);
}

//---------------------------------------------------------------------------
Expand Down Expand Up @@ -407,11 +439,12 @@ void kdtree_2d(int dims[2])

//---------------------------------------------------------------------------
template <typename T>
void single_domain_point_query_3d(int dims[3])
void single_domain_point_query_3d(int dims[3], int fuzz = -1)
{
// Initialization.
T *coords[3];
T *coords[3], *query_coords[3];
make_coords_3d(coords, dims);
make_coords_3d(query_coords, dims, fuzz);
int npts = dims[0] * dims[1] * dims[2];

// NOTE: We don't need the topology.
Expand All @@ -427,9 +460,9 @@ void single_domain_point_query_3d(int dims[3])
for(int i = 0; i < npts; i++)
{
double pt[3];
pt[0] = coords[0][i];
pt[1] = coords[1][i];
pt[2] = coords[2][i];
pt[0] = query_coords[0][i];
pt[1] = query_coords[1][i];
pt[2] = query_coords[2][i];

Q.add(domain0, pt);
}
Expand All @@ -451,6 +484,7 @@ void single_domain_point_query_3d(int dims[3])
EXPECT_EQ(res[badIdx], Q.NotFound);

free_coords_3d(coords);
free_coords_3d(query_coords);
}

//---------------------------------------------------------------------------
Expand Down Expand Up @@ -503,9 +537,27 @@ void single_domain_point_query_2d(int dims[2])
//-----------------------------------------------------------------------------
TEST(conduit_blueprint_mesh_query, kdtree_3d)
{
// Large enough to trigger accelerated search.
int dims[3]={30,30,30};
kdtree_3d<double>(dims);
kdtree_3d<float>(dims);

// 2D in 3D space - large enough that accelerated search should kick in.
// This checks that we can locate points in a plane in 3D using the kdtree.
// For some of these, fuzz the query points in one of the dimensions to test
// the kdtree bounding box fix.
for(int fuzz = 0; fuzz < 2; fuzz++)
{
int dims1[] = {1, 150, 150};
kdtree_3d<double>(dims1, fuzz > 0 ? 0 : -1);
kdtree_3d<float>(dims1, fuzz > 0 ? 0 : -1);
int dims2[] = {150, 1, 150};
kdtree_3d<double>(dims2, fuzz > 0 ? 1 : -1);
kdtree_3d<float>(dims2, fuzz > 0 ? 1 : -1);
int dims3[] = {150, 150, 1};
kdtree_3d<double>(dims3, fuzz > 0 ? 2 : -1);
kdtree_3d<float>(dims3, fuzz > 0 ? 2 : -1);
}
}

//-----------------------------------------------------------------------------
Expand All @@ -527,6 +579,21 @@ TEST(conduit_blueprint_mesh_query, point_query_3d)

single_domain_point_query_3d<double>(dims);
single_domain_point_query_3d<float>(dims);

// 2D in 3D space - large enough that accelerated search should kick in.
// This checks that we can locate points in a plane in 3D using the kdtree.
for(int fuzz = 0; fuzz < 2; fuzz++)
{
int dims1[] = {1, 150, 150};
single_domain_point_query_3d<double>(dims1, fuzz > 0 ? 0 : -1);
single_domain_point_query_3d<float>(dims1, fuzz > 0 ? 0 : -1);
int dims2[] = {150, 1, 150};
single_domain_point_query_3d<double>(dims2, fuzz > 0 ? 1 : -1);
single_domain_point_query_3d<float>(dims2, fuzz > 0 ? 1 : -1);
int dims3[] = {150, 150, 1};
single_domain_point_query_3d<double>(dims3, fuzz > 0 ? 2 : -1);
single_domain_point_query_3d<float>(dims3, fuzz > 0 ? 2 : -1);
}
}

//-----------------------------------------------------------------------------
Expand Down

0 comments on commit c8ea7af

Please sign in to comment.