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;
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) ;
130 extrafields.push_back({name, length, extra.size()}) ;
131 extra.resize(extra.size()+length) ;
132 return std::get<2>(extrafields.back());
146 d=dd ; npt=nnpt ; box=bbox ; Time=T ;
148 for (
int i=0 ; i<
d ; i++)
149 dx[i]=((box[1][i]-box[0][i])/
double(npt[i])) ;
152 printf(
"Window and cutoff: %g %g \n", w, cutoff) ;
161 if (CGPtemp !=
nullptr)
delete CGPtemp ; }
169 vector <CGPoint> * CGPtemp = nullptr ;
183 int get_id(
string nm) ;
184 struct Field * get_field(
string nm) ;
185 Pass set_flags (vector <string> s) ;
188 int set_field_struct() ;
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() ;
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);
209 int idx_FastFirst2SlowFirst (
int n) ;
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)) ; } ;
220 int pass_2 (
bool usetimeavg=
false) ;
224 int compute_fluc_vel (
bool usetimeavg=
false) ;
225 int compute_fluc_rot (
bool usetimeavg=
false) ;
226 bool hasvelfluct=
false, hasrotfluct=false ;
233 int mean_time(
bool temporary=
false) ;
234 int write_vtk(
string sout) ;
235 int write_netCDF (
string sout) ;
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) ;
241 std::pair<size_t, uint8_t*> write_numpy_buffer (
int id,
bool squeeze) ;
249 setWindow<W>(w) ;
return 0 ; }
281 printf(
"Unknown window, check Coarsing::setWindow") ;
283 cutoff = Window->
cutoff() ;
284 printf(
"Window and cutoff: %g %g \n", w, cutoff) ;
292 cutoff = Window->
cutoff() ;
293 printf(
"Window and cutoff: %g %g \n", w, cutoff) ;
296 for (
size_t i=0 ; i<per.size() ; i++)
309 const static int idvel=get_id(
"VAVG") ;
311 auto clip = [&](
int a,
int maxd){
if (
a<0)
return(0) ;
else if (
a>=maxd)
return (maxd-1) ;
else return (
a) ; } ;
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];
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];
331 std::vector<int> lin_idx(n_corners);
332 for (
int c = 0; c < n_corners; ++c) {
334 for (
int dd = 0; dd < D; ++dd)
335 idx += nptcum[dd] * pts[c][dd];
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);
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)
359 val = (*CGPtemp)[lin_idx[i]].fields[0][idvel + j] ;
361 val = CGP[lin_idx[i]].fields[0][idvel + j];
363 result[j] += val * weights[i];
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
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
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