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 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 #ifndef LMP_VECTOR_LIGGGHTS_H
43 #define LMP_VECTOR_LIGGGHTS_H
44
45 #include <cmath>
46 #include "lammps.h"
47
48 namespace LAMMPS_NS {
49
50 //================================================
51 //SOME VERY SIMPLE VECTOR OPERATIONS
52 //================================================
53
54 template<typename T>
vectorConstruct3D(T * v,T x,T y,T z)55 inline void vectorConstruct3D(T *v,T x, T y, T z)
56 {
57 v[0] = x;
58 v[1] = y;
59 v[2] = z;
60 }
61
vectorNormalize3D(double * v)62 inline void vectorNormalize3D(double *v)
63 {
64 const double norm = ::sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
65 const double invnorm = (norm == 0.) ? 0. : 1./norm;
66 v[0] *= invnorm;
67 v[1] *= invnorm;
68 v[2] *= invnorm;
69 }
70
vectorMag3D(const double * v)71 inline double vectorMag3D(const double *v)
72 {
73 return ( ::sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]) );
74 }
75
vectorMag3DSquared(const double * v)76 inline double vectorMag3DSquared(const double *v)
77 {
78 return ( v[0]*v[0]+v[1]*v[1]+v[2]*v[2] );
79 }
80
vectorMag4D(const double * v)81 inline double vectorMag4D(const double *v)
82 {
83 return ( ::sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]+v[3]*v[3]) );
84 }
85
pointDistanceSquared(const double * point1,const double * point2)86 inline double pointDistanceSquared(const double *point1, const double *point2)
87 {
88 return
89 (point1[0]-point2[0]) * (point1[0]-point2[0]) +
90 (point1[1]-point2[1]) * (point1[1]-point2[1]) +
91 (point1[2]-point2[2]) * (point1[2]-point2[2]);
92 }
93
pointDistance(const double * point1,const double * point2)94 inline double pointDistance(const double *point1, const double *point2)
95 {
96 return ::sqrt(pointDistanceSquared(point1, point2));
97 }
98
vectorMag4DSquared(const double * v)99 inline double vectorMag4DSquared(const double *v)
100 {
101 return ( v[0]*v[0]+v[1]*v[1]+v[2]*v[2]+v[3]*v[3] );
102 }
103
vectorDot3D(const double * v1,const double * v2)104 inline double vectorDot3D(const double *v1, const double *v2)
105 {
106 return (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]);
107 }
108
vectorDot2D(const double * v1,const double * v2)109 inline double vectorDot2D(const double *v1, const double *v2)
110 {
111 return (v1[0]*v2[0]+v1[1]*v2[1]);
112 }
113
114 template<typename T>
vectorCopy2D(const T * from,T * to)115 inline void vectorCopy2D(const T *from, T *to)
116 {
117 to[0]=from[0];
118 to[1]=from[1];
119 }
120
vectorFlip3D(double * v)121 inline void vectorFlip3D(double *v)
122 {
123 v[0]=-v[0];
124 v[1]=-v[1];
125 v[2]=-v[2];
126 }
127
128 template<typename T>
vectorCopyN(const T * from,T * to,int N)129 inline void vectorCopyN(const T *from, T *to, int N)
130 {
131 for(int i = 0; i < N; ++i)
132 to[i] = from[i];
133 }
134
135 template<typename T>
vectorCopy3D(const T * from,T * to)136 inline void vectorCopy3D(const T *from, T *to)
137 {
138 to[0]=from[0];
139 to[1]=from[1];
140 to[2]=from[2];
141 }
142
vectorRoundN(double * vec,int N)143 inline void vectorRoundN(double *vec, int N)
144 {
145 for(int i = 0; i < N; ++i)
146 vec[i] = static_cast<double>(round(vec[i]));
147 }
148
vectorAbs3D(double * v)149 inline void vectorAbs3D(double *v)
150 {
151 if(v[0] < 0) v[0] = -v[0];
152 if(v[1] < 0) v[1] = -v[1];
153 if(v[2] < 0) v[2] = -v[2];
154 }
155
156 template<typename T>
vectorMin3D(const T * v,int & dim)157 inline T vectorMin3D(const T *v,int &dim)
158 {
159 if(v[0] < v[1] && v[0] < v[2])
160 {
161 dim = 0;
162 return v[0];
163 }
164
165 if(v[1] < v[2])
166 {
167 dim = 1;
168 return v[1];
169 }
170
171 dim = 2;
172 return v[2];
173 }
174
175 template<typename T>
vectorMin3D(const T * v)176 inline T vectorMin3D(const T *v)
177 {
178 if(v[0] < v[1] && v[0] < v[2])
179 return v[0];
180
181 if(v[1] < v[2])
182 return v[1];
183
184 return v[2];
185 }
186
187 template<typename T>
vectorMax3D(const T * v)188 inline T vectorMax3D(const T *v)
189 {
190 if(v[0] > v[1] && v[0] > v[2])
191 return v[0];
192
193 if(v[1] > v[2])
194 return v[1];
195
196 return v[2];
197 }
198
199 template<typename T>
vectorMaxN(const T * v,const int n)200 inline T vectorMaxN(const T *v, const int n)
201 {
202 T max = v[0];
203 for (int i=1;i<n;++i)
204 max = max > v[i] ? max : v[i];
205 return max;
206 }
207
208 template<typename T>
vectorMinN(const T * v,const int n)209 inline T vectorMinN(const T *v, const int n)
210 {
211 T min = v[0];
212 for (int i=1;i<n;++i)
213 min = min < v[i] ? min : v[i];
214 return min;
215 }
216
vectorComponentMin3D(double const * v1,double const * v2,double * min)217 inline void vectorComponentMin3D(double const*v1,double const*v2,double *min)
218 {
219 if(v1[0] > v2[0])
220 min[0] = v2[0];
221 else
222 min[0] = v1[0];
223
224 if(v1[1] > v2[1])
225 min[1] = v2[1];
226 else
227 min[1] = v1[1];
228
229 if(v1[2] > v2[2])
230 min[2] = v2[2];
231 else
232 min[2] = v1[2];
233 }
234
vectorComponentMax3D(double const * v1,double const * v2,double * max)235 inline void vectorComponentMax3D(double const*v1,double const*v2,double *max)
236 {
237 if(v1[0] > v2[0])
238 max[0] = v1[0];
239 else
240 max[0] = v2[0];
241
242 if(v1[1] > v2[1])
243 max[1] = v1[1];
244 else
245 max[1] = v2[1];
246
247 if(v1[2] > v2[2])
248 max[2] = v1[2];
249 else
250 max[2] = v2[2];
251 }
252
vectorCopy4D(const double * from,double * to)253 inline void vectorCopy4D(const double *from, double *to)
254 {
255 to[0]=from[0];
256 to[1]=from[1];
257 to[2]=from[2];
258 to[3]=from[3];
259 }
260
vectorScalarMultN(int n,double * v,double s)261 inline void vectorScalarMultN(int n,double *v, double s)
262 {
263 for(int i = 0; i < n; ++i)
264 v[i] = s*v[i];
265 }
266
vectorScalarMultN(int n,int * v,double s)267 inline void vectorScalarMultN(int n,int *v, double s)
268 {
269 for(int i = 0; i < n; ++i)
270 v[i] = static_cast<int>(static_cast<double>(s)*v[i]);
271 }
272
vectorScalarMult3D(double * v,double s)273 inline void vectorScalarMult3D(double *v, double s)
274 {
275 v[0]=s*v[0];
276 v[1]=s*v[1];
277 v[2]=s*v[2];
278 }
279
vectorScalarMult3D(const double * v,double s,double * result)280 inline void vectorScalarMult3D(const double *v, double s, double *result)
281 {
282 result[0]=s*v[0];
283 result[1]=s*v[1];
284 result[2]=s*v[2];
285 }
286
vectorScalarDiv3D(double * v,double s)287 inline void vectorScalarDiv3D(double *v, double s)
288 {
289 const double sinv = 1./s;
290 v[0]=sinv*v[0];
291 v[1]=sinv*v[1];
292 v[2]=sinv*v[2];
293 }
294
vectorComponentDiv3D(const double * nom,const double * denom,double * result)295 inline void vectorComponentDiv3D(const double *nom, const double *denom,double *result)
296 {
297 result[0]=nom[0]/denom[0];
298 result[1]=nom[1]/denom[1];
299 result[2]=nom[2]/denom[2];
300 }
301
vectorScalarAdd3D(double * v,double s)302 inline void vectorScalarAdd3D(double *v, double s)
303 {
304 v[0]+=s;
305 v[1]+=s;
306 v[2]+=s;
307 }
308
vectorScalarSubtract3D(double * v,double s)309 inline void vectorScalarSubtract3D(double *v, double s)
310 {
311 v[0]-=s;
312 v[1]-=s;
313 v[2]-=s;
314 }
315
vectorNegate3D(const double * const v,double * const result)316 inline void vectorNegate3D(const double * const v, double * const result)
317 {
318 result[0]=-v[0];
319 result[1]=-v[1];
320 result[2]=-v[2];
321 }
322
vectorNegate3D(double * const v)323 inline void vectorNegate3D(double * const v)
324 {
325 v[0]=-v[0];
326 v[1]=-v[1];
327 v[2]=-v[2];
328 }
329
vectorScalarDiv3D(double * v,double s,double * result)330 inline void vectorScalarDiv3D(double *v, double s, double *result)
331 {
332 const double sinv = 1./s;
333 result[0]=sinv*v[0];
334 result[1]=sinv*v[1];
335 result[2]=sinv*v[2];
336 }
337
vectorAdd3D(const double * v1,const double * v2,double * result)338 inline void vectorAdd3D(const double *v1, const double *v2, double *result)
339 {
340 result[0]=v1[0]+v2[0];
341 result[1]=v1[1]+v2[1];
342 result[2]=v1[2]+v2[2];
343 }
344
vectorAdd4D(double * v1,const double * v2)345 inline void vectorAdd4D(double *v1, const double *v2)
346 {
347 v1[0]+=v2[0];
348 v1[1]+=v2[1];
349 v1[2]+=v2[2];
350 v1[3]+=v2[3];
351 }
352
353 template<typename T>
vectorAddN(T * v1,const T * v2,int n)354 inline void vectorAddN(T *v1, const T *v2, int n)
355 {
356 for(int i = 0; i < n; ++i)
357 v1[i] += v2[i];
358 }
359
vectorAddMultiple3D(const double * v1,double v2factor,const double * v2,double * result)360 inline void vectorAddMultiple3D(const double *v1, double v2factor, const double *v2, double *result)
361 {
362 result[0]=v1[0]+v2factor*v2[0];
363 result[1]=v1[1]+v2factor*v2[1];
364 result[2]=v1[2]+v2factor*v2[2];
365 }
366
vectorSubtract4D(const double * v1,const double * v2,double * result)367 inline void vectorSubtract4D(const double *v1,const double *v2, double *result)
368 {
369 result[0]=v1[0]-v2[0];
370 result[1]=v1[1]-v2[1];
371 result[2]=v1[2]-v2[2];
372 result[3]=v1[3]-v2[3];
373 }
374
vectorSubtract3D(const double * v1,const double * v2,double * result)375 inline void vectorSubtract3D(const double *v1,const double *v2, double *result)
376 {
377 result[0]=v1[0]-v2[0];
378 result[1]=v1[1]-v2[1];
379 result[2]=v1[2]-v2[2];
380 }
381
vectorSubtract2D(const double * v1,const double * v2,double * result)382 inline void vectorSubtract2D(const double *v1,const double *v2, double *result)
383 {
384 result[0]=v1[0]-v2[0];
385 result[1]=v1[1]-v2[1];
386 }
387
388 template<typename T>
vectorMultiN(const T * v1,const T * v2,T * result,const int n)389 inline void vectorMultiN(const T *v1, const T *v2, T *result, const int n)
390 {
391 for(int i = 0; i < n; ++i)
392 result[i] = v1[i] * v2[i];
393 }
394
395 template<typename T>
vectorComponentDivN(const T * v1,const T * v2,T * result,const int n)396 inline void vectorComponentDivN(const T *v1, const T *v2, T *result, const int n)
397 {
398 for(int i = 0; i < n; ++i)
399 result[i] = v1[i] / v2[i];
400 }
401
vectorCross3D(const double * v1,const double * v2,double * result)402 inline void vectorCross3D(const double *v1,const double *v2, double *result)
403 {
404 result[0]=v1[1]*v2[2]-v1[2]*v2[1];
405 result[1]=v1[2]*v2[0]-v1[0]*v2[2];
406 result[2]=v1[0]*v2[1]-v1[1]*v2[0];
407 }
408
vectorCrossMag3D(const double * v1,const double * v2)409 inline double vectorCrossMag3D(const double *v1,const double *v2)
410 {
411 double res[3];
412 res[0]=v1[1]*v2[2]-v1[2]*v2[1];
413 res[1]=v1[2]*v2[0]-v1[0]*v2[2];
414 res[2]=v1[0]*v2[1]-v1[1]*v2[0];
415 return vectorMag3D(res);
416 }
417
triangleArea(const double * const v1,const double * const v2,const double * const v3,const double * const n)418 inline double triangleArea(const double * const v1, const double * const v2, const double * const v3, const double * const n)
419 {
420 // formula: |((v1-v3) x (v2 - v3)).n|
421 return fabs(0.5*(
422 ((v1[1]-v3[1])*(v2[2]-v3[2])-(v1[2]-v3[2])*(v2[1]-v3[1]))*n[0] +
423 ((v1[2]-v3[2])*(v2[0]-v3[0])-(v1[0]-v3[0])*(v2[2]-v3[2]))*n[1] +
424 ((v1[0]-v3[0])*(v2[1]-v3[1])-(v1[1]-v3[1])*(v2[0]-v3[0]))*n[2]
425 ));
426 }
427
428 template<typename T>
vectorProject3D(const T * v,const T * on,T * result)429 inline void vectorProject3D(const T *v, const T *on, T *result)
430 {
431 T norm[3] = {on[0], on[1], on[2]};
432 vectorNormalize3D(norm);
433 vectorScalarMult3D(norm,vectorDot3D(v,norm),result);
434 }
435
436 template<typename T>
vectorZeroize3D(T * v)437 inline void vectorZeroize3D(T *v)
438 {
439 v[0]=0;
440 v[1]=0;
441 v[2]=0;
442 }
443
444 template<typename T>
vectorZeroize4D(T * v)445 inline void vectorZeroize4D(T *v)
446 {
447 v[0]=0;
448 v[1]=0;
449 v[2]=0;
450 v[3]=0;
451 }
452
453 template<typename T>
vectorZeroizeN(T * v,const int n)454 inline void vectorZeroizeN(T *v,const int n)
455 {
456 for(int i = 0; i < n; ++i)
457 v[i]=0;
458 }
459
460 template<typename T>
vectorInitialize3D(T * v,T init)461 inline void vectorInitialize3D(T *v,T init)
462 {
463 v[0]=init;
464 v[1]=init;
465 v[2]=init;
466 }
467
468 template<typename T>
vectorInitializeN(T * v,const int n,const T init)469 inline void vectorInitializeN(T *v,const int n,const T init)
470 {
471 for(int i = 0; i < n; ++i)
472 v[i]=init;
473 }
474
475 template<typename T>
vectorSumN(const T * const v,const unsigned int n)476 inline T vectorSumN(const T * const v, const unsigned int n)
477 {
478 T sum = 0;
479 for (unsigned int i = 0; i < n; ++i)
480 sum += v[i];
481 return sum;
482 }
483
quatIdentity4D(double * q)484 inline void quatIdentity4D(double *q)
485 {
486 q[0]=1.;
487 q[1]=0.;
488 q[2]=0.;
489 q[3]=0.;
490 }
491
quatNormalize4D(double * q)492 inline void quatNormalize4D(double *q)
493 {
494 const double norm = vectorMag4D(q);
495 const double invnorm = (norm == 0.) ? 0. : 1./norm;
496 q[0] *= invnorm;
497 q[1] *= invnorm;
498 q[2] *= invnorm;
499 q[3] *= invnorm;
500 }
501
isIdentityQuat4D(const double * q)502 inline bool isIdentityQuat4D(const double *q)
503 {
504 return
505 (
506 q[0] == 1. &&
507 q[1] == 0. &&
508 q[2] == 0. &&
509 q[3] == 0.
510 );
511 }
512
quatMult4D(const double * const q,const double * const p,double * const result)513 inline void quatMult4D(const double * const q, const double * const p, double * const result)
514 {
515 result[0] = q[0]*p[0] - q[1]*p[1] - q[2]*p[2] - q[3]*p[3];
516 result[1] = q[0]*p[1] + q[1]*p[0] + q[2]*p[3] - q[3]*p[2];
517 result[2] = q[0]*p[2] - q[1]*p[3] + q[2]*p[0] + q[3]*p[1];
518 result[3] = q[0]*p[3] + q[1]*p[2] - q[2]*p[1] + q[3]*p[0];
519 }
520
quatMult4D(double * const q,const double * const p)521 inline void quatMult4D(double * const q, const double * const p)
522 {
523 double tmp[4];
524 quatMult4D(q, p, tmp);
525 vectorCopy4D(tmp, q);
526 }
527
quatInverse4D(double * q,double * result)528 inline void quatInverse4D(double *q, double *result)
529 {
530 double invNorm = 1.0/vectorMag4DSquared(q);
531 result[0] = q[0]*invNorm;
532 result[1] = -q[1]*invNorm;
533 result[2] = -q[2]*invNorm;
534 result[3] = -q[3]*invNorm;
535 }
536
phiToQuat(const double phi,const double * const axis,double * const q)537 inline void phiToQuat(const double phi, const double * const axis, double * const q)
538 {
539 q[0] = ::cos(phi*0.5);
540 const double sinPhi = ::sin(phi*0.5);
541 q[1] = axis[0]*sinPhi;
542 q[2] = axis[1]*sinPhi;
543 q[3] = axis[2]*sinPhi;
544 }
545
normalize_bary(double * v)546 inline void normalize_bary(double *v)
547 {
548 const double mag = v[0]+v[1]+v[2];
549 v[0]/=mag;
550 v[1]/=mag;
551 v[2]/=mag;
552 }
553
vectorToBuf3D(double const * vec,double * buf,int & m)554 inline void vectorToBuf3D(double const*vec,double *buf,int &m)
555 {
556 buf[m++] = vec[0];
557 buf[m++] = vec[1];
558 buf[m++] = vec[2];
559 }
560
bufToVector3D(double * vec,double const * buf,int & m)561 inline void bufToVector3D(double *vec,double const*buf,int &m)
562 {
563 vec[0] = buf[m++];
564 vec[1] = buf[m++];
565 vec[2] = buf[m++];
566 }
567
vectorToBuf4D(double const * vec,double * buf,int & m)568 inline void vectorToBuf4D(double const*vec,double *buf,int &m)
569 {
570 buf[m++] = vec[0];
571 buf[m++] = vec[1];
572 buf[m++] = vec[2];
573 buf[m++] = vec[3];
574 }
575
bufToVector4D(double * vec,double const * buf,int & m)576 inline void bufToVector4D(double *vec,double const*buf,int &m)
577 {
578 vec[0] = buf[m++];
579 vec[1] = buf[m++];
580 vec[2] = buf[m++];
581 vec[3] = buf[m++];
582 }
583
vectorToBufN(double const * vec,double * buf,int & m,int nvalues)584 inline void vectorToBufN(double const*vec,double *buf,int &m, int nvalues)
585 {
586 for(int i = 0; i < nvalues; ++i)
587 buf[m++] = vec[i];
588 }
589
bufToVectorN(double * vec,double const * buf,int & m,int nvalues)590 inline void bufToVectorN(double *vec,double const*buf,int &m, int nvalues)
591 {
592 for(int i = 0; i < nvalues; ++i)
593 vec[i] = buf[m++];
594 }
595
printVec2D(FILE * out,const char * name,double const * vec)596 inline void printVec2D(FILE *out, const char *name, double const *vec)
597 {
598 fprintf(out," vector %s: %e %e\n",name,vec[0],vec[1]);
599 }
600
printVec3D(FILE * out,const char * name,double const * vec)601 inline void printVec3D(FILE *out, const char *name, double const*vec)
602 {
603 fprintf(out," vector %s: %e %e %e\n",name,vec[0],vec[1],vec[2]);
604 }
605
printVec3D(FILE * out,const char * name,int const * vec)606 inline void printVec3D(FILE *out, const char *name, int const*vec)
607 {
608 fprintf(out," vector %s: %d %d %d\n",name,vec[0],vec[1],vec[2]);
609 }
610
printVec4D(FILE * out,const char * name,double const * vec)611 inline void printVec4D(FILE *out, const char *name, double const*vec)
612 {
613 fprintf(out," vector %s: %e %e %e %e\n",name,vec[0],vec[1],vec[2],vec[3]);
614 }
615
printVecN(FILE * out,const char * name,double const * vec,int n)616 inline void printVecN(FILE *out, const char *name, double const *vec, int n)
617 {
618 if(name) fprintf(out," vector %s:",name);
619 for(int i = 0; i < n; ++i)
620 fprintf(out,"%f ",vec[i]);
621 fprintf(out,"\n");
622 }
623
printVecN(FILE * out,const char * name,int const * vec,int n)624 inline void printVecN(FILE *out, const char *name, int const *vec, int n)
625 {
626 if(name) fprintf(out," vector %s:",name);
627 for(int i = 0; i < n; ++i)
628 fprintf(out,"%d ",vec[i]);
629 fprintf(out,"\n");
630 }
631
printMat33(FILE * out,const char * name,double const * const * mat)632 inline void printMat33(FILE *out, const char *name, double const* const* mat)
633 {
634 fprintf(out," matrix %s: %f %f %f\n",name,mat[0][0],mat[0][1],mat[0][2]);
635 fprintf(out," %s: %f %f %f\n",name,mat[1][0],mat[1][1],mat[1][2]);
636 fprintf(out," %s: %f %f %f\n",name,mat[2][0],mat[2][1],mat[2][2]);
637 }
638
639 }
640
641 #endif
642