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_VOLUME_MESH_I_H
43 #define LMP_VOLUME_MESH_I_H
44
45 #define NTRY_MC_VOLUME_MESH_I_H 30000
46 #define NITER_MC_VOLUME_MESH_I_H 5
47 #define TOLERANCE_MC_VOLUME_MESH_I_H 0.05
48
49 /* ----------------------------------------------------------------------
50 constructor(s), destructor
51 ------------------------------------------------------------------------- */
52
53 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
VolumeMesh(LAMMPS * lmp)54 VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::VolumeMesh(LAMMPS *lmp)
55 : TrackingMesh<NUM_NODES>(lmp),
56
57 // TODO should keep volMeshSubdomain up-to-date more often for insertion faces
58 volMesh_(*this->prop().template addGlobalProperty < ScalarContainer<double> > ("volMesh", "comm_none","frame_trans_rot_invariant","restart_no",3)),
59
60 vol_ (*this->prop().template addElementProperty< ScalarContainer<double> > ("vol", "comm_exchange_borders","frame_trans_rot_invariant", "restart_no",3)),
61 volAcc_ (*this->prop().template addElementProperty< ScalarContainer<double> > ("volAcc","comm_exchange_borders","frame_trans_rot_invariant", "restart_no",3)),
62
63 faceNodes_ (*this->prop().template addElementProperty< MultiVectorContainer<int,NUM_FACES,NUM_NODES_PER_FACE> > ("faceNodes ","comm_exchange_borders","frame_invariant", "restart_no")),
64 faceNormals_ (*this->prop().template addElementProperty< MultiVectorContainer<double,NUM_FACES,3> > ("faceNormals ","comm_none","frame_scale_trans_invariant", "restart_no")),
65 isBoundaryFace_ (*this->prop().template addElementProperty< VectorContainer<bool,NUM_FACES> > ("isBoundaryFace","comm_exchange_borders","frame_invariant", "restart_no")),
66
67 nNeighs_ (*this->prop().template addElementProperty< ScalarContainer<int> > ("nNeighs", "comm_exchange_borders","frame_invariant", "restart_no")),
68 neighElems_ (*this->prop().template addElementProperty< VectorContainer<int,NUM_FACES> > ("neighElems", "comm_exchange_borders","frame_invariant", "restart_no"))
69 {
70
71 volMesh_.add(0.);
72 volMesh_.add(0.);
73 volMesh_.add(0.);
74 volMesh_.add(0.);
75
76 this->error->all(FLERR,"have to add neigh topology with comm_exchange_borders");
77 }
78
79 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
~VolumeMesh()80 VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::~VolumeMesh()
81 {}
82
83 /* ----------------------------------------------------------------------
84 add / delete element
85 ------------------------------------------------------------------------- */
86
87 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
addElement(double ** nodeToAdd)88 bool VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::addElement(double **nodeToAdd)
89 {
90 TrackingMesh<NUM_NODES>::addElement(nodeToAdd,-1);
91 this->error->one(FLERR,"TODO line #");
92 this->error->one(FLERR,"TODO auto remove dupl");
93
94 calcVolPropertiesOfNewElement();
95
96 return true;
97 }
98
99 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
deleteElement(int n)100 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::deleteElement(int n)
101 {
102 TrackingMesh<NUM_NODES>::deleteElement(n);
103 }
104
105 /* ----------------------------------------------------------------------
106 calculate properties when adding new element
107 ------------------------------------------------------------------------- */
108
109 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
calcVolPropertiesOfNewElement()110 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::calcVolPropertiesOfNewElement()
111 {
112
113 int f0[3],f1[3],f2[3],f3[3];
114
115 int n = MultiNodeMesh<NUM_NODES>::node_.size()-1;
116
117 // flip node 3 if on the wrong side
118 checkOrientation(n);
119
120 // assign nodes to faces
121 vectorConstruct3D(f0,0,1,2);
122 vectorConstruct3D(f1,0,3,1);
123 vectorConstruct3D(f2,0,2,3);
124 vectorConstruct3D(f1,1,3,2);
125 faceNodes_.set(n,0,f0);
126 faceNodes_.set(n,1,f1);
127 faceNodes_.set(n,2,f2);
128 faceNodes_.set(n,3,f3);
129
130 // calc face normals
131 calcFaceNormals(n);
132
133 // calc volume
134 double vol_elem = calcVol(n);
135 volMesh_(0) += vol_elem;
136 vol_(n) = vol_elem;
137 volAcc_(n) = vol_elem;
138 if(n > 0) volAcc_(n) += volAcc_(n-1);
139 }
140
141 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
checkOrientation(int n)142 inline void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::checkOrientation(int n)
143 {
144 double v01[3],v02[3],v03[3],tmp[3];
145 double **node = this->node_.begin()[n];
146
147 vectorSubtract3D(node[1],node[0],v01);
148 vectorSubtract3D(node[2],node[0],v02);
149 vectorSubtract3D(node[3],node[0],v03);
150
151 vectorCross3D(v01,v02,tmp);
152
153 // if wrong orientation, switch node 0 and 1
154 if(vectorDot3D(tmp,v03) > 0.)
155 {
156 vectorCopy3D(node[0],tmp);
157 vectorCopy3D(node[1],node[0]);
158 vectorCopy3D(tmp,node[1]);
159 }
160 }
161
162 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
calcFaceNormals(int n)163 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::calcFaceNormals(int n)
164 {
165 double v01[3],v02[3],fnormal[3];
166 double **node = this->node_.begin()[n];
167 int **facenodes = this->faceNodes_(n);
168
169 for(int iFace = 0; iFace < NUM_FACES; iFace++)
170 {
171 vectorSubtract3D(node[facenodes[iFace][1]],node[facenodes[iFace][0]],v01);
172 vectorSubtract3D(node[facenodes[iFace][2]],node[facenodes[iFace][0]],v02);
173 vectorCross3D(v01,v02,fnormal);
174 vectorNormalize3D(fnormal);
175
176 faceNormals_.set(n,iFace,fnormal);
177 }
178 }
179
180 /* ----------------------------------------------------------------------
181 recalculate properties on setup (on start and during simulation)
182 ------------------------------------------------------------------------- */
183
184 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
refreshOwned(int setupFlag)185 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::refreshOwned(int setupFlag)
186 {
187 TrackingMesh<NUM_NODES>::refreshOwned(setupFlag);
188 // (re)calculate all properties for owned elements
189
190 recalcLocalVolProperties();
191 }
192 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
refreshGhosts(int setupFlag)193 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::refreshGhosts(int setupFlag)
194 {
195 TrackingMesh<NUM_NODES>::refreshGhosts(setupFlag);
196
197 recalcGhostVolProperties();
198 }
199
200 /* ----------------------------------------------------------------------
201 recalculate properties of local elements
202 ------------------------------------------------------------------------- */
203
204 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
recalcLocalVolProperties()205 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::recalcLocalVolProperties()
206 {
207
208 // volMeshGlobal [volMesh_(0)] and volMeshOwned [volMesh_(1)]
209 // calculated here
210
211 volMesh_(0) = 0.;
212 volMesh_(1) = 0.;
213
214 int nlocal = this->sizeLocal();
215
216 for(int i = 0; i < nlocal; i++)
217 {
218 calcFaceNormals(i);
219
220 vol(i) = calcVol(i);
221 volAcc(i) = vol(i);
222 if(i > 0) volAcc(i) += volAcc(i-1);
223
224 // add to local volume
225 volMesh_(1) += vol(i);
226
227 }
228
229 // mesh vol must be summed up
230 MPI_Sum_Scalar(volMesh_(1),volMesh_(0),this->world);
231
232 }
233
234 /* ----------------------------------------------------------------------
235 recalculate properties of ghost elements
236 ------------------------------------------------------------------------- */
237
238 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
recalcGhostVolProperties()239 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::recalcGhostVolProperties()
240 {
241 double pos[3];
242 int n_succ, n_iter;
243 int nlocal = this->sizeLocal();
244 int nall = this->sizeLocal()+this->sizeGhost();
245
246 // volMeshGhost [volMesh_(2)] and volMeshSubdomain [volMesh_(3)]
247 // calculated here
248
249 // accumulated vol includes owned and ghosts
250 volMesh_(2) = 0.;
251 for(int i = nlocal; i < nall; i++)
252 {
253 calcFaceNormals(i);
254
255 vol(i) = calcVol(i);
256 volAcc(i) = vol(i);
257 if(i > 0) volAcc(i) += volAcc(i-1);
258
259 // add to ghost area
260 volMesh_(2) += vol(i);
261 }
262
263 // calc area of owned and ghost elements in my subdomain
264
265 volMesh_(3) = 0.;
266 double volCheck = 0.;
267
268 if(this->isInsertionMesh())
269 {
270 n_succ = 0;
271 n_iter = 0;
272
273 // iterate long enough so MC has the desired tolerance
274 while( (n_iter < NITER_MC_VOLUME_MESH_I_H) &&
275 (fabs((volCheck-volMeshGlobal()))/volMeshGlobal() > TOLERANCE_MC_VOLUME_MESH_I_H) )
276 {
277 // only generate random positions if I have any mesh elements
278 if(nall)
279 {
280 for(int i = 0; i < NTRY_MC_VOLUME_MESH_I_H; i++)
281 {
282 // pick a random position on owned or ghost element
283 if((generateRandomOwnedGhost(pos) >= 0) && (this->domain->is_in_subdomain(pos)))
284 n_succ++;
285 }
286 }
287 n_iter++;
288 volMesh_(3) = static_cast<double>(n_succ)/static_cast<double>(NTRY_MC_VOLUME_MESH_I_H*n_iter) * (volMeshOwned()+volMeshGhost());
289
290 MPI_Sum_Scalar(volMesh_(3),volCheck,this->world);
291 }
292
293 if(fabs((volCheck-volMeshGlobal()))/volMeshGlobal() > TOLERANCE_MC_VOLUME_MESH_I_H)
294 this->error->all(FLERR,"Local mesh volume calculation failed, try boosting NITER_MC_VOLUME_MESH_I_H");
295
296 // correct so sum of all owned vols is equal to global area
297 volMesh_(3) *= volMeshGlobal()/volCheck;
298 }
299
300 }
301
302 /* ----------------------------------------------------------------------
303 generate a random Element by volAcc
304 ------------------------------------------------------------------------- */
305
306 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
randomOwnedGhostElement()307 int VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::randomOwnedGhostElement()
308 {
309
310 if(!this->isInsertionMesh()) this->error->one(FLERR,"Illegal call for non-insertion mesh");
311 double r = this->random_->uniform() * (volMeshOwned()+volMeshGhost());
312 int nall = this->sizeLocal()+this->sizeGhost()-1;
313 return searchElementByVolAcc(r,0,nall);
314 }
315
316 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
searchElementByVolAcc(double vol,int lo,int hi)317 int VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::searchElementByVolAcc(double vol,int lo, int hi)
318 {
319 if( (lo < 1 || vol > volAcc(lo-1)) && (vol <= volAcc(lo)) )
320 return lo;
321 if( (hi < 1 || vol > volAcc(hi-1)) && (vol <= volAcc(hi)) )
322 return hi;
323
324 int mid = static_cast<int>((lo+hi)/2);
325 if(vol > volAcc(mid))
326 return searchElementByVolAcc(vol,mid,hi);
327 else
328 return searchElementByVolAcc(vol,lo,mid);
329 }
330
331 /* ----------------------------------------------------------------------
332 build neighlist, generate mesh topology
333 ------------------------------------------------------------------------- */
334
335 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
buildNeighbours()336 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::buildNeighbours()
337 {
338 int iFace,jFace;
339
340 // iterate over all elems, over ghosts as well
341 int nall = this->sizeLocal()+this->sizeGhost();
342
343 // inititalize neigh topology - reset to default, ~n
344
345 int neighs[NUM_FACES];
346 bool isb[NUM_FACES];
347 for(int i=0;i<NUM_FACES;i++)
348 {
349 neighs[i] = -1;
350 isb[i] = true;
351 }
352
353 for(int i = 0; i < nall; i++)
354 {
355 nNeighs_.set(i,0);
356 neighElems_.set(i,neighs);
357 isBoundaryFace_.set(i,isb);
358 }
359
360 // build neigh topology, ~n*n/2
361 for(int i = 0; i < nall; i++)
362 {
363 for(int j = i+1; j < nall; j++)
364 {
365
366 if(0 == this->nSharedNodes(i,j)) continue;
367
368 if(shareFace(i,j,iFace,jFace))
369 {
370 neighElems_(i)[nNeighs_(i)] = this->id(i);
371 neighElems_(j)[nNeighs_(j)] = this->id(j);
372 nNeighs_(i)++;
373 nNeighs_(j)++;
374 isBoundaryFace_(i)[iFace] = false;
375 isBoundaryFace_(j)[jFace] = false;
376 }
377 }
378 }
379 }
380
381 /* ----------------------------------------------------------------------
382 isInside etc
383 ------------------------------------------------------------------------- */
384
385 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
isInside(double * p)386 bool VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::isInside(double *p)
387 {
388 // check subdomain
389 if(!this->domain->is_in_subdomain(p)) return false;
390
391 // check bbox
392 if(!this->bbox_.isInside(p)) return false;
393
394 int nall = this->size();
395
396 // brute force
397 for(int i = 0; i < nall; i++)
398 if(isInside(i,p)) return true;
399
400 return false;
401 }
402
403 /* ----------------------------------------------------------------------
404 move, rotate, scale mesh
405 ------------------------------------------------------------------------- */
406
407 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
move(const double * const vecTotal,const double * const vecIncremental)408 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::move(const double * const vecTotal, const double * const vecIncremental)
409 {
410 TrackingMesh<NUM_NODES>::move(vecTotal,vecIncremental);
411 }
412
413 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
move(const double * const vecIncremental)414 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::move(const double * const vecIncremental)
415 {
416 TrackingMesh<NUM_NODES>::move(vecIncremental);
417 }
418
419 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
scale(double factor)420 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::scale(double factor)
421 {
422 TrackingMesh<NUM_NODES>::scale(factor);
423 }
424
425 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
rotate(const double * const totalQ,const double * const dQ,const double * const totalDispl,const double * const dDisp)426 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::rotate(const double * const totalQ, const double * const dQ, const double * const totalDispl, const double * const dDisp)
427 {
428 TrackingMesh<NUM_NODES>::rotate(totalQ,dQ,totalDispl,dDisp);
429
430 }
431
432 template<int NUM_NODES,int NUM_FACES,int NUM_NODES_PER_FACE>
rotate(const double * const dQ,const double * const dDispl)433 void VolumeMesh<NUM_NODES,NUM_FACES,NUM_NODES_PER_FACE>::rotate(const double * const dQ, const double * const dDispl)
434 {
435 TrackingMesh<NUM_NODES>::rotate(dQ,dDispl);
436 }
437
438 #endif
439