NGSolve  4.9
fem/tscalarfe.hpp
00001 #ifndef FILE_TSCALARFE
00002 #define FILE_TSCALARFE
00003 
00004 /*********************************************************************/
00005 /* File:   tscalarfe.hpp                                             */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   25. Mar. 2000                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngfem
00011 {
00012 
00013 
00017   template <int DIM>
00018   class DShapeElement
00019   {
00020     FlatVec<DIM> data;
00021   public:
00023     DShapeElement (FlatVec<DIM> adata) : data(adata) { ; }
00024 
00026     void operator= (AutoDiff<DIM> ad) 
00027     { 
00028       for (int i = 0; i < DIM; i++) 
00029         data(i) = ad.DValue(i); 
00030     }
00031   };
00032 
00033 
00037   template <int DIM>
00038   class DShapeAssign
00039   {
00040     FlatMatrixFixWidth<DIM> dshape;
00041   public:
00043     DShapeAssign (FlatMatrixFixWidth<DIM> mat) : dshape(mat) { ; }
00044 
00046     DShapeElement<DIM> operator[] (int i) const
00047     { return DShapeElement<DIM> (dshape.Row(i)); }
00048 
00050     const DShapeAssign Addr (int i) const
00051     { return DShapeAssign (dshape.Rows(i, dshape.Height())); }
00052   };
00053 
00054 
00055 
00056   
00057 
00061   class EvaluateShapeElement
00062   {
00063     double coef;
00064     double * sum;
00065   public:
00067     EvaluateShapeElement (double acoef, double * asum)
00068       : coef(acoef), sum(asum) { ; }
00069 
00071     void operator= (double val) 
00072     {
00073       *sum += coef * val;
00074     }
00075   };
00076 
00077 
00081   class EvaluateShape
00082   {
00083     const double * coefs;
00084     double * sum;
00085   public:
00087     EvaluateShape (FlatVector<> acoefs, double * asum)
00088       : coefs(&acoefs(0)), sum(asum) { ; }
00089     
00091     EvaluateShape (const double * acoefs, double * asum)
00092       : coefs(acoefs), sum(asum) { ; }
00093 
00095     EvaluateShapeElement operator[] (int i) const
00096     { return EvaluateShapeElement (coefs[i], sum); }
00097 
00099     const EvaluateShape Addr (int i) const
00100     { return EvaluateShape (coefs+i, sum); } 
00101   };
00102 
00103 
00104 
00105   
00107   class EvaluateShapeTransElement
00108   {
00109     double & data;
00110     const double & fac;
00111   public:
00113     EvaluateShapeTransElement (double & adata, const double & afac) 
00114       : data(adata), fac(afac) { ; }
00115 
00117     void operator= (double ad) 
00118     { 
00119       data += ad * fac;
00120     }
00121   };
00122 
00123   class EvaluateShapeTrans
00124   {
00125     double * coefs;
00126     const double & fac;
00127   public:
00128     EvaluateShapeTrans (FlatVector<> acoefs, const double & afac)
00129       : coefs(&acoefs(0)), fac(afac) { ; }
00130 
00131     EvaluateShapeTrans (double * acoefs, const double & afac)
00132       : coefs(acoefs), fac(afac) { ; }
00133 
00134     EvaluateShapeTransElement operator[] (int i) const
00135     { return EvaluateShapeTransElement (coefs[i], fac); }
00136 
00137     const EvaluateShapeTrans Addr (int i) const
00138     { return EvaluateShapeTrans (coefs+i, fac); }
00139   };
00140 
00141 
00142 
00143 
00144 
00145 
00146 
00147 
00148   
00149   template <int DIM>
00150   class EvaluateDShapeElement
00151   {
00152     double data;
00153     Vec<DIM> & sum;
00154   public:
00155     EvaluateDShapeElement (double adata, Vec<DIM> & asum) : data(adata), sum(asum) { ; }
00156     void operator= (AutoDiff<DIM> ad) 
00157     { 
00158       for (int i = 0; i < DIM; i++) 
00159         sum(i) += ad.DValue(i) * data;
00160     }
00161   };
00162 
00163   template <int DIM>
00164   class EvaluateDShape
00165   {
00166     double * coefs;
00167     Vec<DIM> & sum;
00168   public:
00169     EvaluateDShape (FlatVector<> acoefs, Vec<DIM> & asum)
00170       : coefs(&acoefs(0)), sum(asum) { ; }
00171 
00172     EvaluateDShape (double * acoefs, Vec<DIM> & asum)
00173       : coefs(acoefs), sum(asum) { ; }
00174 
00175     EvaluateDShapeElement<DIM> operator[] (int i) const
00176     { return EvaluateDShapeElement<DIM> (coefs[i], sum); }
00177 
00178     const EvaluateDShape Addr (int i) const
00179     { return EvaluateDShape (coefs+i, sum); }
00180   };
00181 
00182 
00183 
00184 
00185 
00186   
00187   
00188   template <int DIM>
00189   class EvaluateDShapeTransElement
00190   {
00191     double & data;
00192     const Vec<DIM> & fac;
00193   public:
00194     EvaluateDShapeTransElement (double & adata, const Vec<DIM> & afac) : data(adata), fac(afac) { ; }
00195     void operator= (AutoDiff<DIM> ad) 
00196     { 
00197       for (int i = 0; i < DIM; i++) 
00198         data += ad.DValue(i) * fac(i);
00199     }
00200   };
00201 
00203   template <int DIM>
00204   class EvaluateDShapeTrans
00205   {
00206     double * coefs;
00207     const Vec<DIM> & fac;
00208   public:
00209     EvaluateDShapeTrans (FlatVector<> acoefs, const Vec<DIM> & afac)
00210       : coefs(&acoefs(0)), fac(afac) { ; }
00211 
00212     EvaluateDShapeTrans (double * acoefs, const Vec<DIM> & afac)
00213       : coefs(acoefs), fac(afac) { ; }
00214 
00215     EvaluateDShapeTransElement<DIM> operator[] (int i) const
00216     { return EvaluateDShapeTransElement<DIM> (coefs[i], fac); }
00217 
00218     const EvaluateDShapeTrans Addr (int i) const
00219     { return EvaluateDShapeTrans (coefs+i, fac); }
00220   };
00221 
00222 
00223 
00224 
00225 
00226 
00227 
00232   template <class FEL, ELEMENT_TYPE ET, int NDOF, int ORDER>
00233   class T_ScalarFiniteElement : public ScalarFiniteElement<ET_trait<ET>::DIM>
00234   {
00235 
00236   public:
00237     
00238   protected:
00239     enum { DIM = ET_trait<ET>::DIM };
00240     
00241     T_ScalarFiniteElement ()
00242       : ScalarFiniteElement<DIM> (ET, NDOF, ORDER) { ; }
00243 
00244     virtual ~T_ScalarFiniteElement() { ; }
00245 
00246   public:
00247 
00248     virtual void CalcShape (const IntegrationPoint & ip, 
00249                             FlatVector<> shape) const
00250     {
00251       // double pt[DIM];
00252       Vec<DIM> pt;
00253       for (int i = 0; i < DIM; i++) pt[i] = ip(i);
00254       FEL::T_CalcShape (&pt(0), shape); 
00255     }
00256 
00257     virtual double
00258     Evaluate (const IntegrationPoint & ip, FlatVector<double> x) const
00259     {
00260       // double pt[DIM];
00261       Vec<DIM> pt;
00262       for (int i = 0; i < DIM; i++) pt[i] = ip(i);
00263 
00264       double sum = 0.0;
00265       EvaluateShape eval(x, &sum); 
00266 
00267       FEL::T_CalcShape (&pt(0), eval); 
00268       return sum;
00269     }  
00270     
00271     virtual void
00272     Evaluate (const IntegrationRule & ir, FlatVector<double> coefs, FlatVector<double> vals) const
00273     {
00274       // double pt[DIM];
00275       Vec<DIM> pt;
00276       for (int i = 0; i < ir.GetNIP(); i++)
00277         {
00278           for (int j = 0; j < DIM; j++) pt[j] = ir[i](j);
00279           
00280           vals(i) = 0.0;
00281           EvaluateShape eval(coefs, &vals(i)); 
00282           FEL::T_CalcShape (&pt(0), eval); 
00283         }
00284     }
00285 
00286 
00287 
00288     static void CalcShapeStat (const IntegrationPoint & ip, 
00289                                FlatVector<> shape)
00290     {
00291       // double pt[DIM];
00292       Vec<DIM> pt;
00293       for (int i = 0; i < DIM; i++) pt[i] = ip(i);
00294       FEL::T_CalcShape (&pt(0), shape); 
00295     }
00296     
00297     virtual void CalcDShape (const IntegrationPoint & ip, 
00298                              FlatMatrixFixWidth<DIM> dshape) const
00299     {
00300       // AutoDiff<DIM> adp[DIM];
00301       Vec<DIM, AutoDiff<DIM> > adp;
00302       for (int i = 0; i < DIM; i++)
00303         adp[i] = AutoDiff<DIM> (ip(i), i);
00304       
00305       DShapeAssign<DIM> ds(dshape); 
00306       FEL::T_CalcShape (&adp(0), ds);
00307     }
00308 
00309     static void CalcDShapeStat (const IntegrationPoint & ip, 
00310                                 FlatMatrixFixWidth<DIM> dshape)
00311     {
00312       // AutoDiff<DIM> adp[DIM];
00313       Vec<DIM, AutoDiff<DIM> > adp;
00314       for (int i = 0; i < DIM; i++)
00315         adp[i] = AutoDiff<DIM> (ip(i), i);
00316       
00317       DShapeAssign<DIM> ds(dshape); 
00318       FEL::T_CalcShape (&adp(0), ds);
00319     }
00320 
00321     virtual void 
00322     CalcMappedDShape (const MappedIntegrationPoint<DIM,DIM> & sip, 
00323                       FlatMatrixFixWidth<DIM> dshape) const
00324     {
00325       // AutoDiff<DIM> adp[DIM];
00326       Vec<DIM, AutoDiff<DIM> > adp;
00327       for (int i = 0; i < DIM; i++)
00328         adp[i].Value() = sip.IP()(i);
00329       
00330       for (int i = 0; i < DIM; i++)
00331         for (int j = 0; j < DIM; j++)
00332           adp[i].DValue(j) = sip.GetJacobianInverse()(i,j);
00333       
00334       DShapeAssign<DIM> ds(dshape); 
00335       FEL::T_CalcShape (&adp(0), ds);
00336     }
00337   };
00338 
00339 
00340 
00341 
00342 
00343 
00344 
00345 
00346 
00347 
00348 
00349 
00350 
00351 
00352 
00353 
00354 
00355 
00356 
00357 
00358 
00359 
00360 
00361 
00362 
00363 
00364 
00365 
00371   template <class FEL, ELEMENT_TYPE ET>
00372   class NGS_DLL_HEADER T_ScalarFiniteElement2 : virtual public ScalarFiniteElement<ET_trait<ET>::DIM>
00373   {
00374   public:
00375     enum { DIM = ET_trait<ET>::DIM };
00376     using ScalarFiniteElement<DIM>::eltype;
00377 
00378     T_ScalarFiniteElement2 () { eltype = ET; }
00379 
00380     virtual void CalcShape (const IntegrationPoint & ip, 
00381                             FlatVector<> shape) const;
00382 
00383     virtual double Evaluate (const IntegrationPoint & ip, FlatVector<double> x) const;
00384     
00385     virtual void Evaluate (const IntegrationRule & ir, FlatVector<double> coefs, FlatVector<double> vals) const;
00386 
00387     virtual void EvaluateTrans (const IntegrationRule & ir, FlatVector<> vals, FlatVector<double> coefs) const;
00388 
00389     virtual void EvaluateGrad (const IntegrationRule & ir, FlatVector<double> coefs, FlatMatrixFixWidth<DIM> vals) const;
00390 
00391     virtual void EvaluateGradTrans (const IntegrationRule & ir, FlatMatrixFixWidth<DIM> vals, FlatVector<double> coefs) const;
00392 
00393     virtual void CalcDShape (const IntegrationPoint & ip, 
00394                              FlatMatrixFixWidth<DIM> dshape) const;
00395 
00396     virtual void 
00397     CalcMappedDShape (const MappedIntegrationPoint<DIM,DIM> & sip, 
00398                       FlatMatrixFixWidth<DIM> dshape) const;
00399 
00400 
00401 
00402   };
00403 
00404   template <typename TFA>
00405   inline void SetZero (TFA & shape, int first, int next)
00406   {
00407     for (int i = first; i < next; i++)
00408       shape[i] = 0.0;
00409   }
00410   
00411   inline void SetZero (EvaluateShape & shape, int first, int next) 
00412   {
00413     ;
00414   }
00415   
00416   inline void SetZero (EvaluateShapeTrans & shape, int first, int next)
00417   {
00418     ;
00419   }
00420 
00421 
00422 }
00423 
00424 
00425 #endif