.. _program_listing_file_src_pcms_transfer_interpolation_base.h:

Program Listing for File interpolation_base.h
=============================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_pcms_transfer_interpolation_base.h>` (``src/pcms/transfer/interpolation_base.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #ifndef PCMS_TRANSFER_INTERPOLATION_BASE_H
   #define PCMS_TRANSFER_INTERPOLATION_BASE_H
   
   #include "pcms/field/evaluator/mls_interpolation.hpp"
   #include "pcms/localization/adj_search.hpp"
   #include "interpolation_helpers.h"
   #include <Omega_h_file.hpp>
   #include "pcms/utility/arrays.h"
   #include <string>
   
   namespace pcms
   {
   class InterpolationBase
   {
   public:
     virtual ~InterpolationBase() = default;
     virtual void eval(
       // TODO: Should these be templated to support different types?
       pcms::Rank1View<double, pcms::HostMemorySpace> source_field,
       pcms::Rank1View<double, pcms::HostMemorySpace> target_field) = 0;
   
     virtual size_t getSourceSize() const = 0;
   
     virtual size_t getTargetSize() const = 0;
   };
   
   class MLSPointCloudInterpolation final : public InterpolationBase
   {
   public:
     template <typename SourceType, typename TargetType>
     MLSPointCloudInterpolation(
       pcms::Rank1View<SourceType, pcms::HostMemorySpace> source_points,
       pcms::Rank1View<TargetType, pcms::HostMemorySpace> target_points, int dim,
       double radius, unsigned min_req_supports = 10, unsigned degree = 3,
       bool adapt_radius = true, double lambda = 0.0, double decay_factor = 5.0)
       : dim_(dim),
         radius_(radius),
         adapt_radius_(adapt_radius),
         degree_(degree),
         min_req_supports_(min_req_supports),
         n_targets_(target_points.size() / dim),
         n_sources_(source_points.size() / dim),
         lambda_(lambda),
         decay_factor_(decay_factor)
     {
       // check if source and target types can fall back to const double
       static_assert(std::is_convertible_v<SourceType, const Omega_h::Real>,
                     "SourceType must be convertible to const double");
       static_assert(std::is_convertible_v<TargetType, const Omega_h::Real>,
                     "TargetType must be convertible to const double");
   
       source_field_ = Omega_h::HostWrite<Omega_h::Real>(
         source_points.size() / dim_, "source field");
       target_field_ = Omega_h::HostWrite<Omega_h::Real>(
         target_points.size() / dim_, "target field");
   
       Omega_h::HostWrite<Omega_h::Real> source_coords_host(source_points.size(),
                                                            "source points");
       Omega_h::HostWrite<Omega_h::Real> target_coords_host(target_points.size(),
                                                            "target points");
       // TODO Remove these copies
       for (int i = 0; i < source_points.size(); ++i) {
         source_coords_host[i] = source_points[i];
       }
       for (int i = 0; i < target_points.size(); ++i) {
         target_coords_host[i] = target_points[i];
       }
       source_coords_ = Omega_h::Reals(source_coords_host);
       target_coords_ = Omega_h::Reals(target_coords_host);
   
       find_supports(min_req_supports_, 3 * min_req_supports_, 100);
     }
   
     void eval(
       pcms::Rank1View<double, pcms::HostMemorySpace> source_field,
       pcms::Rank1View<double, pcms::HostMemorySpace> target_field) override;
   
     SupportResults getSupports() { return supports_; }
     size_t getSourceSize() const override { return source_coords_.size() / dim_; }
     size_t getTargetSize() const override { return target_coords_.size() / dim_; }
   
   private:
     int dim_;           
     double radius_;     
     bool adapt_radius_; 
     unsigned
       degree_; 
     unsigned
       min_req_supports_; 
     double lambda_; 
     double decay_factor_; 
     // InterpolationType interpolation_type_;
     Omega_h::LO n_sources_ = 0;    
     Omega_h::LO n_targets_ = 0;    
     Omega_h::Reals source_coords_; 
     Omega_h::Reals target_coords_; 
     SupportResults supports_; 
     Omega_h::HostWrite<Omega_h::Real> target_field_; 
     Omega_h::HostWrite<Omega_h::Real> source_field_; 
     void fill_support_structure(Omega_h::Write<Omega_h::Real> radii2_l,
                                 Omega_h::Write<Omega_h::LO> num_supports);
   
     void distance_based_pointcloud_search(
       Omega_h::Write<Omega_h::Real> radii2_l,
       Omega_h::Write<Omega_h::LO> num_supports) const;
   
     void find_supports(unsigned min_req_supports = 10,
                        unsigned max_allowed_supports = 30,
                        unsigned max_count = 100);
   };
   
   class MLSMeshInterpolation final : public InterpolationBase
   {
   
   public:
     void eval(
       pcms::Rank1View<double, pcms::HostMemorySpace> source_field,
       pcms::Rank1View<double, pcms::HostMemorySpace> target_field) override;
   
     MLSMeshInterpolation(Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh,
                          double radius, unsigned min_req_supports = 10,
                          unsigned degree = 3, bool adapt_radius = true,
                          double lambda = 0.0, double decay_factor = 5.0);
   
     MLSMeshInterpolation(Omega_h::Mesh& source_mesh, double radius,
                          unsigned min_req_supports = 10, unsigned degree = 3,
                          bool adapt_radius = true, double lambda = 0.0,
                          double decay_factor = 5.0);
   
     size_t getSourceSize() const override;
     size_t getTargetSize() const override;
   
     SupportResults getSupports() { return supports_; }
   
   private:
     double radius_; 
     double lambda_; 
     double decay_factor_; 
     bool adapt_radius_; 
     bool single_mesh_ = false; 
     unsigned
       degree_; 
     unsigned
       min_req_supports_; 
     // InterpolationType interpolation_type_;
   
     Omega_h::Mesh& source_mesh_; 
     // TODO: handle what to do with this when only 1 mesh is provided
     Omega_h::Mesh& target_mesh_;   
     Omega_h::Reals source_coords_; 
     Omega_h::Reals target_coords_; 
     SupportResults supports_; 
     Omega_h::HostWrite<Omega_h::Real> target_field_; 
     Omega_h::HostWrite<Omega_h::Real> source_field_; 
     void find_supports(unsigned min_req_supports = 10,
                        unsigned max_allowed_supports = 30);
   };
   } // namespace pcms
   
   #endif // PCMS_TRANSFER_INTERPOLATION_BASE_H