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