NGSolve  4.9
fem/hdivhofe.hpp
00001 #ifndef FILE_HDIVHOFE_
00002 #define FILE_HDIVHOFE_ 
00003 
00004 /*********************************************************************/
00005 /* File:   hdivhofe.hpp                                              */
00006 /* Author: A. Becirovic, S. Zaglmayr, J. Schoeberl                   */
00007 /* Date:   15. Feb. 2003                                             */
00008 /*********************************************************************/
00009 
00010 
00011 #include "thdivfe.hpp"
00012 
00013 namespace ngfem
00014 {
00015   
00016 
00017 #undef HDIV_OLD
00018 
00019   
00020 
00021 
00023 template <int DIM>
00024 class HDivHighOrderFiniteElement : virtual public HDivFiniteElement<DIM>
00025 {
00026 protected:
00027   int vnums[8];
00028  
00029   INT<3> order_inner;
00030   INT<2> order_face[6];  // 3D only
00031   int order_edge[12];   // 2D only
00032 
00033   // bool augmented;
00034 
00035   bool discontinuous;
00036   bool ho_div_free;
00037   bool only_ho_div;
00038 
00039   typedef IntegratedLegendreMonomialExt T_ORTHOPOL;
00040 
00041 
00042 public:
00043   HDivHighOrderFiniteElement () { ; }
00044   HDivHighOrderFiniteElement (ELEMENT_TYPE aeltype);
00045 
00046   void SetVertexNumbers (FlatArray<int> & avnums);
00047   void SetOrderEdge(FlatArray<int> & oe);
00048   void SetOrderFace (FlatArray<int> & of);
00049   void SetOrderInner (int oi);
00050   void SetOrderFace (FlatArray<INT<2> > & of);
00051   void SetOrderInner (INT<3> oi); 
00052 
00053   void SetDiscontinuous (bool disc) { discontinuous = disc; };  
00054   void SetHODivFree (bool aho_div_free) { ho_div_free = aho_div_free; only_ho_div = only_ho_div && !ho_div_free;};  
00055   void SetOnlyHODiv (bool aonly_ho_div) { only_ho_div = aonly_ho_div; ho_div_free = ho_div_free && !only_ho_div;};  
00056 
00057   virtual void ComputeNDof () = 0;
00058   
00059   int EdgeOrientation (int enr) const
00060   {
00061     const EDGE * edges = ElementTopology::GetEdges (this->eltype);
00062     return (vnums[edges[enr][1]] > vnums[edges[enr][0]]) ? 1 : -1;
00063   }
00064   
00065   virtual void GetFacetDofs(int i, Array<int> & dnums) const
00066   { *testout  << " GetFacetDofs for nothing " << endl; dnums.SetSize(0);}; 
00067 };
00068 
00069 
00070 template <int D>
00071 class HDivHighOrderNormalFiniteElement : public HDivNormalFiniteElement<D>
00072 {
00073 
00074   //public:
00075   // enum { DIM = D };
00076 
00077 protected:
00078   int vnums[4];
00079   INT<2> order_inner;
00080 
00081   int ned; // number of edges in element
00082   int nv; // number of vertices in element
00083 
00084 
00085   bool augmented;
00086 
00087 public:
00089   HDivHighOrderNormalFiniteElement (ELEMENT_TYPE aeltype);
00090 
00091 
00092   void SetVertexNumbers (FlatArray<int> & avnums);
00093 
00094   void SetOrderInner (int oi);
00095   void SetOrderInner (INT<2> oi);
00096 
00097   virtual void ComputeNDof () = 0;
00098 
00099   int EdgeOrientation (int enr) const
00100   {
00101     const EDGE * edges = ElementTopology::GetEdges (this->eltype);
00102     return (vnums[edges[enr][1]] > vnums[edges[enr][0]]) ? 1 : -1;
00103   }
00104   
00105 };
00106 
00107 
00108 template <class T_ORTHOPOL = TrigExtensionMonomial>
00109 class HDivHighOrderNormalSegm : public HDivHighOrderNormalFiniteElement<1>
00110 {
00111 public:
00112 
00113   HDivHighOrderNormalSegm (int aorder);
00114   virtual void ComputeNDof();
00115 
00117   virtual void CalcShape (const IntegrationPoint & ip,
00118                            FlatVector<> shape) const;
00119 
00120 };
00121 
00122 template <class T_ORTHOPOL = TrigExtensionMonomial>
00123 class HDivHighOrderNormalTrig : public HDivHighOrderNormalFiniteElement<2>
00124 {
00125 public:
00126 
00127   HDivHighOrderNormalTrig (int aorder);
00128   virtual void ComputeNDof();
00129 
00131   virtual void CalcShape (const IntegrationPoint & ip,
00132                            FlatVector<> shape) const;
00133 
00134 };
00135 
00136 template <class T_ORTHOPOL = TrigExtensionMonomial>
00137 class HDivHighOrderNormalQuad : public HDivHighOrderNormalFiniteElement<2>
00138 {
00139 public:
00140 
00141   HDivHighOrderNormalQuad (int aorder);
00142   virtual void ComputeNDof();
00143 
00145   virtual void CalcShape (const IntegrationPoint & ip,
00146                            FlatVector<> shape) const;
00147 
00148 };
00149 
00150 
00151 
00152 
00153 template <ELEMENT_TYPE ET> class HDivHighOrderFE;
00154 
00155 
00156 template <ELEMENT_TYPE ET>
00157 class T_HDivHighOrderFiniteElement 
00158   : public HDivHighOrderFiniteElement<ET_trait<ET>::DIM>,
00159     public T_HDivFiniteElement< HDivHighOrderFE<ET>, ET>
00160     
00161 {
00162 protected:
00163   enum { DIM = ET_trait<ET>::DIM };
00164   
00165   using HDivFiniteElement<DIM>::ndof;
00166   using HDivFiniteElement<DIM>::order;
00167   using HDivFiniteElement<DIM>::eltype;
00168   // using HDivFiniteElement<DIM>::dimspace;
00169 
00170   using HDivHighOrderFiniteElement<DIM>::order_edge;
00171   using HDivHighOrderFiniteElement<DIM>::order_face;
00172   using HDivHighOrderFiniteElement<DIM>::order_inner;
00173   using HDivHighOrderFiniteElement<DIM>::ho_div_free;
00174   using HDivHighOrderFiniteElement<DIM>::only_ho_div;
00175   using HDivHighOrderFiniteElement<DIM>::discontinuous;
00176 
00177 
00178   using HDivHighOrderFiniteElement<DIM>::vnums;
00179   
00180   
00181 public:
00182   T_HDivHighOrderFiniteElement () 
00183     : HDivHighOrderFiniteElement<DIM> (ET)
00184   {
00185     for (int i = 0; i < ET_trait<ET>::N_VERTEX; i++)
00186       vnums[i] = i;
00187     // dimspace = DIM;
00188     eltype = ET;
00189   }
00190 
00191   T_HDivHighOrderFiniteElement (int aorder) 
00192   {
00193     if (DIM == 2)
00194       for (int i = 0; i < ET_trait<ET>::N_EDGE; i++)
00195         order_edge[i] = aorder;
00196     else
00197       for (int i=0; i < ET_trait<ET>::N_FACE; i++) 
00198         order_face[i] = INT<2> (aorder,aorder); 
00199     
00200     order_inner = INT<3> (aorder,aorder,aorder);
00201 
00202     for (int i = 0; i < ET_trait<ET>::N_VERTEX; i++)
00203       vnums[i] = i;
00204     // dimspace = DIM;
00205     eltype = ET;
00206   }
00207 
00208   virtual void ComputeNDof();
00209   virtual void GetInternalDofs (Array<int> & idofs) const;
00210 
00211 
00212 };
00213 
00214 
00215 
00216 
00217 template <>
00218 class HDivHighOrderFE<ET_TRIG> : public T_HDivHighOrderFiniteElement<ET_TRIG>
00219 {
00220 public:
00221   HDivHighOrderFE (int aorder);
00222 
00224   template<typename Tx, typename TFA>  
00225   void T_CalcShape (Tx hx[2], TFA & shape) const; 
00226 
00227   virtual void GetFacetDofs(int i, Array<int> & dnums) const; 
00228 };
00229 
00231 
00232 template <>
00233 class HDivHighOrderFE<ET_QUAD> : public T_HDivHighOrderFiniteElement<ET_QUAD>
00234 {
00235 public:
00236   HDivHighOrderFE (int aorder);
00237 
00238   template<typename Tx, typename TFA>  
00239   void T_CalcShape (Tx hx[2], TFA & shape) const; 
00240 
00241   virtual void GetFacetDofs(int i, Array<int> & dnums) const; 
00242 };
00243 
00244 
00245 
00246 // template <class T_ORTHOPOL = TrigExtensionMonomial>
00247 template<> 
00248 class HDivHighOrderFE<ET_TET> : public T_HDivHighOrderFiniteElement<ET_TET>
00249 {
00250    typedef TetShapesInnerLegendre T_INNERSHAPES;
00251    typedef TetShapesFaceLegendre T_FACESHAPES; 
00252 public:
00253   HDivHighOrderFE () { ; }
00254   HDivHighOrderFE (int aorder);
00255   // virtual void ComputeNDof();
00256   // virtual void GetInternalDofs (Array<int> & idofs) const;
00257   
00259   template<typename Tx, typename TFA>  
00260   void T_CalcShape (Tx hx[], TFA & shape) const; 
00261 
00262   virtual void GetFacetDofs(int i, Array<int> & dnums) const; 
00263 };
00264 
00265 
00266 // template <class T_ORTHOPOL = TrigExtensionMonomial>
00267 template<>
00268 class HDivHighOrderFE<ET_PRISM> : public T_HDivHighOrderFiniteElement<ET_PRISM>
00269 {
00270   typedef TrigShapesInnerLegendre T_TRIGFACESHAPES;
00271 public:
00272   HDivHighOrderFE () { ; }
00273   HDivHighOrderFE (int aorder);
00274   // virtual void ComputeNDof();
00275   // virtual void GetInternalDofs (Array<int> & idofs) const;
00276 
00277   /*
00279   virtual void CalcShape (const IntegrationPoint & ip,
00280                           FlatMatrixFixWidth<3> shape) const;
00281 
00283   virtual void CalcDivShape (const IntegrationPoint & ip,
00284                              FlatVector<> shape) const;
00285   */
00286 
00288   template<typename Tx, typename TFA>  
00289   void T_CalcShape (Tx hx[], TFA & shape) const; 
00290 
00291 
00293   //void CalcNumDivShape( const IntegrationPoint & ip,
00294   //                    FlatVector<> divshape) const;
00295 
00296   virtual void GetFacetDofs(int i, Array<int> & dnums) const; 
00297 };
00298 
00299 template<> 
00300 class HDivHighOrderFE<ET_HEX> : public HDivHighOrderFiniteElement<3>
00301 {
00302 public:
00303 
00304   HDivHighOrderFE (int aorder);
00305   virtual void ComputeNDof();
00306   virtual void GetInternalDofs (Array<int> & idofs) const;
00307   
00308 
00310   virtual void CalcShape (const IntegrationPoint & ip,
00311                           FlatMatrixFixWidth<3> shape) const;
00312 
00314  virtual void CalcDivShape (const IntegrationPoint & ip,
00315                              FlatVector<> shape) const;
00317   //void CalcNumDivShape( const IntegrationPoint & ip,
00318   //                    FlatVector<> divshape) const;
00319  virtual void GetFacetDofs(int i, Array<int> & dnums) const; 
00320 
00321 };
00322 
00323 }
00324 
00325 
00326 
00327 #endif
00328 
00329 
00330