mesh bounding box bug and enhancement of points_not_found in LagrangeInterpolator::interpolate
Issue #1132
new
- when creating bounding box of mesh0 in line 167, donot init the min coordinate with a big valus, so the min coordinate will be always 0, the bbox is wrong. Should add these code:
// Create bounding box of mesh0
std::vector<double> x_min_max(2*gdim0);
std::vector<double> coordinates = mesh0.coordinates();
// add ------------------------------
#include <cfloat>
for (std::size_t i = 0; i < gdim0; ++i)
{
x_min_max[i] = dbl_max;
}
// add ------------------------------
for (std::size_t i = 0; i < gdim0; ++i)
{
for (auto it = coordinates.begin() + i; it < coordinates.end(); it += gdim0)
{
x_min_max[i] = std::min(x_min_max[i], *it);
x_min_max[gdim0 + i] = std::max(x_min_max[gdim0 + i], *it);
}
}
2. line 218 insert all not found points in the points_not_found
vector, then MPI::all_to_all these points and search all these points, which is no need. Can check these points if in bounding box before MPI, to reduce communication time and for_loop times. Can merge line 201 to line 247 like these:
// Search this process first for all coordinates in u's local mesh
std::vector<double> points_not_found;
// Get number of MPI processes
std::size_t num_processes = MPI::size(mpi_comm);
// Remaining interpolation points must be found through MPI
// communication. Check first using bounding boxes which process
// may own the points
std::vector<std::vector<double>> potential_points(num_processes);
for (const auto &map_it : coords_to_dofs)
{
// Place interpolation point in x
std::copy(map_it.first.begin(), map_it.first.end(), x.begin());
try
{ // Store values when point is found
u0.eval(_values, _x);
std::vector<std::size_t> dofs = map_it.second;
for (const auto &d : map_it.second)
local_u_vector[d] = values[dof_component_map[d]];
}
catch (std::exception &e)
{
// If not found then it must be searched on the other processes
points_not_found.insert(points_not_found.end(), x.begin(), x.end());
// Find potential owners
for (std::size_t p = 0; p < num_processes; p++)
{
if (p == MPI::rank(mpi_comm))
continue;
// Check if in bounding box
if (in_bounding_box(x, bounding_boxes[p], 1e-12))
{
potential_points[p].insert(potential_points[p].end(),
x.begin(), x.end());
}
}
}
}
May I have the pull request permission to add my fix code to the project?