tricubic_cell.h

// -*-c++-*-
#ifndef STORMM_TRICUBIC_CELL_H
#define STORMM_TRICUBIC_CELL_H
 
#include <vector>
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "Constants/behavior.h"
#include "Constants/scaling.h"
#include "DataTypes/common_types.h"
#include "DataTypes/stormm_vector_types.h"
#include "Parsing/parse.h"
#include "Reporting/error_format.h"
#include "math_enumerators.h"
#include "matrix_ops.h"
 
namespace stormm {
namespace stmath {
 
using card::Hybrid;
using card::HybridTargetLevel;
using constants::CartesianDimension;
using constants::UnitCellAxis;
using data_types::isFloatingPointScalarType;

class TricubicStencil {
public:

  TricubicStencil(Interpolant kind_in = Interpolant::SMOOTHNESS);

  TricubicStencil(const TricubicStencil &original) = default;
  TricubicStencil(TricubicStencil &&original) = default;
  TricubicStencil& operator=(const TricubicStencil &original) = default;
  TricubicStencil& operator=(TricubicStencil &&original) = default;

  Interpolant getKind() const;
  
  const double* data(HybridTargetLevel tier = HybridTargetLevel::HOST) const;
  
#ifdef STORMM_USE_HPC
#  ifdef STORMM_USE_CUDA
  const double* getDeviceViewToHostData() const;
#  endif

  void upload();

  void download();
#endif

  std::vector<double> exportMatrix() const;
  
private:
  Interpolant kind;          
  Hybrid<double> transform;  
};

template <typename T> class TricubicCell {
public:

  TricubicCell();
 
  TricubicCell(const TricubicStencil weights_matrix, const std::vector<double> &bounds,
               const std::vector<T> &f_in, const std::vector<T> &dx_in,
               const std::vector<T> &dy_in, const std::vector<T> &dz_in,
               const std::vector<T> &dxx_in, const std::vector<T> &dxy_in,
               const std::vector<T> &dxz_in, const std::vector<T> &dyy_in,
               const std::vector<T> &dyz_in, const std::vector<T> &dxxx_in,
               const std::vector<T> &dxxy_in, const std::vector<T> &dxxz_in,
               const std::vector<T> &dxyy_in, const std::vector<T> &dxyz_in);
 
  TricubicCell(const TricubicStencil weights_matrix, const std::vector<double> &bounds,
               const std::vector<T> &f_in, const std::vector<T> &dx_in,
               const std::vector<T> &dy_in, const std::vector<T> &dz_in,
               const std::vector<T> &dxy_in, const std::vector<T> &dxz_in,
               const std::vector<T> &dyz_in, const std::vector<T> &dxyz_in);

  T getCoefficient(int i, int j, int k) const;

  std::vector<T> getCoefficients() const;
  
  void setCoefficient(T value, int i, int j, int k);

  T getData(FunctionLevel kind, int i, int j, int k) const;

  void setData(T value, FunctionLevel kind, int i, int j, int k);

  T getCellOrigin(CartesianDimension dim) const;

  T getCellLength(CartesianDimension dim) const;

  void fractionalPosition(T x, T y, T z, T *a_frac, T *b_frac, T *c_frac,
                          ExceptionResponse policy = ExceptionResponse::SILENT,
                          const char* caller = nullptr) const;

  T evaluate(T x, T y, T z) const;

  template <typename T3> T3 derivative(T x, T y, T z) const;

  void secondDerivative(T x, T y, T z, T* result) const;
  void secondDerivative(T x, T y, T z, std::vector<T> *result) const;
  std::vector<T> secondDerivative(T x, T y, T z) const;
  
  void thirdDerivative(T x, T y, T z, T* result) const;
  void thirdDerivative(T x, T y, T z, std::vector<T> *result) const;
  std::vector<T> thirdDerivative(T x, T y, T z) const;
 
private:
  T coefficients[64];  
 
  // The following arrays store their series of values in "Fortran" order: (X0, Y0, Z0),
  // (X1, Y0, Z0), (X0, Y1, Z0), (X1, Y1, Z0), (X0, Y0, Z1), ...
  T f[8];      
  T dx[8];     
  T dy[8];     
  T dz[8];     
  T dxx[8];    
  T dxy[8];    
  T dxz[8];    
  T dyy[8];    
  T dyz[8];    
  T dxxx[8];   
  T dxxy[8];   
  T dxxz[8];   
  T dxyy[8];   
  T dxyz[8];   
 
  // The grid cell boundaries are stored in double precision for accuracy considerations.
  double origin_x;  
  double origin_y;  
  double origin_z;  

  double umat[9];
  
  double invu[9];
};

double3 stencilInternalOffset(UnitCellAxis face_normal, int a_index, int b_index, int c_index,
                              const double* invu);
 
double3 stencilInternalOffset(int a_index, int b_index, int c_index, const double* invu);

void incorporateStencilOrigin(const int95_t orig_x, const int95_t orig_y, const int95_t orig_z,
                              const double3 pt_xyz, double scale_factor, int95_t *point_x,
                              int95_t *point_y, int95_t *point_z);
 
template <typename T>
void incorporateStencilOrigin(T orig_x, T orig_y, T orig_z, const double3 pt_xyz,
                              double scale_factor, T *point_x, T *point_y, T *point_z);

template <typename T>
void fvStencilCoordinates(T orig_x, T orig_y, T orig_z, double scale_factor,
                          UnitCellAxis face_normal, int a_index, int b_index, int c_index,
                          const double* invu, T *point_x, T *point_y, T *point_z);
 
template <typename T>
void fvStencilCoordinates(T orig_x, T orig_y, T orig_z, double scale_factor, int a_index,
                          int b_index, int c_index, const double* invu, T *point_x,
                          T *point_y, T *point_z);
 
template <typename T>
void fvStencilCoordinates(T orig_x, T orig_y, T orig_z, UnitCellAxis face_normal, int a_index,
                          int b_index, int c_index, const double* invu, T *point_x, T *point_y,
                          T *point_z);
 
template <typename T>
void fvStencilCoordinates(T orig_x, T orig_y, T orig_z, int a_index, int b_index, int c_index,
                          const double* invu, T *point_x, T *point_y, T *point_z);
 
} // namespace stmath
} // namespace stormm
 
#include "tricubic_cell.tpp"
 
#endif