mesh bounding box bug and enhancement of points_not_found in LagrangeInterpolator::interpolate

Issue #1132 new
旭丁 created an issue
  1. 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?

Comments (0)

  1. Log in to comment