NGSolve
4.9
|
00001 #ifndef FILE_TOPOLOGY 00002 #define FILE_TOPOLOGY 00003 00004 /*********************************************************************/ 00005 /* File: topology.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 25. Mar. 2000 */ 00008 /*********************************************************************/ 00009 00010 namespace ngfem 00011 { 00012 00013 /* 00014 Toplogy of reference elements 00015 */ 00016 00017 00022 enum ELEMENT_TYPE { ET_POINT = 0, ET_SEGM = 1, 00023 ET_TRIG = 10, ET_QUAD = 11, 00024 ET_TET = 20, ET_PYRAMID = 21, ET_PRISM = 22, ET_HEX = 24 }; 00025 00026 00027 00028 00029 00042 enum NODE_TYPE { NT_VERTEX = 0, NT_EDGE = 1, NT_FACE = 2, NT_CELL = 3 }; 00043 00044 inline void operator++(NODE_TYPE & nt, int) { nt = NODE_TYPE(nt+1); } 00045 00046 00048 typedef double POINT3D[3]; 00049 00051 typedef int EDGE[2]; 00052 00054 typedef int FACE[4]; 00055 00057 typedef double NORMAL[3]; 00058 00059 00061 class NGS_DLL_HEADER ElementTopology 00062 { 00063 public: 00065 static const char * GetElementName (ELEMENT_TYPE et); 00066 00068 static int GetSpaceDim (ELEMENT_TYPE et) 00069 { 00070 switch (et) 00071 { 00072 case ET_POINT: return 0; 00073 case ET_SEGM: return 1; 00074 case ET_TRIG: return 2; 00075 case ET_QUAD: return 2; 00076 case ET_TET: return 3; 00077 case ET_PYRAMID: return 3; 00078 case ET_PRISM: return 3; 00079 case ET_HEX: return 3; 00080 } 00081 return 0; 00082 } 00083 00085 static int GetNVertices (ELEMENT_TYPE et) 00086 { 00087 switch (et) 00088 { 00089 case ET_POINT: return 1; 00090 case ET_SEGM: return 2; 00091 case ET_TRIG: return 3; 00092 case ET_QUAD: return 4; 00093 case ET_TET: return 4; 00094 case ET_PYRAMID: return 5; 00095 case ET_PRISM: return 6; 00096 case ET_HEX: return 8; 00097 } 00098 return 0; 00099 } 00100 00102 static int GetNEdges (ELEMENT_TYPE et) 00103 { 00104 switch (et) 00105 { 00106 case ET_POINT: return 0; 00107 case ET_SEGM: return 1; 00108 case ET_TRIG: return 3; 00109 case ET_QUAD: return 4; 00110 case ET_TET: return 6; 00111 case ET_PYRAMID: return 8; 00112 case ET_PRISM: return 9; 00113 case ET_HEX: return 12; 00114 } 00115 return 0; 00116 } 00117 00118 00120 static int GetNFaces (ELEMENT_TYPE et) 00121 { 00122 switch (et) 00123 { 00124 case ET_POINT: return 0; 00125 case ET_SEGM: return 0; 00126 case ET_TRIG: return 1; 00127 case ET_QUAD: return 1; 00128 case ET_TET: return 4; 00129 case ET_PYRAMID: return 5; 00130 case ET_PRISM: return 5; 00131 case ET_HEX: return 6; 00132 } 00133 return 0; 00134 } 00135 00136 00138 static ELEMENT_TYPE GetFaceType (ELEMENT_TYPE et, int k) 00139 { 00140 switch (et) 00141 { 00142 case ET_TRIG: return ET_TRIG; 00143 case ET_QUAD: return ET_QUAD; 00144 case ET_TET: return ET_TRIG; 00145 case ET_PYRAMID: return (k<4 ? ET_TRIG : ET_QUAD); 00146 case ET_PRISM: return (k<2 ? ET_TRIG : ET_QUAD); 00147 case ET_HEX: return ET_QUAD; 00148 default: 00149 cerr << "*** error in GetFacetType: Unhandled Elementtype" << endl; 00150 return ET_SEGM; 00151 } 00152 } 00153 00154 00155 00156 00158 static int GetNNodes (ELEMENT_TYPE et, NODE_TYPE nt) 00159 { 00160 static const int nn_point[] = { 1, 0, 0, 0 }; 00161 static const int nn_segm[] = { 2, 1, 0, 0 }; 00162 static const int nn_trig[] = { 3, 3, 1, 0 }; 00163 static const int nn_quad[] = { 4, 4, 1, 0 }; 00164 static const int nn_tet[] = { 4, 6, 4, 1 }; 00165 static const int nn_pyramid[] = { 5, 8, 5, 1 }; 00166 static const int nn_prism[] = { 6, 9, 5, 1 }; 00167 static const int nn_hex[] = { 8, 12, 6, 1 }; 00168 switch (et) 00169 { 00170 case ET_POINT: return nn_point[nt]; 00171 case ET_SEGM: return nn_segm[nt]; 00172 case ET_TRIG: return nn_trig[nt]; 00173 case ET_QUAD: return nn_quad[nt]; 00174 case ET_TET: return nn_tet[nt]; 00175 case ET_PYRAMID: return nn_pyramid[nt]; 00176 case ET_PRISM: return nn_prism[nt]; 00177 case ET_HEX: return nn_hex[nt]; 00178 } 00179 return 0; 00180 } 00181 00182 00183 00184 00185 00186 00187 00189 static int GetNFacets (ELEMENT_TYPE et) 00190 { 00191 switch (et) 00192 { 00193 case ET_POINT: return 0; 00194 case ET_SEGM: return 2; 00195 case ET_TRIG: return 3; 00196 case ET_QUAD: return 4; 00197 case ET_TET: return 4; 00198 case ET_PYRAMID: return 5; 00199 case ET_PRISM: return 5; 00200 case ET_HEX: return 6; 00201 } 00202 return 0; 00203 } 00204 00205 00206 00207 00209 static ELEMENT_TYPE GetFacetType (ELEMENT_TYPE et, int k) 00210 { 00211 switch (et) 00212 { 00213 case ET_SEGM: return ET_POINT; 00214 case ET_TRIG: return ET_SEGM; 00215 case ET_QUAD: return ET_SEGM; 00216 case ET_TET: return ET_TRIG; 00217 case ET_PYRAMID: return (k<4 ? ET_TRIG : ET_QUAD); 00218 case ET_PRISM: return (k<2 ? ET_TRIG : ET_QUAD); 00219 case ET_HEX: return ET_QUAD; 00220 default: 00221 throw Exception ("undefined element type in ElementTopology::GetFacetType()\n"); 00222 } 00223 } 00224 00226 static const POINT3D * GetVertices (ELEMENT_TYPE et); 00227 00229 static const EDGE * GetEdges (ELEMENT_TYPE et) 00230 { 00231 static const int segm_edges[1][2] = 00232 { { 0, 1 }}; 00233 00234 00235 static const int trig_edges[3][2] = 00236 { { 2, 0 }, 00237 { 1, 2 }, 00238 { 0, 1 }}; 00239 00240 static const int quad_edges[4][2] = 00241 { { 0, 1 }, 00242 { 2, 3 }, 00243 { 3, 0 }, 00244 { 1, 2 }}; 00245 00246 static const int tet_edges[6][2] = 00247 { { 3, 0 }, 00248 { 3, 1 }, 00249 { 3, 2 }, 00250 { 0, 1 }, 00251 { 0, 2 }, 00252 { 1, 2 }}; 00253 00254 static const int prism_edges[9][2] = 00255 { { 2, 0 }, 00256 { 0, 1 }, 00257 { 2, 1 }, 00258 { 5, 3 }, 00259 { 3, 4 }, 00260 { 5, 4 }, 00261 { 2, 5 }, 00262 { 0, 3 }, 00263 { 1, 4 }}; 00264 00265 static const int pyramid_edges[8][2] = 00266 { { 0, 1 }, 00267 { 1, 2 }, 00268 { 0, 3 }, 00269 { 3, 2 }, 00270 { 0, 4 }, 00271 { 1, 4 }, 00272 { 2, 4 }, 00273 { 3, 4 }}; 00274 00275 static const int hex_edges[12][2] = 00276 { 00277 { 0, 1 }, 00278 { 2, 3 }, 00279 { 3, 0 }, 00280 { 1, 2 }, 00281 { 4, 5 }, 00282 { 6, 7 }, 00283 { 7, 4 }, 00284 { 5, 6 }, 00285 { 0, 4 }, 00286 { 1, 5 }, 00287 { 2, 6 }, 00288 { 3, 7 }, 00289 }; 00290 00291 switch (et) 00292 { 00293 case ET_SEGM: return segm_edges; 00294 case ET_TRIG: return trig_edges; 00295 case ET_QUAD: return quad_edges; 00296 case ET_TET: return tet_edges; 00297 case ET_PYRAMID: return pyramid_edges; 00298 case ET_PRISM: return prism_edges; 00299 case ET_HEX: return hex_edges; 00300 default: 00301 break; 00302 } 00303 cerr << "Ng_GetEdges, illegal element type " << et << endl; 00304 return 0; 00305 } 00306 00308 static const FACE * GetFaces (ELEMENT_TYPE et) 00309 { 00310 static int tet_faces[4][4] = 00311 { { 3, 1, 2, -1 }, 00312 { 3, 2, 0, -1 }, 00313 { 3, 0, 1, -1 }, 00314 { 0, 2, 1, -1 } }; // all faces point into interior! 00315 00316 static int prism_faces[5][4] = 00317 { 00318 { 0, 2, 1, -1 }, 00319 { 3, 4, 5, -1 }, 00320 { 2, 0, 3, 5 }, 00321 { 0, 1, 4, 3 }, 00322 { 1, 2, 5, 4 } 00323 }; 00324 00325 static int pyramid_faces[5][4] = 00326 { 00327 { 0, 1, 4, -1 }, 00328 { 1, 2, 4, -1 }, 00329 { 2, 3, 4, -1 }, 00330 { 3, 0, 4, -1 }, 00331 { 0, 1, 2, 3 } // points into interior! 00332 }; 00333 00334 static int hex_faces[6][4] = 00335 { 00336 { 0, 3, 2, 1 }, 00337 { 4, 5, 6, 7 }, 00338 { 0, 1, 5, 4 }, 00339 { 1, 2, 6, 5 }, 00340 { 2, 3, 7, 6 }, 00341 { 3, 0, 4, 7 } 00342 }; 00343 00344 static int trig_faces[1][4] = 00345 { 00346 { 0, 1, 2, -1 }, 00347 }; 00348 00349 static int quad_faces[1][4] = 00350 { 00351 { 0, 1, 2, 3 }, 00352 }; 00353 00354 switch (et) 00355 { 00356 case ET_TET: return tet_faces; 00357 case ET_PRISM: return prism_faces; 00358 case ET_PYRAMID: return pyramid_faces; 00359 case ET_HEX: return hex_faces; 00360 00361 case ET_TRIG: return trig_faces; 00362 case ET_QUAD: return quad_faces; 00363 00364 case ET_SEGM: 00365 default: 00366 break; 00367 } 00368 00369 cerr << "Ng_GetFaces, illegal element type " << et << endl; 00370 return 0; 00371 } 00372 00374 static NORMAL * GetNormals(ELEMENT_TYPE et); 00375 00376 template <int D> 00377 static FlatVector<Vec<D> > GetNormals(ELEMENT_TYPE et); 00378 00379 00381 static int GetEdgeNr (ELEMENT_TYPE et, int v1, int v2); 00382 00384 static int GetFaceNr (ELEMENT_TYPE et, int v1, int v2, int v3); 00385 }; 00386 00387 00388 00396 class NGS_DLL_HEADER Node 00397 { 00398 NODE_TYPE nt; 00399 int nodenr; 00400 00401 public: 00403 Node () { ; } 00404 00406 Node (NODE_TYPE ant, int anodenr) 00407 : nt(ant), nodenr(anodenr) { ; } 00408 00410 Node (const Node & n2) 00411 { nt = n2.nt; nodenr = n2.nodenr; } 00412 00414 NODE_TYPE GetType () const { return nt; } 00415 00417 int GetNr() const { return nodenr; } 00418 }; 00419 00420 inline int CalcNodeId (ELEMENT_TYPE et, const Node & node) 00421 { 00422 switch (et) 00423 { 00424 case ET_TRIG: 00425 { 00426 static const int nodebase[] = { 0, 3, 6 }; 00427 return nodebase[node.GetType()] + node.GetNr(); 00428 } 00429 00430 case ET_TET: 00431 { 00432 static const int nodebase[] = { 0, 4, 10, 14 }; 00433 return nodebase[node.GetType()] + node.GetNr(); 00434 } 00435 default: 00436 throw Exception (string ("CalcNodeId not implemented for element ") + 00437 ElementTopology::GetElementName (et)); 00438 } 00439 } 00440 00441 inline Node CalcNodeFromId (ELEMENT_TYPE et, int nodeid) 00442 { 00443 switch (et) 00444 { 00445 case ET_TRIG: 00446 { 00447 static const int nodetypes[] = { 0, 0, 0, 1, 1, 1, 2 }; 00448 static const int nodenrs[] = { 0, 1, 2, 0, 1, 2, 0 }; 00449 return Node (NODE_TYPE(nodetypes[nodeid]), nodenrs[nodeid]); 00450 } 00451 case ET_TET: 00452 { 00453 static const int nodetypes[] = { 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3 }; 00454 static const int nodenrs[] = { 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 0 }; 00455 return Node (NODE_TYPE(nodetypes[nodeid]), nodenrs[nodeid]); 00456 } 00457 00458 default: 00459 throw Exception (string ("CalcNodeFromId not implemented for element ") + 00460 ElementTopology::GetElementName (et)); 00461 } 00462 } 00463 00464 ostream & operator<< (ostream & ost, const Node & node); 00465 00471 class NodeSet 00472 { 00473 int set; 00474 public: 00475 00476 NodeSet (NODE_TYPE nt1) 00477 { 00478 set = 1 << nt1; 00479 } 00480 00481 NodeSet (NODE_TYPE nt1, NODE_TYPE nt2) 00482 { 00483 set = (1 << nt1) + (1 << nt2); 00484 } 00485 00486 NodeSet (NODE_TYPE nt1, NODE_TYPE nt2, NODE_TYPE nt3) 00487 { 00488 set = (1 << nt1) + (1 << nt2) + (1 << nt3); 00489 } 00490 00491 NodeSet (NODE_TYPE nt1, NODE_TYPE nt2, NODE_TYPE nt3, NODE_TYPE nt4) 00492 { 00493 set = (1 << nt1) + (1 << nt2) + (1 << nt3) + (1 << nt4); 00494 } 00495 operator int () const { return set; } 00496 }; 00497 00498 00499 00500 00501 00502 class TopologicElement 00503 { 00504 ELEMENT_TYPE et; 00505 ArrayMem<Node, 27> nodes; 00506 00507 public: 00508 void SetElementType (ELEMENT_TYPE aet) { et = aet; } 00509 void Clear() { nodes.SetSize(0); } 00510 void AddNode (const Node & node) { nodes.Append (node); } 00511 00512 00513 ELEMENT_TYPE GetType () const { return et; } 00514 int GetNNodes () const { return nodes.Size(); } 00515 const Node & GetNode (int i) const { return nodes[i]; } 00516 00517 Node GlobalNode (const Node & locnode) 00518 { return Node (locnode.GetType(), nodes[CalcNodeId (et, locnode)].GetNr() ); } 00519 }; 00520 00521 ostream & operator<< (ostream & ost, const TopologicElement & etop); 00522 00523 00524 00525 00526 template <int ET> class ET_trait { }; 00527 00528 template<> class ET_trait<ET_POINT> 00529 { 00530 public: 00531 enum { DIM = 0 }; 00532 enum { N_VERTEX = 1 }; 00533 enum { N_EDGE = 0 }; 00534 enum { N_FACE = 0 }; 00535 enum { N_FACET = 0 }; 00536 00537 static int PolDimension (int order) { return 1; } 00538 static int PolBubbleDimension (INT<1> order) { return 0; } 00539 00540 static ELEMENT_TYPE FaceType(int i) { return ET_POINT; } // dummy 00541 00542 00543 00544 static INT<2> GetEdge (int /* i */); 00545 /* 00546 { 00547 static const int edges[][2] = 00548 { { 0, 1 } }; 00549 return INT<2> (edges[0][0], edges[0][1]); 00550 } 00551 */ 00552 00553 template <typename TVN> 00554 static INT<2> GetEdgeSort (int i, const TVN & vnums); 00555 /* 00556 { 00557 INT<2> e = GetEdge (i); 00558 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 00559 return e; 00560 } 00561 */ 00562 00563 static INT<4> GetFace (int /* i */ ); 00564 /* 00565 { 00566 return INT<4> (-1, -1, -1, -1); 00567 } 00568 */ 00569 template <typename TVN> 00570 static INT<4> GetFaceSort (int /* i */ , const TVN & vnums); 00571 /* 00572 { 00573 return GetFace(0); 00574 } 00575 */ 00576 /* 00577 template <typename TVN> 00578 static int GetClassNr (const TVN & vnums) 00579 { 00580 int classnr = 0; 00581 int sort[3] = { 0, 1 }; 00582 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 1; } 00583 return classnr; 00584 } 00585 00586 template <typename TVN> 00587 static int GetFacetClassNr (int facet, const TVN & vnums) 00588 { 00589 return facet; 00590 } 00591 00592 */ 00593 }; 00594 00595 00596 00597 00598 00599 template<> class ET_trait<ET_SEGM> 00600 { 00601 public: 00602 enum { DIM = 1 }; 00603 enum { N_VERTEX = 2 }; 00604 enum { N_EDGE = 1 }; 00605 enum { N_FACE = 0 }; 00606 enum { N_FACET = 2 }; 00607 00608 static int PolDimension (int order) { return order+1; } 00609 static int PolBubbleDimension (INT<1> order) { return order[0]-1; } 00610 00611 static ELEMENT_TYPE FaceType(int i) { return ET_TRIG; } 00612 00613 static INT<2> GetEdge (int /* i */) 00614 { 00615 static const int edges[][2] = 00616 { { 0, 1 } }; 00617 return INT<2> (edges[0][0], edges[0][1]); 00618 } 00619 00620 template <typename TVN> 00621 static INT<2> GetEdgeSort (int i, const TVN & vnums) 00622 { 00623 INT<2> e = GetEdge (i); 00624 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 00625 return e; 00626 } 00627 00628 00629 static INT<4> GetFace (int /* i */ ) 00630 { 00631 return INT<4> (-1, -1, -1, -1); 00632 } 00633 00634 template <typename TVN> 00635 static INT<4> GetFaceSort (int /* i */ , const TVN & vnums) 00636 { 00637 return GetFace(0); 00638 } 00639 00640 template <typename TVN> 00641 static int GetClassNr (const TVN & vnums) 00642 { 00643 int classnr = 0; 00644 int sort[3] = { 0, 1 }; 00645 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 1; } 00646 return classnr; 00647 } 00648 00649 template <typename TVN> 00650 static int GetFacetClassNr (int facet, const TVN & vnums) 00651 { 00652 return facet; 00653 } 00654 }; 00655 00656 template<> class ET_trait<ET_TRIG> 00657 { 00658 public: 00659 enum { DIM = 2 }; 00660 enum { N_VERTEX = 3 }; 00661 enum { N_EDGE = 3 }; 00662 enum { N_FACE = 1 }; 00663 enum { N_FACET = 3 }; 00664 00665 // static int PolDimension (int order) { return (order+1)*(order+2)/2; } 00666 static int PolDimension (INT<2> order) { return (order[0]+1)*(order[0]+2)/2; } 00667 static int PolBubbleDimension (INT<2> order) { return (order[0]-1)*(order[0]-2)/2; } 00668 00669 static ELEMENT_TYPE FaceType(int i) { return ET_TRIG; } 00670 00671 template <typename Tx, typename Tlam> 00672 static void CalcLambda (const Tx & x, Tlam & lam) 00673 { lam[0] = x[0]; lam[1] = x[1], lam[2] = 1-x[0]-x[1]; } 00674 00675 00676 static INT<2> GetEdge (int i) 00677 { 00678 static const int edges[][2] = 00679 { { 2, 0 }, 00680 { 1, 2 }, 00681 { 0, 1 } }; 00682 return INT<2> (edges[i][0], edges[i][1]); 00683 } 00684 00685 template <typename TVN> 00686 static INT<2> GetEdgeSort (int i, const TVN & vnums) 00687 { 00688 INT<2> e = GetEdge (i); 00689 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 00690 return e; 00691 } 00692 00693 00694 static INT<4> GetFace (int /* i */ ) 00695 { 00696 static const int face[] = { 0, 1, 2, -1 }; 00697 00698 return INT<4> (face[0], face[1], face[2], -1); 00699 } 00700 00701 template <typename TVN> 00702 static INT<4> GetFaceSort (int /* i */ , const TVN & vnums) 00703 { 00704 INT<4> f = GetFace (0); 00705 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 00706 if(vnums[f[1]] > vnums[f[2]]) swap(f[1],f[2]); 00707 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 00708 00709 return f; 00710 } 00711 00712 00713 template <typename TVN> 00714 static int GetClassNr (const TVN & vnums) 00715 { 00716 int classnr = 0; 00717 int sort[3] = { 0, 1, 2 }; 00718 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 1; } 00719 if (vnums[sort[1]] > vnums[sort[2]]) { Swap (sort[1], sort[2]); classnr += 2; } 00720 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 2; } // tricky ! 00721 return classnr; 00722 } 00723 00724 template <typename TVN> 00725 static int GetFacetClassNr (int facet, const TVN & vnums) 00726 { 00727 int sort[3] = { 0, 1, 2 }; 00728 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); } 00729 if (vnums[sort[1]] > vnums[sort[2]]) { Swap (sort[1], sort[2]); } 00730 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); } 00731 00732 static int f2vop[] = { 1, 0, 2 }; 00733 int vop = f2vop[facet]; 00734 for (int i = 0; i < 3; i++) 00735 if (vop == sort[i]) 00736 return i; 00737 return -1; // not possible 00738 00739 // return GetClassNr (vnums) * 3 + facet; 00740 } 00741 }; 00742 00743 template<> class ET_trait<ET_QUAD> 00744 { 00745 public: 00746 enum { DIM = 2 }; 00747 enum { N_VERTEX = 4 }; 00748 enum { N_EDGE = 4 }; 00749 enum { N_FACE = 1 }; 00750 enum { N_FACET = 4 }; 00751 00752 00753 // static int PolDimension (int order) { return (order+1)*(order+1); } 00754 static int PolDimension (INT<2> order) { return (order[0]+1)*(order[1]+1); } 00755 static int PolBubbleDimension (INT<2> order) { return (order[0]-1)*(order[1]-1); } 00756 00757 00758 static ELEMENT_TYPE FaceType(int i) { return ET_QUAD; } 00759 00760 00761 00762 00763 static INT<2> GetEdge (int i) 00764 { 00765 static const int edges[][2] = 00766 { { 0, 1 }, 00767 { 2, 3 }, 00768 { 3, 0 }, 00769 { 1, 2 }}; 00770 return INT<2> (edges[i][0], edges[i][1]); 00771 } 00772 00773 template <typename TVN> 00774 static INT<2> GetEdgeSort (int i, const TVN & vnums) 00775 { 00776 INT<2> e = GetEdge (i); 00777 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 00778 return e; 00779 } 00780 00781 00782 static INT<4> GetFace (int /* i */ ) 00783 { 00784 static const int face[] = 00785 { 0, 1, 2, 3 }; 00786 00787 return INT<4> (face[0], face[1], face[2], face[3]); 00788 } 00789 00790 template <typename TVN> 00791 static INT<4> GetFaceSort (int /* i */ , const TVN & vnums) 00792 { 00793 INT<4> f = GetFace (0); 00794 00795 int fmax = 0; 00796 for (int j=1; j<4; j++) 00797 if (vnums[j] < vnums[fmax]) fmax = j; 00798 00799 int f1 = (fmax+3)%4; 00800 int f2 = (fmax+1)%4; 00801 int fop = (fmax+2)%4; 00802 00803 if(vnums[f2]<vnums[f1]) swap(f1,f2); // fmax > f1 > f2 00804 00805 f[0] = fmax; 00806 f[1] = f1; 00807 f[2] = fop; 00808 f[3] = f2; 00809 00810 return f; 00811 } 00812 00813 00814 template <typename TVN> 00815 static int GetClassNr (const TVN & vnums) 00816 { 00817 int classnr = 0; 00818 00819 int sort[4] = { 0, 1, 2, 3 }; 00820 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 1; } 00821 if (vnums[sort[2]] > vnums[sort[3]]) { Swap (sort[2], sort[3]); classnr += 2; } 00822 if (vnums[sort[0]] > vnums[sort[2]]) { Swap (sort[0], sort[2]); classnr += 4; } 00823 if (vnums[sort[1]] > vnums[sort[3]]) { Swap (sort[1], sort[3]); classnr += 8; } 00824 if (vnums[sort[1]] > vnums[sort[2]]) { Swap (sort[1], sort[2]); classnr += 16; } 00825 00826 return classnr; 00827 } 00828 00829 00830 template <typename TVN> 00831 static int GetFacetClassNr (int facet, const TVN & vnums) 00832 { 00833 return GetClassNr (vnums) * 4 + facet; 00834 } 00835 }; 00836 00837 00838 template<> class ET_trait<ET_TET> 00839 { 00840 public: 00841 enum { DIM = 3 }; 00842 enum { N_VERTEX = 4 }; 00843 enum { N_EDGE = 6 }; 00844 enum { N_FACE = 4 }; 00845 enum { N_FACET = 4 }; 00846 00847 static int PolDimension (int order) { return (order+1)*(order+2)*(order+3)/6; } 00848 static int PolBubbleDimension (INT<3> p) { return (p[0] <= 3) ? 0 : (p[0]-1)*(p[0]-2)*(p[0]-3)/6; } 00849 00850 static ELEMENT_TYPE FaceType(int i) { return ET_TRIG; } 00851 00852 00853 static INT<2> GetEdge (int i) 00854 { 00855 static const int edges[6][2] = 00856 { { 3, 0 }, 00857 { 3, 1 }, 00858 { 3, 2 }, 00859 { 0, 1 }, 00860 { 0, 2 }, 00861 { 1, 2 }}; 00862 return INT<2> (edges[i][0], edges[i][1]); 00863 } 00864 00865 template <typename TVN> 00866 static INT<2> GetEdgeSort (int i, const TVN & vnums) 00867 { 00868 INT<2> e = GetEdge (i); 00869 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 00870 return e; 00871 } 00872 00873 00874 static INT<4> GetFace (int i ) 00875 { 00876 static const int faces[][4] = 00877 { { 3, 1, 2, -1 }, 00878 { 3, 2, 0, -1 }, 00879 { 3, 0, 1, -1 }, 00880 { 0, 2, 1, -1 } }; 00881 00882 return INT<4> (faces[i][0], faces[i][1], faces[i][2], -1); 00883 } 00884 00885 template <typename TVN> 00886 static INT<4> GetFaceSort (int i, const TVN & vnums) 00887 { 00888 INT<4> f = GetFace (i); 00889 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 00890 if(vnums[f[1]] > vnums[f[2]]) swap(f[1],f[2]); 00891 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 00892 00893 return f; 00894 } 00895 00896 00897 template <typename TVN> 00898 static int GetClassNr (const TVN & vnums) 00899 { 00900 int classnr = 0; 00901 00902 int sort[4] = { 0, 1, 2, 3 }; 00903 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 1; } 00904 if (vnums[sort[2]] > vnums[sort[3]]) { Swap (sort[2], sort[3]); classnr += 2; } 00905 if (vnums[sort[0]] > vnums[sort[2]]) { Swap (sort[0], sort[2]); classnr += 4; } 00906 if (vnums[sort[1]] > vnums[sort[3]]) { Swap (sort[1], sort[3]); classnr += 8; } 00907 if (vnums[sort[1]] > vnums[sort[2]]) { Swap (sort[1], sort[2]); classnr += 16; } 00908 00909 return classnr; 00910 } 00911 00912 template <typename TVN> 00913 static int GetFacetClassNr (int facet, const TVN & vnums) 00914 { 00915 int sort[4] = { 0, 1, 2, 3 }; 00916 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); } 00917 if (vnums[sort[2]] > vnums[sort[3]]) { Swap (sort[2], sort[3]); } 00918 if (vnums[sort[0]] > vnums[sort[2]]) { Swap (sort[0], sort[2]); } 00919 if (vnums[sort[1]] > vnums[sort[3]]) { Swap (sort[1], sort[3]); } 00920 if (vnums[sort[1]] > vnums[sort[2]]) { Swap (sort[1], sort[2]); } 00921 00922 for (int i = 0; i < 4; i++) 00923 if (facet == sort[i]) 00924 return i; 00925 return -1; // not possible 00926 00927 // return GetClassNr (vnums) * 4 + facet; 00928 } 00929 00930 00931 }; 00932 00933 template<> class ET_trait<ET_PRISM> 00934 { 00935 public: 00936 enum { DIM = 3 }; 00937 enum { N_VERTEX = 6 }; 00938 enum { N_EDGE = 9 }; 00939 enum { N_FACE = 5 }; 00940 enum { N_FACET = 5 }; 00941 00942 static int PolDimension (int order) { return (order+1)*(order+2)*(order+1)/2; } 00943 static int PolBubbleDimension (INT<3> p) { return (p[0] <= 2) ? 0 : (p[0]-1)*(p[0]-2)*(p[2]-1)/2; } 00944 00945 static ELEMENT_TYPE FaceType(int i) { return (i < 2) ? ET_TRIG : ET_QUAD; } 00946 00947 00948 static INT<2> GetEdge (int i) 00949 { 00950 static const int edges[][2] = 00951 { { 2, 0 }, 00952 { 0, 1 }, 00953 { 2, 1 }, 00954 { 5, 3 }, 00955 { 3, 4 }, 00956 { 5, 4 }, 00957 { 2, 5 }, 00958 { 0, 3 }, 00959 { 1, 4 }}; 00960 return INT<2> (edges[i][0], edges[i][1]); 00961 } 00962 00963 template <typename TVN> 00964 static INT<2> GetEdgeSort (int i, const TVN & vnums) 00965 { 00966 INT<2> e = GetEdge (i); 00967 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 00968 return e; 00969 } 00970 00971 00972 static INT<4> GetFace (int i ) 00973 { 00974 static const int faces[][4] = 00975 { 00976 { 0, 2, 1, -1 }, 00977 { 3, 4, 5, -1 }, 00978 { 2, 0, 3, 5 }, 00979 { 0, 1, 4, 3 }, 00980 { 1, 2, 5, 4 } 00981 }; 00982 00983 return INT<4> (faces[i][0], faces[i][1], faces[i][2], faces[i][3]); 00984 } 00985 00986 template <typename TVN> 00987 static INT<4> GetFaceSort (int i, const TVN & vnums) 00988 { 00989 INT<4> f = GetFace (i); 00990 if (i < 2) 00991 { 00992 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 00993 if(vnums[f[1]] > vnums[f[2]]) swap(f[1],f[2]); 00994 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 00995 return f; 00996 } 00997 else 00998 { 00999 int fmax = 0; 01000 for (int j=1; j<4; j++) 01001 if (vnums[f[j]] < vnums[f[fmax]]) fmax = j; 01002 01003 int f1 = (fmax+3)%4; 01004 int f2 = (fmax+1)%4; 01005 int fop = (fmax+2)%4; 01006 01007 if(vnums[f[f2]]<vnums[f[f1]]) swap(f1,f2); // fmax > f1 > f2 01008 01009 return INT<4> (f[fmax], f[f1], f[fop], f[f2]); 01010 } 01011 } 01012 01013 template <typename TVN> 01014 static int GetClassNr (const TVN & vnums) 01015 { 01016 int classnr = 0; 01017 int sort[3] = { 0, 1, 2 }; 01018 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 1; } 01019 if (vnums[sort[1]] > vnums[sort[2]]) { Swap (sort[1], sort[2]); classnr += 2; } 01020 if (vnums[sort[0]] > vnums[sort[1]]) { Swap (sort[0], sort[1]); classnr += 2; } // tricky ! 01021 return classnr; 01022 } 01023 01024 template <typename TVN> 01025 static int GetFacetClassNr (int facet, const TVN & vnums) 01026 { 01027 return 0; 01028 } 01029 01030 }; 01031 01032 template<> class ET_trait<ET_PYRAMID> 01033 { 01034 public: 01035 enum { DIM = 3 }; 01036 enum { N_VERTEX = 5 }; 01037 enum { N_EDGE = 8 }; 01038 enum { N_FACE = 5 }; 01039 enum { N_FACET = 5 }; 01040 01041 static ELEMENT_TYPE FaceType(int i) { return (i < 4) ? ET_TRIG : ET_QUAD; } 01042 01043 static int PolDimension (int order) { return (order+2)*(order+1)*(2*order+3) / 6; } 01044 static int PolBubbleDimension (INT<3> p) { return (p[0] <= 2) ? 0 : (p[0]-1)*(p[0]-2)*(2*p[0]-3)/6; } 01045 01046 static INT<2> GetEdge (int i) 01047 { 01048 static const int edges[][2] = 01049 { { 0, 1 }, 01050 { 1, 2 }, 01051 { 0, 3 }, 01052 { 3, 2 }, 01053 { 0, 4 }, 01054 { 1, 4 }, 01055 { 2, 4 }, 01056 { 3, 4 }}; 01057 01058 return INT<2> (edges[i][0], edges[i][1]); 01059 } 01060 01061 template <typename TVN> 01062 static INT<2> GetEdgeSort (int i, const TVN & vnums) 01063 { 01064 INT<2> e = GetEdge (i); 01065 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 01066 return e; 01067 } 01068 01069 01070 01071 static INT<4> GetFace (int i ) 01072 { 01073 static const int faces[][4] = 01074 { 01075 { 0, 1, 4, -1 }, 01076 { 1, 2, 4, -1 }, 01077 { 2, 3, 4, -1 }, 01078 { 3, 0, 4, -1 }, 01079 { 0, 3, 2, 1 } 01080 }; 01081 01082 return INT<4> (faces[i][0], faces[i][1], faces[i][2], faces[i][3]); 01083 } 01084 01085 template <typename TVN> 01086 static INT<4> GetFaceSort (int i, const TVN & vnums) 01087 { 01088 INT<4> f = GetFace (i); 01089 if (i < 4) 01090 { 01091 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 01092 if(vnums[f[1]] > vnums[f[2]]) swap(f[1],f[2]); 01093 if(vnums[f[0]] > vnums[f[1]]) swap(f[0],f[1]); 01094 return f; 01095 } 01096 else 01097 { 01098 int fmax = 0; 01099 for (int j=1; j<4; j++) 01100 if (vnums[f[j]] < vnums[f[fmax]]) fmax = j; 01101 01102 int f1 = (fmax+3)%4; 01103 int f2 = (fmax+1)%4; 01104 int fop = (fmax+2)%4; 01105 01106 if(vnums[f[f2]]<vnums[f[f1]]) swap(f1,f2); // fmax > f1 > f2 01107 01108 return INT<4> (f[fmax], f[f1], f[fop], f[f2]); 01109 } 01110 } 01111 01112 template <typename TVN> 01113 static int GetClassNr (const TVN & vnums) 01114 { 01115 return 0; 01116 } 01117 01118 template <typename TVN> 01119 static int GetFacetClassNr (int facet, const TVN & vnums) 01120 { 01121 return 0; 01122 } 01123 01124 }; 01125 01126 01127 01128 template<> class ET_trait<ET_HEX> 01129 { 01130 public: 01131 enum { DIM = 3 }; 01132 enum { N_VERTEX = 8 }; 01133 enum { N_EDGE = 12 }; 01134 enum { N_FACE = 6 }; 01135 enum { N_FACET = 6 }; 01136 01137 static ELEMENT_TYPE FaceType(int i) { return ET_QUAD; } 01138 01139 static inline int PolDimension (int order) { return (order+1)*(order+2)*(order+1)/2; } 01140 static inline int PolBubbleDimension (INT<3> p) { return (p[0] <= 1) ? 0 : (p[0]-1)*(p[1]-1)*(p[2]-1); } 01141 01142 static INT<2> GetEdge (int i) 01143 { 01144 static const int edges[][2] = 01145 { 01146 { 0, 1 }, 01147 { 2, 3 }, 01148 { 3, 0 }, 01149 { 1, 2 }, 01150 { 4, 5 }, 01151 { 6, 7 }, 01152 { 7, 4 }, 01153 { 5, 6 }, 01154 { 0, 4 }, 01155 { 1, 5 }, 01156 { 2, 6 }, 01157 { 3, 7 }, 01158 }; 01159 return INT<2> (edges[i][0], edges[i][1]); 01160 } 01161 01162 template <typename TVN> 01163 static INT<2> GetEdgeSort (int i, const TVN & vnums) 01164 { 01165 INT<2> e = GetEdge (i); 01166 if (vnums[e[0]] > vnums[e[1]]) swap (e[0], e[1]); 01167 return e; 01168 } 01169 01170 01171 static INT<4> GetFace (int i ) 01172 { 01173 static const int faces[][4] = 01174 { 01175 { 0, 3, 2, 1 }, 01176 { 4, 5, 6, 7 }, 01177 { 0, 1, 5, 4 }, 01178 { 1, 2, 6, 5 }, 01179 { 2, 3, 7, 6 }, 01180 { 3, 0, 4, 7 } 01181 }; 01182 01183 return INT<4> (faces[i][0], faces[i][1], faces[i][2], faces[i][3]); 01184 } 01185 01186 01187 template <typename TVN> 01188 static INT<4> GetFaceSort (int i, const TVN & vnums) 01189 { 01190 INT<4> f = GetFace (i); 01191 01192 int fmax = 0; 01193 for (int j=1; j<4; j++) 01194 if (vnums[f[j]] < vnums[f[fmax]]) fmax = j; 01195 01196 int f1 = (fmax+3)%4; 01197 int f2 = (fmax+1)%4; 01198 int fop = (fmax+2)%4; 01199 01200 if(vnums[f[f2]]<vnums[f[f1]]) swap(f1,f2); // fmax > f1 > f2 01201 01202 return INT<4> (f[fmax], f[f1], f[fop], f[f2]); 01203 } 01204 01205 01206 template <typename TVN> 01207 static int GetClassNr (const TVN & vnums) 01208 { 01209 return 0; 01210 } 01211 01212 template <typename TVN> 01213 static int GetFacetClassNr (int facet, const TVN & vnums) 01214 { 01215 return 0; 01216 } 01217 01218 }; 01219 01220 01221 inline int PolBubbleDimension (ELEMENT_TYPE ET, const INT<1> & order) 01222 { 01223 switch (ET) 01224 { 01225 case ET_SEGM: return ET_trait<ET_SEGM>::PolBubbleDimension (order); 01226 default: return 0; 01227 } 01228 } 01229 01230 inline int PolBubbleDimension (ELEMENT_TYPE ET, const INT<2> & order) 01231 { 01232 switch (ET) 01233 { 01234 case ET_TRIG: return ET_trait<ET_TRIG>::PolBubbleDimension (order); 01235 case ET_QUAD: return ET_trait<ET_QUAD>::PolBubbleDimension (order); 01236 default: return 0; 01237 } 01238 } 01239 01240 inline int PolBubbleDimension (ELEMENT_TYPE ET, const INT<3> & order) 01241 { 01242 switch (ET) 01243 { 01244 case ET_TET: return ET_trait<ET_TET>::PolBubbleDimension (order); 01245 case ET_PRISM: return ET_trait<ET_PRISM>::PolBubbleDimension (order); 01246 case ET_PYRAMID: return ET_trait<ET_PYRAMID>::PolBubbleDimension (order); 01247 case ET_HEX: return ET_trait<ET_HEX>::PolBubbleDimension (order); 01248 default: return 0; 01249 } 01250 } 01251 01252 01253 } 01254 01255 #endif