1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35
36 Christoph Kloss (DCS Computing GmbH, Linz)
37 Christoph Kloss (JKU Linz)
38 Philippe Seil (JKU Linz)
39 Andreas Aigner (DCS Computing GmbH, Linz)
40
41 Copyright 2012- DCS Computing GmbH, Linz
42 Copyright 2009-2012 JKU Linz
43 ------------------------------------------------------------------------- */
44
45 #ifndef LMP_GENERAL_CONTAINER_I_H
46 #define LMP_GENERAL_CONTAINER_I_H
47
48 /* ----------------------------------------------------------------------
49 constructors
50 ------------------------------------------------------------------------- */
51
52 template<typename T, int NUM_VEC, int LEN_VEC>
GeneralContainer(const char * _id)53 GeneralContainer<T,NUM_VEC,LEN_VEC>::GeneralContainer(const char *_id)
54 : ContainerBase(_id),
55 numElem_(0),
56 maxElem_(GROW_CONTAINER()),
57 defaultValue_(0)
58 {
59 create<T>(arr_,GROW_CONTAINER(),NUM_VEC,LEN_VEC);
60 }
61
62 template<typename T, int NUM_VEC, int LEN_VEC>
GeneralContainer(const char * _id,const char * _comm,const char * _ref,const char * _restart,int _scalePower)63 GeneralContainer<T,NUM_VEC,LEN_VEC>::GeneralContainer(const char *_id, const char *_comm, const char *_ref, const char *_restart, int _scalePower)
64 : ContainerBase(_id, _comm, _ref, _restart, _scalePower),
65 numElem_(0),
66 maxElem_(GROW_CONTAINER()),
67 defaultValue_(0)
68 {
69 create<T>(arr_,GROW_CONTAINER(),NUM_VEC,LEN_VEC);
70 }
71
72 template<typename T, int NUM_VEC, int LEN_VEC>
GeneralContainer(GeneralContainer<T,NUM_VEC,LEN_VEC> const & orig)73 GeneralContainer<T,NUM_VEC,LEN_VEC>::GeneralContainer(GeneralContainer<T,NUM_VEC,LEN_VEC> const &orig)
74 : ContainerBase(orig),
75 numElem_(orig.numElem_),
76 maxElem_(orig.numElem_),
77 defaultValue_(orig.defaultValue_)
78 {
79 create<T>(arr_,maxElem_,NUM_VEC,LEN_VEC);
80 for(int i=0;i<maxElem_;i++)
81 for(int ii=0;ii<NUM_VEC;ii++)
82 for(int jj=0;jj<LEN_VEC;jj++)
83 arr_[i][ii][jj] = orig.arr_[i][ii][jj];
84 }
85
86 /* ----------------------------------------------------------------------
87 destructor
88 ------------------------------------------------------------------------- */
89
90 template<typename T, int NUM_VEC, int LEN_VEC>
~GeneralContainer()91 GeneralContainer<T,NUM_VEC,LEN_VEC>::~GeneralContainer()
92 {
93 destroy<T>(arr_);
94 }
95
96 /* ----------------------------------------------------------------------
97 check if data is of type double
98 ------------------------------------------------------------------------- */
99
100 template<typename T, int NUM_VEC, int LEN_VEC>
isDoubleData()101 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::isDoubleData()
102 {
103 // partial templatization does not work
104 // std::is_same<T,double>::value is from C++11
105 // this is work-around
106
107 if(sizeof(T) == sizeof(double))
108 return true;
109 else
110 return false;
111 }
112
113 template<typename T, int NUM_VEC, int LEN_VEC>
isIntData()114 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::isIntData()
115 {
116 // partial templatization does not work
117 // std::is_same<T,double>::value is from C++11
118 // this is work-around
119
120 if(sizeof(T) == sizeof(int))
121 return true;
122 else
123 return false;
124 }
125
126 /* ----------------------------------------------------------------------
127 add element(s)
128 ------------------------------------------------------------------------- */
129
130 template<typename T, int NUM_VEC, int LEN_VEC>
subtract(GeneralContainer<T,NUM_VEC,LEN_VEC> const & A,GeneralContainer<T,NUM_VEC,LEN_VEC> const & minusB)131 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::subtract(GeneralContainer<T,NUM_VEC,LEN_VEC> const &A,
132 GeneralContainer<T,NUM_VEC,LEN_VEC> const &minusB)
133 {
134 int len = size();
135 int lenA = A.size();
136 int lenB = minusB.size();
137
138 if(lenA != lenB)
139 return false;
140
141 if(len < lenA)
142 addUninitialized(lenA-len);
143
144 for(int i = 0; i < len; i++)
145 for(int j = 0; j < NUM_VEC; j++)
146 for(int k = 0; k < LEN_VEC; k++)
147 arr_[i][j][k] = A(i)[j][k] - minusB(i)[j][k];
148
149 return true;
150 }
151
152 /* ----------------------------------------------------------------------
153 add element(s)
154 ------------------------------------------------------------------------- */
155
156 template<typename T, int NUM_VEC, int LEN_VEC>
add(T ** elem)157 void GeneralContainer<T,NUM_VEC,LEN_VEC>::add(T** elem)
158 {
159 if(numElem_ == maxElem_)
160 {
161 grow<T>(arr_,maxElem_+GROW_CONTAINER(),NUM_VEC,LEN_VEC);
162 maxElem_ += GROW_CONTAINER();
163 }
164 for(int i=0;i<NUM_VEC;i++)
165 for(int j=0;j<LEN_VEC;j++)
166 arr_[numElem_][i][j] = elem[i][j];
167 numElem_++;
168 }
169
170 template<typename T, int NUM_VEC, int LEN_VEC>
addZero()171 void GeneralContainer<T,NUM_VEC,LEN_VEC>::addZero()
172 {
173 if(numElem_ == maxElem_)
174 {
175 grow<T>(arr_,maxElem_+GROW_CONTAINER(),NUM_VEC,LEN_VEC);
176 maxElem_ += GROW_CONTAINER();
177 }
178 for(int i=0;i<NUM_VEC;i++)
179 for(int j=0;j<LEN_VEC;j++)
180 arr_[numElem_][i][j] = static_cast<T>(0);
181 numElem_++;
182 }
183
184 template<typename T, int NUM_VEC, int LEN_VEC>
addUninitialized(int n)185 void GeneralContainer<T,NUM_VEC,LEN_VEC>::addUninitialized(int n)
186 {
187 numElem_ += n;
188 if(numElem_ >= maxElem_)
189 {
190 T init_val = static_cast<T>(0);
191 grow(arr_,numElem_+GROW_CONTAINER(),NUM_VEC,LEN_VEC);
192 for(int i = numElem_; i < numElem_+GROW_CONTAINER(); i++)
193 for(int j=0;j<NUM_VEC;j++)
194 for(int k=0;k<LEN_VEC;k++)
195 arr_[i][j][k] = init_val;
196 maxElem_ = numElem_ + GROW_CONTAINER();
197 }
198 }
199
200 /* ----------------------------------------------------------------------
201 delete an element
202 ------------------------------------------------------------------------- */
203
204 template<typename T, int NUM_VEC, int LEN_VEC>
del(int n)205 void GeneralContainer<T,NUM_VEC,LEN_VEC>::del(int n)
206 {
207
208 numElem_--;
209 if(numElem_ == n) return;
210 for(int i=0;i<NUM_VEC;i++)
211 for(int j=0;j<LEN_VEC;j++)
212 arr_[n][i][j] = arr_[numElem_][i][j];
213 }
214
215 /* ----------------------------------------------------------------------
216 copy element data
217 ------------------------------------------------------------------------- */
218
219 template<typename T, int NUM_VEC, int LEN_VEC>
copy(int from,int to)220 void GeneralContainer<T,NUM_VEC,LEN_VEC>::copy(int from,int to)
221 {
222 for(int i=0;i<NUM_VEC;i++)
223 for(int j=0;j<LEN_VEC;j++)
224 arr_[to][i][j] = arr_[from][i][j];
225 }
226
227 /* ----------------------------------------------------------------------
228 delete an element
229 ------------------------------------------------------------------------- */
230
231 template<typename T, int NUM_VEC, int LEN_VEC>
delForward(int n,bool scale,bool translate,bool rotate)232 void GeneralContainer<T,NUM_VEC,LEN_VEC>::delForward(int n,bool scale,bool translate,bool rotate)
233 {
234 // do only delete property if it is a forward comm property
235 if(!decidePackUnpackOperation(OPERATION_COMM_FORWARD, scale, translate, rotate))
236 return;
237
238 numElem_--;
239 if(numElem_ == n) return;
240 for(int i=0;i<NUM_VEC;i++)
241 for(int j=0;j<LEN_VEC;j++)
242 arr_[n][i][j] = arr_[numElem_][i][j];
243 }
244
245 /* ----------------------------------------------------------------------
246 clear reverse properties, i.e. reset all of them to 0
247 ------------------------------------------------------------------------- */
248
249 template<typename T, int NUM_VEC, int LEN_VEC>
clearReverse(bool scale,bool translate,bool rotate)250 void GeneralContainer<T,NUM_VEC,LEN_VEC>::clearReverse(bool scale,bool translate,bool rotate)
251 {
252 // do only reset property if it is a reverse comm property
253 if(!decidePackUnpackOperation(OPERATION_COMM_REVERSE, scale, translate, rotate))
254 return;
255
256 int len = size();
257 for(int i = 0; i < len; i++)
258 for(int j = 0; j < NUM_VEC; j++)
259 for(int k = 0; k < LEN_VEC; k++)
260 arr_[i][j][k] = 0.;
261 }
262
263 /* ----------------------------------------------------------------------
264 delete an element if restart
265 ------------------------------------------------------------------------- */
266
267 template<typename T, int NUM_VEC, int LEN_VEC>
delRestart(int n,bool scale,bool translate,bool rotate)268 void GeneralContainer<T,NUM_VEC,LEN_VEC>::delRestart(int n,bool scale,bool translate,bool rotate)
269 {
270 // do only delete property if it is a restart property
271 if(!decidePackUnpackOperation(OPERATION_RESTART, scale, translate, rotate))
272 return;
273
274 numElem_--;
275 if(numElem_ == n) return;
276 for(int i=0;i<NUM_VEC;i++)
277 for(int j=0;j<LEN_VEC;j++)
278 arr_[n][i][j] = arr_[numElem_][i][j];
279 }
280
281 /* ----------------------------------------------------------------------
282 delete all elements if restart
283 ------------------------------------------------------------------------- */
284
285 template<typename T, int NUM_VEC, int LEN_VEC>
delRestart(bool scale,bool translate,bool rotate)286 void GeneralContainer<T,NUM_VEC,LEN_VEC>::delRestart(bool scale,bool translate,bool rotate)
287 {
288 // do only delete property if it is a restart property
289 if(!decidePackUnpackOperation(OPERATION_RESTART, scale, translate, rotate))
290 return;
291
292 numElem_ = 0;
293 }
294
295 /* ----------------------------------------------------------------------
296 get an element
297 ------------------------------------------------------------------------- */
298
299 template<typename T, int NUM_VEC, int LEN_VEC>
get(int n,T ** elem)300 void GeneralContainer<T,NUM_VEC,LEN_VEC>::get(int n, T** elem)
301 {
302 for(int i=0;i<NUM_VEC;i++)
303 for(int j=0;j<LEN_VEC;j++)
304 elem[i][j] = arr_[n][i][j];
305 }
306
307 /* ----------------------------------------------------------------------
308 operator()
309 ------------------------------------------------------------------------- */
310
311 template<typename T, int NUM_VEC, int LEN_VEC>
operator()312 T**& GeneralContainer<T,NUM_VEC,LEN_VEC>::operator() (int n)
313 {
314 return arr_[n];
315 }
316
317 template<typename T, int NUM_VEC, int LEN_VEC>
operator()318 T** const& GeneralContainer<T,NUM_VEC,LEN_VEC>::operator() (int n) const
319 {
320 return arr_[n];
321 }
322
323 /* ----------------------------------------------------------------------
324 set all data by copy from other container
325 ------------------------------------------------------------------------- */
326
327 template<typename T, int NUM_VEC, int LEN_VEC>
setFromContainer(ContainerBase * cont)328 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::setFromContainer(ContainerBase *cont)
329 {
330 GeneralContainer<T,NUM_VEC,LEN_VEC> *gcont = static_cast<GeneralContainer<T,NUM_VEC,LEN_VEC>* >(cont);
331
332 if(size() != gcont->size() || nVec() != gcont->nVec() || lenVec() != gcont->lenVec())
333 return false;
334
335 int len = size();
336 for(int n = 0; n < len; n++)
337 for(int i=0;i<NUM_VEC;i++)
338 for(int j=0;j<LEN_VEC;j++)
339 {
340 arr_[n][i][j] = gcont->arr_[n][i][j];
341
342 }
343
344 return true;
345 }
346
347 /* ----------------------------------------------------------------------
348 average from other container
349 ------------------------------------------------------------------------- */
350
351 template<typename T, int NUM_VEC, int LEN_VEC>
calcAvgFromContainer()352 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::calcAvgFromContainer()
353 {
354
355 GeneralContainer<T,NUM_VEC,LEN_VEC> *gcont = static_cast<GeneralContainer<T,NUM_VEC,LEN_VEC>* >(container_statistics_raw_data_);
356 GeneralContainer<T,1,1> *gscale = dynamic_cast<GeneralContainer<T,1,1>* >(container_statistics_scale_data_);
357 GeneralContainer<T,1,1> *gscaleAvg = dynamic_cast<GeneralContainer<T,1,1>* >(container_statistics_scale_average_data_);
358
359 // source has to be defined
360 if (!gcont)
361 return false;
362
363 // only use if identical dimensions
364 if(size() != gcont->size() || nVec() != gcont->nVec() || lenVec() != gcont->lenVec())
365 return false;
366
367 const int len = size();
368
369 T epsilon = std::numeric_limits<T>::epsilon();
370
371 if (enable_favre_)
372 {
373 for(int n = 0; n < len; n++)
374 {
375 const double scale = (*gscaleAvg)(n)[0][0] < epsilon ? 0.0 : (*gscale)(n)[0][0]/(*gscaleAvg)(n)[0][0];
376 for(int i=0;i<NUM_VEC;i++)
377 for(int j=0;j<LEN_VEC;j++)
378 {
379 const T contribution = gcont->arr_[n][i][j];
380
381 if(MathExtraLiggghts::abs(arr_[n][i][j]) < epsilon)
382 arr_[n][i][j] = contribution;
383 else
384 arr_[n][i][j] = (1.-weighting_factor_*scale)*arr_[n][i][j] +
385 weighting_factor_*scale*contribution;
386 }
387 }
388 }
389 else
390 {
391 for(int n = 0; n < len; n++)
392 for(int i=0;i<NUM_VEC;i++)
393 for(int j=0;j<LEN_VEC;j++)
394 {
395 const T contribution = gcont->arr_[n][i][j];
396
397 if(MathExtraLiggghts::abs(arr_[n][i][j]) < epsilon)
398 arr_[n][i][j] = contribution;
399 else
400 arr_[n][i][j] = (1.-weighting_factor_)*arr_[n][i][j] +
401 weighting_factor_*contribution;
402 }
403 }
404 return true;
405 }
406
407 /* ----------------------------------------------------------------------
408 mean square from other container
409 ------------------------------------------------------------------------- */
410
411 template<typename T, int NUM_VEC, int LEN_VEC>
calcMeanSquareFromContainer()412 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::calcMeanSquareFromContainer()
413 {
414
415 GeneralContainer<T,NUM_VEC,LEN_VEC> *gcont = static_cast<GeneralContainer<T,NUM_VEC,LEN_VEC>* >(container_statistics_raw_data_);
416 GeneralContainer<T,1,1> *gscale = dynamic_cast<GeneralContainer<T,1,1>* >(container_statistics_scale_data_);
417 GeneralContainer<T,1,1> *gscaleAvg = dynamic_cast<GeneralContainer<T,1,1>* >(container_statistics_scale_average_data_);
418
419 // at least source has to be defined
420 if (!gcont)
421 return false;
422
423 // only copy if identical
424 if(size() != gcont->size() || nVec() != gcont->nVec() || lenVec() != gcont->lenVec())
425 return false;
426
427 const int len = size();
428
429 T epsilon = std::numeric_limits<T>::epsilon();
430
431 if (enable_favre_)
432 {
433 for(int n = 0; n < len; n++)
434 {
435 const double scale = (*gscaleAvg)(n)[0][0] < epsilon ? 0.0 : (*gscale)(n)[0][0]/(*gscaleAvg)(n)[0][0];
436 for(int i=0;i<NUM_VEC;i++)
437 for(int j=0;j<LEN_VEC;j++)
438 {
439 const T contribution = gcont->arr_[n][i][j];
440
441 if(MathExtraLiggghts::abs(arr_[n][i][j]) < epsilon)
442 arr_[n][i][j] = contribution*contribution;
443 else
444 arr_[n][i][j] = (1.-weighting_factor_*scale)*arr_[n][i][j] +
445 weighting_factor_*scale*contribution*contribution;
446 }
447
448 }
449 }
450 else
451 {
452 for(int n = 0; n < len; n++)
453 for(int i=0;i<NUM_VEC;i++)
454 for(int j=0;j<LEN_VEC;j++)
455 {
456 const T contribution = gcont->arr_[n][i][j];
457
458 if(MathExtraLiggghts::abs(arr_[n][i][j]) < epsilon)
459 arr_[n][i][j] = contribution*contribution;
460 else
461 arr_[n][i][j] = (1.-weighting_factor_)*arr_[n][i][j] +
462 weighting_factor_*contribution*contribution;
463 }
464
465 }
466
467 return true;
468 }
469
470 // This is the averaging for the scaling arrays
471 template<typename T, int NUM_VEC, int LEN_VEC>
calcSumFromContainer()472 bool GeneralContainer<T,NUM_VEC,LEN_VEC>::calcSumFromContainer()
473 {
474
475 GeneralContainer<T,NUM_VEC,LEN_VEC> *gcont = static_cast<GeneralContainer<T,NUM_VEC,LEN_VEC>* >(container_statistics_raw_data_);
476
477 // at least source has to be defined
478 if (!gcont)
479 return false;
480
481 // only copy if identical
482 if(size() != gcont->size() || nVec() != gcont->nVec() || lenVec() != gcont->lenVec())
483 return false;
484
485 const int len = size();
486 for(int n = 0; n < len; n++)
487 for(int i=0;i<NUM_VEC;i++)
488 for(int j=0;j<LEN_VEC;j++)
489 {
490 arr_[n][i][j] = (1.-weighting_factor_)*arr_[n][i][j] +
491 weighting_factor_*(*gcont)(n)[i][j];
492
493 if (arr_[n][i][j] < std::numeric_limits<T>::epsilon())
494 arr_[n][i][j] = 0;
495 }
496
497 return true;
498 }
499
500 /* ---------------------------------------------------------------------- */
501
502 template<typename T, int NUM_VEC, int LEN_VEC>
setToDefault(int n)503 void GeneralContainer<T,NUM_VEC,LEN_VEC>::setToDefault(int n)
504 {
505
506 for(int i = 0; i < NUM_VEC; i++)
507 for(int j = 0; j < LEN_VEC; j++)
508 arr_[n][i][j] = defaultValue_;
509 }
510
511 template<typename T, int NUM_VEC, int LEN_VEC>
set(int n,T ** elem)512 void GeneralContainer<T,NUM_VEC,LEN_VEC>::set(int n, T** elem)
513 {
514 for(int i = 0; i < NUM_VEC; i++)
515 for(int j = 0; j < LEN_VEC; j++)
516 arr_[n][i][j] = elem[i][j];
517 }
518
519 template<typename T, int NUM_VEC, int LEN_VEC>
set(int n,int m,T * elem)520 void GeneralContainer<T,NUM_VEC,LEN_VEC>::set(int n, int m, T* elem)
521 {
522 for(int j = 0; j < LEN_VEC; j++)
523 arr_[n][m][j] = elem[j];
524 }
525
526 template<typename T, int NUM_VEC, int LEN_VEC>
setAll(T def)527 void GeneralContainer<T,NUM_VEC,LEN_VEC>::setAll(T def)
528 {
529 int len = size();
530 for(int n = 0; n < len; n++)
531 for(int i = 0; i < NUM_VEC; i++)
532 for(int j = 0; j < LEN_VEC; j++)
533 arr_[n][i][j] = def;
534 }
535
536 template<typename T, int NUM_VEC, int LEN_VEC>
setAll(int to,T def)537 void GeneralContainer<T,NUM_VEC,LEN_VEC>::setAll(int to,T def)
538 {
539 int len = std::min(to,size());
540 for(int n = 0; n < len; n++)
541 for(int i = 0; i < NUM_VEC; i++)
542 for(int j = 0; j < LEN_VEC; j++)
543 arr_[n][i][j] = def;
544 }
545
546 template<typename T, int NUM_VEC, int LEN_VEC>
begin()547 T*** GeneralContainer<T,NUM_VEC,LEN_VEC>::begin()
548 {
549 return arr_;
550 }
551
552 template<typename T, int NUM_VEC, int LEN_VEC>
begin_slow_dirty()553 void* GeneralContainer<T,NUM_VEC,LEN_VEC>::begin_slow_dirty()
554 {
555 return (void*) arr_;
556 }
557
558 template<typename T, int NUM_VEC, int LEN_VEC>
getElemSize()559 int GeneralContainer<T,NUM_VEC,LEN_VEC>::getElemSize()
560 {
561 return NUM_VEC*LEN_VEC*sizeof(T);
562 }
563
564 /* ----------------------------------------------------------------------
565 min,max
566 ------------------------------------------------------------------------- */
567
568 template<typename T, int NUM_VEC, int LEN_VEC>
max_scalar()569 T GeneralContainer<T,NUM_VEC,LEN_VEC>::max_scalar()
570 {
571 T max = arr_[0][0][0];
572
573 int len = size();
574 for(int i = 0; i < len; i++)
575 for(int j = 0; j < NUM_VEC; j++)
576 for(int k = 0; k < LEN_VEC; k++)
577 if(arr_[i][j][k] > max)
578 max = arr_[i][j][k];
579
580 return max;
581 }
582
583 template<typename T, int NUM_VEC, int LEN_VEC>
min_scalar()584 T GeneralContainer<T,NUM_VEC,LEN_VEC>::min_scalar()
585 {
586 T min = arr_[0][0][0];
587
588 int len = size();
589 for(int i = 0; i < len; i++)
590 for(int j = 0; j < NUM_VEC; j++)
591 for(int k = 0; k < LEN_VEC; k++)
592 if(arr_[i][j][k] < min)
593 min = arr_[i][j][k];
594
595 return min;
596 }
597
598 /* ----------------------------------------------------------------------
599 translate, rotate, scale
600 ------------------------------------------------------------------------- */
601
602 template<typename T, int NUM_VEC, int LEN_VEC>
scale(double factor)603 void GeneralContainer<T,NUM_VEC,LEN_VEC>::scale(double factor)
604 {
605 if(isScaleInvariant()) return;
606
607 double factorApplied = 1.;
608 for(int i = 0; i < scalePower_; i++)
609 factorApplied *= factor;
610
611 int len = size();
612 for(int i = 0; i < len; i++)
613 for(int j = 0; j < NUM_VEC;j++)
614 for(int k = 0; k < LEN_VEC; k++)
615 arr_[i][j][k] *= factorApplied;
616 }
617
618 template<typename T, int NUM_VEC, int LEN_VEC>
move(const double * const delta)619 void GeneralContainer<T,NUM_VEC,LEN_VEC>::move(const double * const delta)
620 {
621 if(isTranslationInvariant()) return;
622
623 int len = size();
624
625 for(int i = 0; i < len; i++)
626 for(int j = 0; j < NUM_VEC; j++)
627 for(int k = 0; k < LEN_VEC; k++)
628 arr_[i][j][k] += delta[k];
629 }
630
631 template<typename T, int NUM_VEC, int LEN_VEC>
moveElement(const int i,const double * const delta)632 void GeneralContainer<T,NUM_VEC,LEN_VEC>::moveElement(const int i, const double * const delta)
633 {
634 if(isTranslationInvariant()) return;
635
636 for(int j = 0; j < NUM_VEC; j++)
637 for(int k = 0; k < LEN_VEC; k++)
638 arr_[i][j][k] += delta[k];
639 }
640
641 template<typename T, int NUM_VEC, int LEN_VEC>
rotate(const double * const dQ)642 void GeneralContainer<T,NUM_VEC,LEN_VEC>::rotate(const double * const dQ)
643 {
644 if(isRotationInvariant()) return;
645
646 // ATTENTION: only correct for 3D vectors
647 int len = size();
648 for(int i = 0; i < len; i++)
649 for(int j = 0; j < NUM_VEC; j++)
650 MathExtraLiggghts::vec_quat_rotate(arr_[i][j],dQ);
651 }
652
653 /* ----------------------------------------------------------------------
654 buffer size for all elements, push / pop for all elements
655 used for global properties
656 ------------------------------------------------------------------------- */
657
658 template<typename T, int NUM_VEC, int LEN_VEC>
bufSize(int operation,bool scale,bool translate,bool rotate)659 int GeneralContainer<T,NUM_VEC,LEN_VEC>::bufSize(int operation,bool scale,bool translate,bool rotate) const
660 {
661 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
662 return 0;
663
664 if(!this->decideCommOperation(operation))
665 return 0;
666
667 return (1 + size()*NUM_VEC*LEN_VEC);
668 }
669
670 template<typename T, int NUM_VEC, int LEN_VEC>
pushToBuffer(double * buf,int operation,bool scale,bool translate,bool rotate)671 int GeneralContainer<T,NUM_VEC,LEN_VEC>::pushToBuffer(double *buf,int operation,bool scale,bool translate, bool rotate)
672 {
673 //TODO throw error if sizeof(T) > sizeof(double)
674
675 int m = 0;
676
677 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
678 return 0;
679
680 int len = size();
681
682 buf[m++] = static_cast<double>(len);
683
684 for(int i = 0; i < len; i++)
685 for(int j = 0; j < NUM_VEC; j++)
686 for(int k = 0; k < LEN_VEC; k++)
687 buf[m++] = static_cast<double>(arr_[i][j][k]);
688
689 return (1 + len*NUM_VEC*LEN_VEC);
690 }
691
692 template<typename T, int NUM_VEC, int LEN_VEC>
popFromBuffer(double * buf,int operation,bool scale,bool translate,bool rotate)693 int GeneralContainer<T,NUM_VEC,LEN_VEC>::popFromBuffer(double *buf,int operation,bool scale,bool translate, bool rotate)
694 {
695 int nNew, m = 0;
696
697 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
698 return 0;
699
700 if(decideCreateNewElements(operation))
701 {
702 T** tmp;
703 create<T>(tmp,NUM_VEC,LEN_VEC);
704
705 nNew = static_cast<int>(buf[m++]);
706
707 for(int i = 0; i < nNew; i++)
708 {
709 for(int j = 0; j < NUM_VEC; j++)
710 for(int k = 0; k < LEN_VEC; k++)
711 tmp[j][k] = static_cast<T>(buf[m++]);
712 add(tmp);
713 }
714
715 destroy<T>(tmp);
716
717 return (1 + nNew*NUM_VEC*LEN_VEC);
718 }
719 else return 0;
720 }
721
722 /* ----------------------------------------------------------------------
723 buffer size for a list of elements, push / pop a list of elements
724 used for borders, fw and rev comm for element properties
725 ------------------------------------------------------------------------- */
726
727 template<typename T, int NUM_VEC, int LEN_VEC>
elemListBufSize(int n,int operation,bool scale,bool translate,bool rotate)728 int GeneralContainer<T,NUM_VEC,LEN_VEC>::elemListBufSize(int n,int operation,bool scale,bool translate,bool rotate)
729 {
730 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
731 return 0;
732
733 if(!this->decideCommOperation(operation))
734 return 0;
735
736 return (n*NUM_VEC*LEN_VEC);
737 }
738
739 template<typename T, int NUM_VEC, int LEN_VEC>
pushElemListToBuffer(int n,int * list,int * wraplist,double * buf,int operation,double * dlo,double * dhi,bool scale,bool translate,bool rotate)740 int GeneralContainer<T,NUM_VEC,LEN_VEC>::pushElemListToBuffer(int n, int *list, int *wraplist, double *buf,int operation, double *dlo, double *dhi, bool scale,bool translate, bool rotate)
741 {
742 int i,m = 0;
743
744 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
745 return 0;
746
747 if(!this->decideCommOperation(operation))
748 return 0;
749
750 for(int ii = 0; ii < n; ii++)
751 {
752 i = list[ii];
753 for(int j = 0; j < NUM_VEC; j++)
754 for(int k = 0; k < LEN_VEC; k++)
755 {
756 buf[m] = static_cast<double>(arr_[i][j][k]);
757 if (wrapPeriodic())
758 {
759 const int wrap = wraplist[ii];
760 if (wrap != IS_GHOST)
761 {
762 if ((k == 0 && wrap == IS_GHOST_WRAP_DIM_0_NEG) ||
763 (k == 1 && wrap == IS_GHOST_WRAP_DIM_1_NEG) ||
764 (k == 2 && wrap == IS_GHOST_WRAP_DIM_2_NEG) )
765 buf[m] -= dhi[k] - dlo[k];
766 else if ((k == 0 && wrap == IS_GHOST_WRAP_DIM_0_POS) ||
767 (k == 1 && wrap == IS_GHOST_WRAP_DIM_1_POS) ||
768 (k == 2 && wrap == IS_GHOST_WRAP_DIM_2_POS) )
769 buf[m] += dhi[k] - dlo[k];
770 }
771 }
772 m++;
773 }
774 }
775
776 return (n*NUM_VEC*LEN_VEC);
777 }
778
779 template<typename T, int NUM_VEC, int LEN_VEC>
popElemListFromBuffer(int first,int n,double * buf,int operation,bool scale,bool translate,bool rotate)780 int GeneralContainer<T,NUM_VEC,LEN_VEC>::popElemListFromBuffer(int first, int n, double *buf,int operation,bool scale,bool translate, bool rotate)
781 {
782 int m = 0;
783
784 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
785 return 0;
786
787 bool pullBuf = decideCommOperation(operation);
788
789 bool createElem = decideCreateNewElements(operation);
790
791 T** tmp;
792 create<T>(tmp,NUM_VEC,LEN_VEC);
793
794 for(int i = first; i < first+n; i++)
795 {
796 for(int j = 0; j < NUM_VEC; j++)
797 for(int k = 0; k < LEN_VEC; k++)
798 (createElem ? tmp[j][k] : arr_[i][j][k]) = (pullBuf ? static_cast<T>(buf[m++]) : static_cast<T>(0));
799
800 if(createElem) add(tmp);
801 }
802
803 destroy<T>(tmp);
804
805 return m;
806 }
807
808 template<typename T, int NUM_VEC, int LEN_VEC>
pushElemListToBufferReverse(int first,int n,double * buf,int operation,bool scale,bool translate,bool rotate)809 int GeneralContainer<T,NUM_VEC,LEN_VEC>::pushElemListToBufferReverse(int first, int n, double *buf,int operation,bool scale,bool translate, bool rotate)
810 {
811 int m = 0;
812
813 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
814 return 0;
815
816 for(int i = first; i < first+n; i++)
817 {
818 for(int j = 0; j < NUM_VEC; j++)
819 for(int k = 0; k < LEN_VEC; k++)
820 buf[m++] = static_cast<double>(arr_[i][j][k]);
821 }
822
823 return (n*NUM_VEC*LEN_VEC);
824 }
825
826 template<typename T, int NUM_VEC, int LEN_VEC>
popElemListFromBufferReverse(int n,int * list,double * buf,int operation,bool scale,bool translate,bool rotate)827 int GeneralContainer<T,NUM_VEC,LEN_VEC>::popElemListFromBufferReverse(int n, int *list,double *buf,int operation,bool scale,bool translate, bool rotate)
828 {
829 int i,m = 0;
830
831 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
832 return 0;
833
834 if(COMM_TYPE_REVERSE == this->communicationType())
835 {
836
837 for(int ii = 0; ii < n; ii++)
838 {
839 i = list[ii];
840 for(int j = 0; j < NUM_VEC; j++)
841 for(int k = 0; k < LEN_VEC; k++)
842 arr_[i][j][k] += static_cast<T>(buf[m++]);
843 }
844 }
845 else if(sizeof(int) == sizeof(T) && COMM_TYPE_REVERSE_BITFIELD == this->communicationType())
846 {
847
848 for(int ii = 0; ii < n; ii++)
849 {
850 i = list[ii];
851 for(int j = 0; j < NUM_VEC; j++)
852 for(int k = 0; k < LEN_VEC; k++)
853 arr_[i][j][k] = (T) (static_cast<int>(arr_[i][j][k]) | static_cast<int>(buf[m++]));
854 }
855 }
856
857 return (n*NUM_VEC*LEN_VEC);
858 }
859
860 /* ----------------------------------------------------------------------
861 buffer size for a single element, push / pop a single element
862 used for exchange of single elements
863 ------------------------------------------------------------------------- */
864
865 template<typename T, int NUM_VEC, int LEN_VEC>
elemBufSize(int operation,bool scale,bool translate,bool rotate)866 int GeneralContainer<T,NUM_VEC,LEN_VEC>::elemBufSize(int operation,bool scale,bool translate,bool rotate)
867 {
868
869 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
870 return 0;
871
872 if(!this->decideCommOperation(operation))
873 return 0;
874
875 return (NUM_VEC*LEN_VEC);
876 }
877
878 template<typename T, int NUM_VEC, int LEN_VEC>
pushElemToBuffer(int i,double * buf,int operation,bool scale,bool translate,bool rotate)879 int GeneralContainer<T,NUM_VEC,LEN_VEC>::pushElemToBuffer(int i, double *buf,int operation,bool scale,bool translate, bool rotate)
880 {
881 int m = 0;
882
883 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
884 return 0;
885
886 if(!this->decideCommOperation(operation))
887 return 0;
888
889 for(int j = 0; j < NUM_VEC; j++)
890 for(int k = 0; k < LEN_VEC; k++)
891 buf[m++] = static_cast<double>(arr_[i][j][k]);
892
893 return m;
894 }
895
896 template<typename T, int NUM_VEC, int LEN_VEC>
popElemFromBuffer(double * buf,int operation,bool scale,bool translate,bool rotate)897 int GeneralContainer<T,NUM_VEC,LEN_VEC>::popElemFromBuffer(double *buf,int operation,bool scale,bool translate, bool rotate)
898 {
899 int m = 0;
900
901 if(!this->decidePackUnpackOperation(operation,scale,translate,rotate))
902 return 0;
903
904 bool pullBuf = decideCommOperation(operation);
905
906 T** tmp;
907 create<T>(tmp,NUM_VEC,LEN_VEC);
908
909 for(int j = 0; j < NUM_VEC; j++)
910 for(int k = 0; k < LEN_VEC; k++)
911 tmp[j][k] = pullBuf ? static_cast<T>(buf[m++]) : static_cast<T>(0);
912
913 add(tmp);
914 destroy<T>(tmp);
915
916 return m;
917 }
918
919 #endif
920