mesh_parameters.h

// -*-c++-*-
#ifndef STORMM_MESH_PARAMETERS_H
#define STORMM_MESH_PARAMETERS_H
 
#include <string>
#include <vector>
#include "copyright.h"
#include "Constants/behavior.h"
#include "DataTypes/common_types.h"
#include "DataTypes/stormm_vector_types.h"
#include "Math/math_enumerators.h"
#include "Math/rounding.h"
#include "Numerics/split_fixed_precision.h"
#include "Topology/atomgraph.h"
#include "Topology/atomgraph_abstracts.h"
#include "Topology/atomgraph_enumerators.h"
#include "Trajectory/coordinateframe.h"
#include "Trajectory/coordinate_series.h"
#include "Trajectory/phasespace.h"
#include "structure_enumerators.h"
 
namespace stormm {
namespace structure {
 
using constants::CartesianDimension;
using constants::UnitCellAxis;
using data_types::getStormmScalarTypeName;
using data_types::isSignedIntegralScalarType;
using data_types::isFloatingPointScalarType;
using data_types::isFloatingPointHpcVectorType;
using stmath::Interpolant;
using stmath::roundUp;
using topology::AtomGraph;
using topology::NonbondedKit;
using topology::UnitCellType;
using trajectory::CoordinateFrame;
using trajectory::CoordinateFrameReader;
using trajectory::CoordinateSeries;
using trajectory::CoordinateSeriesReader;
using trajectory::determineUnitCellTypeByShape;
  
constexpr int default_mesh_scaling_bits = 40;

constexpr int mesh_nonoverflow_bits = 46;

constexpr int max_mesh_definition_bits = 77;
  
struct MeshParamKit {

  MeshParamKit(int na_in, int nb_in, int nc_in, int95_t orig_x_in, int95_t orig_y_in,
               int95_t orig_z_in, double scale_in, double inv_scale_in, int scale_bits_in,
               const double* umat_in, const double* invu_in, const double* full_umat_in,
               const double* full_invu_in, const double* widths_in, const int95_t* fp_invu_in,
               double max_span_in, BoundaryCondition bounds_in, Interpolant stencil_kind_in,
               UnitCellType unit_cell_in);

  MeshParamKit(const MeshParamKit &original) = default;
  MeshParamKit(MeshParamKit &&original) = default;
  
  const int na;                    
  const int nb;                    
  const int nc;                    
  const int95_t orig_x;            
  const int95_t orig_y;            
  const int95_t orig_z;            
  const double scale;              
  const double inv_scale;          
  const float scale_f;             
  const float inv_scale_f;         
  const int scale_bits;            
  const double umat[9];            
  const double invu[9];            
  const double full_umat[9];       
  const double full_invu[9];       
  const double widths[3];          
  const int95_t fp_invu[9];        
  const double max_span;           
  const BoundaryCondition bounds;  
  const Interpolant stencil_kind;  
  const UnitCellType unit_cell;    
};
  
class MeshParameters {
public:

  MeshParameters(int na_in, int nb_in, int nc_in, double origin_x_in, double origin_y_in,
                 double origin_z_in, const std::vector<double> &element_vectors,
                 int scale_bits_in = default_mesh_scaling_bits,
                 Interpolant stencil_kind_in = Interpolant::SMOOTHNESS);
 
  MeshParameters();
 
  MeshParameters(int na_in, int nb_in, int nc_in, double origin_x_in, double origin_y_in,
                 double origin_z_in, double element_x, double element_y, double element_z,
                 int scale_bits_in = default_mesh_scaling_bits,
                 Interpolant stencil_kind_in = Interpolant::SMOOTHNESS);
 
  MeshParameters(int na_in, int nb_in, int nc_in, double origin_x_in, double origin_y_in,
                 double origin_z_in, double element_width,
                 int scale_bits_in = default_mesh_scaling_bits,
                 Interpolant stencil_kind_in = Interpolant::SMOOTHNESS);

  MeshParameters(const MeshParameters &original) = default;
  MeshParameters(MeshParameters &&original) = default;
  MeshParameters& operator=(const MeshParameters &other) = default;
  MeshParameters& operator=(MeshParameters &&other) = default;
  
  int getAxisElementCount(UnitCellAxis dim) const;
  int getAxisElementCount(CartesianDimension dim) const;
  
  template <typename Tcoord> std::vector<Tcoord> getMeshOrigin() const;
  double getMeshOrigin(CartesianDimension dim) const;

  template <typename T3> T3 getMeshOriginAsTuple() const;

  std::vector<int95_t> getMeshOriginAsFP() const;
  int95_t getMeshOriginAsFP(CartesianDimension dim) const;

  UnitCellType getMeshCellType() const;

  BoundaryCondition getBoundaryConditions() const;

  Interpolant getStencilKind() const;
  
  template <typename Tcoord> std::vector<Tcoord> getMeshElementVector(UnitCellAxis dim) const;
  template <typename Tcoord>
  std::vector<Tcoord> getMeshElementVector(CartesianDimension dim) const;

  template <typename T3> T3 getMeshElementVectorAsTuple(UnitCellAxis dim) const;
  template <typename T3> T3 getMeshElementVectorAsTuple(CartesianDimension dim) const;
  
  std::vector<int95_t> getMeshElementVectorAsFP(UnitCellAxis dim) const;
  std::vector<int95_t> getMeshElementVectorAsFP(CartesianDimension dim) const;

  template <typename Tcoord> std::vector<Tcoord> getMeshTransform() const;

  template <typename Tcoord> std::vector<Tcoord> getMeshInverseTransform() const;

  std::vector<int95_t> getMeshInverseTransformAsFP() const;

  int getScalingBits() const;
  
  double getScalingFactor() const;

  double getInverseScalingFactor() const;

  std::vector<int95_t> getAxisCoordinates(UnitCellAxis mesh_axis,
                                          CartesianDimension cart_axis) const;

  double getMaximumSpan() const;

  std::string printDimensions() const;
  
  MeshParamKit data() const;

  const MeshParameters* getSelfPointer() const;
  
  void setMeshDimension(int n_in, UnitCellAxis mesh_axis);
  void setMeshDimension(const std::vector<int> &n_in);
  
  void setOrigin(double v, CartesianDimension cart_axis);
  void setOrigin(int95_t v, CartesianDimension cart_axis);
  void setOrigin(const std::vector<double> &v);
  void setOrigin(const std::vector<int95_t> &v);

  void setScalingBits(int scale_bits_in);

  void defineElement(const std::vector<double> &element_vectors);

  void setBoundaryCondition(BoundaryCondition boundary_in);

  void setStencilKind(Interpolant stencil_kind_in);
  
private:
  int na;                       
  int nb;                       
  int nc;                       
  int95_t origin_x;             
  int95_t origin_y;             
  int95_t origin_z;             
  int scale_bits;               
  double scale_factor;          
  double inverse_scale_factor;  
  UnitCellType unit_cell;       
  BoundaryCondition boundary;   
  Interpolant stencil_kind;     

  double element_umat[9];

  float sp_element_umat[9];

  double element_invu[9];

  float sp_element_invu[9];

  double full_umat[9];

  double full_invu[9];
  
  double widths[3];

  float sp_widths[3];

  int95_t fp_element_invu[9];

  double maximum_span;

  void validateMeshDimensions() const;

  void validateFixedPrecisionBits() const;

  double maximumSpan() const;
};

MeshParameters getMeasurements(const AtomGraph *ag, const CoordinateFrame *cf, double padding,
                               double spacing, int scale_bits_in);
 
MeshParameters getMeasurements(const AtomGraph *ag, const CoordinateFrame *cf,
                               const std::vector<double> &mesh_bounds, double spacing,
                               int scale_bits_in);
 
MeshParameters getMeasurements(const AtomGraph *ag, const CoordinateFrame *cf, double padding,
                               const std::vector<double> &spacing, int scale_bits_in);
 
template <typename Tcoord>
MeshParameters getMeasurements(const AtomGraph &ag, const CoordinateSeries<Tcoord> &cs,
                               double padding, const std::vector<double> &spacing,
                               int scale_bits_in);
 
MeshParameters getMeasurements(const std::vector<double> &mesh_bounds,
                               const std::vector<double> &spacing, int scale_bits_in);
  
} // namespace structure
} // namespace stormm
 
#include "mesh_parameters.tpp"
 
#endif