interpolate_3rd_order.h

Go to the documentation of this file.

/*
This file is part of the Ristra portage project.
Please see the license file at the root of this repository, or at:
    https://github.com/laristra/portage/blob/master/LICENSE
*/

#ifndef SRC_INTERPOLATE_INTERPOLATE_3RD_ORDER_H_
#define SRC_INTERPOLATE_INTERPOLATE_3RD_ORDER_H_

#include <cassert>
#include <stdexcept>
#include <algorithm>
#include <string>
#include <iostream>
#include <utility>
#include <vector>

#include "portage/support/portage.h"
#include "portage/interpolate/quadfit.h"
#include "portage/driver/parts.h"

namespace Portage {

template<int D, Entity_kind on_what,
         typename SourceMeshType,
         typename TargetMeshType,
         typename SourceStateType,
         typename TargetStateType = SourceStateType>
class Interpolate_3rdOrder {

 public:
  Interpolate_3rdOrder(SourceMeshType const & source_mesh,
                       TargetMeshType const & target_mesh,
                       SourceStateType const & source_state,
                       NumericTolerances_t num_tols) :
      source_mesh_(source_mesh),
      target_mesh_(target_mesh),
      source_state_(source_state),
      interp_var_name_("VariableNameNotSet"),
      source_vals_(nullptr),
      num_tols_(num_tols) {}

  //  Interpolate_3rdOrder(const Interpolate_3rdOrder &) = delete;

  Interpolate_3rdOrder & operator = (const Interpolate_3rdOrder &) = delete;

  ~Interpolate_3rdOrder() = default;



  void set_interpolation_variable(std::string const & interp_var_name,
                                  Limiter_type limiter_type = NOLIMITER,
                                  Boundary_Limiter_type boundary_limiter_type = BND_NOLIMITER,
                                  Portage::vector<Vector<D>>* gradients = nullptr) {
    interp_var_name_ = interp_var_name;

    // Extract the field data from the statemanager

    source_state_.mesh_get_data(on_what, interp_var_name, &source_vals_);

    // Compute the limited quadfits for the field

    Limited_Quadfit<D, on_what, SourceMeshType, SourceStateType>
        limqfit(source_mesh_, source_state_, interp_var_name, 
                limiter_type, boundary_limiter_type);


    int nentities = source_mesh_.end(on_what)-source_mesh_.begin(on_what);
    quadfits_.resize(nentities);

    // call transform functor to take the values of the variable on
    // the cells and compute a "limited" quadfit of the field on the
    // cells (for transform definition, see portage.h)

    // Even though we defined Portage::transform (to be
    // thrust::transform or boost::transform) in portage.h, the
    // compiler is not able to disambiguate this call and is getting
    // confused. So we will explicitly state that this is Portage::transform

    Portage::transform(source_mesh_.begin(on_what), source_mesh_.end(on_what),
                       quadfits_.begin(), limqfit);
  }

  double operator() (int const targetCellID,
                     std::vector<Weights_t> const & sources_and_weights) const {
    // not implemented for all types - see specialization for cells and nodes

    std::cerr << "Interpolation operator not implemented for this entity type"
              << std::endl;
    return 0.;
  }

  constexpr static int order = 3;

 private:
  SourceMeshType const & source_mesh_;
  TargetMeshType const & target_mesh_;
  SourceStateType const & source_state_;
  std::string interp_var_name_;
  double const * source_vals_;
  NumericTolerances_t num_tols_;

  // Portage::vector is generalization of std::vector and
  // Wonton::Vector<D*(D+3)/2> is a geometric vector
  Portage::vector<Vector<D*(D+3)/2>> quadfits_;
};





template<int D,
         typename SourceMeshType,
         typename TargetMeshType,
         typename SourceStateType,
         typename TargetStateType>
class Interpolate_3rdOrder<
  D, Entity_kind::CELL,
  SourceMeshType, TargetMeshType,
  SourceStateType, TargetStateType> {

  // useful aliases
  using Parts = PartPair<
    D, SourceMeshType, SourceStateType,
    TargetMeshType, TargetStateType
  >;

 public:
  Interpolate_3rdOrder(SourceMeshType const & source_mesh,
                       TargetMeshType const & target_mesh,
                       SourceStateType const & source_state,
                       NumericTolerances_t num_tols,
                       const Parts* const parts = nullptr) :
      source_mesh_(source_mesh),
      target_mesh_(target_mesh),
      source_state_(source_state),
      interp_var_name_("VariableNameNotSet"),
      source_vals_(nullptr),
      num_tols_(num_tols),
      parts_(parts) {}



  void set_interpolation_variable(std::string const & interp_var_name,
                                  Limiter_type limiter_type = NOLIMITER,
                                  Boundary_Limiter_type boundary_limiter_type = BND_NOLIMITER,
                                  Portage::vector<Vector<D>>* gradients = nullptr) {

    interp_var_name_ = interp_var_name;

    // Extract the field data from the statemanager

    source_state_.mesh_get_data(Entity_kind::CELL, interp_var_name, &source_vals_);

    // Compute the limited quadfits for the field

    Limited_Quadfit<D, Entity_kind::CELL, SourceMeshType, SourceStateType>
        limqfit(source_mesh_, source_state_, interp_var_name_, limiter_type, boundary_limiter_type);

    int nentities = source_mesh_.end(Entity_kind::CELL)-source_mesh_.begin(Entity_kind::CELL);
    quadfits_.resize(nentities);

    // call transform functor to take the values of the variable on
    // the cells and compute a "limited" quadfit of the field on the
    // cells (for transform definition, see portage.h)

    // Even though we defined Portage::transform (to be
    // thrust::transform or boost::transform) in portage.h, the
    // compiler is not able to disambiguate this call and is getting
    // confused. So we will explicitly state that this is Portage::transform

    Portage::transform(source_mesh_.begin(Entity_kind::CELL), source_mesh_.end(Entity_kind::CELL),
                       quadfits_.begin(), limqfit);
  }

  //  Interpolate_3rdOrder(const Interpolate_3rdOrder &) = delete;

  Interpolate_3rdOrder & operator = (const Interpolate_3rdOrder &) = delete;

  ~Interpolate_3rdOrder() = default;


  double operator() (int const targetCellID,
                     std::vector<Weights_t> const & sources_and_weights) const {

    int nsrccells = sources_and_weights.size();
    if (!nsrccells) {
#ifdef DEBUG
      std::cerr << "WARNING: No source cells contribute to target cell." <<
        std::endl;
#endif
      return 0.0;
    }

    double totalval = 0.0;

    // contribution of the source cell is its field value weighted by
    // its "weight" (in this case, its 0th moment/area/volume)


    for (int j = 0; j < nsrccells; ++j) {
      int srccell = sources_and_weights[j].entityID;
      // int N = D*(D+3)/2;
      std::vector<double> xsect_weights = sources_and_weights[j].weights;
      double xsect_volume = xsect_weights[0];

      if (xsect_volume <= num_tols_.min_absolute_volume)
        continue;  // no intersection

      Point<D> srccell_centroid;
      source_mesh_.cell_centroid(srccell, &srccell_centroid);

      Point<D> xsect_centroid;
      for (int i = 0; i < D; ++i)
        xsect_centroid[i] = xsect_weights[1+i]/xsect_volume;  // (1st moment)/vol

      Vector<D*(D+3)/2> quadfit = quadfits_[srccell];
      Vector<D> vec = xsect_centroid - srccell_centroid;
      Vector<D*(D+3)/2> dvec;
      for (int j = 0; j < D; ++j) {
        int j1 = D;
        dvec[j] = vec[j];
        // Add the quadratic terms
        for (int k = 0; k <= j; ++k) {
          dvec[j1] = dvec[k]*dvec[j];
          j1 += 1;
          //dvec[D+j+k] = dvec[k]*dvec[j];
        }
      }
      double val = source_vals_[srccell] + dot(quadfit,dvec);
      val *= xsect_volume;
      totalval += val;
    }

    // Normalize the value by sum of all the 0th weights (which is the
    // same as the total volume of the source cell)

    totalval /= target_mesh_.cell_volume(targetCellID);

    return totalval;
  }

  constexpr static int order = 3;

 private:
  SourceMeshType const & source_mesh_;
  TargetMeshType const & target_mesh_;
  SourceStateType const & source_state_;
  std::string interp_var_name_;
  double const * source_vals_;
  NumericTolerances_t num_tols_;
  Parts const* parts_;

  // Portage::vector is generalization of std::vector and
  // Wonton::Vector<D> is a geometric vector
  Portage::vector<Vector<D*(D+3)/2>> quadfits_;
};






template<int D,
         typename SourceMeshType,
         typename TargetMeshType,
         typename SourceStateType,
         typename TargetStateType>
class Interpolate_3rdOrder<
  D, Entity_kind::NODE,
  SourceMeshType, TargetMeshType,
  SourceStateType, TargetStateType> {

 public:
  Interpolate_3rdOrder(SourceMeshType const & source_mesh,
                       TargetMeshType const & target_mesh,
                       SourceStateType const & source_state,
                       NumericTolerances_t num_tols) :
      source_mesh_(source_mesh),
      target_mesh_(target_mesh),
      source_state_(source_state),
      interp_var_name_("VariableNameNotSet"),
      source_vals_(nullptr),
      num_tols_(num_tols) {}

  //  Interpolate_3rdOrder(const Interpolate_3rdOrder &) = delete;

  Interpolate_3rdOrder & operator = (const Interpolate_3rdOrder &) = delete;

  ~Interpolate_3rdOrder() = default;



  void set_interpolation_variable(std::string const & interp_var_name,
                                  Limiter_type limiter_type = NOLIMITER,
                                  Boundary_Limiter_type boundary_limiter_type = BND_NOLIMITER,
                                  Portage::vector<Vector<D>>* gradients = nullptr) {

    interp_var_name_ = interp_var_name;

    // Extract the field data from the statemanager

    source_state_.mesh_get_data(Entity_kind::NODE, interp_var_name, &source_vals_);

    // Compute the limited quadfits for the field

    Limited_Quadfit<D, Entity_kind::NODE, SourceMeshType, SourceStateType>
        limqfit(source_mesh_, source_state_, interp_var_name, limiter_type, boundary_limiter_type);

    int nentities = source_mesh_.end(Entity_kind::NODE)-source_mesh_.begin(Entity_kind::NODE);
    quadfits_.resize(nentities);

    // call transform functor to take the values of the variable on
    // the cells and compute a "limited" quadfit of the field on the
    // cells (for transform definition, see portage.h)

    // Even though we defined Portage::transform (to be
    // thrust::transform or boost::transform) in portage.h, the
    // compiler is not able to disambiguate this call and is getting
    // confused. So we will explicitly state that this is Portage::transform

    Portage::transform(source_mesh_.begin(Entity_kind::NODE), source_mesh_.end(Entity_kind::NODE),
                       quadfits_.begin(), limqfit);
  }


  double operator() (const int targetNodeID,
                     std::vector<Weights_t> const & sources_and_weights) const {

    int nsrcnodes = sources_and_weights.size();
    if (!nsrcnodes) {
#ifdef DEBUG
      std::cerr << "WARNING: No source nodes contribute to target node." <<
        std::endl;
#endif
      return 0.0;
    }

    double totalval = 0.0;

    // contribution of the source cell is its field value weighted by
    // its "weight" (in this case, its 0th moment/area/volume)


    for (int j = 0; j < nsrcnodes; ++j) {
      int srcnode = sources_and_weights[j].entityID;
      std::vector<double> xsect_weights = sources_and_weights[j].weights;
      double xsect_volume = xsect_weights[0];

      if (xsect_volume <= num_tols_.min_absolute_volume)
        continue;  // no intersection

      // note: here we are getting the node coord, not the centroid of
      // the dual cell
      Point<D> srcnode_coord;
      source_mesh_.node_get_coordinates(srcnode, &srcnode_coord);

      Point<D> xsect_centroid;
      for (int i = 0; i < D; ++i)
        // (1st moment)/(vol)
        xsect_centroid[i] = xsect_weights[1+i]/xsect_volume;

      Vector<D*(D+3)/2> quadfit = quadfits_[srcnode];
      Vector<D> vec = xsect_centroid - srcnode_coord;
      Vector<D*(D+3)/2> dvec;
      for (int j = 0; j < D; ++j) {
        int j1 = D;
        dvec[j] = vec[j];
        // Add the quadratic terms
        for (int k = 0; k <= j; ++k) {
          dvec[j1] = dvec[k]*dvec[j];
          j1 += 1;
          // dvec[D+j+k] = dvec[k]*dvec[j];
        }
      }
      double val = source_vals_[srcnode] + dot(quadfit,dvec);
      val *= xsect_volume;
      totalval += val;
    }

    // Normalize the value by volume of the target dual cell

    totalval /= target_mesh_.dual_cell_volume(targetNodeID);

    return totalval;
  }

  constexpr static int order = 3;

 private:
  SourceMeshType const & source_mesh_;
  TargetMeshType const & target_mesh_;
  SourceStateType const & source_state_;
  std::string interp_var_name_;
  double const * source_vals_;
  NumericTolerances_t num_tols_;

  // Portage::vector is generalization of std::vector and
  // Wonton::Vector<D> is a geometric vector
  Portage::vector<Vector<D*(D+3)/2>> quadfits_;
};
}  // namespace Portage

#endif  // SRC_INTERPOLATE_INTERPOLATE_3RD_ORDER_H_