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 Christoph Kloss (DCS Computing GmbH, Linz)
39 Christoph Kloss (JKU Linz)
40 Philippe Seil (JKU Linz)
41
42 Copyright 2012- DCS Computing GmbH, Linz
43 Copyright 2009-2012 JKU Linz
44 ------------------------------------------------------------------------- */
45
46 #ifndef LMP_MULTI_NODE_MESH_I_H
47 #define LMP_MULTI_NODE_MESH_I_H
48
49 /* ----------------------------------------------------------------------
50 consturctors
51 ------------------------------------------------------------------------- */
52
53 template<int NUM_NODES>
MultiNodeMesh(LAMMPS * lmp)54 MultiNodeMesh<NUM_NODES>::MultiNodeMesh(LAMMPS *lmp)
55 : AbstractMesh(lmp),
56 node_("node"),
57 node_orig_(0),
58 nodesLastRe_("nodesLastRe"),
59 center_("center"),
60 rBound_("rBound"),
61 random_(new RanPark(lmp,"179424799")), // big prime #
62 mesh_id_(0),
63 precision_(EPSILON_PRECISION),
64 min_feature_length_(-1.),
65 element_exclusion_list_(0),
66 autoRemoveDuplicates_(false),
67 nMove_(0),
68 nScale_(0),
69 nTranslate_(0),
70 nRotate_(0),
71 store_vel(0),
72 store_omega(0),
73 step_store_vel(0),
74 step_store_omega(0),
75 stepLastReset_(-1)
76 {
77 vectorZeroize3D(global_vel);
78 quatIdentity4D(global_quaternion);
79 quatIdentity4D(prev_quaternion);
80 center_.setWrapPeriodic(true);
81 node_.setWrapPeriodic(true);
82 }
83
84 /* ----------------------------------------------------------------------
85 destructor
86 ------------------------------------------------------------------------- */
87
88 template<int NUM_NODES>
~MultiNodeMesh()89 MultiNodeMesh<NUM_NODES>::~MultiNodeMesh()
90 {
91 if(node_orig_) delete node_orig_;
92 delete random_;
93 if(mesh_id_) delete []mesh_id_;
94 }
95
96 /* ----------------------------------------------------------------------
97 set ID and mesh accuracy (latter used for mesh topology)
98 set healing options
99 ------------------------------------------------------------------------- */
100
101 template<int NUM_NODES>
setMeshID(const char * _id)102 void MultiNodeMesh<NUM_NODES>::setMeshID(const char *_id)
103 {
104 if(mesh_id_) delete []mesh_id_;
105 mesh_id_ = new char[strlen(_id)+1];
106 strcpy(mesh_id_,_id);
107 }
108
109 template<int NUM_NODES>
setPrecision(double _precision)110 void MultiNodeMesh<NUM_NODES>::setPrecision(double _precision)
111 {
112 precision_ = _precision;
113 }
114
115 template<int NUM_NODES>
setMinFeatureLength(double _min_feature_length)116 void MultiNodeMesh<NUM_NODES>::setMinFeatureLength(double _min_feature_length)
117 {
118 min_feature_length_ = _min_feature_length;
119 }
120
121 template<int NUM_NODES>
setElementExclusionList(FILE * _file)122 void MultiNodeMesh<NUM_NODES>::setElementExclusionList(FILE *_file)
123 {
124 element_exclusion_list_ = _file;
125 }
126
127 template<int NUM_NODES>
autoRemoveDuplicates()128 void MultiNodeMesh<NUM_NODES>::autoRemoveDuplicates()
129 {
130 autoRemoveDuplicates_ = true;
131 }
132
133 /* ----------------------------------------------------------------------
134 add an element - only called at mesh construction
135 i.e. only used to construct local elements
136 ------------------------------------------------------------------------- */
137
138 template<int NUM_NODES>
addElement(double ** nodeToAdd)139 bool MultiNodeMesh<NUM_NODES>::addElement(double **nodeToAdd)
140 {
141
142 double avg[3];
143
144 if(nodesAreEqual(nodeToAdd[0],nodeToAdd[1]) || nodesAreEqual(nodeToAdd[1],nodeToAdd[2]) ||
145 nodesAreEqual(nodeToAdd[0],nodeToAdd[2]) )
146 return false;
147
148 // add node
149 node_.add(nodeToAdd);
150
151 int n = sizeLocal();
152
153 // calculate center
154 vectorZeroize3D(avg);
155 for(int i = 0; i < NUM_NODES; i++)
156 vectorAdd3D(nodeToAdd[i],avg,avg);
157 vectorScalarDiv3D(avg,static_cast<double>(NUM_NODES));
158 center_.add(avg);
159
160 // extend bbox
161 this->extendToElem(n);
162
163 // calculate rounding radius
164 double rb = 0.;
165 double vec[3];
166 for(int i = 0; i < NUM_NODES; i++)
167 {
168 vectorSubtract3D(center_(n),node_(n)[i],vec);
169 rb = std::max(rb,vectorMag3D(vec));
170 }
171 rBound_.add(rb);
172
173 if(autoRemoveDuplicates_)
174 {
175 for(int i = 0; i < n; i++)
176 {
177
178 if(this->nSharedNodes(i,n) == NUM_NODES)
179 {
180
181 node_.del(n);
182 center_.del(n);
183 rBound_.del(n);
184 return false;
185 }
186 }
187 }
188
189 return true;
190 }
191
192 /* ----------------------------------------------------------------------
193 delete an element - may delete an owned or ghost
194 ------------------------------------------------------------------------- */
195
196 template<int NUM_NODES>
deleteElement(int n)197 void MultiNodeMesh<NUM_NODES>::deleteElement(int n)
198 {
199 node_.del(n);
200 if(node_orig_) node_orig_->del(n);
201 center_.del(n);
202 rBound_.del(n);
203
204 // do not re-calc bbox here
205 }
206
207 /* ----------------------------------------------------------------------
208 recalculate properties on setup (on start and during simulation)
209 ------------------------------------------------------------------------- */
210
211 template<int NUM_NODES>
refreshOwned(int setupFlag)212 void MultiNodeMesh<NUM_NODES>::refreshOwned(int setupFlag)
213 {
214 int ilo = 0, ihi = sizeLocal();
215
216 if(isDeforming())
217 updateCenterRbound(ilo,ihi);
218
219 storeNodePosRebuild();
220
221 if(node_orig_ && setupFlag)
222 storeNodePosOrig(ilo,ihi);
223
224 // nothing more to do here, necessary initialitation done in addElement()
225 }
226
227 template<int NUM_NODES>
refreshGhosts(int setupFlag)228 void MultiNodeMesh<NUM_NODES>::refreshGhosts(int setupFlag)
229 {
230 int ilo = sizeLocal(), ihi = sizeLocal()+sizeGhost();
231
232 if(isDeforming())
233 updateCenterRbound(ilo,ihi);
234
235 if(node_orig_ && setupFlag)
236 storeNodePosOrig(ilo,ihi);
237
238 // nothing more to do here, necessary initialitation done in addElement()
239 }
240
241 /* ----------------------------------------------------------------------
242 comparison
243 ------------------------------------------------------------------------- */
244
245 template<int NUM_NODES>
nodesAreEqual(int iElem,int iNode,int jElem,int jNode)246 bool MultiNodeMesh<NUM_NODES>::nodesAreEqual(int iElem, int iNode, int jElem, int jNode)
247 {
248 for(int i=0;i<3;i++)
249 if(!MathExtraLiggghts::compDouble(node_(iElem)[iNode][i],node_(jElem)[jNode][i],precision_))
250 return false;
251 return true;
252 }
253
254 template<int NUM_NODES>
nodesAreEqual(double * nodeToCheck1,double * nodeToCheck2)255 bool MultiNodeMesh<NUM_NODES>::nodesAreEqual(double *nodeToCheck1,double *nodeToCheck2)
256 {
257 for(int i=0;i<3;i++)
258 if(!MathExtraLiggghts::compDouble(nodeToCheck1[i],nodeToCheck2[i],precision_))
259 return false;
260 return true;
261 }
262
263 template<int NUM_NODES>
containsNode(int iElem,double * nodeToCheck)264 int MultiNodeMesh<NUM_NODES>::containsNode(int iElem, double *nodeToCheck)
265 {
266 for(int iNode = 0; iNode < NUM_NODES; iNode++)
267 {
268 if(MathExtraLiggghts::compDouble(node_(iElem)[iNode][0],nodeToCheck[0],precision_) &&
269 MathExtraLiggghts::compDouble(node_(iElem)[iNode][1],nodeToCheck[1],precision_) &&
270 MathExtraLiggghts::compDouble(node_(iElem)[iNode][2],nodeToCheck[2],precision_))
271 return iNode;
272 }
273 return -1;
274 }
275
276 /* ----------------------------------------------------------------------
277 return if elemens share node, returns lowest iNode and corresponding jNode
278 ------------------------------------------------------------------------- */
279
280 template<int NUM_NODES>
share2Nodes(int iElem,int jElem,int & iNode1,int & jNode1,int & iNode2,int & jNode2)281 bool MultiNodeMesh<NUM_NODES>::share2Nodes(int iElem, int jElem,
282 int &iNode1, int &jNode1, int &iNode2, int &jNode2)
283 {
284 // broad phase
285 double dist[3], radsum;
286 int nShared = 0;
287 vectorSubtract3D(center_(iElem),center_(jElem),dist);
288 radsum = rBound_(iElem) + rBound_(jElem);
289
290 if(vectorMag3DSquared(dist) > radsum*radsum)
291 {
292 iNode1 = jNode1 = iNode2 = jNode2 = -1;
293
294 return false;
295 }
296
297 // narrow phase
298 for(int i=0;i<NUM_NODES;i++){
299 for(int j=0;j<NUM_NODES;j++){
300 if(MultiNodeMesh<NUM_NODES>::nodesAreEqual(iElem,i,jElem,j)){
301 if(0 == nShared)
302 {
303 iNode1 = i;
304 jNode1 = j;
305 }
306 else
307 {
308 iNode2 = i;
309 jNode2 = j;
310
311 return true;
312 }
313 nShared++;
314 }
315 }
316 }
317
318 iNode1 = jNode1 = iNode2 = jNode2 = -1;
319
320 return false;
321 }
322
323 /* ----------------------------------------------------------------------
324 return the number of shared nodes
325 ------------------------------------------------------------------------- */
326
327 template<int NUM_NODES>
nSharedNodes(int iElem,int jElem)328 int MultiNodeMesh<NUM_NODES>::nSharedNodes(int iElem, int jElem)
329 {
330 double dist[3], radsum;
331 int nShared = 0;
332
333 // broad phase
334 vectorSubtract3D(center_(iElem),center_(jElem),dist);
335 radsum = rBound_(iElem) + rBound_(jElem);
336 if(vectorMag3DSquared(dist) > radsum*radsum)
337 return 0;
338
339 // narrow phase
340 for(int i=0;i<NUM_NODES;i++)
341 for(int j=0;j<NUM_NODES;j++)
342 if(MultiNodeMesh<NUM_NODES>::nodesAreEqual(iElem,i,jElem,j))
343 nShared++;
344
345 return nShared;
346 }
347
348 /* ----------------------------------------------------------------------
349 register and unregister mesh movement
350 on registration, return bool staing if this is first mover on this mesh
351 ------------------------------------------------------------------------- */
352
353 template<int NUM_NODES>
registerMove(bool _scale,bool _translate,bool _rotate)354 bool MultiNodeMesh<NUM_NODES>::registerMove(bool _scale, bool _translate, bool _rotate)
355 {
356 bool isFirst = true;
357 if(nMove_ > 0)
358 isFirst = false;
359
360 nMove_ ++;
361 if(_scale) nScale_++;
362 if(_translate) nTranslate_++;
363 if(_rotate) nRotate_++;
364
365 if(isFirst)
366 {
367 int nall = sizeLocal()+sizeGhost();
368
369 double **tmp;
370 this->memory->template create<double>(tmp,NUM_NODES,3,"MultiNodeMesh:tmp");
371
372 if(node_orig_ || (0 == nall && 0 == sizeGlobal()))
373 error->one(FLERR,"Illegal situation in MultiNodeMesh<NUM_NODES>::registerMove");
374
375 node_orig_ = new MultiVectorContainer<double,NUM_NODES,3>("node_orig");
376 node_orig_->setWrapPeriodic(true);
377 for(int i = 0; i < nall; i++)
378 {
379 for(int j = 0; j < NUM_NODES; j++)
380 vectorCopy3D(node_(i)[j],tmp[j]);
381
382 node_orig_->add(tmp);
383 }
384
385 this->memory->template destroy<double>(tmp);
386 }
387
388 return isFirst;
389 }
390
391 template<int NUM_NODES>
unregisterMove(bool _scale,bool _translate,bool _rotate)392 void MultiNodeMesh<NUM_NODES>::unregisterMove(bool _scale, bool _translate, bool _rotate)
393 {
394 nMove_ --;
395 if(_scale) nScale_--;
396 if(_translate) nTranslate_--;
397 if(_rotate) nRotate_--;
398
399 bool del = true;
400 if(nMove_ > 0)
401 del = false;
402
403 if(del)
404 {
405 delete node_orig_;
406 node_orig_ = NULL;
407 }
408 }
409
410 /* ----------------------------------------------------------------------
411 store current node position as original node position for use by moving mesh
412 ------------------------------------------------------------------------- */
413
414 template<int NUM_NODES>
storeNodePosOrig(int ilo,int ihi)415 void MultiNodeMesh<NUM_NODES>::storeNodePosOrig(int ilo, int ihi)
416 {
417 if(!node_orig_)
418 error->one(FLERR,"Internal error in MultiNodeMesh<NUM_NODES>::storeNodePosOrig");
419
420 int nall = this->sizeLocal()+this->sizeGhost();
421 int capacity = this->node_orig_->capacity();
422 if(capacity < nall)
423 this->node_orig_->addUninitialized(nall - capacity);
424
425 for(int i = ilo; i < ihi; i++)
426 for(int j = 0; j < NUM_NODES; j++)
427 {
428 vectorCopy3D(node_(i)[j],node_orig(i)[j]);
429
430 }
431 }
432
433 /* ----------------------------------------------------------------------
434 reset mesh nodes to original position, done before movements are added
435 ------------------------------------------------------------------------- */
436
437 template<int NUM_NODES>
resetToOrig()438 bool MultiNodeMesh<NUM_NODES>::resetToOrig()
439 {
440 if(!node_orig_)
441 error->all(FLERR,"Internal error in MultiNodeMesh<NUM_NODES>::resetNodesToOrig");
442
443 int ntimestep = update->ntimestep;
444
445 if(stepLastReset_ < ntimestep)
446 {
447 int nall = sizeLocal() + sizeGhost();
448 stepLastReset_ = ntimestep;
449 for(int i = 0; i < nall; i++)
450 for(int j = 0; j < NUM_NODES; j++)
451 vectorCopy3D(node_orig(i)[j],node_(i)[j]);
452
453 return true;
454 }
455 return false;
456 }
457
458 /* ----------------------------------------------------------------------
459 move mesh by amount vecTotal, starting from original position
460 ------------------------------------------------------------------------- */
461
462 template<int NUM_NODES>
move(const double * const vecTotal,const double * const vecIncremental)463 void MultiNodeMesh<NUM_NODES>::move(const double * const vecTotal, const double * const vecIncremental)
464 {
465 if(!isTranslating())
466 this->error->all(FLERR,"Illegal call, need to register movement first");
467
468 const int n = sizeLocal() + sizeGhost();
469
470 resetToOrig();
471
472 for(int i = 0; i < n; i++)
473 {
474 vectorZeroize3D(center_(i));
475
476 for(int j = 0; j < NUM_NODES; j++)
477 {
478 vectorAdd3D(node_(i)[j],vecTotal,node_(i)[j]);
479 vectorAdd3D(node_(i)[j],center_(i),center_(i));
480 }
481 vectorScalarDiv3D(center_(i),static_cast<double>(NUM_NODES));
482 }
483
484 if (store_vel)
485 {
486 if (step_store_vel != update->ntimestep)
487 {
488 step_store_vel = update->ntimestep;
489 vectorZeroize3D(global_vel);
490 }
491 vectorAddMultiple3D(global_vel, 1.0/update->dt, vecIncremental, global_vel);
492 }
493
494 updateGlobalBoundingBox();
495 }
496
497 /* ----------------------------------------------------------------------
498 move mesh incrementally by amount vecIncremental
499 ------------------------------------------------------------------------- */
500
501 template<int NUM_NODES>
move(const double * const vecIncremental)502 void MultiNodeMesh<NUM_NODES>::move(const double * const vecIncremental)
503 {
504
505 int n = sizeLocal() + sizeGhost();
506
507 for(int i = 0; i < n; i++)
508 {
509 for(int j = 0; j < NUM_NODES; j++)
510 vectorAdd3D(node_(i)[j],vecIncremental,node_(i)[j]);
511
512 vectorAdd3D(center_(i),vecIncremental,center_(i));
513 }
514
515 if (store_vel)
516 {
517 if (step_store_vel != update->ntimestep)
518 {
519 step_store_vel = update->ntimestep;
520 vectorZeroize3D(global_vel);
521 }
522 vectorAddMultiple3D(global_vel, 1.0/update->dt, vecIncremental, global_vel);
523 }
524
525 updateGlobalBoundingBox();
526 }
527 /* ----------------------------------------------------------------------
528 move mesh incrementally by amount vecIncremental
529 ------------------------------------------------------------------------- */
530
531 template<int NUM_NODES>
moveElement(const int i,const double * const vecIncremental)532 void MultiNodeMesh<NUM_NODES>::moveElement(const int i, const double * const vecIncremental)
533 {
534 for(int j = 0; j < NUM_NODES; j++)
535 vectorAdd3D(node_(i)[j],vecIncremental,node_(i)[j]);
536
537 vectorAdd3D(center_(i),vecIncremental,center_(i));
538
539 extendToElem(bbox_,i);
540 }
541
542 /* ----------------------------------------------------------------------
543 rotate mesh interface, takes both total angle and dAngle in rad
544 assumes axis stays the same over time
545 ------------------------------------------------------------------------- */
546
547 template<int NUM_NODES>
rotate(const double totalAngle,const double dAngle,const double * const axis,const double * const p)548 void MultiNodeMesh<NUM_NODES>::rotate(const double totalAngle, const double dAngle, const double * const axis, const double * const p)
549 {
550 double totalQ[4],dQ[4], axisNorm[3], origin[3];
551
552 // rotates around axis through p
553
554 // normalize axis
555 vectorCopy3D(axis,axisNorm);
556 vectorScalarDiv3D(axisNorm,vectorMag3D(axisNorm));
557
558 // quat for total rotation from original position
559 totalQ[0] = cos(totalAngle*0.5);
560 for(int i=0;i<3;i++)
561 totalQ[i+1] = axis[i]*sin(totalAngle*0.5);
562
563 // quat for rotation since last time-step
564 dQ[0] = cos(dAngle*0.5);
565 for(int i = 0; i < 3; i++)
566 dQ[i+1] = axis[i]*sin(dAngle*0.5);
567
568 vectorCopy3D(p,origin);
569
570 // apply rotation around center axis + displacement
571 // = rotation around axis through p
572 rotate(totalQ,dQ,origin);
573 }
574
575 template<int NUM_NODES>
rotate(const double * const totalQ,const double * const dQ,const double * const origin)576 void MultiNodeMesh<NUM_NODES>::rotate(const double * const totalQ, const double * const dQ, const double * const origin)
577 {
578 if(!isRotating())
579 this->error->all(FLERR,"Illegal call, need to register movement first");
580
581 resetToOrig();
582
583 int n = sizeLocal() + sizeGhost();
584
585 bool trans = vectorMag3DSquared(origin) > 0.;
586
587 // perform total rotation for data in this class
588 for(int i = 0; i < n; i++)
589 {
590 vectorZeroize3D(center_(i));
591
592 for(int j = 0; j < NUM_NODES; j++)
593 {
594 if(trans) vectorSubtract3D(node_(i)[j],origin,node_(i)[j]);
595 MathExtraLiggghts::vec_quat_rotate(node_(i)[j], totalQ, node_(i)[j]);
596 if(trans) vectorAdd3D(node_(i)[j],origin,node_(i)[j]);
597 vectorAdd3D(node_(i)[j],center_(i),center_(i));
598 }
599 vectorScalarDiv3D(center_(i),static_cast<double>(NUM_NODES));
600 }
601
602 if (store_omega)
603 {
604 if (step_store_omega != update->ntimestep)
605 {
606 step_store_omega = update->ntimestep;
607 vectorCopy4D(global_quaternion, prev_quaternion);
608 }
609 vectorCopy4D(totalQ, global_quaternion);
610 }
611
612 updateGlobalBoundingBox();
613 }
614
615 /* ----------------------------------------------------------------------
616 rotate mesh interface, takes only dAngle
617 ------------------------------------------------------------------------- */
618
619 template<int NUM_NODES>
rotate(const double dAngle,const double * const axis,const double * const p)620 void MultiNodeMesh<NUM_NODES>::rotate(const double dAngle, const double * const axis, const double * const p)
621 {
622 double dQ[4], axisNorm[3], origin[3];
623
624 // rotates around axis through p
625
626 // normalize axis
627 vectorCopy3D(axis,axisNorm);
628 vectorScalarDiv3D(axisNorm,vectorMag3D(axisNorm));
629
630 // quat for rotation since last time-step
631 dQ[0] = cos(dAngle*0.5);
632 for(int i = 0; i < 3; i++)
633 dQ[i+1] = axisNorm[i]*sin(dAngle*0.5);
634
635 vectorCopy3D(p,origin);
636
637 // apply rotation around center axis + displacement
638 // = rotation around axis through p
639 rotate(dQ,origin);
640 }
641
642 template<int NUM_NODES>
rotate(const double * const dQ,const double * const origin)643 void MultiNodeMesh<NUM_NODES>::rotate(const double * const dQ, const double * const origin)
644 {
645
646 int n = sizeLocal() + sizeGhost();
647
648 bool trans = vectorMag3DSquared(origin) > 0.;
649
650 // perform total rotation for data in this class
651
652 for(int i = 0; i < n; i++)
653 {
654 vectorZeroize3D(center_(i));
655 for(int j = 0; j < NUM_NODES; j++)
656 {
657 if(trans) vectorSubtract3D(node_(i)[j],origin,node_(i)[j]);
658 MathExtraLiggghts::vec_quat_rotate(node_(i)[j], dQ,node_(i)[j]);
659 if(trans) vectorAdd3D(node_(i)[j],origin,node_(i)[j]);
660 vectorAdd3D(node_(i)[j],center_(i),center_(i));
661 }
662 vectorScalarDiv3D(center_(i),static_cast<double>(NUM_NODES));
663 }
664
665 if (store_omega)
666 {
667 if (step_store_omega != update->ntimestep)
668 {
669 step_store_omega = update->ntimestep;
670 vectorCopy4D(global_quaternion, prev_quaternion);
671 }
672 quatMult4D(global_quaternion, dQ);
673 }
674
675 updateGlobalBoundingBox();
676 }
677
678 /* ----------------------------------------------------------------------
679 scale mesh
680 ------------------------------------------------------------------------- */
681
682 template<int NUM_NODES>
scale(double factor)683 void MultiNodeMesh<NUM_NODES>::scale(double factor)
684 {
685
686 int n = sizeLocal() + sizeGhost();
687
688 for(int i = 0; i < n; i++)
689 {
690 vectorZeroize3D(center_(i));
691 for(int j = 0; j < NUM_NODES; j++)
692 {
693 node_(i)[j][0] *= factor;
694 node_(i)[j][1] *= factor;
695 node_(i)[j][2] *= factor;
696 vectorAdd3D(node_(i)[j],center_(i),center_(i));
697 }
698 vectorScalarDiv3D(center_(i),static_cast<double>(NUM_NODES));
699
700 // calculate rounding radius
701 double rb = 0.;
702 double vec[3];
703 for(int j = 0; j < NUM_NODES; j++)
704 {
705 vectorSubtract3D(center_(i),node_(i)[j],vec);
706 rb = std::max(rb,vectorMag3D(vec));
707 }
708 rBound_(i) = rb;
709 }
710
711 updateGlobalBoundingBox();
712 }
713
714 /* ----------------------------------------------------------------------
715 update center and rbound from node data
716 ------------------------------------------------------------------------- */
717
718 template<int NUM_NODES>
updateCenterRbound(int ilo,int ihi)719 void MultiNodeMesh<NUM_NODES>::updateCenterRbound(int ilo, int ihi)
720 {
721 for(int i = ilo; i < ihi; i++)
722 {
723 vectorZeroize3D(center_(i));
724 for(int j = 0; j < NUM_NODES; j++)
725 vectorAdd3D(node_(i)[j],center_(i),center_(i));
726 vectorScalarDiv3D(center_(i),static_cast<double>(NUM_NODES));
727
728 // calculate rounding radius
729 double rb = 0.;
730 double vec[3];
731 for(int j = 0; j < NUM_NODES; j++)
732 {
733 vectorSubtract3D(center_(i),node_(i)[j],vec);
734 rb = std::max(rb,vectorMag3D(vec));
735 }
736 rBound_(i) = rb;
737 }
738
739 updateGlobalBoundingBox();
740 }
741
742 /* ----------------------------------------------------------------------
743 bounding box funtions
744 ------------------------------------------------------------------------- */
745
746 template<int NUM_NODES>
getElementBoundingBoxOnSubdomain(int const n)747 BoundingBox MultiNodeMesh<NUM_NODES>::getElementBoundingBoxOnSubdomain(int const n)
748 {
749 BoundingBox ret;
750 extendToElem(ret,n);
751 ret.shrinkToSubbox(this->domain->sublo,this->domain->subhi);
752 return ret;
753 }
754
755 template<int NUM_NODES>
getGlobalBoundingBox()756 BoundingBox MultiNodeMesh<NUM_NODES>::getGlobalBoundingBox() const
757 {
758 return bbox_;
759 }
760
761 template<int NUM_NODES>
updateGlobalBoundingBox()762 void MultiNodeMesh<NUM_NODES>::updateGlobalBoundingBox()
763 {
764 bbox_.reset();
765
766 int n = sizeLocal();
767
768 for(int i = 0; i < n; i++)
769 extendToElem(bbox_,i);
770 bbox_.extendToParallel(this->world);
771 }
772
773 template<int NUM_NODES>
extendToElem(int const nElem)774 void MultiNodeMesh<NUM_NODES>::extendToElem(int const nElem)
775 {
776 for(int i = 0; i < NUM_NODES; ++i)
777 bbox_.extendToContain(node_(nElem)[i]);
778 }
779
780 template<int NUM_NODES>
extendToElem(BoundingBox & box,int const nElem)781 void MultiNodeMesh<NUM_NODES>::extendToElem(BoundingBox &box, int const nElem)
782 {
783 for(int i = 0; i < NUM_NODES; ++i)
784 box.extendToContain(node_(nElem)[i]);
785 }
786
787 /* ----------------------------------------------------------------------
788 decide if any node has moved far enough to trigger re-build
789 ------------------------------------------------------------------------- */
790
791 template<int NUM_NODES>
decideRebuild()792 bool MultiNodeMesh<NUM_NODES>::decideRebuild()
793 {
794 // just return for non-moving mesh
795 if(!isMoving() && !isDeforming()) return false;
796
797 double ***node = node_.begin();
798 double ***old = nodesLastRe_.begin();
799 int flag = 0;
800 int nlocal = sizeLocal();
801 double triggersq = 0.25*this->neighbor->skin*this->neighbor->skin;
802
803 if(nlocal != nodesLastRe_.size())
804 this->error->one(FLERR,"Internal error in MultiNodeMesh::decide_rebuild()");
805
806 for(int iTri = 0; iTri < nlocal; iTri++)
807 {
808 for(int iNode = 0; iNode < NUM_NODES; iNode++)
809 {
810 double deltaX[3];
811 vectorSubtract3D(node[iTri][iNode],old[iTri][iNode],deltaX);
812 double distSq = deltaX[0]*deltaX[0] + deltaX[1]*deltaX[1] + deltaX[2]*deltaX[2];
813 if(distSq > triggersq){
814
815 flag = 1;
816 }
817 }
818 if (flag) break;
819 }
820
821 // allreduce result
822 MPI_Max_Scalar(flag,this->world);
823
824 if(flag) return true;
825 else return false;
826 }
827
828 /* ----------------------------------------------------------------------
829 store node pos at last re-build
830 ------------------------------------------------------------------------- */
831
832 template<int NUM_NODES>
storeNodePosRebuild()833 void MultiNodeMesh<NUM_NODES>::storeNodePosRebuild()
834 {
835 // just return for non-moving mesh
836 if(!isMoving() && !isDeforming()) return;
837
838 int nlocal = sizeLocal();
839 double ***node = node_.begin();
840
841 nodesLastRe_.clearContainer();
842 for(int i = 0; i < nlocal; i++)
843 nodesLastRe_.add(node[i]);
844 }
845 /* ----------------------------------------------------------------------
846 calculate simple center of mass, NOT weighted with element area
847 ------------------------------------------------------------------------- */
848
849 template<int NUM_NODES>
center_of_mass(double * _com)850 void MultiNodeMesh<NUM_NODES>::center_of_mass(double *_com)
851 {
852 int nlocal = sizeLocal();
853 int nprocs = this->comm->nprocs;
854 vectorZeroize3D(_com);
855
856 for(int i = 0; i < nlocal; i++)
857 vectorAdd3D(_com,center_(i),_com);
858
859 vectorScalarDiv3D(_com,static_cast<double>(nlocal));
860
861 //printVec3D(this->screen,"_com on one proc",_com);
862
863 if(1 < nprocs)
864 {
865 double result[4];
866 vectorCopy3D(_com,result);
867 result[3] = static_cast<double>(nlocal);
868
869 double *result_all;
870 int size_all = MPI_Allgather_Vector(result,4,result_all,this->world);
871
872 if(size_all != 4*nprocs)
873 this->error->one(FLERR,"internal error");
874
875 vectorZeroize3D(_com);
876 double com_weighted[3];
877 double weightsum = 0.;
878 for(int iproc = 0; iproc < nprocs; iproc++)
879 {
880 vectorScalarMult3D(&result_all[iproc*4],static_cast<double>(result_all[iproc*4+3]),com_weighted);
881 weightsum += static_cast<double>(result_all[iproc*4+3]);
882 vectorAdd3D(_com,com_weighted,_com);
883 //printVec3D(this->screen,"com_weighted",com_weighted);
884 //fprintf(this->screen,"weightsum %f\n",weightsum);
885 }
886 vectorScalarDiv3D(_com,weightsum);
887
888 delete []result_all;
889 }
890
891 }
892
893 template<int NUM_NODES>
get_global_vel(double * vel)894 void MultiNodeMesh<NUM_NODES>::get_global_vel(double *vel)
895 {
896 if (!store_vel)
897 return;
898 if (step_store_vel != update->ntimestep)
899 {
900 step_store_vel = update->ntimestep;
901 vectorZeroize3D(global_vel);
902 }
903 vectorCopy3D(global_vel, vel);
904 }
905
906 template<int NUM_NODES>
get_global_omega(double * omega)907 void MultiNodeMesh<NUM_NODES>::get_global_omega(double *omega)
908 {
909 if (!store_omega)
910 return;
911 if (step_store_omega != update->ntimestep)
912 {
913 step_store_omega = update->ntimestep;
914 vectorCopy4D(global_quaternion, prev_quaternion);
915 }
916 double dQ[4];
917 double invPrevQ[4];
918 quatInverse4D(prev_quaternion, invPrevQ);
919 quatMult4D(invPrevQ, global_quaternion, dQ);
920 dQ[0] = fmax(-1.0, fmin(1.0, dQ[0]));
921 const double dAngle = 2.0*acos(dQ[0]);
922 if (fabs(dAngle) > 1e-12)
923 {
924 const double SinHalfdAngle = sin(dAngle*0.5);
925 const double multi = dAngle/(SinHalfdAngle*update->dt);
926 vectorScalarMult3D(&(dQ[1]), multi, omega);
927 }
928 else
929 vectorZeroize3D(omega);
930 }
931
932 #endif
933