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 #include <stdlib.h>
43 #include <string.h>
44 #include <cmath>
45 #include <algorithm>
46 #include "region_mesh_tet.h"
47 #include "lammps.h"
48 #include "bounding_box.h"
49 #include "tri_mesh.h"
50 #include "memory.h"
51 #include "error.h"
52 #include "domain.h"
53 #include "vector_liggghts.h"
54 #include "mpi_liggghts.h"
55 #include "math_extra_liggghts.h"
56 #include "input_mesh_tet.h"
57
58 // include last to ensure correct macros
59 #include "domain_definitions.h"
60
61 #define DELTA_TET 1000
62
63 using namespace LAMMPS_NS;
64
65 /* ---------------------------------------------------------------------- */
66
RegTetMesh(LAMMPS * lmp,int narg,char ** arg)67 RegTetMesh::RegTetMesh(LAMMPS *lmp, int narg, char **arg) :
68 Region(lmp, narg, arg),
69 bounding_box_mesh(*new BoundingBox(BIG,-BIG,BIG,-BIG,BIG,-BIG)),
70 neighList(*new RegionNeighborList<interpolate_no>(lmp)),
71 tri_mesh(*new TriMesh(lmp))
72 {
73 if(narg < 14) error->all(FLERR,"Illegal region mesh/tet command");
74 options(narg-14,&arg[14]);
75
76 if(scaleflag) error->all(FLERR,"Lattice scaling not implemented for region mesh/tet, please use 'units box'");
77
78 //TODO parse all_in yes
79 //TODO disallow !interior and all_in yes
80
81 if(strcmp(arg[2],"file"))
82 error->all(FLERR,"Illegal region mesh/tet command, expecting keyword 'scale'");
83 char *filename = arg[3];
84
85 if(strcmp(arg[4],"scale"))
86 error->all(FLERR,"Illegal region mesh/tet command, expecting keyword 'scale'");
87 scale_fact = atof(arg[5]);
88 if(strcmp(arg[6],"move"))
89 error->all(FLERR,"Illegal region mesh/tet command, expecting keyword 'move'");
90 off_fact[0] = atof(arg[7]);
91 off_fact[1] = atof(arg[8]);
92 off_fact[2] = atof(arg[9]);
93 if(strcmp(arg[10],"rotate"))
94 error->all(FLERR,"Illegal region mesh/tet command, expecting keyword 'rotate'");
95 rot_angle[0] = atof(arg[11]);
96 rot_angle[1] = atof(arg[12]);
97 rot_angle[2] = atof(arg[13]);
98
99 node = NULL;
100 center = NULL;
101 rbound = NULL;
102 rbound_max = 0.;
103
104 n_face_neighs = NULL;
105 face_neighs = NULL;
106 n_face_neighs_node = NULL;
107
108 n_node_neighs = NULL;
109 node_neighs = NULL;
110
111 n_surfaces = NULL;
112 surfaces = NULL;
113
114 volume = NULL;
115 acc_volume = NULL;
116 nTet = 0;
117 nTetMax = 0;
118 total_volume = 0.;
119
120 // manage input
121 InputMeshTet *my_input = new InputMeshTet(lmp, 0, NULL);
122 my_input->meshtetfile(filename,this,true);
123 delete my_input;
124
125 // extent of region and mesh
126
127 set_extent_mesh();
128
129 if (interior) {
130 bboxflag = 1;
131 set_extent_region();
132 tri_mesh.useAsInsertionMesh(false);
133 build_neighs();
134 build_surface();
135 tri_mesh.initialSetup();
136 } else bboxflag = 0;
137
138 cmax = 1;
139 contact = new Contact[cmax];
140 }
141
142 /* ---------------------------------------------------------------------- */
143
~RegTetMesh()144 RegTetMesh::~RegTetMesh()
145 {
146 delete [] contact;
147
148 delete &bounding_box_mesh;
149 delete &neighList;
150 delete &tri_mesh;
151
152 memory->destroy(node);
153 memory->destroy(center);
154 memory->sfree(volume);
155 memory->sfree(acc_volume);
156 }
157
158 /* ----------------------------------------------------------------------
159 inside = 1 if x,y,z is inside or on surface
160 inside = 0 if x,y,z is outside and not on surface
161 ------------------------------------------------------------------------- */
162
inside(double x,double y,double z)163 int RegTetMesh::inside(double x, double y, double z)
164 {
165 double pos[3];
166 pos[0] = x; pos[1] = y; pos[2] = z;
167
168 // check subdomain
169 if(!domain->is_in_subdomain(pos)) return 0;
170
171 // check bbox, only if exists
172 if(bboxflag)
173 {
174 if(pos[0] < extent_xlo || pos[0] > extent_xhi) return 0;
175 if(pos[1] < extent_ylo || pos[1] > extent_yhi) return 0;
176 if(pos[2] < extent_zlo || pos[2] > extent_zhi) return 0;
177 }
178
179 // brute force naive search
180 int inside_mesh = 0;
181 for(int i = 0; i < nTet; i++)
182 {
183 inside_mesh = inside_mesh + is_inside_tet(i,pos);
184 if(inside_mesh) break;
185 }
186
187 //fprintf(screen,"checking pos %f %f %f, result %d; ntet %d\n",x,y,z,inside_mesh,nTet);
188
189 return inside_mesh;
190 }
191
192 /* ---------------------------------------------------------------------- */
193
surface_interior(double * x,double cutoff)194 int RegTetMesh::surface_interior(double *x, double cutoff)
195 {
196 error->one(FLERR,"This feature is not available for tet mesh regions");
197 return 0;
198 }
199
200 /* ---------------------------------------------------------------------- */
201
surface_exterior(double * x,double cutoff)202 int RegTetMesh::surface_exterior(double *x, double cutoff)
203 {
204 error->one(FLERR,"This feature is not available for tet mesh regions");
205 return 0;
206 }
207
208 /* ---------------------------------------------------------------------- */
209
generate_random(double * pos,bool subdomain_flag)210 void RegTetMesh::generate_random(double *pos,bool subdomain_flag)
211 {
212 if(!interior) error->all(FLERR,"Impossible to generate random points on tet mesh region with side = out");
213
214 int ntry = 0;
215
216 do
217 {
218 mesh_randpos(pos);
219 ntry++;
220 }
221 while(ntry < 10000 && (!subdomain_flag || !domain->is_in_subdomain(pos)));
222
223 if(10000 == ntry)
224 error->one(FLERR,"internal error");
225 }
226
227 /* ---------------------------------------------------------------------- */
228
229 // generates a random point within the region and has a min distance from surface
230 // i.e. generate random point in region "shrunk" by cut
generate_random_shrinkby_cut(double * pos,double cut,bool subdomain_flag)231 void RegTetMesh::generate_random_shrinkby_cut(double *pos,double cut,bool subdomain_flag)
232 {
233 int ntry = 0;
234 bool is_near_surface = false;
235 int barysign = -1;
236
237 for(int i = 0; i < nTet; i++)
238 {
239
240 }
241
242 do
243 {
244 ntry++;
245
246 int iTetChosen = mesh_randpos(pos);
247 is_near_surface = false;
248 double delta[3];
249
250 // check all surfaces of chosen tet
251 for(int is = 0; is < n_surfaces[iTetChosen]; is++)
252 {
253 int iSurf = surfaces[iTetChosen][is];
254
255 if(tri_mesh.resolveTriSphereContact(-1,iSurf,cut,pos,delta,barysign) < 0)
256 {
257 is_near_surface = true;
258 break;
259 }
260
261 }
262
263 if(is_near_surface)
264 continue;
265
266 // check all surfaces of face neigh tets
267 for(int iFaceNeigh = 0; iFaceNeigh < n_face_neighs[iTetChosen]; iFaceNeigh++)
268 {
269 for(int is = 0; is < n_surfaces[face_neighs[iTetChosen][iFaceNeigh]]; is++)
270 {
271 int iSurf = surfaces[face_neighs[iTetChosen][iFaceNeigh]][is];
272
273 if(tri_mesh.resolveTriSphereContact(-1,iSurf,cut,pos,delta,barysign) < 0)
274 {
275 is_near_surface = true;
276 break;
277 }
278
279 }
280 }
281 if(is_near_surface)
282 continue;
283
284 // check all surfaces of node neigh tets
285 for(int iNodeNeigh = 0; iNodeNeigh < n_node_neighs[iTetChosen]; iNodeNeigh++)
286 {
287
288 for(int is = 0; is < n_surfaces[node_neighs[iTetChosen][iNodeNeigh]]; is++)
289 {
290 int iSurf = surfaces[node_neighs[iTetChosen][iNodeNeigh]][is];
291
292 if(tri_mesh.resolveTriSphereContact(-1,iSurf,cut,pos,delta,barysign) < 0)
293 {
294 is_near_surface = true;
295 break;
296 }
297
298 }
299 }
300 }
301 // pos has to be within region, and within cut of region surface
302 while(ntry < 10000 && (is_near_surface || (!subdomain_flag || !domain->is_in_subdomain(pos))));
303
304 if(10000 == ntry)
305 {
306 error->one(FLERR,"internal error");
307 }
308 }
309
310 /* ---------------------------------------------------------------------- */
311
add_tet(double ** n)312 void RegTetMesh::add_tet(double **n)
313 {
314 double ctr[3];
315
316 if(nTet == nTetMax) grow_arrays();
317
318 vectorZeroize3D(ctr);
319 for(int i=0;i<4;i++)
320 {
321 vectorCopy3D(n[i],node[nTet][i]);
322 vectorAdd3D(ctr,node[nTet][i],ctr);
323 }
324 vectorScalarDiv3D(ctr,4.);
325 vectorCopy3D(ctr,center[nTet]);
326
327 double vol = volume_of_tet(nTet);
328 if(vol < 0.)
329 {
330 // flip nodes 0 and 3
331 double node0[3];
332 vectorCopy3D(node[nTet][0],node0);
333 vectorCopy3D(node[nTet][3],node[nTet][0]);
334 vectorCopy3D(node0,node[nTet][3]);
335 }
336
337 vol = volume_of_tet(nTet);
338 if(vol < 0.) error->all(FLERR,"Fatal error: RegTetMesh::add_tet: vol < 0");
339
340 volume[nTet] = vol;
341 total_volume += volume[nTet];
342 acc_volume[nTet] = volume[nTet];
343 if(nTet > 0) acc_volume[nTet] += acc_volume[nTet-1];
344 nTet++;
345 }
346
347 /* ---------------------------------------------------------------------- */
348
build_neighs()349 void RegTetMesh::build_neighs()
350 {
351 neighList.reset();
352
353 for(int i = 0; i < nTet; i++)
354 {
355 vectorZeroize3D(center[i]);
356
357 for(int j=0;j<4;j++)
358 vectorAdd3D(node[i][j],center[i],center[i]);
359 vectorScalarDiv3D(center[i],4.);
360
361 double rb = 0., vec[3];
362 for(int j = 0; j < 4; j++)
363 {
364 vectorSubtract3D(center[i],node[i][j],vec);
365 rb = std::max(rb,vectorMag3D(vec));
366 }
367 rbound[i] = rb;
368 if(rb > rbound_max)
369 rbound_max = rb;
370
371 }
372
373 neighList.setBoundingBox(bounding_box_mesh, rbound_max);
374
375 for(int i = 0; i < nTet; i++)
376 {
377 n_face_neighs[i] = 0;
378 n_node_neighs[i] = 0;
379 vectorZeroizeN(n_face_neighs_node[i],4);
380 }
381
382 bool badMesh = false;
383 for(int i = 0; i < nTet; i++)
384 {
385 std::vector<int> overlaps;
386 neighList.hasOverlapWith(center[i], rbound[i],overlaps);
387
388 for(size_t icontainer = 0; icontainer < overlaps.size(); icontainer++)
389 {
390 int iOverlap = overlaps[icontainer];
391
392 if(iOverlap < 0 || iOverlap >= nTet)
393 this->error->one(FLERR,"Mesh error: internal error");
394
395 int nNodesEqual = 0;
396 std::vector<int> iNodesInvolved, iOverlapNodesInvolved;
397
398 for(int iNode = 0; iNode < 4; iNode++)
399 {
400 for(int iOverlapNode = 0; iOverlapNode < 4; iOverlapNode++)
401 {
402 if(nodesAreEqual(node[i][iNode],node[iOverlap][iOverlapNode],1.e-12))
403 {
404 iNodesInvolved.push_back(iNode);
405 iOverlapNodesInvolved.push_back(iOverlapNode);
406 nNodesEqual++;
407 }
408 }
409 }
410
411 if(0 == nNodesEqual)
412 continue;
413
414 else if(3 > nNodesEqual)
415 {
416
417 if(100 == n_node_neighs[i])
418 badMesh = true;
419 else
420 {
421 node_neighs[i][n_node_neighs[i]] = iOverlap;
422 n_node_neighs[i]++;
423 }
424
425 if(100 == n_node_neighs[iOverlap])
426 badMesh = true;
427 else
428 {
429 node_neighs[iOverlap][n_node_neighs[iOverlap]] = i;
430 n_node_neighs[iOverlap]++;
431 }
432 }
433
434 else if(3 == nNodesEqual)
435 {
436
437 face_neighs[i][n_face_neighs[i]++] = iOverlap;
438 face_neighs[iOverlap][n_face_neighs[iOverlap]++] = i;
439
440 n_face_neighs_node[i][iNodesInvolved[0]]++;
441 n_face_neighs_node[i][iNodesInvolved[1]]++;
442 n_face_neighs_node[i][iNodesInvolved[2]]++;
443 n_face_neighs_node[iOverlap][iOverlapNodesInvolved[0]]++;
444 n_face_neighs_node[iOverlap][iOverlapNodesInvolved[1]]++;
445 n_face_neighs_node[iOverlap][iOverlapNodesInvolved[2]]++;
446
447 //TODO: build up the tri mesh here, link the two meshes
448 //TODO: worst case: broad phase (wie auf zettel skizziert)
449 }
450 else
451 {
452
453 char errstr[256];
454
455 sprintf(errstr,"duplicate elements %d and %d in tet for region %s",i,iOverlap,id);
456 error->one(FLERR,errstr);
457 }
458 }
459
460 neighList.insert(center[i], rbound[i],i);
461 }
462 if (badMesh)
463 error->warningAll(FLERR,"Region mesh/tet: too many node neighbors, mesh is of bad quality; 'all_in yes' might not work correctly");
464
465 for(int i = 0; i < nTet; i++)
466 {
467
468 }
469 }
470
471 /* ---------------------------------------------------------------------- */
472
build_surface()473 void RegTetMesh::build_surface()
474 {
475
476 double **nodeTmp = create<double>(nodeTmp,3,3);
477
478 int n_surf_elems = 0;
479 for(int i = 0; i < nTet; i++)
480 {
481 n_surfaces[i] = 0;
482
483 int dummy;
484
485 // 4 neighs
486 if(3 == n_face_neighs_node[i][0] && 3 == n_face_neighs_node[i][1] &&
487 3 == n_face_neighs_node[i][2] && 3 == n_face_neighs_node[i][3])
488 {
489 if(! (4 == n_face_neighs[i]))
490 error->one(FLERR,"assertion failed");
491 continue;
492 }
493
494 // 3 neighs
495 if(3 == n_face_neighs_node[i][0] || 3 == n_face_neighs_node[i][1] ||
496 3 == n_face_neighs_node[i][2] || 3 == n_face_neighs_node[i][3])
497 {
498 if(!(3 == n_face_neighs[i]))
499 error->one(FLERR,"assertion failed");
500
501 int which;
502 MathExtraLiggghts::max(n_face_neighs_node[i],4,which);
503 for(int k=0;k<3;k++){
504 nodeTmp[0][k] = node[i][(which+1)%4][k];
505 nodeTmp[1][k] = node[i][(which+2)%4][k];
506 nodeTmp[2][k] = node[i][(which+3)%4][k];
507 }
508 tri_mesh.addElement(nodeTmp,-1);
509 surfaces[i][n_surfaces[i]] = n_surf_elems;
510 n_surfaces[i]++;
511 n_surf_elems++;
512 }
513 // 0 neighs
514 else if(0 == MathExtraLiggghts::max(n_face_neighs_node[i],4,dummy))
515 {
516 if(!(0 == n_face_neighs[i]))
517 error->one(FLERR,"assertion failed");
518
519 // add 0 1 2 , 1 2 3, 2 3 0, 3 0 1
520 for(int iAdd = 0; iAdd < 4; iAdd++)
521 {
522 for(int k=0;k<3;k++){
523 nodeTmp[0][k] = node[i][(iAdd+0)%4][k];
524 nodeTmp[1][k] = node[i][(iAdd+1)%4][k];
525 nodeTmp[2][k] = node[i][(iAdd+2)%4][k];
526 }
527 tri_mesh.addElement(nodeTmp,-1);
528 surfaces[i][n_surfaces[i]] = n_surf_elems;
529 n_surfaces[i]++;
530 n_surf_elems++;
531 }
532 }
533 // 2 neighs
534 else if( (2 == MathExtraLiggghts::max(n_face_neighs_node[i],4,dummy)) &&
535 (1 == MathExtraLiggghts::min(n_face_neighs_node[i],4,dummy)) &&
536 (6 == vectorSumN(n_face_neighs_node[i],4))
537 )
538 {
539 if(! (2 == n_face_neighs[i]))
540 error->one(FLERR,"assertion failed");
541
542 int which_hi_1;
543 MathExtraLiggghts::max(n_face_neighs_node[i],4,which_hi_1);
544 int which_lo_1;
545 MathExtraLiggghts::min(n_face_neighs_node[i],4,which_lo_1);
546 int which_2 = (which_hi_1+1)%4;
547 if(which_2 == which_lo_1)
548 which_2 = (which_lo_1+1)%4;
549 int which_3 = 6 - (which_hi_1+which_lo_1+which_2);
550
551 int which_lo_2,which_hi_2;
552 if(n_face_neighs_node[i][which_3] > n_face_neighs_node[i][which_2])
553 {
554 which_hi_2 = which_3;
555 which_lo_2 = which_2;
556 }
557 else
558 {
559 which_hi_2 = which_2;
560 which_lo_2 = which_3;
561 }
562
563 // add lo1 lo2 hi1 and lo1 lo2 hi2
564
565 for(int k=0;k<3;k++){
566 nodeTmp[0][k] = node[i][which_lo_1][k];
567 nodeTmp[1][k] = node[i][which_lo_2][k];
568 nodeTmp[2][k] = node[i][which_hi_1][k];
569 }
570 tri_mesh.addElement(nodeTmp,-1);
571 surfaces[i][n_surfaces[i]] = n_surf_elems;
572 n_surfaces[i]++;
573 n_surf_elems++;
574 for(int k=0;k<3;k++){
575 nodeTmp[0][k] = node[i][which_lo_1][k];
576 nodeTmp[1][k] = node[i][which_lo_2][k];
577 nodeTmp[2][k] = node[i][which_hi_2][k];
578 }
579 tri_mesh.addElement(nodeTmp,-1);
580 surfaces[i][n_surfaces[i]] = n_surf_elems;
581 n_surfaces[i]++;
582 n_surf_elems++;
583 }
584 // 1 neighs
585 else if(0 == MathExtraLiggghts::min(n_face_neighs_node[i],4,dummy))
586 {
587 if(! (1 == n_face_neighs[i]))
588 error->one(FLERR,"assertion failed");
589
590 int which_lo = dummy;
591
592 for(int iAdd = 0; iAdd < 3; iAdd++)
593 {
594 for(int k=0;k<3;k++){
595 nodeTmp[0][k] = node[i][(4+which_lo-2+iAdd)%4][k];
596 nodeTmp[1][k] = node[i][(4+which_lo-1+iAdd)%4][k];
597 nodeTmp[2][k] = node[i][(4+which_lo+0+iAdd)%4][k];
598 }
599 tri_mesh.addElement(nodeTmp,-1);
600 surfaces[i][n_surfaces[i]] = n_surf_elems;
601 n_surfaces[i]++;
602 n_surf_elems++;
603 }
604
605 }
606 else
607 error->one(FLERR,"unrecognized");
608 }
609 destroy<double>(nodeTmp);
610 }
611
612 /* ---------------------------------------------------------------------- */
613
nodesAreEqual(double * nodeToCheck1,double * nodeToCheck2,double precision)614 bool RegTetMesh::nodesAreEqual(double *nodeToCheck1,double *nodeToCheck2,double precision)
615 {
616 for(int i=0;i<3;i++)
617 if(!MathExtraLiggghts::compDouble(nodeToCheck1[i],nodeToCheck2[i],precision))
618 return false;
619 return true;
620 }
621
622 /* ---------------------------------------------------------------------- */
623
grow_arrays()624 void RegTetMesh::grow_arrays()
625 {
626 nTetMax += DELTA_TET;
627 node = (double***)(memory->grow(node,nTetMax, 4 , 3, "vtk_tet_node"));
628 center = (double**)(memory->grow(center,nTetMax, 3, "vtk_tet_center"));
629 rbound = (double*)(memory->grow(rbound,nTetMax, "vtk_tet_rbound"));
630
631 n_face_neighs = (int*)(memory->grow(n_face_neighs,nTetMax, "vtk_tet_n_face_neighs"));
632 face_neighs = (int**)(memory->grow(face_neighs,nTetMax,4,"vtk_tet_face_neighs"));
633 n_face_neighs_node = (int**)(memory->grow(n_face_neighs_node,nTetMax,4, "vtk_tet_n_face_neighs_node"));
634
635 n_node_neighs = (int*)(memory->grow(n_node_neighs,nTetMax, "vtk_tet_n_node_neighs"));
636 node_neighs = (int**)(memory->grow(node_neighs,nTetMax,100,"vtk_tet_node_neighs"));
637
638 n_surfaces = (int*)(memory->grow(n_surfaces,nTetMax, "vtk_tet_n_surfaces"));
639 surfaces = (int**)(memory->grow(surfaces,nTetMax,4,"vtk_tet_surfaces"));
640
641 volume = (double*)(memory->srealloc(volume,nTetMax*sizeof(double),"vtk_tet_volume"));
642 acc_volume = (double*)(memory->srealloc(acc_volume,nTetMax*sizeof(double),"vtk_tet_acc_volume"));
643 }
644
645 /* ---------------------------------------------------------------------- */
646
n_tet()647 int RegTetMesh::n_tet()
648 {
649 return nTet;
650 }
651
652 /* ---------------------------------------------------------------------- */
653
total_vol()654 double RegTetMesh::total_vol()
655 {
656 return total_volume;
657 }
658
659 /* ---------------------------------------------------------------------- */
660
tet_vol(int i)661 double RegTetMesh::tet_vol(int i)
662 {
663 return volume[i];
664 }
665
666 /* ---------------------------------------------------------------------- */
667
tet_acc_vol(int i)668 double RegTetMesh::tet_acc_vol(int i)
669 {
670 return acc_volume[i];
671 }
672
673 /* ---------------------------------------------------------------------- */
674
volume_of_tet(int iTet)675 inline double RegTetMesh::volume_of_tet(int iTet)
676 {
677 return volume_of_tet(node[iTet][0],node[iTet][1],node[iTet][2],node[iTet][3]);
678 }
679
680 /* ---------------------------------------------------------------------- */
681
is_inside_tet(int iTet,double * pos)682 inline int RegTetMesh::is_inside_tet(int iTet,double *pos)
683 {
684 double vol1,vol2,vol3,vol4;
685
686 vol1 = volume_of_tet(node[iTet][0], node[iTet][1], node[iTet][2], pos );
687 vol2 = volume_of_tet(node[iTet][0], node[iTet][1], pos, node[iTet][3]);
688 vol3 = volume_of_tet(node[iTet][0], pos, node[iTet][2], node[iTet][3]);
689 vol4 = volume_of_tet(pos , node[iTet][1], node[iTet][2], node[iTet][3]);
690
691 if(vol1 > 0. && vol2 > 0. && vol3 > 0. && vol4 > 0.) return 1;
692 return 0;
693 }
694
695 /* ---------------------------------------------------------------------- */
696
volume_of_tet(double * v0,double * v1,double * v2,double * v3)697 double RegTetMesh::volume_of_tet(double* v0, double* v1, double* v2, double* v3)
698 {
699 double A[3];
700 A[0] = v3[0] - v1[0];
701 A[1] = v3[1] - v1[1];
702 A[2] = v3[2] - v1[2];
703
704 double B[3];
705 B[0] = v2[0] - v1[0];
706 B[1] = v2[1] - v1[1];
707 B[2] = v2[2] - v1[2];
708
709 double C[3];
710 C[0] = v0[0] - v1[0];
711 C[1] = v0[1] - v1[1];
712 C[2] = v0[2] - v1[2];
713
714 // cross product A x B
715 double cp[] = {
716 A[1]*B[2] - A[2]*B[1],
717 -A[0]*B[2] + A[2]*B[0],
718 A[0]*B[1] - A[1]*B[0]
719 };
720
721 // dot with C
722 double volume = cp[0] * C[0] + cp[1] * C[1] + cp[2] * C[2];
723 volume /= 6.;
724 return volume;
725 }
726
727 /* ---------------------------------------------------------------------- */
728
set_extent_mesh()729 inline void RegTetMesh::set_extent_mesh()
730 {
731 for(int i = 0; i < nTet; i++)
732 for(int j=0;j<4;j++)
733 bounding_box_mesh.extendToContain(node[i][j]);
734 }
735
736 /* ---------------------------------------------------------------------- */
737
set_extent_region()738 inline void RegTetMesh::set_extent_region()
739 {
740 extent_xlo = extent_ylo = extent_zlo = BIG;
741 extent_xhi = extent_yhi = extent_zhi = -BIG;
742
743 for(int i = 0; i < nTet; i++)
744 {
745 for(int j=0;j<4;j++)
746 {
747 if(node[i][j][0] < extent_xlo) extent_xlo = node[i][j][0];
748 if(node[i][j][1] < extent_ylo) extent_ylo = node[i][j][1];
749 if(node[i][j][2] < extent_zlo) extent_zlo = node[i][j][2];
750
751 if(node[i][j][0] > extent_xhi) extent_xhi = node[i][j][0];
752 if(node[i][j][1] > extent_yhi) extent_yhi = node[i][j][1];
753 if(node[i][j][2] > extent_zhi) extent_zhi = node[i][j][2];
754 }
755 }
756 }
757
758 /* ---------------------------------------------------------------------- */
759
mesh_randpos(double * pos)760 inline int RegTetMesh::mesh_randpos(double *pos)
761 {
762 int iTriChosen = tet_rand_tri();
763 tet_randpos(iTriChosen,pos);
764 if(pos[0] == 0. && pos[1] == 0. && pos[2] == 0.)
765 error->one(FLERR,"illegal RegTetMesh::mesh_randpos");
766
767 return iTriChosen;
768 }
769
770 /* ---------------------------------------------------------------------- */
771
tet_rand_tri()772 inline int RegTetMesh::tet_rand_tri()
773 {
774 //SIMPLISTIC
775 /*
776 double rd = total_volume * random->uniform();
777 int chosen = 0;
778 while (rd > acc_volume[chosen] && chosen < nTet-1) chosen++;
779 return chosen;
780 */
781
782 // FAST
783
784 double rd = total_volume * random->uniform();
785
786 int i = nTet/2;
787 int imin = 0;
788 int imax = nTet -1;
789 int ntry = 0;
790
791 // binary search
792 do
793 {
794 if(rd < acc_volume[i] && ((i == 0) ? true : (rd > acc_volume[i-1])))
795 return i;
796
797 if(imax == imin)
798 error->one(FLERR,"internal error");
799
800 // must go up
801 if(rd > acc_volume[i])
802 {
803 imin = i;
804 i = (imax+imin) / 2;
805 if (i == imin) i++;
806 }
807 // must go down
808 else if (rd < acc_volume[i-1])
809 {
810 imax = i;
811 i = (imax+imin) / 2;
812 if (i == imin) i++;
813 if (i == imax && i > 0) i--;
814 }
815 else
816 error->one(FLERR,"internal error");
817
818 ntry++;
819 }
820 while(ntry < 10000);
821
822 error->one(FLERR,"internal error");
823 return 0;
824 }
825
826 /* ---------------------------------------------------------------------- */
827
volume_mc(int n_test,bool cutflag,double cut,double & vol_global,double & vol_local)828 void RegTetMesh::volume_mc(int n_test,bool cutflag,double cut,double &vol_global,double &vol_local)
829 {
830 double pos[3], volume_in_local = 0., vol_in_local_all;
831
832 //error->all(FLERR,"end");
833
834 //NO TODO: implementation for cutflag = true
835
836 if(total_volume == 0.) error->all(FLERR,"mesh/tet region has zero volume, cannot continue");
837
838 vol_global = total_volume;
839
840 for(int iTet = 0; iTet < nTet; iTet++)
841 {
842
843 for(int iNode = 0; iNode < 5; iNode++)
844 {
845
846 double weight = (iNode<4) ? (volume[iTet]*0.1) : (volume[iTet]*0.6);
847
848 if(iNode<4)
849 vectorCopy3D(node[iTet][iNode],pos);
850 else
851 vectorCopy3D(center[iTet],pos);
852
853 if(!domain->is_in_domain(pos))
854 error->one(FLERR,"mesh point outside simulation domain");
855
856 // check if point is in subdomain
857 if(domain->is_in_subdomain(pos))
858 volume_in_local += weight;
859 }
860 }
861
862 MPI_Sum_Scalar(volume_in_local,vol_in_local_all,world);
863 if(vol_in_local_all < 1e-13)
864 error->all(FLERR,"Unable to calculate region volume - are you operating on a 2d region?");
865
866 // return calculated values
867 vol_local = volume_in_local;
868
869 // sum of local volumes will not be equal to global volume because of
870 // different random generator states - correct this now
871 vol_local *= (vol_global/vol_in_local_all);
872
873 }
874