NGSolve  4.9
fem/l2hofe.hpp
00001 #ifndef FILE_L2HOFE
00002 #define FILE_L2HOFE
00003 
00004 /*********************************************************************/
00005 /* File:   l2hofe.hpp                                                */
00006 /* Author: Start                                                     */
00007 /* Date:   6. Feb. 2003                                              */
00008 /*********************************************************************/
00009 
00010 
00011 #include "tscalarfe.hpp"
00012 #include "precomp.hpp"
00013 
00014 
00015 namespace ngfem
00016 {
00017 
00022   template<int D>
00023   class L2HighOrderFiniteElement : virtual public ScalarFiniteElement<D>
00024   {
00025   protected:
00026 
00027     enum { DIM = D };
00028     int vnums[8];  
00029     INT<DIM> order_inner; 
00030 
00031     using ScalarFiniteElement<DIM>::ndof;
00032     using ScalarFiniteElement<DIM>::order;
00033     using ScalarFiniteElement<DIM>::eltype;
00034 
00035 
00036   public:
00038     template <typename TA>
00039     void SetVertexNumbers (const TA & avnums)
00040     { for (int i = 0; i < avnums.Size(); i++) vnums[i] = avnums[i]; }
00041 
00042     void SetVertexNumber (int nr, int vnum) { vnums[nr] = vnum; }
00043 
00045     // void SetOrder (int p)  { order_inner = p; }
00046 
00048     void SetOrder (INT<DIM> p) { order_inner = p; }
00049 
00051     virtual void ComputeNDof () = 0; 
00052   
00053     virtual void GetInternalDofs (Array<int> & idofs) const; 
00054 
00055     virtual void PrecomputeTrace () = 0;
00056 
00057     void CalcTraceMatrix (int facet, FlatMatrix<> & trace) const;
00058 
00059     virtual void GetTrace (int facet, FlatVector<> coefs, FlatVector<> fcoefs) const
00060     {
00061       Matrix<> trace(fcoefs.Size(), coefs.Size());
00062       CalcTraceMatrix(facet, trace);
00063       fcoefs = trace * coefs;
00064     }
00065 
00066     virtual void GetTraceTrans (int facet, FlatVector<> fcoefs, FlatVector<> coefs) const
00067     {
00068       Matrix<> trace(fcoefs.Size(), coefs.Size());
00069       CalcTraceMatrix(facet, trace);
00070       coefs = Trans (trace) * fcoefs;
00071     }
00072   };
00073 
00074 
00075 
00080   // template <ELEMENT_TYPE ET> class L2HighOrderFE;
00081 
00082 
00086   template <ELEMENT_TYPE ET>
00087   class T_L2HighOrderFiniteElement : 
00088     public L2HighOrderFiniteElement<ET_trait<ET>::DIM>,
00089     public ET_trait<ET> 
00090   {
00091   protected:
00092     enum { DIM = ET_trait<ET>::DIM };
00093 
00094     using ScalarFiniteElement<DIM>::ndof;
00095     using ScalarFiniteElement<DIM>::order;
00096     using ScalarFiniteElement<DIM>::eltype;
00097 
00098     using L2HighOrderFiniteElement<DIM>::vnums;
00099     using L2HighOrderFiniteElement<DIM>::order_inner;
00100 
00101     using ET_trait<ET>::N_VERTEX;
00102     using ET_trait<ET>::N_EDGE;
00103     using ET_trait<ET>::N_FACE;
00104     using ET_trait<ET>::FaceType;
00105     using ET_trait<ET>::GetEdgeSort;
00106     using ET_trait<ET>::GetFaceSort;
00107     using ET_trait<ET>::PolDimension;
00108 
00109 
00110   public:
00111 
00112     T_L2HighOrderFiniteElement () 
00113     { eltype = ET; }
00114 
00115     T_L2HighOrderFiniteElement (int aorder) 
00116     {
00117       for (int i = 0; i < ET_trait<ET>::N_VERTEX; i++) vnums[i] = i;
00118       eltype = ET;
00119 
00120       order = aorder;
00121       order_inner = aorder;
00122       ndof = PolDimension (aorder);
00123     }
00124 
00125     virtual void ComputeNDof();
00126   };
00127 
00128 
00129 
00130 
00131 
00132 
00133   template <ELEMENT_TYPE ET> class L2HighOrderFE_Shape;
00134 
00135   template <int DIM>
00136   class PrecomputedScalShapes
00137   {
00138   public:
00139     Matrix<> shapes;
00140     Matrix<> dshapes;
00141     
00142     PrecomputedScalShapes (int nip, int ndof)
00143       : shapes(nip, ndof), dshapes (DIM*nip, ndof)
00144     { ; }
00145 
00146   };
00147 
00148 
00149   template <ELEMENT_TYPE ET, template <ELEMENT_TYPE ET> class SHAPES = L2HighOrderFE_Shape> 
00150   class NGS_DLL_HEADER L2HighOrderFE :  public T_L2HighOrderFiniteElement<ET>,
00151                                         public T_ScalarFiniteElement2<SHAPES<ET>, ET >
00152   { 
00153   protected:
00154     using T_L2HighOrderFiniteElement<ET>::ndof;
00155     using T_L2HighOrderFiniteElement<ET>::order;
00156     using T_L2HighOrderFiniteElement<ET>::vnums;
00157 
00158     enum { DIM = ET_trait<ET>::DIM };
00159     typedef PrecomputedShapesContainer<PrecomputedScalShapes<DIM> > TPRECOMP;
00160     static TPRECOMP precomp;
00161 
00162     typedef HashTable<INT<2>, Matrix<>*> TPRECOMP_TRACE;
00163     static TPRECOMP_TRACE precomp_trace;
00164   public:
00165     L2HighOrderFE () { ; }
00166 
00167     L2HighOrderFE (int aorder)
00168       : T_L2HighOrderFiniteElement<ET> (aorder) { ; }
00169 
00170 
00171     virtual void PrecomputeTrace ();
00172     virtual void PrecomputeShapes (const IntegrationRule & ir);
00173 
00174     virtual void Evaluate (const IntegrationRule & ir, FlatVector<double> coefs, FlatVector<double> vals) const;
00175 
00176     virtual void EvaluateGrad (const IntegrationRule & ir, FlatVector<> coefs, FlatMatrixFixWidth<DIM> values) const;
00177 
00178     virtual void EvaluateGradTrans (const IntegrationRule & ir, FlatMatrixFixWidth<DIM> values, FlatVector<> coefs) const;
00179 
00180     virtual void GetTrace (int facet, FlatVector<> coefs, FlatVector<> fcoefs) const;
00181 
00182     virtual void GetTraceTrans (int facet, FlatVector<> fcoefs, FlatVector<> coefs) const;
00183   };
00184 
00185 
00186 
00187 
00188 
00189 
00190 
00191 
00192 
00193 
00194 
00198   template <>
00199   class L2HighOrderFE_Shape<ET_SEGM> : public L2HighOrderFE<ET_SEGM>
00200   {
00201   public:
00202     template<typename Tx, typename TFA>  
00203     void T_CalcShape (Tx x[], TFA & shape) const
00204     {
00205       Tx lam[2] = { x[0], 1-x[0] };
00206       INT<2> e = GetEdgeSort (0, vnums);
00207       LegendrePolynomial (order, lam[e[1]]-lam[e[0]], shape);
00208     }
00209   };
00210 
00211 
00215   template <> 
00216   class L2HighOrderFE_Shape<ET_TRIG> : public L2HighOrderFE<ET_TRIG>
00217   {
00218   public:
00219     template<typename Tx, typename TFA>  
00220     void T_CalcShape (Tx x[], TFA & shape) const
00221     {
00222       Tx lam[3] = { x[0], x[1], 1-x[0]-x[1] };
00223       INT<4> f = GetFaceSort (0, vnums);
00224       int p = order_inner[0];
00225       DubinerBasis::Eval (p, lam[f[0]], lam[f[1]], shape);
00226     }
00227   };
00228   
00229     
00233   template <> 
00234   class L2HighOrderFE_Shape<ET_QUAD> : public L2HighOrderFE<ET_QUAD>
00235   {
00236   public:
00237     template<typename Tx, typename TFA>  
00238     void T_CalcShape (Tx hx[], TFA & shape) const
00239     {
00240       Tx x = hx[0], y = hx[1];
00241       
00242       Tx sigma[4] = {(1-x)+(1-y),x+(1-y),x+y,(1-x)+y};  
00243       
00244       INT<4> f = GetFaceSort (0, vnums);  
00245       
00246       Tx xi = sigma[f[0]]-sigma[f[1]]; 
00247       Tx eta = sigma[f[0]]-sigma[f[3]]; 
00248       
00249       int n = max(order_inner[0],order_inner[1]);
00250       ArrayMem<Tx, 20> polx(n+1), poly(n+1);
00251       
00252       LegendrePolynomial (n, xi, polx);
00253       LegendrePolynomial (n, eta, poly);
00254       
00255       for (int i = 0, ii = 0; i <= order_inner[0]; i++)
00256         for (int j = 0; j <= order_inner[1]; j++)
00257           shape[ii++] = polx[i] * poly[j];
00258     }
00259   };
00260 
00261 
00265   template <> 
00266   class L2HighOrderFE_Shape<ET_TET> : public L2HighOrderFE<ET_TET>
00267   {
00268   public:
00269     template<typename Tx, typename TFA>  
00270     void T_CalcShape (Tx x[], TFA & shape) const
00271     {
00272       Tx lami[4] = { x[0], x[1], x[2], 1-x[0]-x[1]-x[2] };
00273       
00274       int sort[4];
00275       for (int i = 0; i < 4; i++) sort[i] = i;
00276     
00277       if (vnums[sort[0]] > vnums[sort[1]]) Swap (sort[0], sort[1]);
00278       if (vnums[sort[2]] > vnums[sort[3]]) Swap (sort[2], sort[3]);
00279       if (vnums[sort[0]] > vnums[sort[2]]) Swap (sort[0], sort[2]);
00280       if (vnums[sort[1]] > vnums[sort[3]]) Swap (sort[1], sort[3]);
00281       if (vnums[sort[1]] > vnums[sort[2]]) Swap (sort[1], sort[2]);
00282 
00283       Tx lamis[4];
00284       for (int i = 0; i < 4; i++)
00285         lamis[i] = lami[sort[i]];
00286 
00287       ArrayMem<Tx, 20> memx(sqr(order+1));
00288       ArrayMem<Tx, 20> memy(sqr(order+1));
00289 
00290       FlatMatrix<Tx> polsx(order+1, &memx[0]);
00291       FlatMatrix<Tx> polsy(order+1, &memy[0]);
00292       VectorMem<10, Tx> polsz(order+1);
00293     
00294       for (int i = 0; i <= order; i++)
00295         JacobiPolynomial (order, 2*lamis[0]-1, 2*i+2, 0, polsx.Row(i));
00296       for (int i = 0; i <= order; i++)
00297         ScaledJacobiPolynomial (order, lamis[1]-lamis[2]-lamis[3], 1-lamis[0], 2*i+1, 0, polsy.Row(i));
00298 
00299       ScaledLegendrePolynomial (order, lamis[2]-lamis[3], lamis[2]+lamis[3], polsz);
00300 
00301       for (int i = 0, ii = 0; i <= order; i++)
00302         for (int j = 0; j <= order-i; j++)
00303           for (int k = 0; k <= order-i-j; k++, ii++)
00304             shape[ii] = polsz(k) * polsy(k, j) * polsx(j+k, i);
00305     }
00306 
00307   };
00308 
00309 
00313   template <> 
00314   class L2HighOrderFE_Shape<ET_PRISM> : public L2HighOrderFE<ET_PRISM>
00315   {
00316   public:
00317     template<typename Tx, typename TFA>  
00318     void T_CalcShape (Tx hx[], TFA & shape) const
00319     {
00320       Tx lami[3] = { hx[0], hx[1], 1-hx[0]-hx[1] };
00321 
00322       int sort[3];
00323       for (int i = 0; i < 3; i++) sort[i] = i;
00324     
00325       if (vnums[sort[0]] > vnums[sort[1]]) Swap (sort[0], sort[1]);
00326       if (vnums[sort[1]] > vnums[sort[2]]) Swap (sort[1], sort[2]);
00327       if (vnums[sort[0]] > vnums[sort[1]]) Swap (sort[0], sort[1]);
00328 
00329       Tx lamis[3];
00330       for (int i = 0; i < 3; i++)
00331         lamis[i] = lami[sort[i]];
00332 
00333       Tx x = lamis[0];
00334       // Tx y = lamis[1];
00335       Tx z = hx[2];
00336 
00337       int p=order_inner[0];
00338       int q=order_inner[1];
00339 
00340       ArrayMem<Tx, 20> memx(sqr(p+1));
00341       FlatMatrix<Tx> polsx(p+1, &memx[0]);
00342 
00343       VectorMem<10, Tx> polsy(p+1);
00344       VectorMem<10, Tx> polsz(q+1);
00345     
00346       for (int i = 0; i <= p; i++)
00347         JacobiPolynomial (p, 2*x-1, 2*i+1, 0, polsx.Row(i));
00348 
00349       ScaledLegendrePolynomial (order, lamis[1]-lamis[2], lamis[1]+lamis[2], polsy);
00350       LegendrePolynomial (order, 2*z-1, polsz);
00351 
00352       int ii = 0;
00353       for (int k = 0; k <= q; k++)
00354         for (int i = 0; i <= p; i++)
00355           for (int j = 0; j <= p-i; j++)
00356             shape[ii++] = polsx(j,i) * polsy(j) * polsz(k);
00357     }
00358 
00359   };
00360 
00361 
00365   template <> 
00366   class L2HighOrderFE_Shape<ET_PYRAMID> : public L2HighOrderFE<ET_PYRAMID>
00367   {
00368   public:
00369     template<typename Tx, typename TFA>  
00370     void T_CalcShape (Tx x[], TFA & shape) const;
00371   };
00372 
00373 
00377   template <> 
00378   class L2HighOrderFE_Shape<ET_HEX> : public L2HighOrderFE<ET_HEX>
00379   {
00380   public:
00381     template<typename Tx, typename TFA>  
00382     void T_CalcShape (Tx x[], TFA & shape) const;
00383   };
00384 }
00385 
00386 #endif