NDDEM
WindowLibrary.h
Go to the documentation of this file.
1 
2 enum class Windows {Rect3D, Sphere3DIntersect, SphereNDIntersect, Lucy3D, Lucy3DFancyInt, Hann3D, RectND, LucyND, LucyND_Periodic} ;
3 
6 class LibBase {
7 public:
8  LibBase(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; }
9  LibBase() {}
10  virtual ~LibBase() {}
11 
12  struct Data * data ;
13  double w, d ;
14 
15  virtual double window (double) = 0 ; //Purely virtual, need to be defined in derived classes
16  virtual std::pair<double,double> window_contact_weight (int p, int q, const v1d & loc) {
17  double rp=distance(p, loc) ;
18  double rq=distance(q, loc) ;
19  double wpqs = window_avg (rp, rq) ;
20  double wpqf = window_int (rp, rq) ;
21  return (make_pair(wpqs, wpqf)) ;
22  }
23  virtual double window_int(double r1, double r2) {return window_avg(r1, r2) ; }
24  virtual double window_avg (double r1, double r2) {return (0.5*(window(r1)+window(r2))) ; }
25  virtual double distance (int id, v1d loc) {double res=0 ; for (int i=0 ; i<d ; i++) res+=(loc[i]-data->pos[i][id])*(loc[i]-data->pos[i][id]) ; return sqrt(res) ; }
26  virtual double cutoff (void) {return 2.5*w ;}
27 } ;
28 //---------------------
29 class LibLucy3D : public LibBase {
30 public :
31  LibLucy3D(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; cst=105./(16*M_PI*w*w*w) ; }
32  LibLucy3D() {} ;
33  double Lucy (double r) {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)) ; }}
34  double window(double r) { return (Lucy(r)) ;}
35  double cutoff (void) {return 2*w ;} // added by benjy to get rid of extra (hopefully useless) data
36  //double window_int(double r1, double r2) {return window_avg(r1, r2) ; }
37  //double window_avg (double r1, double r2) {return (0.5*(Lucy(r1)+Lucy(r2))) ; }
38  protected:
39  double cst ;
40 };
41 //--------------------
42 class LibLucy3DFancyInt : public LibLucy3D {
43 public :
44  LibLucy3DFancyInt (struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; cst=105./(16*M_PI*w*w*w) ;}
45  std::pair<double,double> window_contact_weight (int p, int q, const v1d & loc)
46  {
47  double res = 0 ;
48  double prev, cur, first ;
49  v1d locp (d,0) ;
50  v1d lpq(d,0) ;
51  for (int i=0 ; i<d ;i++)
52  {
53  lpq[i] = data->pos[i][q] - data->pos[i][p] ;
54  locp[i] = data->pos[i][p] ;
55  }
56 
57  //double length = 0 ; for (int i=0 ; i<d ;i++) length += lpq[i]*lpq[i] ; length = sqrt(length)/Nsteps ;
58 
59  first=prev=cur=Lucy(distance(p,loc)) ;
60  for (int i=0 ; i<Nsteps ; i++)
61  {
62  for (int j=0; j<d ; j++)
63  locp[j] += (1./Nsteps) * lpq[j] ;
64  cur=Lucy(distancevec(locp,loc)) ;
65  res += (cur+prev)/(2*Nsteps) ;
66  //printf("%g | %g | %g %g %g\n", res, cur, distance(q,locp), Lucy(distance(p,loc)), Lucy(distance(q,loc))) ;
67  prev=cur ;
68  }
69  //printf("%g %g \n", res, LibLucy3D::window_int(distance(p,loc),distance(q,loc))) ;
70  return (make_pair(window_avg(first,cur),res)) ;
71  }
72  double distancevec (v1d l1, v1d loc) {double res=0 ; for (int i=0 ; i<d ; i++) res+=(loc[i]-l1[i])*(loc[i]-l1[i]) ; return sqrt(res) ;}
73  void set_integrationsteps (int steps) {if (steps<1) printf("Less than 1 step is not meaningful.") ; if (steps>1000) printf("You've chosen a very large number of integration steps, you may want to reconsider") ; Nsteps=steps ; }
74 private:
75  int Nsteps = 10 ;
76 };
77 
78 
79 
80 //---------------------
81 class LibRect3D : public LibBase {
82 public:
83  LibRect3D(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; cst= 1./(4./3.*M_PI*w*w*w) ; }
84  double cst ;
85  double window (double r) {if (r>=w) return 0 ; else return cst ;}
86 };
87 //---------------------
88 class LibRectND : public LibBase {
89 public:
90  LibRectND(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; }
91  double window (double r) {if (r>=w) return 0 ; else {double a =1 ; for (int b=0 ; b<d ; b++,a*=w) {} return 1/a ;}}
92 };
93 //---------------------
94 class LibSphere3DIntersect : public LibBase {
95 public:
96  LibSphere3DIntersect(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; cst=1./(4./3.*M_PI*w*w*w) ; }
97  double cst ;
98  double result ;
99  double distance(int id, v1d loc)
100  {
101  double dst=0 ;
102  for (int i=0 ; i<d ; i++)
103  dst+=(loc[i]-data->pos[i][id])*(loc[i]-data->pos[i][id]) ;
104  dst= sqrt(dst) ;
105  double r = data->radius[id] ;
106  if (dst>w+r) return 0 ;
107  else if (dst<=fabs(w-r))
108  {
109  if (w>r)
110  return cst ;
111  else
112  return (1/(4./3. * M_PI * r*r*r)) ;
113  }
114  else
115  {
116  double vol = M_PI * (w+r-dst) * (w+r-dst) * (dst*dst + 2*dst*r - 3*r*r + 2*dst*w + 6*r*w - 3*w*w) / (12.*dst) ;
117  return (vol * cst / (4./3.*M_PI*r*r*r)) ;
118  }
119 
120  }
121  double cutoff (void) {
122  double maxr=0 ;
123  if (data->N==0) return 2*w ;
124  for (int i=0 ; i<data -> N ; i++)
125  if (maxr<data->radius[i])
126  maxr=data->radius[i] ;
127  return maxr+w ;
128  }
129  double window (double r) {return r ; } // The value calculated by the distance measurement is the one returned, a bit of a hack there ...
130  double windowreal (double r) {if (r>=w) return 0 ; else return cst ;}
131  virtual std::pair<double,double> window_contact_weight (int p, int q, const v1d & loc)
132  {
133  double rp=LibBase::distance(p, loc) ;
134  double rq=LibBase::distance(q, loc) ;
135  double wpqs = window_avg (rp, rq) ;
136  double wpqf = 0 ;
137 
138  v1d locp (d,0) ;
139  v1d lpq(d,0) ;
140  double normp=0, normlpq=0 ;
141  double b = 0 ;
142  for (int i=0 ; i<d ;i++)
143  {
144  lpq[i] = data->pos[i][q] - data->pos[i][p] ;
145  locp[i] = data->pos[i][p]-loc[i] ;
146  normp += locp[i]*locp[i] ;
147  normlpq += lpq[i]*lpq[i] ;
148  b += locp[i]*lpq[i] ;
149  }
150  b*=2 ;
151 
152  double Delta = b*b-4*normlpq*(normp-w*w) ;
153  if (Delta<=0)
154  wpqs = 0 ;
155  else
156  {
157  double alpha1 = (-b - sqrt(Delta))/(2*normlpq) ;
158  double alpha2 = (-b + sqrt(Delta))/(2*normlpq) ;
159 
160  if (alpha1<0 && alpha2>0 && alpha2<1)
161  wpqf = cst * alpha2 ;
162  else if (alpha1<0 && alpha2>1)
163  wpqf = cst ;
164  else if (alpha1>0 && alpha1<1 && alpha2>1)
165  wpqf = cst * (1-alpha1) ;
166  else if (alpha1>0 && alpha1<1 && alpha2>0 && alpha2<1) //this is a VERY small window ... ...
167  wpqf = cst*(alpha2-alpha1) ;
168  else
169  wpqf = 0 ;
170 
171  }
172  return (make_pair(wpqs,wpqf)) ;
173  }
174  virtual double window_avg (double r1, double r2) {return (0.5*(windowreal(r1)+windowreal(r2))) ; }
175 };
176 //------------------------------------------
177 class LibSphereNDIntersect : public LibBase {
178 public:
179  LibSphereNDIntersect(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; cst=1./(nballvolume()* pow(w, d)) ; }
180  double nballvolume () {return pow(M_PI,d/2.)/tgamma(d/2.+1) ; }
181  double ncapvolume (double r, double h)
182  {
183  if (h<0 || h>2*r) return 0 ;
184  if (h<=r)
185  {
186  //printf("%g %g\n", 1/2. * nballvolume() * pow(r, d) * boost::math::ibeta ((d+1)/2., 1/2., (2*r*h-h*h)/(r*r)), M_PI*h*h/3*(3*r-h)) ;
187  return 1/2. * nballvolume() * pow(r, d) * boost::math::ibeta ((d+1)/2., 1/2., (2*r*h-h*h)/(r*r)) ;
188  }
189  else
190  {
191  //printf("%g %g %g %g\n", h, r, nballvolume() * pow(r, d) - 1/2. * nballvolume() * pow(r, d) * (boost::math::ibeta((d+1)/2., 1/2., (2*r*(2*r-h)-(2*r-h)*(2*r-h))/(r*r))),
192  // 4/3.*M_PI*r*r*r - M_PI*(2*r-h)*(2*r-h)/3*(3*r-(2*r-h))) ;
193  return nballvolume() * pow(r, d) - 1/2. * nballvolume() * pow(r, d) * boost::math::ibeta((d+1)/2., 1/2., (2*r*(2*r-h)-(2*r-h)*(2*r-h))/(r*r)) ; // Completement normalised beta function (1-I_x(a,b))
194  }
195  }
196  double cst ;
197  double result ;
198  double distance(int id, v1d loc)
199  {
200  double r=0 ;
201  for (int i=0 ; i<d ; i++)
202  r+=(loc[i]-data->pos[i][id])*(loc[i]-data->pos[i][id]) ;
203  r = sqrt(r) ;
204  double r2 = data->radius[id] ;
205  double &r1 = w ;
206  if (r>r1+r2) return 0 ;
207  else if (r<=fabs(r1-r2))
208  {
209  if (r1>r2)
210  return cst ;
211  else
212  return 1./(nballvolume()*pow(r2, d)) ;
213  }
214  else
215  {
216  double h1 = (r2+r1-r)*(r2-r1+r)/(2*r) ;
217  double h2 = (r1+r2-r)*(r1-r2+r)/(2*r) ;
218 
219  double vol = ncapvolume(r1,h1) + ncapvolume(r2,h2) ;
220  return vol * cst / (nballvolume() * pow(r2,d)) ;
221  }
222  }
223  double cutoff (void) {
224  double maxr=0 ;
225  if (data->N==0) return 2*w ;
226  for (int i=0 ; i<data -> N ; i++)
227  if (maxr<data->radius[i])
228  maxr=data->radius[i] ;
229  return maxr+w ;
230  }
231  double window (double r) {return r ; } // The value calculated by the distance measurement is the one returned, a bit of a hack there ...
232  double windowreal (double r) {if (r>=w) return 0 ; else return cst ;}
233  virtual std::pair<double,double> window_contact_weight (int p, int q, const v1d & loc)
234  {
235  double rp=LibBase::distance(p, loc) ;
236  double rq=LibBase::distance(q, loc) ;
237  double wpqs = window_avg (rp, rq) ;
238  double wpqf = 0 ;
239 
240  v1d locp (d,0) ;
241  v1d lpq(d,0) ;
242  double normp=0, normlpq=0 ;
243  double b = 0 ;
244  for (int i=0 ; i<d ;i++)
245  {
246  lpq[i] = data->pos[i][q] - data->pos[i][p] ;
247  locp[i] = data->pos[i][p]-loc[i] ;
248  normp += locp[i]*locp[i] ;
249  normlpq += lpq[i]*lpq[i] ;
250  b += locp[i]*lpq[i] ;
251  }
252  b*=2 ;
253 
254  double Delta = b*b-4*normlpq*(normp-w*w) ;
255  if (Delta<=0)
256  wpqs = 0 ;
257  else
258  {
259  double alpha1 = (-b - sqrt(Delta))/(2*normlpq) ;
260  double alpha2 = (-b + sqrt(Delta))/(2*normlpq) ;
261 
262  if (alpha1<0 && alpha2>0 && alpha2<1)
263  wpqf = cst * alpha2 ;
264  else if (alpha1<0 && alpha2>1)
265  wpqf = cst ;
266  else if (alpha1>0 && alpha1<1 && alpha2>1)
267  wpqf = cst * (1-alpha1) ;
268  else if (alpha1>0 && alpha1<1 && alpha2>0 && alpha2<1) //this is a VERY small window ... ...
269  wpqf = cst*(alpha2-alpha1) ;
270  else
271  wpqf = 0 ;
272 
273  }
274  return (make_pair(wpqs,wpqf)) ;
275  }
276  virtual double window_avg (double r1, double r2) {return (0.5*(windowreal(r1)+windowreal(r2))) ; }
277 };
278 //---------------------
279 class LibHann3D : public LibBase {
280 public:
281  LibHann3D(struct Data * D, double ww, double dd) { data=D; w=ww ; d=dd ; cst= 3*M_PI/(2.*(M_PI*M_PI-6)*w*w*w) ; }
282  double cst ;
283  double window (double r) {if (r>=w) return 0 ; else return cst*pow(cos(M_PI*r/(2.*w)),2) ;}
284 } ;
285 //---------------------
286 class LibLucyND : public LibBase {
287 public:
288  LibLucyND(struct Data * D, double ww, double dd)
289  {
290  data=D; w=ww ; d=dd ;
291  double Vol=pow(M_PI,d/2.)/(tgamma(d/2.+1)) ; // N-ball volume
292  scale = Vol * d * (-3./(d+4) + 8./(d+3) - 6./(d+2) + 1./d) ;
293  scale = 1/scale ;
294  scale = scale / (pow(w, d)) ;
295  printf("Lucy function scaling : %f \n", scale) ;
296  }
297  double scale ;
298  double Lucy (double r) {if (r>=w) return 0 ; else {double f=r/w ; return (scale*(-3*f*f*f*f + 8*f*f*f - 6*f*f +1)) ; }}
299  double window(double r) override {return (Lucy(r)) ;}
300 } ;
301 //---------------------
302 class LibLucyND_Periodic : public LibLucyND
303 {
304  public :
305  LibLucyND_Periodic(struct Data * D, double ww, double dd, int periodic, vector<int>boxes, vector<double> deltas)
306  : LibLucyND(D,ww,dd-__builtin_popcount(periodic))
307  {
308  printf("%g ", scale) ;
309  maskperiodic = 0 ;
310  for (size_t i =0 ; i<boxes.size() ; i++)
311  if ((periodic&(1<<i)) && boxes[i]==1)
312  {
313  maskperiodic |= (1<<i) ;
314  scale /= deltas[i] ;
315  }
316  printf("NB: do not use periodic_atoms with LibLucyND_Periodic! %g \n", scale) ;
317  }
318 
319  double distance (int id, v1d loc) override
320  {
321  double res=0 ;
322  for (int i=0 ; i<d ; i++)
323  if ((maskperiodic&(1<<i))==0)
324  res+=(loc[i]-data->pos[i][id])*(loc[i]-data->pos[i][id]) ;
325  return sqrt(res) ;
326  }
327 
328  int maskperiodic ;
329 };
sqrt
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Definition: ArrayCwiseUnaryOps.h:187
cos
EIGEN_DEVICE_FUNC const CosReturnType cos() const
Definition: ArrayCwiseUnaryOps.h:237
Windows
Windows
Definition: WindowLibrary.h:2
Windows::Lucy3DFancyInt
@ Lucy3DFancyInt
Windows::LucyND_Periodic
@ LucyND_Periodic
Windows::SphereNDIntersect
@ SphereNDIntersect
Windows::Sphere3DIntersect
@ Sphere3DIntersect
LibBase
A window base class that needs to be specialised to a specific CG window.
Definition: WindowLibrary.h:6
LibBase::window_contact_weight
virtual std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:16
LibBase::~LibBase
virtual ~LibBase()
Definition: WindowLibrary.h:10
LibBase::window
virtual double window(double)=0
LibBase::distance
virtual double distance(int id, v1d loc)
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:25
LibBase::window_int
virtual double window_int(double r1, double r2)
Definition: WindowLibrary.h:23
LibBase::window_avg
virtual double window_avg(double r1, double r2)
Definition: WindowLibrary.h:24
LibBase::data
struct Data * data
Definition: WindowLibrary.h:12
LibBase::d
double d
Definition: WindowLibrary.h:13
LibBase::LibBase
LibBase()
Definition: WindowLibrary.h:9
LibBase::cutoff
virtual double cutoff(void)
Definition: WindowLibrary.h:26
LibBase::LibBase
LibBase(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:8
LibBase::w
double w
Definition: WindowLibrary.h:13
LibHann3D
Definition: WindowLibrary.h:279
LibHann3D::window
double window(double r)
Definition: WindowLibrary.h:283
LibHann3D::LibHann3D
LibHann3D(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:281
LibHann3D::cst
double cst
Definition: WindowLibrary.h:282
LibLucy3DFancyInt
Definition: WindowLibrary.h:42
LibLucy3DFancyInt::LibLucy3DFancyInt
LibLucy3DFancyInt(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:44
LibLucy3DFancyInt::Nsteps
int Nsteps
Definition: WindowLibrary.h:75
LibLucy3DFancyInt::window_contact_weight
std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:45
LibLucy3DFancyInt::distancevec
double distancevec(v1d l1, v1d loc)
Definition: WindowLibrary.h:72
LibLucy3DFancyInt::set_integrationsteps
void set_integrationsteps(int steps)
Definition: WindowLibrary.h:73
LibLucy3D
Definition: WindowLibrary.h:29
LibLucy3D::window
double window(double r)
Definition: WindowLibrary.h:34
LibLucy3D::cutoff
double cutoff(void)
Definition: WindowLibrary.h:35
LibLucy3D::cst
double cst
Definition: WindowLibrary.h:39
LibLucy3D::LibLucy3D
LibLucy3D(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:31
LibLucy3D::Lucy
double Lucy(double r)
Definition: WindowLibrary.h:33
LibLucy3D::LibLucy3D
LibLucy3D()
Definition: WindowLibrary.h:32
LibLucyND_Periodic
Definition: WindowLibrary.h:303
LibLucyND_Periodic::LibLucyND_Periodic
LibLucyND_Periodic(struct Data *D, double ww, double dd, int periodic, vector< int >boxes, vector< double > deltas)
Definition: WindowLibrary.h:305
LibLucyND_Periodic::maskperiodic
int maskperiodic
Definition: WindowLibrary.h:328
LibLucyND_Periodic::distance
double distance(int id, v1d loc) override
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:319
LibLucyND
Definition: WindowLibrary.h:286
LibLucyND::scale
double scale
Definition: WindowLibrary.h:297
LibLucyND::Lucy
double Lucy(double r)
Definition: WindowLibrary.h:298
LibLucyND::window
double window(double r) override
Definition: WindowLibrary.h:299
LibLucyND::LibLucyND
LibLucyND(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:288
LibRect3D
Definition: WindowLibrary.h:81
LibRect3D::window
double window(double r)
Definition: WindowLibrary.h:85
LibRect3D::cst
double cst
Definition: WindowLibrary.h:84
LibRect3D::LibRect3D
LibRect3D(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:83
LibRectND
Definition: WindowLibrary.h:88
LibRectND::LibRectND
LibRectND(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:90
LibRectND::window
double window(double r)
Definition: WindowLibrary.h:91
LibSphere3DIntersect
Definition: WindowLibrary.h:94
LibSphere3DIntersect::cst
double cst
Definition: WindowLibrary.h:97
LibSphere3DIntersect::window_avg
virtual double window_avg(double r1, double r2)
Definition: WindowLibrary.h:174
LibSphere3DIntersect::LibSphere3DIntersect
LibSphere3DIntersect(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:96
LibSphere3DIntersect::windowreal
double windowreal(double r)
Definition: WindowLibrary.h:130
LibSphere3DIntersect::window
double window(double r)
Definition: WindowLibrary.h:129
LibSphere3DIntersect::cutoff
double cutoff(void)
Definition: WindowLibrary.h:121
LibSphere3DIntersect::window_contact_weight
virtual std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:131
LibSphere3DIntersect::result
double result
Definition: WindowLibrary.h:98
LibSphere3DIntersect::distance
double distance(int id, v1d loc)
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:99
LibSphereNDIntersect
Definition: WindowLibrary.h:177
LibSphereNDIntersect::nballvolume
double nballvolume()
Definition: WindowLibrary.h:180
LibSphereNDIntersect::cutoff
double cutoff(void)
Definition: WindowLibrary.h:223
LibSphereNDIntersect::LibSphereNDIntersect
LibSphereNDIntersect(struct Data *D, double ww, double dd)
Definition: WindowLibrary.h:179
LibSphereNDIntersect::window_contact_weight
virtual std::pair< double, double > window_contact_weight(int p, int q, const v1d &loc)
Definition: WindowLibrary.h:233
LibSphereNDIntersect::windowreal
double windowreal(double r)
Definition: WindowLibrary.h:232
LibSphereNDIntersect::window_avg
virtual double window_avg(double r1, double r2)
Definition: WindowLibrary.h:276
LibSphereNDIntersect::cst
double cst
Definition: WindowLibrary.h:196
LibSphereNDIntersect::distance
double distance(int id, v1d loc)
function for mixed particle id / vector informations.
Definition: WindowLibrary.h:198
LibSphereNDIntersect::window
double window(double r)
Definition: WindowLibrary.h:231
LibSphereNDIntersect::result
double result
Definition: WindowLibrary.h:197
LibSphereNDIntersect::ncapvolume
double ncapvolume(double r, double h)
Definition: WindowLibrary.h:181
Data::radius
double * radius
Particle radius.
Definition: Coarsing.h:97
Data::N
int N
Number of particles.
Definition: Coarsing.h:96
Data::pos
vector< double * > pos
Particle positions.
Definition: Coarsing.h:100
lpq
@ lpq
Definition: Typedefs.h:19
radius
@ radius
Definition: Typedefs.h:19
v1d
vector< double > v1d
Definition: Typedefs.h:9
N
int N
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:514
Eigen::internal::first
EIGEN_CONSTEXPR Index first(const T &x) EIGEN_NOEXCEPT
Definition: IndexedViewHelper.h:81
a
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1181
Data
Data structure handling point data and contact data.
Definition: Coarsing.h:93
Generated by doxygen 1.9.1