NGSolve  4.9
fem/elementtopology.hpp
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