Skip to content

Commit

Permalink
Merge pull request #1228 from LLNL/bugfix/whitlock/partition_point_me…
Browse files Browse the repository at this point in the history 
…rge_fix

Fix point_merge in partitioner.
  • Loading branch information
@BradWhitlock
BradWhitlock committed Jan 11, 2024
2 parents 4fceb68 + ef5a29c commit 07f7a1b
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 112 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ and this project aspires to adhere to [Semantic Versioning](https://semver.org/s
- 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.
- The `conduit::blueprint::mpi::mesh::generate_partition_field()` function was not treating polyhedral topologies correctly, leading to unusable partitioning fields.
- The point merging algorithm in the Blueprint partitioner was corrected so it should no longer produce occasional duplicate points when merging coordsets.

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

Expand Down
141 changes: 29 additions & 112 deletions src/libs/blueprint/conduit_blueprint_mesh_partition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4689,9 +4689,6 @@ class point_merge
const std::vector<coord_system> &systems, index_t dimension,
double tolerance);

void truncate_merge(const std::vector<Node> &coordsets,
const std::vector<coord_system> &systems, index_t dimension, double tolerance);

static void xyz_to_rtp(double x, double y, double z, double &out_r, double &out_t, double &out_p);
// TODO
static void xyz_to_rz (double x, double y, double z, double &out_r, double &out_z);
Expand Down Expand Up @@ -4902,14 +4899,19 @@ using vec2 = vector<double,2>;
using vec3 = vector<double,3>;

/**
@brief A spatial search structure used to merge points within a given tolerance
@brief A spatial search structure used to merge points within a given tolerance.

@note This class differs from conduit::blueprint::mesh::utils::kdtree in that it
is designed to be constructed on the fly as points are inserted into it
and it contains the sorted data, which lends itself better to use with
multiple coordsets. The two classes were introduced independently and in
the future, it may be better to consolidate them.
*/
template<typename VectorType, typename DataType>
class kdtree
{
private:
using Float = typename VectorType::value_type;
// using IndexType = conduit_index_t;
public:
constexpr static auto dimension = std::tuple_size<typename VectorType::data_type>::value;
using vector_type = VectorType;
Expand Down Expand Up @@ -5121,22 +5123,7 @@ class kdtree
{
const bool left_contains = current->left->bb.contains(point, tolerance);
const bool right_contains = current->right->bb.contains(point, tolerance);
if(!left_contains && !right_contains)
{
// ERROR! This shouldn't happen, the tree must've been built improperly
retval = nullptr;
}
else if(left_contains)
{
// Traverse left
retval = find_point(current->left, depth+1, point, tolerance);
}
else if(right_contains)
{
// Traverse right
retval = find_point(current->right, depth+1, point, tolerance);
}
else // (left_contains && right_contains)
if(left_contains && right_contains)
{
// Rare, but possible due to tolerance.
// Check if the left side has the point without tolerance
Expand All @@ -5152,6 +5139,21 @@ class kdtree
: find_point(current->left, depth+1, point, tolerance);
}
}
else if(left_contains)
{
// Traverse left
retval = find_point(current->left, depth+1, point, tolerance);
}
else if(right_contains)
{
// Traverse right
retval = find_point(current->right, depth+1, point, tolerance);
}
else //if(!left_contains && !right_contains)
{
// ERROR! This shouldn't happen, the tree must've been built improperly
retval = nullptr;
}
}
else
{
Expand Down Expand Up @@ -5699,9 +5701,7 @@ point_merge::merge_data(const std::vector<Node> &coordsets,
const std::vector<coord_system> &systems, index_t dimension, double tolerance)
{
#define USE_SPATIAL_SEARCH_MERGE
#if defined(USE_TRUNCATE_PRECISION_MERGE)
truncate_merge(coordsets, systems, dimension, tolerance);
#elif defined(USE_SPATIAL_SEARCH_MERGE)
#if defined(USE_SPATIAL_SEARCH_MERGE)
spatial_search_merge(coordsets, systems, dimension, tolerance);
#else
simple_merge_data(coordsets, systems, dimension, tolerance);
Expand Down Expand Up @@ -5999,7 +5999,7 @@ point_merge::spatial_search_merge(const std::vector<Node> &coordsets,
const auto &coordset = coordsets[i];

// To be invoked on every coordinate
const auto merge = [&](float64 *p, index_t) {
const auto merge = [&](const float64 *p, index_t) {
vec3 key;
key.v[0] = p[0]; key.v[1] = p[1]; key.v[2] = p[2];
const auto potential_id = new_coords.size() / dimension;
Expand All @@ -6022,10 +6022,11 @@ point_merge::spatial_search_merge(const std::vector<Node> &coordsets,
}
};

const auto translate_merge = [&](float64 *p, index_t d) {
const auto translate_merge = [&](const float64 *p, index_t d) {
float64 tp[3];
translate_system(systems[i], coord_system::cartesian,
p[0], p[1], p[2], p[0], p[1], p[2]);
merge(p, d);
p[0], p[1], p[2], tp[0], tp[1], tp[2]);
merge(tp, d);
};

// Invoke the proper lambda on each coordinate
Expand All @@ -6045,90 +6046,6 @@ point_merge::spatial_search_merge(const std::vector<Node> &coordsets,
<< ", nodes in tree " << point_records.nodes() << std::endl);
}

//-----------------------------------------------------------------------------
void
point_merge::truncate_merge(const std::vector<Node> &coordsets,
const std::vector<coord_system> &systems, index_t dimension, double tolerance)
{
PM_DEBUG_PRINT("Truncate merging!" << std::endl);
// Determine what to scale each value by
// TODO: Be dynamic
(void)tolerance;
double scale = 0.;
{
auto decimal_places = 4u;
static const std::array<double, 7u> lookup = {
1.,
(2u << 4),
(2u << 7),
(2u << 10),
(2u << 14),
(2u << 17),
(2u << 20)
};
if(decimal_places < lookup.size())
{
scale = lookup[decimal_places];
}
else
{
scale = lookup[6];
}
}

/*index_t size = */reserve_vectors(coordsets, dimension);

// Iterate each of the coordinate sets
using fp_type = int64;
using tup = std::tuple<fp_type, fp_type, fp_type>;
std::map<tup, index_t> point_records;

for(size_t i = 0u; i < coordsets.size(); i++)
{
const auto &coordset = coordsets[i];

// To be invoked on every coordinate
const auto merge = [&](float64 *p, index_t) {
tup key = std::make_tuple(
static_cast<fp_type>(std::round(p[0] * scale)),
static_cast<fp_type>(std::round(p[1] * scale)),
static_cast<fp_type>(std::round(p[2] * scale)));
auto res = point_records.insert({key, {}});
if(res.second)
{
const index_t id = (index_t)(new_coords.size() / dimension);
res.first->second = id;
old_to_new_ids[i].push_back(id);
for(index_t j = 0; j < dimension; j++)
{
new_coords.push_back(p[j]);
}
}
else
{
old_to_new_ids[i].push_back(res.first->second);
}
};

const auto translate_merge = [&](float64 *p, index_t d) {
translate_system(systems[i], out_system,
p[0], p[1], p[2], p[0], p[1], p[2]);
merge(p, d);
};

// Invoke the proper lambda on each coordinate
if(systems[i] != out_system
&& systems[i] != coord_system::logical)
{
iterate_coordinates(coordset, translate_merge);
}
else
{
iterate_coordinates(coordset, merge);
}
}
}

//-----------------------------------------------------------------------------
void
point_merge::xyz_to_rtp(double x, double y, double z, double &out_r, double &out_t, double &out_p)
Expand Down

0 comments on commit 07f7a1b

Please sign in to comment.