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_TET_MESH_I_H
43 #define LMP_TET_MESH_I_H
44 
45   /* ----------------------------------------------------------------------
46    calculate volume of a tet
47   ------------------------------------------------------------------------- */
48 
calcVol(int n)49   inline double TetMesh::calcVol(int n)
50   {
51      return calcTetVol(node_(n)[0],node_(n)[1],node_(n)[2],node_(n)[3]);
52   }
53 
calcTetVol(double * v0,double * v1,double * v2,double * v3)54   inline double TetMesh::calcTetVol(double* v0, double* v1, double* v2, double* v3)
55   {
56      double A[3];
57      A[0] = v3[0] - v1[0];
58      A[1] = v3[1] - v1[1];
59      A[2] = v3[2] - v1[2];
60 
61      double B[3];
62      B[0] = v2[0] - v1[0];
63      B[1] = v2[1] - v1[1];
64      B[2] = v2[2] - v1[2];
65 
66      double C[3];
67      C[0] = v0[0] - v1[0];
68      C[1] = v0[1] - v1[1];
69      C[2] = v0[2] - v1[2];
70 
71      // cross product A x B
72      double cp[] = {
73         A[1]*B[2] - A[2]*B[1],
74        -A[0]*B[2] + A[2]*B[0],
75         A[0]*B[1] - A[1]*B[0]
76      };
77 
78      // dot with C
79      double volume = cp[0] * C[0] + cp[1] * C[1] + cp[2] * C[2];
80      volume /= 6.;
81      return volume;
82   }
83 
84   /* ----------------------------------------------------------------------
85    check if point if inside tet
86   ------------------------------------------------------------------------- */
87 
isInside(int iTet,double * pos)88   inline bool TetMesh::isInside(int iTet,double *pos)
89   {
90       double ***node = node_.begin();
91       double vol1,vol2,vol3,vol4;
92 
93       vol1 = calcTetVol(node[iTet][0], node[iTet][1], node[iTet][2], pos          );
94       if(vol1 < 0.) return false;
95       vol2 = calcTetVol(node[iTet][0], node[iTet][1], pos,           node[iTet][3]);
96       if(vol2 < 0.) return false;
97       vol3 = calcTetVol(node[iTet][0], pos,           node[iTet][2], node[iTet][3]);
98       if(vol3 < 0.) return false;
99       vol4 = calcTetVol(pos          , node[iTet][1], node[iTet][2], node[iTet][3]);
100       if(vol4 < 0.) return false;
101 
102       return true;
103   }
104 
105   /* ----------------------------------------------------------------------
106    generates a random point on the surface that lies within my subbox
107   ------------------------------------------------------------------------- */
108 
generateRandomSubbox(double * pos)109   inline int TetMesh::generateRandomSubbox(double *pos)
110   {
111       int index;
112       do
113       {
114           index = generateRandomOwnedGhost(pos);
115       }
116       while(!domain->is_in_subdomain(pos));
117 
118       return index;
119   }
120 
121   /* ----------------------------------------------------------------------
122    generates a random point on the surface that lies within my subbox
123    additionally, point is within a distance of delta from active edges
124    so can generate a point within a distance delta from boundary on planar meshes
125    coarse approch: test only against active edges of chosen triangle and its neigbors
126   ------------------------------------------------------------------------- */
127 
generateRandomSubboxWithin(double * pos,double delta)128   inline int TetMesh::generateRandomSubboxWithin(double *pos,double delta)
129   {
130       // TODO
131       error->all(FLERR,"all_in 'yes' not yet implemented");
132       return 0;
133   }
134 
135   /* ----------------------------------------------------------------------
136    generates a random point on the surface that lies on an owned or ghost element
137   ------------------------------------------------------------------------- */
138 
generateRandomOwnedGhost(double * pos)139   inline int TetMesh::generateRandomOwnedGhost(double *pos)
140   {
141     double s,t,u,tmp,bary_coo[4];
142     int nTet = sizeLocal() + sizeGhost();
143 
144     // step 1 - choose triangle
145     int chosen = randomOwnedGhostElement();
146 
147     if(chosen >= nTet || chosen < 0)
148     {
149 
150         error->one(FLERR,"TriMesh::generate_random error");
151         return -1;
152     }
153 
154     s = random_->uniform();
155     t = random_->uniform();
156     u = random_->uniform();
157 
158     if(s+t > 1.)
159     {
160         s = 1.-s;
161         t = 1.-t;
162     }
163     if(t+u > 1.)
164     {
165         tmp = u;
166         u = 1.-s-t;
167         t = 1.-tmp;
168     }
169     else if(s+t+u > 1.)
170     {
171         tmp = u;
172         u = s+t+u-1.;
173         s = 1.-t-tmp;
174     }
175 
176     bary_coo[0] = 1.-s-t-u;
177     bary_coo[1] = s;
178     bary_coo[2] = t;
179     bary_coo[3] = u;
180 
181     baryToCart(chosen,bary_coo,pos);
182 
183     return chosen;
184   }
185 
baryToCart(int iTet,double * bary_coo,double * pos)186   inline void TetMesh::baryToCart(int iTet,double *bary_coo,double *pos)
187   {
188     double ***node = node_.begin();
189     for(int i=0;i<3;i++)
190        pos[i] = bary_coo[0] * node[iTet][0][i] +
191                 bary_coo[1] * node[iTet][1][i] +
192                 bary_coo[2] * node[iTet][2][i] +
193                 bary_coo[3] * node[iTet][3][i];
194   }
195 
196   /* ----------------------------------------------------------------------
197    check if two elements share a face
198   ------------------------------------------------------------------------- */
199 
shareFace(int i,int j,int & iFace,int & jFace)200   bool TetMesh::shareFace(int i, int j, int &iFace, int &jFace)
201   {
202 
203       //TODO: detect faces of volumes, facenormals etc
204       this->error->all(FLERR,"END");
205       return false;
206   }
207 
208 #endif
209