Coarsing.h

Go to the documentation of this file.

 
#ifndef COARSING_H
#define COARSING_H
 
#include <cstdlib>
#include <cstdio>
#include <vector>
#include "Typedefs.h"
#include <boost/random.hpp>
#include <fstream>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/crc.hpp>
#include <map>
#include <variant>
 
#ifdef NETCDF
#include <netcdf.h>
#endif
 
#ifdef NRRDIO
#include "../NrrdIO-1.11.0-src/NrrdIO.h"
#endif
 
#ifdef MATLAB
#include "mat.h"
#endif
 
#include "Eigen/Dense"
#include "Eigen/Geometry"
 
 
using namespace std ;
 
double Volume (int d , double R) ; 
 
enum TensorOrder {NONE=-1, SCALAR=0, VECTOR=1, TENSOR=2} ;
enum FieldType {Defined, Particle, Fluctuation, Contact} ;
enum AverageType {None, Final, Intermediate, Both, Pre5, IntermediateAndPre5} ;
enum Pass {Pass1=1, Pass2=2, Pass3=4, Pass4=8, Pass5=16,
           VelFluct=256, RotFluct=512} ;
 
inline Pass operator|(Pass a, Pass b){return static_cast<Pass>(static_cast<int>(a) | static_cast<int>(b));}
//inline Pass operator|=(Pass a, const Pass b){return static_cast<Pass>(static_cast<int>(a) | static_cast<int>(b));} // Not working for some reason
inline bool operator& (Pass a, Pass b) { return (static_cast<int>(a) & static_cast<int>(b)) ; }
 
//=========================================================
class CGPoint
{
public :
    CGPoint(int dd, v1d loc) {d=dd ; location=loc ; }
 
    v2d fields ;    
    v2d subfields ; 
    v1d location ;  
    //Useful things
    vector <int> neighbors ; 
    int d ; 
} ;
//-------------------------
struct Field {
 uint64_t flag ; 
 string name ; 
 TensorOrder type ; 
 FieldType ftype ; 
 Pass passlevel ; 
 int datalocation = -1; 
};
//-------------------------
struct Data {
public:
    Data () : radius(nullptr), mass(nullptr), Imom(nullptr), id1(nullptr), id2(nullptr) {}
int N = 0 ; 
double * radius ; 
double * mass ; 
double *Imom ; 
vector <double *> pos ; 
vector <double *> vel ;  
vector <double *> omega ; 
vector <double *> orient ; 
 
// Superquadrics
vector <double *> superquadric ; 
 
v2d vel_fluc ; 
v2d rot_fluc ; 
 
int Ncf ; 
double * id1 ; 
double * id2 ; 
vector <double *> pospq ; 
vector <double *> lpq ;  
vector <double *> fpq; 
vector <double *> mpq; 
vector <double *> mqp ; 
 
// Exta fields if needed
vector<double*> extra ;
vector<std::tuple<std::string, int, int>> extrafields ;
 
// Some useful functions
int Nnonper=-1 ; 
int random_test (int N, int Ncf, int d, v2d box ) ; 
int compute_lpq (int d) ; 
int periodic_atoms (int d, v2d bounds, int pbc, v1d Delta, bool omegainclude) ; 
int clean_periodic_atoms () {if (Nnonper==-1) printf("ERR: must call periodic_atoms before cleaning the periodic_atoms\n") ; else N=Nnonper ; return 0 ; } 
bool check_field_availability(string name) ;
int add_extra_field (int length, std::string name)
{
  extrafields.push_back({name, length, extra.size()}) ;
  extra.resize(extra.size()+length) ;
  return std::get<2>(extrafields.back());
}
} ;
//------------------------------------------------------
#include "WindowLibrary.h"
 
//=========================================================
class Coarsing
{
public :
    Coarsing (int dd, v1i nnpt, v2d bbox, int T) : cT(0), flags(0)
    {
        d=dd ; npt=nnpt ; box=bbox ; Time=T ;
        dx.resize(d, 0) ;
        for (int i=0 ; i<d ; i++)
          dx[i]=((box[1][i]-box[0][i])/double(npt[i])) ;
        double w= (*std::min_element(dx.begin(),dx.end())*2) ; // w automatically set
        cutoff=2.5*w ; //TODO
        printf("Window and cutoff: %g %g \n", w, cutoff) ;
        //for (int i=0 ; i<d ; i++)
        // printf("%d %d %g %g %g|", d, npt[i], box[1][i], box[0][i], dx[i]) ; fflush(stdout) ;
        grid_generate() ;
        //grid_neighbour() ;
        set_field_struct() ;
        Window = new LibLucy3D( &data, w, d) ;
    }
    ~Coarsing() { if (Window != nullptr) delete Window ;
                  if (CGPtemp != nullptr) delete CGPtemp ; }
 
    int d ; 
    int Npt; 
    int Time; 
    int cT ; 
    double cutoff ; 
    vector <CGPoint> CGP ; 
    vector <CGPoint> * CGPtemp = nullptr ; 
    vector <int> npt; 
    vector <int> nptcum ; 
    v1d dx ; 
    v2d box ; 
    LibBase * Window = nullptr ; 
    std::variant<Eigen::Matrix3d, Eigen::Quaternion<double>> original_orientation = Eigen::Quaternion<double>(0,0,0,1) ;     
 
    // Fields variable and function
    vector <struct Field > FIELDS ; 
    unsigned int flags ; 
    vector <int> Fidx ;  
    vector <std::pair<Field *, int>> Fcols ; 
 
    //vector <string> Fields, Fname ; ///< Flagged field names
    //vector <TensorOrder> Ftype ; ///< Flagged field types
    int get_id(string nm) ; 
    struct Field * get_field(string nm) ; 
    Pass set_flags (vector <string> s) ; 
 
    // Subfield variables and functions
    vector <struct Field> SUBFIELDS ; 
    vector <std::pair<struct Field *, uint64_t>> subflags ;
    vector <std::pair<Field *, std::vector<std::pair<Field*, int>> >> SFcols ; 
    uint64_t sf_mask_eigen ;
    int grid_getsubfields() ;
    struct Field * get_subfield(string nm) ; 
    
    // Grid functions
    int set_field_struct() ; 
    int add_extra_field(string name, TensorOrder order, FieldType type) ; 
    int setWindow (Windows win, double w, vector <bool> per ={}, vector<int> boxes = {}, vector<double> deltas = {}, vector<vector<double>> bounds={{}}) ; 
    template <Windows W> int setWindow () ; 
    template <Windows W> int setWindow (double w) ; 
    template <Windows W> int setWindow (vector<vector<double>> &) ; 
    template <Windows W> int setWindow (double w, vector<bool> per, vector<int> boxes, vector<double> deltas) ; 
    int grid_generate() ; 
    int grid_neighbour() ; 
    std::map<std::string, size_t> grid_setfields() ; 
    int grid_getfields() ; 
    v2d get_bounds() ; 
    CGPoint * reverseloop (string type) ; 
    int find_closest (int id) ; 
    int find_closest_pq (int id) ; 
    v1d interpolate_vel(int id, bool usetimeavg=false) { return interpolate_vel_nearest (id, usetimeavg) ; } 
    v1d interpolate_rot(int id, bool usetimeavg=false) { return interpolate_rot_nearest (id, usetimeavg) ; } 
    v1d interpolate_vel_nearest (int id, bool usetimeavg=false) ; 
    v1d interpolate_rot_nearest (int id, bool usetimeavg=false) ; 
    v1d interpolate_vel_trilinear (int id, bool usetimeavg) ; 
    template <int D> v1d interpolate_vel_multilinear (int id, bool usetimeavg);
    template <typename T> void add_rotated_quat(double * p, double weight, Eigen::Quaternion<double> q, T original) ; 
    
    int idx_FastFirst2SlowFirst (int n) ; 
 
    // Windowing functions
    //double window(double r) {Lucy(r) ; }
    //double window_int (v1d r1, v1d lpq, v1d x) {printf("Numerical integration of wpqf not implemented\n") ; } ///< Numerical integration: not implemented
    //double window_int(double r1, double r2) {return window_avg(r1, r2) ; } ///< Overload to avoid integration ...
    //double window_avg (double r1, double r2) {return (0.5*(Lucy(r1)+Lucy(r2))) ; }
    //double Lucy (double r) {static double cst=105./(16*M_PI*w*w*w) ; if (r>=w) return 0 ; else {double f=r/w ; return (cst*(-3*f*f*f*f + 8*f*f*f - 6*f*f +1)) ; }}
    double normdiff (v1d a, v1d b) {double res=0 ; for (int i=0 ; i<d ; i++) res+=(a[i]-b[i])*(a[i]-b[i]) ; return (sqrt(res)) ; } ; 
    // Coarse graining functions
    int pass_1 () ; 
    int pass_2 (bool usetimeavg=false) ; 
    int pass_3 () ; 
    int pass_4 () ; 
    int pass_5 () ; 
    int compute_fluc_vel (bool usetimeavg=false) ; 
    int compute_fluc_rot (bool usetimeavg=false) ; 
    bool hasvelfluct=false, hasrotfluct=false ;
 
    // Data handling functions
 
    struct Data data ; 
 
    // Time and output handling
    int mean_time(bool temporary=false) ; 
    int write_vtk(string sout) ; 
    int write_netCDF (string sout) ; 
    int write_NrrdIO (string path) ; 
    int write_matlab (string path, bool squeeze = false) ; 
    int write_numpy (string path, bool squeeze = false) ; 
    int write_numpy_npy (string path, bool squeeze) ;
    std::pair<size_t, uint8_t*> write_numpy_locbuffer (bool squeeze) {return write_numpy_buffer(-2, squeeze) ; }
    std::pair<size_t, uint8_t*> write_numpy_buffer (int id, bool squeeze) ;
} ;
 
 
//-------------------------------------------------------
template <Windows W>
int Coarsing::setWindow ()
{ double w= (*std::min_element(dx.begin(),dx.end())*1) ; // w automatically set
  setWindow<W>(w) ; return 0 ; }
//-------------------------------------------------------
template <Windows W>
int Coarsing::setWindow (double w)
{
  static_assert(W != Windows::LucyND_Periodic) ;
  switch (W) {
      case Windows::Rect3D :
         Window=new LibRect3D (&data, w, d) ;
        break ;
      case Windows::Sphere3DIntersect :
        Window=new LibSphere3DIntersect (&data, w, d) ;
        break ;
      case Windows::Sphere3DIntersect_MonteCarlo:
        Window= new LibSphere3DIntersect_MonteCarlo(&data, w, d) ; 
        break ; 
      case Windows::SphereNDIntersect :
        Window=new LibSphereNDIntersect (&data, w, d) ;
        break ;
      case Windows::Lucy3D :
        Window=new LibLucy3D (&data, w, d) ;
        break ;
      case Windows::Lucy3DFancyInt :
        Window=new LibLucy3DFancyInt (&data, w, d) ;
        break ;
      case Windows::Hann3D :
        Window=new LibHann3D (&data, w, d) ;
        break ;
      case Windows::RectND :
        Window=new LibRectND (&data, w, d) ;
        break ;
      case Windows::LucyND :
        Window=new LibLucyND (&data, w, d) ;
        break ;
      default:
        printf("Unknown window, check Coarsing::setWindow #2") ;
  }
  cutoff = Window->cutoff() ;
  printf("Window and cutoff: %g %g \n", w, cutoff) ;
  return 0 ;
}
//-------------------------------------------------------
template <Windows W>
int Coarsing::setWindow (double w, vector<bool> per, vector<int> boxes, vector<double> deltas)
{
  static_assert(W == Windows::LucyND_Periodic) ;
  cutoff = Window->cutoff() ;
  printf("Window and cutoff: %g %g \n", w, cutoff) ;
 
  int p = 0 ;
  for (size_t i=0 ; i<per.size() ; i++)
    if (per[i])
        p |= (1<<i) ;
 
  Window = new LibLucyND_Periodic (&data,w,d,p,boxes,deltas) ;
return 0 ;
}
//-----------------------------------------------------------
template <Windows W>
int Coarsing::setWindow (vector<vector<double>> & bounds)
{
  static_assert(W == Windows::RVE) ;
  double scale = 1 ;
  for (size_t i =0 ; i<bounds.size() ; i++)
      scale /= (bounds[i][1] - bounds[i][0]) ;
 
  Window = new LibRVE (scale) ;
return 0 ;
}
 
//-----------------------------------------------------------------------------------------
template <int D>
v1d Coarsing::interpolate_vel_multilinear (int id, bool usetimeavg)
{
int pts[1<<D][D] ;
const static int idvel=get_id("VAVG") ;
 
auto clip = [&](double a, int maxd){if (a<0.) return(0) ; else if (a>=maxd) return (maxd-1) ; else return (static_cast<int>(a)) ; } ;
 
 
// Determine the floor and ceil indices in each dimension
std::vector<int> i0(D), i1(D);
for (int dd = 0; dd < D; ++dd) {
    double val = (data.pos[dd][id] - CGP[0].location[dd]) / dx[dd];
    i0[dd] = clip(std::floor(val), npt[dd]);
    i1[dd] = clip(std::ceil(val), npt[dd]);
    if (npt[dd]==1) continue ; 
    if (i0[dd]==i1[dd])
    {
      if (i0[dd]==0) i1[dd]++ ;
      else if (i0[dd]==npt[dd]-1) i0[dd]-- ; 
      else i1[dd] ++ ; // If we are exactly on a point, just push one of the corner. The weights should be 0 anyway for this point.  
    }
}
 
// Construct all corner indices of the cube
int n_corners = 1<<D ; 
for (int c = 0; c < n_corners; ++c) {
    for (int dd = 0; dd < D; ++dd) {
        pts[c][dd] = ((c >> dd) & 1) ? i1[dd] : i0[dd];
    }
}
 
// Convert multi-dim index to linear index
std::vector<int> lin_idx(n_corners);
for (int c = 0; c < n_corners; ++c) {
    int idx = 0;
    for (int dd = 0; dd < D; ++dd)
        idx += nptcum[dd] * pts[c][dd];
    lin_idx[c] = idx;
}
 
std::vector<double> weights(n_corners, 1.0);
for (int i = 0; i < n_corners; ++i) {
    for (int j = 0; j < D; ++j) {
        double x0 = CGP[lin_idx[i]].location[j];
        double t = (data.pos[j][id] - x0) / dx[j];
        if (((i >> j) & 1) == 0)
            weights[i] *= npt[j]==1?0.5:(1.0 - t);
        else
            weights[i] *= npt[j]==1?0.5:(1.0 + t);
    }
}
 
// Interpolate field
std::vector<double> result(D, 0.0);
for (int i = 0; i < n_corners; ++i)
    for (int j = 0; j < D; ++j) 
    {
      double val ; 
      if (usetimeavg)
        val = (*CGPtemp)[lin_idx[i]].fields[0][idvel + j] ; 
      else 
        val = CGP[lin_idx[i]].fields[0][idvel + j];
      result[j] += val * weights[i];
    }
 
return result ; 
}
 
 
#endif