35 #include <boost/random.hpp>
37 #include <boost/math/special_functions/factorials.hpp>
38 #include <boost/math/special_functions/beta.hpp>
39 #include <boost/crc.hpp>
47 #include "../NrrdIO-1.11.0-src/NrrdIO.h"
56 double Volume (
int d ,
double R) ;
89 int datalocation = -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) ;
129 extrafields.push_back({name, length, extra.size()}) ;
130 extra.resize(extra.size()+length) ;
131 return std::get<2>(extrafields.back());
145 d=dd ; npt=nnpt ; box=bbox ; Time=T ;
147 for (
int i=0 ; i<
d ; i++)
148 dx[i]=((box[1][i]-box[0][i])/
double(npt[i])) ;
151 printf(
"Window and cutoff: %g %g \n", w, cutoff) ;
160 if (CGPtemp !=
nullptr)
delete CGPtemp ; }
168 vector <CGPoint> * CGPtemp = nullptr ;
182 int get_id(
string nm) ;
183 struct Field * get_field(
string nm) ;
184 Pass set_flags (vector <string> s) ;
187 int set_field_struct() ;
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() ;
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);
208 int idx_FastFirst2SlowFirst (
int n) ;
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)) ; } ;
219 int pass_2 (
bool usetimeavg=
false) ;
223 int compute_fluc_vel (
bool usetimeavg=
false) ;
224 int compute_fluc_rot (
bool usetimeavg=
false) ;
225 bool hasvelfluct=
false, hasrotfluct=false ;
232 int mean_time(
bool temporary=
false) ;
233 int write_vtk(
string sout) ;
234 int write_netCDF (
string sout) ;
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) ;
240 std::pair<size_t, uint8_t*> write_numpy_buffer (
int id,
bool squeeze) ;
248 setWindow<W>(w) ;
return 0 ; }
280 printf(
"Unknown window, check Coarsing::setWindow") ;
282 cutoff = Window->
cutoff() ;
283 printf(
"Window and cutoff: %g %g \n", w, cutoff) ;
291 cutoff = Window->
cutoff() ;
292 printf(
"Window and cutoff: %g %g \n", w, cutoff) ;
295 for (
size_t i=0 ; i<per.size() ; i++)
308 const static int idvel=get_id(
"VAVG") ;
310 auto clip = [&](
int a,
int maxd){
if (
a<0)
return(0) ;
else if (
a>=maxd)
return (maxd-1) ;
else return (
a) ; } ;
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]);
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];
330 std::vector<int> lin_idx(n_corners);
331 for (
int c = 0; c < n_corners; ++c) {
333 for (
int dd = 0; dd < D; ++dd)
334 idx += nptcum[dd] * pts[c][dd];
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);
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)
358 val = (*CGPtemp)[lin_idx[i]].fields[0][idvel + j] ;
360 val = CGP[lin_idx[i]].fields[0][idvel + j];
362 result[j] += val * weights[i];
Windows
Definition: WindowLibrary.h:2
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
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