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 vector <double *> orient ;
104 
107 
108 int Ncf ;
109 double * id1 ;
110 double * id2 ;
111 vector <double *> pospq ;
112 vector <double *> lpq ;
113 vector <double *> fpq;
114 vector <double *> mpq;
115 vector <double *> mqp ;
116 
117 // Exta fields if needed
118 vector<double*> extra ;
119 vector<std::tuple<std::string, int, int>> extrafields ;
120 
121 // Some useful functions
122 int Nnonper=-1 ;
123 int random_test (int N, int Ncf, int d, v2d box ) ;
124 int compute_lpq (int d) ;
125 int periodic_atoms (int d, v2d bounds, int pbc, v1d Delta, bool omegainclude) ;
126 int clean_periodic_atoms () {if (Nnonper==-1) printf("ERR: must call periodic_atoms before cleaning the periodic_atoms\n") ; else N=Nnonper ; return 0 ; }
127 bool check_field_availability(string name) ;
128 int add_extra_field (int length, std::string name)
129 {
130  extrafields.push_back({name, length, extra.size()}) ;
131  extra.resize(extra.size()+length) ;
132  return std::get<2>(extrafields.back());
133 }
134 } ;
135 //------------------------------------------------------
136 #include "WindowLibrary.h"
137 
138 //=========================================================
141 class Coarsing
142 {
143 public :
144  Coarsing (int dd, v1i nnpt, v2d bbox, int T) : cT(0), flags(0)
145  {
146  d=dd ; npt=nnpt ; box=bbox ; Time=T ;
147  dx.resize(d, 0) ;
148  for (int i=0 ; i<d ; i++)
149  dx[i]=((box[1][i]-box[0][i])/double(npt[i])) ;
150  double w= (*std::min_element(dx.begin(),dx.end())*2) ; // w automatically set
151  cutoff=2.5*w ; //TODO
152  printf("Window and cutoff: %g %g \n", w, cutoff) ;
153  //for (int i=0 ; i<d ; i++)
154  // printf("%d %d %g %g %g|", d, npt[i], box[1][i], box[0][i], dx[i]) ; fflush(stdout) ;
155  grid_generate() ;
156  //grid_neighbour() ;
157  set_field_struct() ;
158  Window = new LibLucy3D( &data, w, d) ;
159  }
160  ~Coarsing() { if (Window != nullptr) delete Window ;
161  if (CGPtemp != nullptr) delete CGPtemp ; }
162 
163  int d ;
164  int Npt;
165  int Time;
166  int cT ;
167  double cutoff ;
168  vector <CGPoint> CGP ;
169  vector <CGPoint> * CGPtemp = nullptr ;
170  vector <int> npt;
171  vector <int> nptcum ;
172  v1d dx ;
173  v2d box ;
174  LibBase * Window = nullptr ;
175 
176 
177  // Fields variable and function
178  unsigned int flags ;
179  vector <string> Fields, Fname ;
180  vector <int> Fidx ;
181  vector <TensorOrder> Ftype ;
182  vector <struct Field > FIELDS ;
183  int get_id(string nm) ;
184  struct Field * get_field(string nm) ;
185  Pass set_flags (vector <string> s) ;
186 
187  // Grid functions
188  int set_field_struct() ;
189  int add_extra_field(string name, TensorOrder order, FieldType type) ;
190  int setWindow (Windows win, double w, vector <bool> per ={}, vector<int> boxes = {}, vector<double> deltas = {}) ;
191  template <Windows W> int setWindow () ;
192  template <Windows W> int setWindow (double w) ;
193  template <Windows W> int setWindow (double w, vector<bool> per, vector<int> boxes, vector<double> deltas) ;
194  int grid_generate() ;
195  int grid_neighbour() ;
196  std::map<std::string, size_t> grid_setfields() ;
197  vector<FieldType> grid_getfields() ;
198  v2d get_bounds() ;
199  CGPoint * reverseloop (string type) ;
200  int find_closest (int id) ;
201  int find_closest_pq (int id) ;
202  v1d interpolate_vel(int id, bool usetimeavg=false) { return interpolate_vel_nearest (id, usetimeavg) ; }
203  v1d interpolate_rot(int id, bool usetimeavg=false) { return interpolate_rot_nearest (id, usetimeavg) ; }
204  v1d interpolate_vel_nearest (int id, bool usetimeavg=false) ;
205  v1d interpolate_rot_nearest (int id, bool usetimeavg=false) ;
206  v1d interpolate_vel_trilinear (int id, bool usetimeavg) ;
207  template <int D> v1d interpolate_vel_multilinear (int id, bool usetimeavg);
208 
209  int idx_FastFirst2SlowFirst (int n) ;
210 
211  // Windowing functions
212  //double window(double r) {Lucy(r) ; }
213  //double window_int (v1d r1, v1d lpq, v1d x) {printf("Numerical integration of wpqf not implemented\n") ; } ///< Numerical integration: not implemented
214  //double window_int(double r1, double r2) {return window_avg(r1, r2) ; } ///< Overload to avoid integration ...
215  //double window_avg (double r1, double r2) {return (0.5*(Lucy(r1)+Lucy(r2))) ; }
216  //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)) ; }}
217  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)) ; } ;
218  // Coarse graining functions
219  int pass_1 () ;
220  int pass_2 (bool usetimeavg=false) ;
221  int pass_3 () ;
222  int pass_4 () ;
223  int pass_5 () ;
224  int compute_fluc_vel (bool usetimeavg=false) ;
225  int compute_fluc_rot (bool usetimeavg=false) ;
226  bool hasvelfluct=false, hasrotfluct=false ;
227 
228  // Data handling functions
229 
230  struct Data data ;
231 
232  // Time and output handling
233  int mean_time(bool temporary=false) ;
234  int write_vtk(string sout) ;
235  int write_netCDF (string sout) ;
236  int write_NrrdIO (string path) ;
237  int write_matlab (string path, bool squeeze = false) ;
238  int write_numpy (string path, bool squeeze = false) ;
239  int write_numpy_npy (string path, bool squeeze) ;
240  std::pair<size_t, uint8_t*> write_numpy_locbuffer (bool squeeze) {return write_numpy_buffer(-2, squeeze) ; }
241  std::pair<size_t, uint8_t*> write_numpy_buffer (int id, bool squeeze) ;
242 } ;
243 
244 
245 //-------------------------------------------------------
246 template <Windows W>
248 { double w= (*std::min_element(dx.begin(),dx.end())*1) ; // w automatically set
249  setWindow<W>(w) ; return 0 ; }
250 //-------------------------------------------------------
251 template <Windows W>
252 int Coarsing::setWindow (double w)
253 {
254  static_assert(W != Windows::LucyND_Periodic) ;
255  switch (W) {
256  case Windows::Rect3D :
257  Window=new LibRect3D (&data, w, d) ;
258  break ;
260  Window=new LibSphere3DIntersect (&data, w, d) ;
261  break ;
263  Window=new LibSphereNDIntersect (&data, w, d) ;
264  break ;
265  case Windows::Lucy3D :
266  Window=new LibLucy3D (&data, w, d) ;
267  break ;
269  Window=new LibLucy3DFancyInt (&data, w, d) ;
270  break ;
271  case Windows::Hann3D :
272  Window=new LibHann3D (&data, w, d) ;
273  break ;
274  case Windows::RectND :
275  Window=new LibRectND (&data, w, d) ;
276  break ;
277  case Windows::LucyND :
278  Window=new LibLucyND (&data, w, d) ;
279  break ;
280  default:
281  printf("Unknown window, check Coarsing::setWindow") ;
282  }
283  cutoff = Window->cutoff() ;
284  printf("Window and cutoff: %g %g \n", w, cutoff) ;
285  return 0 ;
286 }
287 //-------------------------------------------------------
288 template <Windows W>
289 int Coarsing::setWindow (double w, vector<bool> per, vector<int> boxes, vector<double> deltas)
290 {
291  static_assert(W == Windows::LucyND_Periodic) ;
292  cutoff = Window->cutoff() ;
293  printf("Window and cutoff: %g %g \n", w, cutoff) ;
294 
295  int p = 0 ;
296  for (size_t i=0 ; i<per.size() ; i++)
297  if (per[i])
298  p |= (1<<i) ;
299 
300  Window = new LibLucyND_Periodic (&data,w,d,p,boxes,deltas) ;
301 return 0 ;
302 }
303 //-----------------------------------------------------------
304 //-----------------------------------------------------------------------------------------
305 template <int D>
307 {
308 int pts[1<<D][D] ;
309 const static int idvel=get_id("VAVG") ;
310 
311 auto clip = [&](int a, int maxd){if (a<0) return(0) ; else if (a>=maxd) return (maxd-1) ; else return (a) ; } ;
312 
313 
314 // Determine the floor and ceil indices in each dimension
315 std::vector<int> i0(D), i1(D);
316 for (int dd = 0; dd < D; ++dd) {
317  double val = (data.pos[dd][id] - CGP[0].location[dd]) / dx[dd];
318  i0[dd] = clip(std::floor(val), npt[dd]);
319  i1[dd] = clip(std::ceil(val), npt[dd]);
320 }
321 
322 // Construct all corner indices of the cube
323 int n_corners = 1<<D ;
324 for (int c = 0; c < n_corners; ++c) {
325  for (int dd = 0; dd < D; ++dd) {
326  pts[c][dd] = ((c >> dd) & 1) ? i1[dd] : i0[dd];
327  }
328 }
329 
330 // Convert multi-dim index to linear index
331 std::vector<int> lin_idx(n_corners);
332 for (int c = 0; c < n_corners; ++c) {
333  int idx = 0;
334  for (int dd = 0; dd < D; ++dd)
335  idx += nptcum[dd] * pts[c][dd];
336  lin_idx[c] = idx;
337 }
338 
339 
340 std::vector<double> weights(n_corners, 1.0);
341 for (int i = 0; i < n_corners; ++i) {
342  for (int j = 0; j < D; ++j) {
343  double x0 = CGP[lin_idx[i]].location[j];
344  double t = (data.pos[j][id] - x0) / dx[j];
345  if (((i >> j) & 1) == 0)
346  weights[i] *= (1.0 - t);
347  else
348  weights[i] *= t;
349  }
350 }
351 
352 // Interpolate field
353 std::vector<double> result(D, 0.0);
354 for (int i = 0; i < n_corners; ++i)
355  for (int j = 0; j < D; ++j)
356  {
357  double val ;
358  if (usetimeavg)
359  val = (*CGPtemp)[lin_idx[i]].fields[0][idvel + j] ;
360  else
361  val = CGP[lin_idx[i]].fields[0][idvel + j];
362 
363  result[j] += val * weights[i];
364  }
365 return result ;
366 }
367 
368 
369 #endif
EIGEN_DEVICE_FUNC const FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
EIGEN_DEVICE_FUNC const CeilReturnType ceil() const
Definition: ArrayCwiseUnaryOps.h:495
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Definition: ArrayCwiseUnaryOps.h:187
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:142
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:217
vector< int > Fidx
Where the fields is referenced in the fields vector in the CGPoint. -1 if not flagged.
Definition: Coarsing.h:180
int Npt
Number of coarse graining points.
Definition: Coarsing.h:164
vector< TensorOrder > Ftype
Flagged field types.
Definition: Coarsing.h:181
vector< double * > omega
Particle angular velocity.
Definition: Coarsing.h:102
double cutoff
CG width, and cutoff.
Definition: Coarsing.h:167
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:1793
int d
Dimension.
Definition: Coarsing.h:79
vector< double * > mpq
Moment of particle 1 on 2.
Definition: Coarsing.h:114
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:111
bool operator&(Pass a, Pass b)
Definition: Coarsing.h:66
double * id2
Index of the second particle in contact.
Definition: Coarsing.h:110
vector< double * > fpq
Force at contact.
Definition: Coarsing.h:113
int Time
Total timesteps.
Definition: Coarsing.h:165
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:144
vector< struct Field > FIELDS
All allowed fields (initialized in grid_getfields)
Definition: Coarsing.h:182
vector< double * > pos
Particle positions.
Definition: Coarsing.h:100
v2d box
CG point location.
Definition: Coarsing.h:173
TensorOrder
Definition: Coarsing.h:58
unsigned int flags
Pointer to the averaging window
Definition: Coarsing.h:178
int clean_periodic_atoms()
Clean periodic atoms.
Definition: Coarsing.h:126
int Ncf
Number of contacts.
Definition: Coarsing.h:108
vector< double * > orient
Particle angular velocity.
Definition: Coarsing.h:103
vector< double * > extra
Definition: Coarsing.h:118
int cT
Current timestep.
Definition: Coarsing.h:166
CGPoint(int dd, v1d loc)
Definition: Coarsing.h:73
int d
Number of dimensions.
Definition: Coarsing.h:163
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:109
FieldType ftype
Mainly used to identified the type of user defined fields.
Definition: Coarsing.h:87
~Coarsing()
Definition: Coarsing.h:160
v2d vel_fluc
Fluctuating velocity. Should not be externally provided but calculated, using the function Coarsing::...
Definition: Coarsing.h:105
std::pair< size_t, uint8_t * > write_numpy_locbuffer(bool squeeze)
Definition: Coarsing.h:240
vector< std::tuple< std::string, int, int > > extrafields
Definition: Coarsing.h:119
FieldType
Definition: Coarsing.h:59
double * Imom
Particle moment of inertia.
Definition: Coarsing.h:99
vector< string > Fields
Definition: Coarsing.h:179
vector< double * > lpq
Branch vector of the contact, use compute_lpq() to populate this.
Definition: Coarsing.h:112
v1d interpolate_vel_multilinear(int id, bool usetimeavg)
Definition: Coarsing.h:306
vector< double * > mqp
Moment of particle 2 on 1.
Definition: Coarsing.h:115
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:203
Pass
Definition: Coarsing.h:61
v1d dx
Distances between CG points.
Definition: Coarsing.h:172
int add_extra_field(int length, std::string name)
Definition: Coarsing.h:128
v2d rot_fluc
Fluctuating angular velocity. Should not be externally provided but calculated, using the function Co...
Definition: Coarsing.h:106
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:202
vector< int > nptcum
Cumulated number of points per dimensions (usefull for quick finding of the closest CG for a grain)
Definition: Coarsing.h:171
vector< double * > vel
Particle velocity.
Definition: Coarsing.h:101
vector< int > npt
Number of points per dimension.
Definition: Coarsing.h:170
int setWindow()
Set the windowing function.
Definition: Coarsing.h:247
Data()
Definition: Coarsing.h:95
vector< CGPoint > CGP
List of Coarse Graining points.
Definition: Coarsing.h:168
@ 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