Program Listing for File interpolation_base.h

↰ Return to documentation for file (src/pcms/transfer/interpolation_base.h)

#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