1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle 2 // 3 // See the LICENSE.txt file in the Gmsh root directory for license information. 4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues. 5 6 #ifndef MTETRAHEDRON_H 7 #define MTETRAHEDRON_H 8 9 #include "MElement.h" 10 11 /* 12 * MTetrahedron 13 * 14 * v 15 * . 16 * ,/ 17 * / 18 * 2 19 * ,/|`\ 20 * ,/ | `\ 21 * ,/ '. `\ 22 * ,/ | `\ 23 * ,/ | `\ 24 * 0-----------'.--------1 --> u 25 * `\. | ,/ 26 * `\. | ,/ 27 * `\. '. ,/ 28 * `\. |/ 29 * `3 30 * `\. 31 * ` w 32 * 33 */ 34 class MTetrahedron : public MElement { 35 protected: 36 MVertex *_v[4]; _getEdgeVertices(const int num,std::vector<MVertex * > & v)37 void _getEdgeVertices(const int num, std::vector<MVertex *> &v) const 38 { 39 v[0] = _v[edges_tetra(num, 0)]; 40 v[1] = _v[edges_tetra(num, 1)]; 41 } _getFaceVertices(const int num,std::vector<MVertex * > & v)42 void _getFaceVertices(const int num, std::vector<MVertex *> &v) const 43 { 44 v[0] = _v[faces_tetra(num, 0)]; 45 v[1] = _v[faces_tetra(num, 1)]; 46 v[2] = _v[faces_tetra(num, 2)]; 47 } 48 49 public: 50 MTetrahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num = 0, 51 int part = 0) MElement(num,part)52 : MElement(num, part) 53 { 54 _v[0] = v0; 55 _v[1] = v1; 56 _v[2] = v2; 57 _v[3] = v3; 58 } 59 MTetrahedron(const std::vector<MVertex *> &v, int num = 0, int part = 0) MElement(num,part)60 : MElement(num, part) 61 { 62 for(int i = 0; i < 4; i++) _v[i] = v[i]; 63 } ~MTetrahedron()64 ~MTetrahedron() {} getDim()65 virtual int getDim() const { return 3; } getNumVertices()66 virtual std::size_t getNumVertices() const { return 4; } getVertex(int num)67 virtual MVertex *getVertex(int num) { return _v[num]; } getVertex(int num)68 virtual const MVertex *getVertex(int num) const { return _v[num]; } setVertex(int num,MVertex * v)69 virtual void setVertex(int num, MVertex *v) { _v[num] = v; } getNumEdges()70 virtual int getNumEdges() const { return 6; } getEdge(int num)71 virtual MEdge getEdge(int num) const 72 { 73 return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]); 74 } numEdge2numVertex(int numEdge,int numVert)75 virtual int numEdge2numVertex(int numEdge, int numVert) const 76 { 77 return edges_tetra(numEdge, numVert); 78 } getNumEdgesRep(bool curved)79 virtual int getNumEdgesRep(bool curved) { return 6; } 80 virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, 81 SVector3 *n); getEdgeVertices(const int num,std::vector<MVertex * > & v)82 virtual void getEdgeVertices(const int num, std::vector<MVertex *> &v) const 83 { 84 v.resize(2); 85 _getEdgeVertices(num, v); 86 } getNumFaces()87 virtual int getNumFaces() { return 4; } getFace(int num)88 virtual MFace getFace(int num) const 89 { 90 return MFace(_v[faces_tetra(num, 0)], _v[faces_tetra(num, 1)], 91 _v[faces_tetra(num, 2)]); 92 } 93 virtual bool getFaceInfo(const MFace &face, int &ithFace, int &sign, 94 int &rot) const; getNumFacesRep(bool curved)95 virtual int getNumFacesRep(bool curved) { return 4; } getFaceRep(bool curved,int num,double * x,double * y,double * z,SVector3 * n)96 virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, 97 SVector3 *n) 98 { 99 _getFaceRep(_v[faces_tetra(num, 0)], _v[faces_tetra(num, 1)], 100 _v[faces_tetra(num, 2)], x, y, z, n); 101 } getFaceVertices(const int num,std::vector<MVertex * > & v)102 virtual void getFaceVertices(const int num, std::vector<MVertex *> &v) const 103 { 104 v.resize(3); 105 _getFaceVertices(num, v); 106 } getType()107 virtual int getType() const { return TYPE_TET; } getTypeForMSH()108 virtual int getTypeForMSH() const { return MSH_TET_4; } getTypeForUNV()109 virtual int getTypeForUNV() const { return 111; } // solid linear tetrahedron getTypeForVTK()110 virtual int getTypeForVTK() const { return 10; } getStringForPOS()111 virtual const char *getStringForPOS() const { return "SS"; } getStringForBDF()112 virtual const char *getStringForBDF() const { return "CTETRA"; } getStringForDIFF()113 virtual const char *getStringForDIFF() const { return "ElmT4n3D"; } getStringForINP()114 virtual const char *getStringForINP() const { return "C3D4"; } getStringForKEY()115 virtual const char *getStringForKEY() const { return "_SOLID"; } getStringForTOCHNOG()116 virtual const char *getStringForTOCHNOG() const { return "-tet4"; } reverse()117 virtual void reverse() 118 { 119 MVertex *tmp = _v[0]; 120 _v[0] = _v[1]; 121 _v[1] = tmp; 122 } getMat(double mat[3][3])123 void getMat(double mat[3][3]) const 124 { 125 mat[0][0] = _v[1]->x() - _v[0]->x(); 126 mat[0][1] = _v[2]->x() - _v[0]->x(); 127 mat[0][2] = _v[3]->x() - _v[0]->x(); 128 mat[1][0] = _v[1]->y() - _v[0]->y(); 129 mat[1][1] = _v[2]->y() - _v[0]->y(); 130 mat[1][2] = _v[3]->y() - _v[0]->y(); 131 mat[2][0] = _v[1]->z() - _v[0]->z(); 132 mat[2][1] = _v[2]->z() - _v[0]->z(); 133 mat[2][2] = _v[3]->z() - _v[0]->z(); 134 } 135 virtual double getVolume(); getVolumeSign()136 virtual int getVolumeSign() { return (getVolume() >= 0) ? 1 : -1; } 137 virtual double gammaShapeMeasure(); 138 virtual double getInnerRadius(); 139 virtual double getOuterRadius(); 140 virtual double etaShapeMeasure(); 141 virtual void xyz2uvw(double xyz[3], double uvw[3]) const; getNode(int num,double & u,double & v,double & w)142 virtual void getNode(int num, double &u, double &v, double &w) const 143 { 144 switch(num) { 145 case 0: 146 u = 0.; 147 v = 0.; 148 w = 0.; 149 break; 150 case 1: 151 u = 1.; 152 v = 0.; 153 w = 0.; 154 break; 155 case 2: 156 u = 0.; 157 v = 1.; 158 w = 0.; 159 break; 160 case 3: 161 u = 0.; 162 v = 0.; 163 w = 1.; 164 break; 165 default: 166 u = 0.; 167 v = 0.; 168 w = 0.; 169 break; 170 } 171 } barycenterUVW()172 virtual SPoint3 barycenterUVW() const { return SPoint3(.25, .25, .25); } isInside(double u,double v,double w)173 virtual bool isInside(double u, double v, double w) const 174 { 175 double tol = getTolerance(); 176 if(u < (-tol) || v < (-tol) || w < (-tol) || u > ((1. + tol) - v - w)) 177 return false; 178 return true; 179 } 180 virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts); 181 virtual SPoint3 circumcenter(); edges_tetra(const int edge,const int vert)182 static int edges_tetra(const int edge, const int vert) 183 { 184 static const int e[6][2] = {{0, 1}, {1, 2}, {2, 0}, {3, 0}, {3, 2}, {3, 1}}; 185 return e[edge][vert]; 186 } faces_tetra(const int face,const int vert)187 static int faces_tetra(const int face, const int vert) 188 { 189 static const int f[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}}; 190 return f[face][vert]; 191 } faces2edge_tetra(const int face,const int edge)192 static int faces2edge_tetra(const int face, const int edge) 193 { 194 // return -iedge - 1 if edge is inverted 195 // iedge + 1 otherwise 196 static const int e[4][3] = { 197 {-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}}; 198 return e[face][edge]; 199 } 200 virtual int numCommonNodesInDualGraph(const MElement *const other) const; getEdgeSolin(int num)201 virtual MEdge getEdgeSolin(int num) 202 { 203 static const int eSolin[6][2] = {{0, 1}, {1, 2}, {2, 0}, 204 {0, 3}, {2, 3}, {1, 3}}; 205 return MEdge(_v[eSolin[num][0]], _v[eSolin[num][1]]); 206 } getFaceSolin(int num)207 virtual MFace getFaceSolin(int num) 208 { 209 static const int fSolin[4][3] = { 210 {0, 1, 2}, {0, 1, 3}, {0, 2, 3}, {1, 2, 3}}; 211 return MFace(_v[fSolin[num][0]], _v[fSolin[num][1]], _v[fSolin[num][2]]); 212 } 213 }; 214 215 /* 216 * MTetrahedron10 217 * 218 * 2 219 * ,/|`\ 220 * ,/ | `\ 221 * ,6 '. `5 222 * ,/ 8 `\ 223 * ,/ | `\ 224 * 0--------4--'.--------1 225 * `\. | ,/ 226 * `\. | ,9 227 * `7. '. ,/ 228 * `\. |/ 229 * `3 230 * 231 */ 232 class MTetrahedron10 : public MTetrahedron { 233 protected: 234 MVertex *_vs[6]; 235 236 public: 237 MTetrahedron10(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, 238 MVertex *v4, MVertex *v5, MVertex *v6, MVertex *v7, 239 MVertex *v8, MVertex *v9, int num = 0, int part = 0) MTetrahedron(v0,v1,v2,v3,num,part)240 : MTetrahedron(v0, v1, v2, v3, num, part) 241 { 242 _vs[0] = v4; 243 _vs[1] = v5; 244 _vs[2] = v6; 245 _vs[3] = v7; 246 _vs[4] = v8; 247 _vs[5] = v9; 248 for(int i = 0; i < 6; i++) _vs[i]->setPolynomialOrder(2); 249 } 250 MTetrahedron10(const std::vector<MVertex *> &v, int num = 0, int part = 0) MTetrahedron(v,num,part)251 : MTetrahedron(v, num, part) 252 { 253 for(int i = 0; i < 6; i++) _vs[i] = v[4 + i]; 254 for(int i = 0; i < 6; i++) _vs[i]->setPolynomialOrder(2); 255 } ~MTetrahedron10()256 ~MTetrahedron10() {} getPolynomialOrder()257 virtual int getPolynomialOrder() const { return 2; } getNumVertices()258 virtual std::size_t getNumVertices() const { return 10; } getVertex(int num)259 virtual MVertex *getVertex(int num) 260 { 261 return num < 4 ? _v[num] : _vs[num - 4]; 262 } getVertex(int num)263 virtual const MVertex *getVertex(int num) const 264 { 265 return num < 4 ? _v[num] : _vs[num - 4]; 266 } setVertex(int num,MVertex * v)267 virtual void setVertex(int num, MVertex *v) 268 { 269 if(num < 4) 270 _v[num] = v; 271 else 272 _vs[num - 4] = v; 273 } getVertexUNV(int num)274 virtual MVertex *getVertexUNV(int num) 275 { 276 static const int map[10] = {0, 4, 1, 5, 2, 6, 7, 9, 8, 3}; 277 return getVertex(map[num]); 278 } getVertexBDF(int num)279 virtual MVertex *getVertexBDF(int num) 280 { 281 static const int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8}; 282 return getVertex(map[num]); 283 } getVertexVTK(int num)284 virtual MVertex *getVertexVTK(int num) 285 { 286 static const int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8}; 287 return getVertex(map[num]); 288 } getVertexDIFF(int num)289 virtual MVertex *getVertexDIFF(int num) { return getVertexBDF(num); } getVertexINP(int num)290 virtual MVertex *getVertexINP(int num) { return getVertexBDF(num); } getVertexKEY(int num)291 virtual MVertex *getVertexKEY(int num) { return getVertexBDF(num); } getNumEdgeVertices()292 virtual int getNumEdgeVertices() const { return 6; } 293 virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, 294 SVector3 *n); 295 virtual int getNumEdgesRep(bool curved); 296 virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, 297 SVector3 *n); 298 virtual int getNumFacesRep(bool curved); getEdgeVertices(const int num,std::vector<MVertex * > & v)299 virtual void getEdgeVertices(const int num, std::vector<MVertex *> &v) const 300 { 301 v.resize(3); 302 MTetrahedron::_getEdgeVertices(num, v); 303 v[2] = _vs[num]; 304 } getFaceVertices(const int num,std::vector<MVertex * > & v)305 virtual void getFaceVertices(const int num, std::vector<MVertex *> &v) const 306 { 307 v.resize(6); 308 MTetrahedron::_getFaceVertices(num, v); 309 static const int f[4][3] = {{2, 1, 0}, {0, 5, 3}, {3, 4, 2}, {5, 1, 4}}; 310 v[3] = _vs[f[num][0]]; 311 v[4] = _vs[f[num][1]]; 312 v[5] = _vs[f[num][2]]; 313 } getTypeForMSH()314 virtual int getTypeForMSH() const { return MSH_TET_10; } getTypeForUNV()315 virtual int getTypeForUNV() const 316 { 317 return 118; 318 } // solid parabolic tetrahedron getTypeForVTK()319 virtual int getTypeForVTK() const { return 24; } getStringForPOS()320 virtual const char *getStringForPOS() const { return "SS2"; } getStringForBDF()321 virtual const char *getStringForBDF() const { return "CTETRA"; } getStringForDIFF()322 virtual const char *getStringForDIFF() const { return "ElmT10n3D"; } getStringForINP()323 virtual const char *getStringForINP() const { return "C3D10"; } getStringForKEY()324 virtual const char *getStringForKEY() const { return "_SOLID"; } getStringForTOCHNOG()325 virtual const char *getStringForTOCHNOG() const { return "-tet10"; } reverse()326 virtual void reverse() 327 { 328 MVertex *tmp; 329 tmp = _v[0]; 330 _v[0] = _v[1]; 331 _v[1] = tmp; 332 tmp = _vs[1]; 333 _vs[1] = _vs[2]; 334 _vs[2] = tmp; 335 tmp = _vs[5]; 336 _vs[5] = _vs[3]; 337 _vs[3] = tmp; 338 } getNode(int num,double & u,double & v,double & w)339 virtual void getNode(int num, double &u, double &v, double &w) const 340 { 341 num < 4 ? MTetrahedron::getNode(num, u, v, w) : 342 MElement::getNode(num, u, v, w); 343 } xyz2uvw(double xyz[3],double uvw[3])344 void xyz2uvw(double xyz[3], double uvw[3]) const 345 { 346 return MElement::xyz2uvw(xyz, uvw); 347 } 348 }; 349 350 /* tet order 3 FIXME: check the plot 351 * 352 * 2 353 * ,/|`\ 354 * ,8 | `7 E = order - 1 355 * ,/ 13 `\ C = 4 + 6*E 356 * ,9 | `6 F = ((order - 1)*(order - 2))/2 357 * ,/ | `\ N = total number of vertices 358 * 0-----4-----'.--5-----1 359 * `\. | ,/ Interior vertex numbers 360 * 11. 12 ,15 for edge 0 <= i <= 5: 4+i*E to 4+(i+1)*E-1 361 * `\. '. 14 for face 0 <= j <= 3: C+j*F to C+(j+1)*F-1 362 * 10\.|/ in volume : C+4*F to N-1 363 * `3 364 * 365 */ 366 367 typedef std::vector<int> IndicesReversed; 368 369 class MTetrahedronN : public MTetrahedron { 370 static std::map<int, IndicesReversed> _order2indicesReversedTet; 371 372 protected: 373 std::vector<MVertex *> _vs; 374 const char _order; 375 376 public: 377 MTetrahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, 378 const std::vector<MVertex *> &v, char order, int num = 0, 379 int part = 0) MTetrahedron(v0,v1,v2,v3,num,part)380 : MTetrahedron(v0, v1, v2, v3, num, part), _vs(v), _order(order) 381 { 382 for(std::size_t i = 0; i < _vs.size(); i++) 383 _vs[i]->setPolynomialOrder(_order); 384 } 385 MTetrahedronN(const std::vector<MVertex *> &v, char order, int num = 0, 386 int part = 0) MTetrahedron(v[0],v[1],v[2],v[3],num,part)387 : MTetrahedron(v[0], v[1], v[2], v[3], num, part), _order(order) 388 { 389 for(std::size_t i = 4; i < v.size(); i++) _vs.push_back(v[i]); 390 for(std::size_t i = 0; i < _vs.size(); i++) 391 _vs[i]->setPolynomialOrder(_order); 392 } ~MTetrahedronN()393 ~MTetrahedronN() {} getPolynomialOrder()394 virtual int getPolynomialOrder() const { return _order; } getNumVertices()395 virtual std::size_t getNumVertices() const { return 4 + _vs.size(); } getVertex(int num)396 virtual MVertex *getVertex(int num) 397 { 398 return num < 4 ? _v[num] : _vs[num - 4]; 399 } getVertex(int num)400 virtual const MVertex *getVertex(int num) const 401 { 402 return num < 4 ? _v[num] : _vs[num - 4]; 403 } setVertex(int num,MVertex * v)404 virtual void setVertex(int num, MVertex *v) 405 { 406 if(num < 4) 407 _v[num] = v; 408 else 409 _vs[num - 4] = v; 410 } getNumEdgeVertices()411 virtual int getNumEdgeVertices() const { return 6 * (_order - 1); } getNumFaceVertices()412 virtual int getNumFaceVertices() const 413 { 414 if(getIsAssimilatedSerendipity()) 415 return 0; 416 else 417 return 4 * ((_order - 1) * (_order - 2)) / 2; 418 } getEdgeVertices(const int num,std::vector<MVertex * > & v)419 virtual void getEdgeVertices(const int num, std::vector<MVertex *> &v) const 420 { 421 v.resize(_order + 1); 422 MTetrahedron::_getEdgeVertices(num, v); 423 int j = 2; 424 const int ie = (num + 1) * (_order - 1); 425 for(int i = num * (_order - 1); i != ie; ++i) v[j++] = _vs[i]; 426 } getFaceVertices(const int num,std::vector<MVertex * > & v)427 virtual void getFaceVertices(const int num, std::vector<MVertex *> &v) const 428 { 429 if(getIsAssimilatedSerendipity()) { v.resize(3 * _order); } 430 else { 431 v.resize((_order + 1) * (_order + 2) / 2); 432 } 433 434 MTetrahedron::_getFaceVertices(num, v); 435 int count = 2; 436 437 int n = _order - 1; 438 for(int i = 0; i < 3; i++) { 439 if(faces2edge_tetra(num, i) > 0) { 440 int edge_num = faces2edge_tetra(num, i) - 1; 441 for(int j = 0; j < n; j++) v[++count] = _vs[n * edge_num + j]; 442 } 443 else { 444 int edge_num = -faces2edge_tetra(num, i) - 1; 445 for(int j = n - 1; j >= 0; j--) v[++count] = _vs[n * edge_num + j]; 446 } 447 } 448 449 if((int)v.size() > count + 1) { 450 int start = 6 * n + num * (n - 1) * n / 2; 451 for(int i = 0; i < (n - 1) * n / 2; i++) { v[++count] = _vs[start + i]; } 452 } 453 } getNumVolumeVertices()454 virtual int getNumVolumeVertices() const 455 { 456 if(getIsAssimilatedSerendipity()) 457 return 0; 458 else 459 return ((_order - 1) * (_order - 2) * (_order - 3)) / 6; 460 } getTypeForMSH()461 virtual int getTypeForMSH() const 462 { 463 // (p+1)*(p+2)*(p+3)/6 464 if(_order == 1 && _vs.size() + 4 == 4) return MSH_TET_4; 465 if(_order == 2 && _vs.size() + 4 == 10) return MSH_TET_10; 466 if(_order == 3 && _vs.size() + 4 == 20) return MSH_TET_20; 467 if(_order == 4 && _vs.size() + 4 == 35) return MSH_TET_35; 468 if(_order == 5 && _vs.size() + 4 == 56) return MSH_TET_56; 469 if(_order == 6 && _vs.size() + 4 == 84) return MSH_TET_84; 470 if(_order == 7 && _vs.size() + 4 == 120) return MSH_TET_120; 471 if(_order == 8 && _vs.size() + 4 == 165) return MSH_TET_165; 472 if(_order == 9 && _vs.size() + 4 == 220) return MSH_TET_220; 473 if(_order == 10 && _vs.size() + 4 == 286) return MSH_TET_286; 474 475 if(_order == 3 && _vs.size() + 4 == 16) return MSH_TET_16; 476 if(_order == 4 && _vs.size() + 4 == 22) return MSH_TET_22; 477 if(_order == 5 && _vs.size() + 4 == 28) return MSH_TET_28; 478 if(_order == 6 && _vs.size() + 4 == 34) return MSH_TET_34; 479 if(_order == 7 && _vs.size() + 4 == 40) return MSH_TET_40; 480 if(_order == 8 && _vs.size() + 4 == 46) return MSH_TET_46; 481 if(_order == 9 && _vs.size() + 4 == 52) return MSH_TET_52; 482 if(_order == 10 && _vs.size() + 4 == 58) return MSH_TET_58; 483 Msg::Error("No MSH type found for P%d tetrahedron with %d nodes", _order, 484 4 + _vs.size()); 485 return 0; 486 } 487 virtual void reverse(); 488 virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, 489 SVector3 *n); 490 virtual int getNumEdgesRep(bool curved); 491 virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, 492 SVector3 *n); 493 virtual int getNumFacesRep(bool curved); getNode(int num,double & u,double & v,double & w)494 virtual void getNode(int num, double &u, double &v, double &w) const 495 { 496 num < 4 ? MTetrahedron::getNode(num, u, v, w) : 497 MElement::getNode(num, u, v, w); 498 } xyz2uvw(double xyz[3],double uvw[3])499 void xyz2uvw(double xyz[3], double uvw[3]) const 500 { 501 return MElement::xyz2uvw(xyz, uvw); 502 } 503 }; 504 505 #endif 506