nearest_point.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/mpi.templates.h>
#include <deal.II/base/mpi_noncontiguous_partitioner.h>
#include <deal.II/base/mpi_remote_point_evaluation.h>
#include <deal.II/base/point.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/types.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/mapping.h>
 
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include <deal.II/matrix_free/fe_point_evaluation.h>
 
#include <deal.II/numerics/rtree.h>
#include <deal.II/numerics/vector_tools_evaluate.h>
 
#include <meltpooldg/level_set/level_set_tools.hpp>
#include <meltpooldg/level_set/nearest_point_data.hpp>
 
#include <boost/geometry/index/rtree.hpp>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <functional>
#include <optional>
#include <string>
#include <vector>
 
 
namespace MeltPoolDG::LevelSet::Tools
{
  template <int dim, typename number>
  class NearestPoint
  {
    using VectorType      = dealii::LinearAlgebra::distributed::Vector<number>;
    using BlockVectorType = dealii::LinearAlgebra::distributed::BlockVector<number>;
 
  public:
    NearestPoint(const dealii::Mapping<dim>    &mapping,
                 const dealii::DoFHandler<dim> &dof_handler_signed_distance,
                 const VectorType              &signed_distance,
                 const BlockVectorType         &normal_vector,
                 dealii::Utilities::MPI::RemotePointEvaluation<dim, dim>   &remote_point_evaluation,
                 const NearestPointData<number>                            &nearest_point_data,
                 std::optional<std::reference_wrapper<dealii::TimerOutput>> timer_output = {});
 
    void
    reinit(const dealii::DoFHandler<dim> *dof_handler_src_in,
           const dealii::DoFHandler<dim> *dof_handler_dst_in = nullptr);
 
    template <int n_components = 1>
    void
    extend_interface_values(VectorType                                &solution_out,
                            const VectorType                          &solution_in,
                            const bool                                 zero_out  = true,
                            const std::function<number(const number)> &operation = {}) const;
 
 
    const std::vector<dealii::Point<dim>> &
    get_projected_points() const;
 
    const std::vector<std::vector<dealii::types::global_dof_index>> &
    get_projected_points_dof_indices() const;
 
    void
    write_to_file(const std::string filename_without_ending) const;
 
    mutable std::vector<dealii::Point<dim>> points_not_found;
 
  private:
    const dealii::Mapping<dim> &mapping;
 
    const dealii::DoFHandler<dim> &dof_handler_ls;
 
    const VectorType &signed_distance;
 
    const BlockVectorType &normal_vector;
 
    const NearestPointData<number> &nearest_point_data;
 
    dealii::Utilities::MPI::RemotePointEvaluation<dim, dim> &remote_point_evaluation;
 
    // Tolerance to be reached for the distance of the projected points to the distance = 0
    // isosurface
    const number tol_distance;
 
    const number narrow_band_threshold;
 
    const number tolerance_normal_vector;
 
    const MPI_Comm mpi_comm;
 
    dealii::ConditionalOStream pcout;
    std::optional<std::reference_wrapper<dealii::TimerOutput>> timer_output;
 
    const dealii::DoFHandler<dim> *dof_handler_src = nullptr;
 
    const dealii::DoFHandler<dim> *dof_handler_dst = nullptr;
 
    // --- Output data: projection results ---
 
    std::vector<dealii::Point<dim>> projected_points_at_interface;
 
    std::vector<std::vector<dealii::types::global_dof_index>> dof_indices;
 
    std::vector<dealii::Point<dim>> stencil;
 
    // --- Internal state ---
 
    int total_points_rpe = 0;
 
    bool is_reinit_called = false;
 
    bool input_vector_is_cut = false;
 
    // --- Internal data for fast path (ArborX-based projection) ---
 
    dealii::Utilities::MPI::NoncontiguousPartitioner partitioner;
 
    std::vector<dealii::types::global_dof_index> locally_owned_surface_indices;
 
    std::vector<dealii::types::global_dof_index> ghost_surface_indices;
 
    std::vector<std::pair<std::pair<unsigned int, unsigned int>, std::vector<dealii::Point<dim>>>>
      surface_cells_and_unit_points;
 
    template <int n_components = 1>
    void
    extend_interface_values_nearest_point_fast(
      VectorType                                &solution_out,
      const VectorType                          &solution_in,
      const bool                                 zero_out  = true,
      const std::function<number(const number)> &operation = {}) const;
 
    void
    clear_cached_data();
 
    void
    register_dof_handlers(const dealii::DoFHandler<dim> *dof_handler_src_in,
                          const dealii::DoFHandler<dim> *dof_handler_dst_in);
 
    void
    collect_narrow_band_support_points();
 
    void
    run_projection_algorithm();
 
    bool
    local_compute_normal_correction(std::vector<dealii::Point<dim>> &y);
 
    bool
    local_compute_normal_and_tangential_correction_coquerelle(std::vector<dealii::Point<dim>> &y);
 
    bool
    local_compute_normal_and_tangential_correction(std::vector<dealii::Point<dim>> &y);
 
    void
    local_compute_nearest_point();
 
    void
    local_compute_nearest_point_fast();
  };
} // namespace MeltPoolDG::LevelSet::Tools