NDDEM
Coarsing.h
Go to the documentation of this file.
1 
28 #ifndef COARSING_H
29 #define COARSING_H
30 
31 #include <cstdlib>
32 #include <cstdio>
33 #include <vector>
34 #include "Typedefs.h"
35 #include <boost/random.hpp>
36 #include <fstream>
37 #include <boost/math/special_functions/factorials.hpp>
38 #include <boost/math/special_functions/beta.hpp>
39 #include <boost/crc.hpp>
40 #include <map>
41 
42 #ifdef NETCDF
43 #include <netcdf.h>
44 #endif
45 
46 #ifdef NRRDIO
47 #include "../NrrdIO-1.11.0-src/NrrdIO.h"
48 #endif
49 
50 #ifdef MATLAB
51 #include "mat.h"
52 #endif
53 
54 using namespace std ;
55 
56 double Volume (int d , double R) ;
57 
58 enum TensorOrder {NONE=-1, SCALAR=0, VECTOR=1, TENSOR=2} ;
61 enum Pass {Pass1=1, Pass2=2, Pass3=4, Pass4=8, Pass5=16,
62  VelFluct=256, RotFluct=512} ;
63 
64 inline Pass operator|(Pass a, Pass b){return static_cast<Pass>(static_cast<int>(a) | static_cast<int>(b));}
65 //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
66 inline bool operator& (Pass a, Pass b) { return (static_cast<int>(a) & static_cast<int>(b)) ; }
67 
68 //=========================================================
70 class CGPoint
71 {
72 public :
73  CGPoint(int dd, v1d loc) {d=dd ; location=loc ; }
74 
77  //Useful things
78  vector <int> neighbors ;
79  int d ;
80 } ;
81 //-------------------------
83 struct Field {
85  string name ;
89  int datalocation = -1;
90 };
91 //-------------------------
93 struct Data {
94 public:
95  Data () : radius(nullptr), mass(nullptr), Imom(nullptr), id1(nullptr), id2(nullptr) {}
96 int N = 0 ;
97 double * radius ;
98 double * mass ;
99 double *Imom ;
100 vector <double *> pos ;
101 vector <double *> vel ;
102 vector <double *> omega ;
103 
106 
107 int Ncf ;
108 double * id1 ;
109 double * id2 ;
110 vector <double *> pospq ;
111 vector <double *> lpq ;
112 vector <double *> fpq;
113 vector <double *> mpq;
114 vector <double *> mqp ;
115 
116 // Exta fields if needed
117 vector<double*> extra ;
118 vector<std::tuple<std::string, int, int>> extrafields ;
119 
120 // Some useful functions
121 int Nnonper=-1 ;
122 int random_test (int N, int Ncf, int d, v2d box ) ;
123 int compute_lpq (int d) ;
124 int periodic_atoms (int d, v2d bounds, int pbc, v1d Delta, bool omegainclude) ;
125 int clean_periodic_atoms () {if (Nnonper==-1) printf("ERR: must call periodic_atoms before cleaning the periodic_atoms\n") ; else N=Nnonper ; return 0 ; }
126 bool check_field_availability(string name) ;
127 int add_extra_field (int length, std::string name)
128 {
129  extrafields.push_back({name, length, extra.size()}) ;
130  extra.resize(extra.size()+length) ;
131  return std::get<2>(extrafields.back());
132 }
133 } ;
134 //------------------------------------------------------
135 #include "WindowLibrary.h"
136 
137 //=========================================================
140 class Coarsing
141 {
142 public :
143  Coarsing (int dd, v1i nnpt, v2d bbox, int T) : cT(0), flags(0)
144  {
145  d=dd ; npt=nnpt ; box=bbox ; Time=T ;
146  dx.resize(d, 0) ;
147  for (int i=0 ; i<d ; i++)
148  dx[i]=((box[1][i]-box[0][i])/double(npt[i])) ;
149  double w= (*std::min_element(dx.begin(),dx.end())*2) ; // w automatically set
150  cutoff=2.5*w ; //TODO
151  printf("Window and cutoff: %g %g \n", w, cutoff) ;
152  //for (int i=0 ; i<d ; i++)
153  // printf("%d %d %g %g %g|", d, npt[i], box[1][i], box[0][i], dx[i]) ; fflush(stdout) ;
154  grid_generate() ;
155  //grid_neighbour() ;
156  set_field_struct() ;
157  Window = new LibLucy3D( &data, w, d) ;
158  }
159  ~Coarsing() { if (Window != nullptr) delete Window ;
160  if (CGPtemp != nullptr) delete CGPtemp ; }
161 
162  int d ;
163  int Npt;
164  int Time;
165  int cT ;
166  double cutoff ;
167  vector <CGPoint> CGP ;
168  vector <CGPoint> * CGPtemp = nullptr ;
169  vector <int> npt;
170  vector <int> nptcum ;
171  v1d dx ;
172  v2d box ;
173  LibBase * Window = nullptr ;
174 
175 
176  // Fields variable and function
177  unsigned int flags ;
178  vector <string> Fields, Fname ;
179  vector <int> Fidx ;
180  vector <TensorOrder> Ftype ;
181  vector <struct Field > FIELDS ;
182  int get_id(string nm) ;
183  struct Field * get_field(string nm) ;
184  Pass set_flags (vector <string> s) ;
185 
186  // Grid functions
187  int set_field_struct() ;
188  int add_extra_field(string name, TensorOrder order, FieldType type) ;
189  int setWindow (Windows win, double w, vector <bool> per ={}, vector<int> boxes = {}, vector<double> deltas = {}) ;
190  template <Windows W> int setWindow () ;
191  template <Windows W> int setWindow (double w) ;
192  template <Windows W> int setWindow (double w, vector<bool> per, vector<int> boxes, vector<double> deltas) ;
193  int grid_generate() ;
194  int grid_neighbour() ;
195  std::map<std::string, size_t> grid_setfields() ;
196  vector<FieldType> grid_getfields() ;
197  v2d get_bounds() ;
198  CGPoint * reverseloop (string type) ;
199  int find_closest (int id) ;
200  int find_closest_pq (int id) ;
201  v1d interpolate_vel(int id, bool usetimeavg=false) { return interpolate_vel_nearest (id, usetimeavg) ; }
202  v1d interpolate_rot(int id, bool usetimeavg=false) { return interpolate_rot_nearest (id, usetimeavg) ; }
203  v1d interpolate_vel_nearest (int id, bool usetimeavg=false) ;
204  v1d interpolate_rot_nearest (int id, bool usetimeavg=false) ;
205  v1d interpolate_vel_trilinear (int id, bool usetimeavg) ;
206  template <int D> v1d interpolate_vel_multilinear (int id, bool usetimeavg);
207 
208  int idx_FastFirst2SlowFirst (int n) ;
209 
210  // Windowing functions
211  //double window(double r) {Lucy(r) ; }
212  //double window_int (v1d r1, v1d lpq, v1d x) {printf("Numerical integration of wpqf not implemented\n") ; } ///< Numerical integration: not implemented
213  //double window_int(double r1, double r2) {return window_avg(r1, r2) ; } ///< Overload to avoid integration ...
214  //double window_avg (double r1, double r2) {return (0.5*(Lucy(r1)+Lucy(r2))) ; }
215  //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)) ; }}
216  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)) ; } ;
217  // Coarse graining functions
218  int pass_1 () ;
219  int pass_2 (bool usetimeavg=false) ;
220  int pass_3 () ;
221  int pass_4 () ;
222  int pass_5 () ;
223  int compute_fluc_vel (bool usetimeavg=false) ;
224  int compute_fluc_rot (bool usetimeavg=false) ;
225  bool hasvelfluct=false, hasrotfluct=false ;
226 
227  // Data handling functions
228 
229  struct Data data ;
230 
231  // Time and output handling
232  int mean_time(bool temporary=false) ;
233  int write_vtk(string sout) ;
234  int write_netCDF (string sout) ;
235  int write_NrrdIO (string path) ;
236  int write_matlab (string path, bool squeeze = false) ;
237  int write_numpy (string path, bool squeeze = false) ;
238  int write_numpy_npy (string path, bool squeeze) ;
239  std::pair<size_t, uint8_t*> write_numpy_locbuffer (bool squeeze) {return write_numpy_buffer(-2, squeeze) ; }
240  std::pair<size_t, uint8_t*> write_numpy_buffer (int id, bool squeeze) ;
241 } ;
242 
243 
244 //-------------------------------------------------------
245 template <Windows W>
247 { double w= (*std::min_element(dx.begin(),dx.end())*1) ; // w automatically set
248  setWindow<W>(w) ; return 0 ; }
249 //-------------------------------------------------------
250 template <Windows W>
251 int Coarsing::setWindow (double w)
252 {
253  static_assert(W != Windows::LucyND_Periodic) ;
254  switch (W) {
255  case Windows::Rect3D :
256  Window=new LibRect3D (&data, w, d) ;
257  break ;
259  Window=new LibSphere3DIntersect (&data, w, d) ;
260  break ;
262  Window=new LibSphereNDIntersect (&data, w, d) ;
263  break ;
264  case Windows::Lucy3D :
265  Window=new LibLucy3D (&data, w, d) ;
266  break ;
268  Window=new LibLucy3DFancyInt (&data, w, d) ;
269  break ;
270  case Windows::Hann3D :
271  Window=new LibHann3D (&data, w, d) ;
272  break ;
273  case Windows::RectND :
274  Window=new LibRectND (&data, w, d) ;
275  break ;
276  case Windows::LucyND :
277  Window=new LibLucyND (&data, w, d) ;
278  break ;
279  default:
280  printf("Unknown window, check Coarsing::setWindow") ;
281  }
282  cutoff = Window->cutoff() ;
283  printf("Window and cutoff: %g %g \n", w, cutoff) ;
284  return 0 ;
285 }
286 //-------------------------------------------------------
287 template <Windows W>
288 int Coarsing::setWindow (double w, vector<bool> per, vector<int> boxes, vector<double> deltas)
289 {
290  static_assert(W == Windows::LucyND_Periodic) ;
291  cutoff = Window->cutoff() ;
292  printf("Window and cutoff: %g %g \n", w, cutoff) ;
293 
294  int p = 0 ;
295  for (size_t i=0 ; i<per.size() ; i++)
296  if (per[i])
297  p |= (1<<i) ;
298 
299  Window = new LibLucyND_Periodic (&data,w,d,p,boxes,deltas) ;
300 return 0 ;
301 }
302 //-----------------------------------------------------------
303 //-----------------------------------------------------------------------------------------
304 template <int D>
306 {
307 int pts[1<<D][D] ;
308 const static int idvel=get_id("VAVG") ;
309 
310 auto clip = [&](int a, int maxd){if (a<0) return(0) ; else if (a>=maxd) return (maxd-1) ; else return (a) ; } ;
311 
312 
313 // Determine the floor and ceil indices in each dimension
314 std::vector<int> i0(D), i1(D);
315 for (int dd = 0; dd < D; ++dd) {
316  double val = (data.pos[dd][id] - CGP[0].location[dd]) / dx[dd];
317  i0[dd] = clip(std::floor(val), npt[dd]);
318  i1[dd] = clip(std::ceil(val), npt[dd]);
319 }
320 
321 // Construct all corner indices of the cube
322 int n_corners = 1<<D ;
323 for (int c = 0; c < n_corners; ++c) {
324  for (int dd = 0; dd < D; ++dd) {
325  pts[c][dd] = ((c >> dd) & 1) ? i1[dd] : i0[dd];
326  }
327 }
328 
329 // Convert multi-dim index to linear index
330 std::vector<int> lin_idx(n_corners);
331 for (int c = 0; c < n_corners; ++c) {
332  int idx = 0;
333  for (int dd = 0; dd < D; ++dd)
334  idx += nptcum[dd] * pts[c][dd];
335  lin_idx[c] = idx;
336 }
337 
338 
339 std::vector<double> weights(n_corners, 1.0);
340 for (int i = 0; i < n_corners; ++i) {
341  for (int j = 0; j < D; ++j) {
342  double x0 = CGP[lin_idx[i]].location[j];
343  double t = (data.pos[j][id] - x0) / dx[j];
344  if (((i >> j) & 1) == 0)
345  weights[i] *= (1.0 - t);
346  else
347  weights[i] *= t;
348  }
349 }
350 
351 // Interpolate field
352 std::vector<double> result(D, 0.0);
353 for (int i = 0; i < n_corners; ++i)
354  for (int j = 0; j < D; ++j)
355  {
356  double val ;
357  if (usetimeavg)
358  val = (*CGPtemp)[lin_idx[i]].fields[0][idvel + j] ;
359  else
360  val = CGP[lin_idx[i]].fields[0][idvel + j];
361 
362  result[j] += val * weights[i];
363  }
364 return result ;
365 }
366 
367 
368 #endif
Windows
Definition: WindowLibrary.h:2
@ Lucy3DFancyInt
@ LucyND_Periodic
@ SphereNDIntersect
@ Sphere3DIntersect
Data computed for a single coarse graining point.
Definition: Coarsing.h:71
Main Coarse graining class.
Definition: Coarsing.h:141
A window base class that needs to be specialised to a specific CG window.
Definition: WindowLibrary.h:6
virtual double cutoff(void)
Definition: WindowLibrary.h:26
Definition: WindowLibrary.h:279
Definition: WindowLibrary.h:42
Definition: WindowLibrary.h:29
Definition: WindowLibrary.h:303
Definition: WindowLibrary.h:286
Definition: WindowLibrary.h:81
Definition: WindowLibrary.h:88
Definition: WindowLibrary.h:94
Definition: WindowLibrary.h:177
double normdiff(v1d a, v1d b)
Definition: Coarsing.h:216
vector< int > Fidx
Where the fields is referenced in the fields vector in the CGPoint. -1 if not flagged.
Definition: Coarsing.h:179
int Npt
Number of coarse graining points.
Definition: Coarsing.h:163
vector< TensorOrder > Ftype
Flagged field types.
Definition: Coarsing.h:180
vector< double * > omega
Particle angular velocity.
Definition: Coarsing.h:102
double cutoff
CG width, and cutoff.
Definition: Coarsing.h:166
AverageType
Definition: Coarsing.h:60
Pass operator|(Pass a, Pass b)
Definition: Coarsing.h:64
double Volume(int d, double R)
Compute a sphere volume in dimension D.
Definition: Coarsing.cpp:1785
int d
Dimension.
Definition: Coarsing.h:79
vector< double * > mpq
Moment of particle 1 on 2.
Definition: Coarsing.h:113
double * radius
Particle radius.
Definition: Coarsing.h:97
TensorOrder type
Tensorial order of the field: SCALAR, VECTOR or TENSOR.
Definition: Coarsing.h:86
Pass passlevel
Identify at which moment the field gets calculated.
Definition: Coarsing.h:88
vector< double * > pospq
Location of contact point.
Definition: Coarsing.h:110
bool operator&(Pass a, Pass b)
Definition: Coarsing.h:66
double * id2
Index of the second particle in contact.
Definition: Coarsing.h:109
vector< double * > fpq
Force at contact.
Definition: Coarsing.h:112
int Time
Total timesteps.
Definition: Coarsing.h:164
double * mass
Particle masses.
Definition: Coarsing.h:98
vector< int > neighbors
All the neighbors of the point given the window. 1st index is the point itself.
Definition: Coarsing.h:78
Coarsing(int dd, v1i nnpt, v2d bbox, int T)
Definition: Coarsing.h:143
vector< struct Field > FIELDS
All allowed fields (initialized in grid_getfields)
Definition: Coarsing.h:181
vector< double * > pos
Particle positions.
Definition: Coarsing.h:100
v2d box
CG point location.
Definition: Coarsing.h:172
TensorOrder
Definition: Coarsing.h:58
unsigned int flags
Pointer to the averaging window
Definition: Coarsing.h:177
int clean_periodic_atoms()
Clean periodic atoms.
Definition: Coarsing.h:125
int Ncf
Number of contacts.
Definition: Coarsing.h:107
vector< double * > extra
Definition: Coarsing.h:117
int cT
Current timestep.
Definition: Coarsing.h:165
CGPoint(int dd, v1d loc)
Definition: Coarsing.h:73
int d
Number of dimensions.
Definition: Coarsing.h:162
v1d location
Location of the coarse graining point.
Definition: Coarsing.h:76
uint64_t flag
Flag for the given field.
Definition: Coarsing.h:84
double * id1
Index of first particle in contact.
Definition: Coarsing.h:108
FieldType ftype
Mainly used to identified the type of user defined fields.
Definition: Coarsing.h:87
~Coarsing()
Definition: Coarsing.h:159
v2d vel_fluc
Fluctuating velocity. Should not be externally provided but calculated, using the function Coarsing::...
Definition: Coarsing.h:104
std::pair< size_t, uint8_t * > write_numpy_locbuffer(bool squeeze)
Definition: Coarsing.h:239
vector< std::tuple< std::string, int, int > > extrafields
Definition: Coarsing.h:118
FieldType
Definition: Coarsing.h:59
double * Imom
Particle moment of inertia.
Definition: Coarsing.h:99
vector< string > Fields
Definition: Coarsing.h:178
vector< double * > lpq
Branch vector of the contact, use compute_lpq() to populate this.
Definition: Coarsing.h:111
v1d interpolate_vel_multilinear(int id, bool usetimeavg)
Definition: Coarsing.h:305
vector< double * > mqp
Moment of particle 2 on 1.
Definition: Coarsing.h:114
string name
Name for the given field.
Definition: Coarsing.h:85
v1d interpolate_rot(int id, bool usetimeavg=false)
Interpolate the angular velocity.
Definition: Coarsing.h:202
Pass
Definition: Coarsing.h:61
v1d dx
Distances between CG points.
Definition: Coarsing.h:171
int add_extra_field(int length, std::string name)
Definition: Coarsing.h:127
v2d rot_fluc
Fluctuating angular velocity. Should not be externally provided but calculated, using the function Co...
Definition: Coarsing.h:105
v2d fields
1st dimension is time, second are fields
Definition: Coarsing.h:75
v1d interpolate_vel(int id, bool usetimeavg=false)
Interpolate the velocity.
Definition: Coarsing.h:201
vector< int > nptcum
Cumulated number of points per dimensions (usefull for quick finding of the closest CG for a grain)
Definition: Coarsing.h:170
vector< double * > vel
Particle velocity.
Definition: Coarsing.h:101
vector< int > npt
Number of points per dimension.
Definition: Coarsing.h:169
int setWindow()
Set the windowing function.
Definition: Coarsing.h:246
Data()
Definition: Coarsing.h:95
vector< CGPoint > CGP
List of Coarse Graining points.
Definition: Coarsing.h:167
@ Final
Definition: Coarsing.h:60
@ None
Definition: Coarsing.h:60
@ Both
Definition: Coarsing.h:60
@ Intermediate
Definition: Coarsing.h:60
@ VECTOR
Definition: Coarsing.h:58
@ SCALAR
Definition: Coarsing.h:58
@ NONE
Definition: Coarsing.h:58
@ TENSOR
Definition: Coarsing.h:58
@ Defined
Definition: Coarsing.h:59
@ Particle
Definition: Coarsing.h:59
@ Fluctuation
Definition: Coarsing.h:59
@ Contact
Definition: Coarsing.h:59
@ Pass1
Definition: Coarsing.h:61
@ Pass3
Definition: Coarsing.h:61
@ Pass4
Definition: Coarsing.h:61
@ Pass2
Definition: Coarsing.h:61
@ VelFluct
Definition: Coarsing.h:62
@ RotFluct
Definition: Coarsing.h:62
@ Pass5
Definition: Coarsing.h:61
vector< vector< double > > v2d
Definition: Typedefs.h:10
vector< int > v1i
Definition: Typedefs.h:18
@ Imom
Definition: Typedefs.h:19
@ id2
Definition: Typedefs.h:19
@ id1
Definition: Typedefs.h:19
@ radius
Definition: Typedefs.h:19
@ mass
Definition: Typedefs.h:19
vector< double > v1d
Definition: Typedefs.h:9
uint d
int N
int write_NrrdIO(string path, int d, vector< vector< float >> &colors)
Writer for NRRD colormaps.
Definition: io.cpp:116
type
The type the bitset is encoded with.
Definition: bitset.hpp:44
Definition: json.hpp:5678
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1181
PUGI_IMPL_FN I min_element(I begin, I end, const Pred &pred)
Definition: pugixml.cpp:7604
unsigned __int64 uint64_t
Definition: stdint.h:136
Data structure handling point data and contact data.
Definition: Coarsing.h:93
Contains Field informations.
Definition: Coarsing.h:83