interpolate_1st_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 PORTAGE_INTERPOLATE_INTERPOLATE_1ST_ORDER_H_
#define PORTAGE_INTERPOLATE_INTERPOLATE_1ST_ORDER_H_


#include <cassert>
#include <string>
#include <iostream>
#include <utility>
#include <vector>
#include <cmath>

// portage includes
#include "portage/intersect/dummy_interface_reconstructor.h"
#include "portage/support/portage.h"
#include "portage/driver/parts.h"

// wonton includes
#include "wonton/support/CoordinateSystem.h"

#ifdef HAVE_TANGRAM
#include "tangram/driver/driver.h"
#include "tangram/driver/CellMatPoly.h"
#include "tangram/support/MatPoly.h"
#endif

namespace Portage {

template<int D,
         Entity_kind on_what,
         typename SourceMeshType,
         typename TargetMeshType,
         typename SourceStateType,
         typename TargetStateType,
         typename T,
         template<class, int, class, class>
           class InterfaceReconstructorType = DummyInterfaceReconstructor,
         class Matpoly_Splitter = void,
         class Matpoly_Clipper = void,
         class CoordSys = Wonton::DefaultCoordSys
         >
class Interpolate_1stOrder {

#ifdef HAVE_TANGRAM
  using InterfaceReconstructor =
      Tangram::Driver<InterfaceReconstructorType, D, SourceMeshType,
                      Matpoly_Splitter, Matpoly_Clipper>;
#endif

 public:

  // Constructor with interface reconstructor
#ifdef HAVE_TANGRAM
  Interpolate_1stOrder(SourceMeshType const & source_mesh,
                       TargetMeshType const & target_mesh,
                       SourceStateType const & source_state,
                       NumericTolerances_t num_tols,
                       std::shared_ptr<InterfaceReconstructor> ir) :
      source_mesh_(source_mesh),
      target_mesh_(target_mesh),
      source_state_(source_state),
      interp_var_name_("VariableNameNotSet"),
      source_vals_(nullptr),
      num_tols_(num_tols),
      interface_reconstructor_(ir) {}
#endif

  Interpolate_1stOrder(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_1stOrder(const Interpolate_1stOrder &) = delete;

  //  Interpolate_1stOrder & operator = (const Interpolate_1stOrder &) = delete;

  ~Interpolate_1stOrder() = default;


  void set_material(int m) {
    matid_ = m;
  }  // set_material


  // Even though 1st order accurate interpolation does not require
  // limiters and only higher order ones do, we need to have the
  // limiter_type as an argument. This is because the templated driver
  // does not know whether its being called with 1st order or higher
  // order interpolators and so all interpolators need to have a
  // uniform interface

  void set_interpolation_variable(std::string const & interp_var_name,
                                  Portage::vector<Wonton::Vector<D>>* gradients = nullptr) {
    interp_var_name_ = interp_var_name;
    field_type_ = source_state_.field_type(Entity_kind::CELL, interp_var_name);
    if (field_type_ == Field_type::MESH_FIELD)
      source_state_.mesh_get_data(Entity_kind::CELL, interp_var_name,
                                  &source_vals_);
    else
      source_state_.mat_get_celldata(interp_var_name, matid_, &source_vals_);
  }  // set_interpolation_variable

  T operator() (int const targetEntityId,
                     std::vector<Weights_t> const & sources_and_weights) const {
    std::cerr << "Interpolation operator not implemented for this entity type"
              << std::endl;
    return T(0.0);
  }
  
  constexpr static int order = 1;

 private:
  SourceMeshType const & source_mesh_;
  TargetMeshType const & target_mesh_;
  SourceStateType const & source_state_;
  std::string interp_var_name_;
  T const * source_vals_;
  int matid_ = 0;
  Field_type field_type_ = Field_type::UNKNOWN_TYPE_FIELD;
  NumericTolerances_t num_tols_;
#ifdef HAVE_TANGRAM
  std::shared_ptr<InterfaceReconstructor> interface_reconstructor_;
#endif
};  // interpolate_1st_order base




template<int D,
         typename SourceMeshType,
         typename TargetMeshType,
         typename SourceStateType,
         typename TargetStateType,
         typename T,
         template<class, int, class, class> class InterfaceReconstructorType,
         class Matpoly_Splitter,
         class Matpoly_Clipper,
         class CoordSys>
class Interpolate_1stOrder<
  D, Entity_kind::CELL,
  SourceMeshType, TargetMeshType,
  SourceStateType, TargetStateType,
  T,
  InterfaceReconstructorType,
  Matpoly_Splitter, Matpoly_Clipper, CoordSys> {

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

#ifdef HAVE_TANGRAM
  using InterfaceReconstructor =
      Tangram::Driver<InterfaceReconstructorType, D, SourceMeshType,
                      Matpoly_Splitter, Matpoly_Clipper>;
#endif

 public:

#ifdef HAVE_TANGRAM
  Interpolate_1stOrder(SourceMeshType const & source_mesh,
                       TargetMeshType const & target_mesh,
                       SourceStateType const & source_state,
                       NumericTolerances_t num_tols,
                       std::shared_ptr<InterfaceReconstructor> ir,
                       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),
      interface_reconstructor_(ir),
      parts_(parts) {}
#endif

  Interpolate_1stOrder(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) {}


  //  Interpolate_1stOrder(const Interpolate_1stOrder &) = delete;

  //  Interpolate_1stOrder & operator = (const Interpolate_1stOrder &) = delete;

  ~Interpolate_1stOrder() = default;



  void set_material(int m) {
    matid_ = m;
  }  // set_material



  // Even though 1st order accurate interpolation does not require
  // limiters and only higher order ones do, we need to have the
  // limiter_type as an argument. This is because the templated driver
  // does not know whether its being called with 1st order or higher
  // order interpolators and so all interpolators need to have a
  // uniform interface

  void set_interpolation_variable(std::string const & interp_var_name,
                                  Portage::vector<Wonton::Vector<D>>* gradients = nullptr) {
    interp_var_name_ = interp_var_name;
    field_type_ = source_state_.field_type(Entity_kind::CELL, interp_var_name);
    if (field_type_ == Field_type::MESH_FIELD)
      source_state_.mesh_get_data(Entity_kind::CELL, interp_var_name,
                                  &source_vals_);
    else
      source_state_.mat_get_celldata(interp_var_name, matid_, &source_vals_);
  }  // set_interpolation_variable


  T operator() (int const targetCellID,
                     std::vector<Weights_t> const & sources_and_weights) const
  {
    int nsrccells = sources_and_weights.size();
    if (!nsrccells) return T(0.0);

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

    T val(0.0);
    double wtsum0 = 0.0;

    int nsummed = 0;
    if (field_type_ == Field_type::MESH_FIELD) {
      for (auto const& wt : sources_and_weights) {
        int srccell = wt.entityID;
        std::vector<double> pair_weights = wt.weights;
        if (fabs(pair_weights[0]) < num_tols_.min_absolute_volume)
          continue;  // skip small intersections
        val += source_vals_[srccell] * pair_weights[0];
        wtsum0 += pair_weights[0];
        nsummed++;
      }
    } else if (field_type_ == Field_type::MULTIMATERIAL_FIELD) {
      for (auto const& wt : sources_and_weights) {
        int srccell = wt.entityID;
        std::vector<double> pair_weights = wt.weights;
        if (fabs(pair_weights[0]) < num_tols_.min_absolute_volume)
          continue;  // skip small intersections
        int matcell = source_state_.cell_index_in_material(srccell, matid_);
        val += source_vals_[matcell] * pair_weights[0];  // 1st order
        wtsum0 += pair_weights[0];
        nsummed++;
      }
    }

    // Normalize the value by the volume of the intersection of the
    // target cells with the source mesh. This will do the right thing
    // for single-material and multi-material remap (conservative and
    // constant preserving) if there is NO mismatch between source and
    // target mesh boundaries. IF THERE IS A MISMATCH, THIS WILL
    // PRESERVE CONSTANT VALUES BUT NOT BE CONSERVATIVE. THEN WE HAVE
    // TO DO A SEMI-LOCAL OR GLOBAL REPAIR.

    // We use the * operator instead of / so as to reduce the number
    // of requirements on generic variable type
    
    if (nsummed)
      val *= (1.0/wtsum0);

    return val;
  }  // operator()
  
  constexpr static int order = 1;

 private:
  SourceMeshType const & source_mesh_;
  TargetMeshType const & target_mesh_;
  SourceStateType const & source_state_;
  std::string interp_var_name_;
  T const * source_vals_;
  int matid_ = 0;
  Field_type field_type_ = Field_type::UNKNOWN_TYPE_FIELD;
  NumericTolerances_t num_tols_;
#ifdef HAVE_TANGRAM
  std::shared_ptr<InterfaceReconstructor> interface_reconstructor_;
#endif
  Parts const* parts_;
};  // interpolate_1st_order specialization for cell





template<int D,
         typename SourceMeshType,
         typename TargetMeshType,
         typename SourceStateType,
         typename TargetStateType,
         typename T,
         template<class, int, class, class> class InterfaceReconstructorType,
         class Matpoly_Splitter,
         class Matpoly_Clipper,
         class CoordSys>
class Interpolate_1stOrder<
  D, Entity_kind::NODE,
  SourceMeshType, TargetMeshType,
  SourceStateType, TargetStateType,
  T,
  InterfaceReconstructorType,
  Matpoly_Splitter, Matpoly_Clipper, CoordSys> {

  // useful aliases
#ifdef HAVE_TANGRAM
  using InterfaceReconstructor =
      Tangram::Driver<InterfaceReconstructorType, D, SourceMeshType,
                      Matpoly_Splitter, Matpoly_Clipper>;
#endif

 public:

  // Constructor with interface reconstructor
#ifdef HAVE_TANGRAM
  Interpolate_1stOrder(SourceMeshType const & source_mesh,
                       TargetMeshType const & target_mesh,
                       SourceStateType const & source_state,
                       NumericTolerances_t num_tols,
                       std::shared_ptr<InterfaceReconstructor> ir) :
      source_mesh_(source_mesh),
      target_mesh_(target_mesh),
      source_state_(source_state),
      interp_var_name_("VariableNameNotSet"),
      source_vals_(nullptr),
      num_tols_(num_tols),
      interface_reconstructor_(ir) {}
#endif

  Interpolate_1stOrder(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_1stOrder(const Interpolate_1stOrder &) = delete;

  //  Interpolate_1stOrder & operator = (const Interpolate_1stOrder &) = delete;

  ~Interpolate_1stOrder() = default;


  void set_material(int m) {
    matid_ = m;
  }  // set_material


  // Even though 1st order accurate interpolation does not require
  // limiters and only higher order ones do, we need to have the
  // limiter_type as an argument. This is because the templated driver
  // does not know whether its being called with 1st order or higher
  // order interpolators and so all interpolators need to have a
  // uniform interface

  void set_interpolation_variable(std::string const & interp_var_name,
                                  Portage::vector<Wonton::Vector<D>>* gradients = nullptr) {
    interp_var_name_ = interp_var_name;
    field_type_ = source_state_.field_type(Entity_kind::NODE, interp_var_name);
    if (field_type_ == Field_type::MESH_FIELD)
      source_state_.mesh_get_data(Entity_kind::NODE, interp_var_name,
                                  &source_vals_);
    else {
      source_state_.mat_get_celldata(interp_var_name, matid_, &source_vals_);
      std::cerr << "Cannot remap NODE-centered multi-material data" << "\n";
    }
  }  // set_interpolation_variable


  T operator() (int const targetNodeID,
                std::vector<Weights_t> const & sources_and_weights) const
  {
    if (field_type_ != Field_type::MESH_FIELD) return T(0.0);

    int nsrcdualcells = sources_and_weights.size();
    if (!nsrcdualcells) return T(0.0);

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

    T val(0.0);
    double wtsum0 = 0.0;
    int nsummed = 0;
    for (auto const& wt : sources_and_weights) {
      int srcnode = wt.entityID;
      std::vector<double> pair_weights = wt.weights;
      if (fabs(pair_weights[0]) < num_tols_.min_absolute_volume)
        continue;  // skip small intersections
      val += source_vals_[srcnode] * pair_weights[0];  // 1st order
      wtsum0 += pair_weights[0];
      nsummed++;
    }

    // Normalize the value by the volume of the intersection of the
    // target cells with the source mesh. This will do the right thing
    // for single-material and multi-material remap (conservative and
    // constant preserving) if there is NO mismatch between source and
    // target mesh boundaries. IF THERE IS A MISMATCH, THIS WILL
    // PRESERVE CONSTANT VALUES BUT NOT BE CONSERVATIVE. THEN WE HAVE
    // TO DO A SEMI-LOCAL OR GLOBAL REPAIR.

    // We use the * operator instead of / so as to reduce the number
    // of requirements on generic variable type
    
    if (nsummed)
      val *= (1.0/wtsum0);

    return val;
  }  // operator()

  constexpr static int order = 1;

 private:
  SourceMeshType const & source_mesh_;
  TargetMeshType const & target_mesh_;
  SourceStateType const & source_state_;
  std::string interp_var_name_;
  T const * source_vals_;
  int matid_ = 0;
  Field_type field_type_ = Field_type::UNKNOWN_TYPE_FIELD;
  NumericTolerances_t num_tols_;
#ifdef HAVE_TANGRAM
  std::shared_ptr<InterfaceReconstructor> interface_reconstructor_;
#endif
};  // interpolate_1st_order specialization for nodes



}  // namespace Portage

#endif  // PORTAGE_INTERPOLATE_INTERPOLATE_1ST_ORDER_H_