NGSolve  4.9
fem/diffop.hpp
00001 #ifndef FILE_DIFFOP
00002 #define FILE_DIFFOP
00003 
00004 /*********************************************************************/
00005 /* File:   diffop.hpp                                                */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   24. Nov. 2009                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngfem
00011 {
00012 
00013 
00014 
00021   template<class DOP>
00022   class NGS_DLL_HEADER DiffOp
00023   {
00024   public:
00025     // enum { DIM_ELEMENT = TDOP::DIM_ELEMENT };
00026     // enum { DIM_SPACE = TDOP::DIM_SPACE };
00027   
00035     template <typename FEL, typename MIP, typename MAT>
00036     static void GenerateMatrix (const FEL & fel, const MIP & mip,
00037                                 MAT & mat, LocalHeap & lh)
00038     {
00039       ;
00040     }
00041 
00046     template <typename FEL, typename MIP, class TVX, class TVY>
00047     static void Apply (const FEL & fel, const MIP & mip,
00048                        const TVX & x, TVY & y,
00049                        LocalHeap & lh)
00050     {
00051       typedef typename TVY::TSCAL TSCAL;
00052 
00053       HeapReset hr(lh);
00054 
00055       FlatMatrixFixHeight<DOP::DIM_DMAT, TSCAL> mat(DOP::DIM*fel.GetNDof(), lh);
00056       DOP::GenerateMatrix (fel, mip, mat, lh);
00057       y = mat * x;
00058     }
00059 
00061     template <typename FEL, class MIR, class TVX, class TVY>
00062     static void ApplyIR (const FEL & fel, const MIR & mir,
00063                          const TVX & x, TVY & y,
00064                          LocalHeap & lh)
00065     {
00066       for (int i = 0; i < mir.Size(); i++)
00067         Apply (fel, mir[i], x, y.Row(i), lh);
00068     }
00069 
00070 
00072     template <typename FEL, typename MIP, class TVX, class TVY>
00073     static void ApplyTrans (const FEL & fel, const MIP & mip,
00074                             const TVX & x, TVY & y,
00075                             LocalHeap & lh) 
00076     {
00077       typedef typename TVY::TSCAL TSCAL;
00078 
00079       HeapReset hr(lh);
00080 
00081       FlatMatrixFixHeight<DOP::DIM_DMAT, TSCAL> mat(DOP::DIM*fel.GetNDof(), lh);
00082       DOP::GenerateMatrix (fel, mip, mat, lh);
00083       y = Trans (mat) * x;
00084     }
00085 
00086 
00088     template <typename FEL, typename MIR, class TVX, class TVY>
00089     static void ApplyTransIR (const FEL & fel, const MIR & mir,
00090                               const TVX & x, TVY & y,
00091                               LocalHeap & lh) 
00092     {
00093       typedef typename TVY::TSCAL TSCAL;
00094 
00095       HeapReset hr(lh);
00096       FlatVector<TSCAL> hy(y.Size(), lh);
00097 
00098       y = 0.0;
00099       for (int i = 0; i < mir.Size(); i++)
00100         {
00101           ApplyTrans (fel, mir[i], x.Row(i), hy, lh);
00102           y += hy;
00103         }
00104     }
00105 
00106 
00107 
00108 
00109 
00111     template <typename MIP, class TVX>
00112     static void Transform (const MIP & mip, TVX & x)
00113     {
00114       ;
00115     }
00116 
00118     template <typename MIP, class TVX>
00119     static void TransformT (const MIP & mip, TVX & x)
00120     {
00121       ;
00122     }
00123 
00124     /*
00126     template <typename FEL, class TVD, class TVY, int D>
00127     static void ApplyGrid (const FEL & fel, 
00128                            const IntegrationRuleTP<D> & ir,
00129                            const TVY & hv, 
00130                            TVD & dvecs, 
00131                            LocalHeap & lh)
00132     {
00133       ;
00134     }
00135 
00137     template <typename FEL, class TVD, class TVY, int D>
00138     static void ApplyTransGrid (const FEL & fel, 
00139                                 const IntegrationRuleTP<D> & ir,
00140                                 const TVD & dvecs, 
00141                                 TVY & hv, LocalHeap & lh)
00142     {
00143       ;
00144     }
00145     */
00146   };
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00159   class DifferentialOperator
00160   {
00161   public:
00163     NGS_DLL_HEADER virtual ~DifferentialOperator ();
00165     virtual int Dim() const = 0;
00167     virtual bool Boundary() const { return false; }
00169     virtual void
00170     CalcMatrix (const FiniteElement & fel,
00171                 const BaseMappedIntegrationPoint & mip,
00172                 FlatMatrix<double> mat, 
00173                 LocalHeap & lh) const 
00174     {
00175       cerr << "DifferentialOperator::CalcMatrix called for base class" << endl;
00176     }
00177 
00178     virtual void
00179     Apply (const FiniteElement & fel,
00180            const BaseMappedIntegrationPoint & mip,
00181            FlatVector<double> x, 
00182            FlatVector<double> flux,
00183            LocalHeap & lh) const
00184     {
00185       FlatMatrix<> mat(Dim(), fel.GetNDof(), lh);
00186       CalcMatrix (fel, mip, mat, lh);
00187       flux = mat * x;
00188     }
00189 
00190     virtual void
00191     Apply (const FiniteElement & fel,
00192            const BaseMappedIntegrationPoint & mip,
00193            FlatVector<Complex> x, 
00194            FlatVector<Complex> flux,
00195            LocalHeap & lh) const
00196     {
00197       FlatMatrix<> mat(Dim(), fel.GetNDof(), lh);
00198       CalcMatrix (fel, mip, mat, lh);
00199       flux = mat * x;
00200     }
00201 
00202     virtual void
00203     Apply (const FiniteElement & fel,
00204            const BaseMappedIntegrationRule & mir,
00205            FlatVector<double> x, 
00206            FlatMatrix<double> flux,
00207            LocalHeap & lh) const
00208     {
00209       for (int i = 0; i < mir.Size(); i++)
00210         Apply (fel, mir[i], x, flux.Row(i), lh);
00211     }
00212 
00213     virtual void
00214     Apply (const FiniteElement & fel,
00215            const BaseMappedIntegrationRule & mir,
00216            FlatVector<Complex> x, 
00217            FlatMatrix<Complex> flux,
00218            LocalHeap & lh) const
00219     {
00220       for (int i = 0; i < mir.Size(); i++)
00221         Apply (fel, mir[i], x, flux.Row(i), lh);
00222     }
00223 
00224 
00225 
00226 
00227     virtual void
00228     ApplyTrans (const FiniteElement & fel,
00229                 const BaseMappedIntegrationPoint & mip,
00230                 FlatVector<double> flux,
00231                 FlatVector<double> x, 
00232                 LocalHeap & lh) const 
00233     {
00234       FlatMatrix<> mat(Dim(), fel.GetNDof(), lh);
00235       CalcMatrix (fel, mip, mat, lh);
00236       flux = mat * x;
00237     }
00238 
00239 
00240     virtual void
00241     ApplyTrans (const FiniteElement & fel,
00242                 const BaseMappedIntegrationPoint & mip,
00243                 FlatVector<Complex> flux,
00244                 FlatVector<Complex> x, 
00245                 LocalHeap & lh) const 
00246     {
00247       FlatMatrix<> mat(Dim(), fel.GetNDof(), lh);
00248       CalcMatrix (fel, mip, mat, lh);
00249       flux = mat * x;
00250     }
00251 
00252 
00253     virtual void
00254     ApplyTrans (const FiniteElement & fel,
00255                 const BaseMappedIntegrationRule & mir,
00256                 FlatMatrix<double> flux,
00257                 FlatVector<double> x, 
00258                 LocalHeap & lh) const 
00259     {
00260       FlatVector<double> hx(x.Size(), lh);
00261       x = 0.0;
00262       for (int i = 0; i < mir.Size(); i++)
00263         {
00264           ApplyTrans (fel, mir[i], flux.Row(i), hx, lh);
00265           x += hx;
00266         }
00267     }
00268 
00269     virtual void
00270     ApplyTrans (const FiniteElement & fel,
00271                 const BaseMappedIntegrationRule & mir,
00272                 FlatMatrix<Complex> flux,
00273                 FlatVector<Complex> x, 
00274                 LocalHeap & lh) const 
00275     {
00276       FlatVector<Complex> hx(x.Size(), lh);
00277       x = 0.0;
00278       for (int i = 0; i < mir.Size(); i++)
00279         {
00280           ApplyTrans (fel, mir[i], flux.Row(i), hx, lh);
00281           x += hx;
00282         }
00283     }
00284   };
00285 
00286 
00287   class BlockDifferentialOperator : public DifferentialOperator
00288   {
00289     const DifferentialOperator & diffop;
00290     int dim;
00291     int comp;
00292   public:
00293     BlockDifferentialOperator (const DifferentialOperator & adiffop, 
00294                                int adim, int acomp = -1)
00295       : diffop(adiffop), dim(adim), comp(acomp) { ; }
00296 
00298     virtual int Dim() const { return dim*diffop.Dim(); }
00299     virtual bool Boundary() const { return diffop.Boundary(); }
00300 
00301 
00302 
00303     virtual void
00304     CalcMatrix (const FiniteElement & fel,
00305                 const BaseMappedIntegrationPoint & mip,
00306                 FlatMatrix<double> mat, 
00307                 LocalHeap & lh) const 
00308     {
00309       FlatMatrix<double> mat1(diffop.Dim(), fel.GetNDof(), lh);
00310       diffop.CalcMatrix (fel, mip, mat1, lh);
00311       mat = 0;
00312       
00313       if (comp == -1)
00314         for (int i = 0; i < mat1.Height(); i++)
00315           for (int j = 0; j < mat1.Width(); j++)
00316             for (int k = 0; k < dim; k++)
00317               mat(dim*i+k, dim*j+k) = mat1(i,j);
00318       else
00319         for (int i = 0; i < mat1.Height(); i++)
00320           for (int j = 0; j < mat1.Width(); j++)
00321             mat(dim*i+comp, dim*j+comp) = mat1(i,j);
00322     }
00323 
00324   };
00325 
00326 
00327 
00331   template <typename DIFFOP>
00332   class T_DifferentialOperator : public DifferentialOperator
00333   {
00334   protected:
00335     enum { DIM_SPACE   = DIFFOP::DIM_SPACE };
00336     enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
00337     enum { DIM_DMAT    = DIFFOP::DIM_DMAT };
00338     enum { DIM         = DIFFOP::DIM };
00339 
00340   public:
00341     virtual int Dim() const { return DIFFOP::DIM_DMAT; }
00342     virtual bool Boundary() const { return int(DIM_SPACE) > int(DIM_ELEMENT); }
00343 
00344     virtual void
00345     CalcMatrix (const FiniteElement & bfel,
00346                 const BaseMappedIntegrationPoint & bmip,
00347                 FlatMatrix<double> mat, 
00348                 LocalHeap & lh) const
00349     {
00350       const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip =
00351         static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip);
00352       DIFFOP::GenerateMatrix (bfel, mip, mat, lh);
00353     }
00354 
00355     virtual void
00356     Apply (const FiniteElement & bfel,
00357            const BaseMappedIntegrationPoint & bmip,
00358            FlatVector<double> x, 
00359            FlatVector<double> flux,
00360            LocalHeap & lh) const
00361     {
00362       const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip =
00363         static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip);
00364       DIFFOP::Apply (bfel, mip, x, flux, lh);
00365     }
00366 
00367 
00368     virtual void
00369     Apply (const FiniteElement & bfel,
00370            const BaseMappedIntegrationRule & bmir,
00371            FlatVector<double> x, 
00372            FlatMatrix<double> flux,
00373            LocalHeap & lh) const
00374     {
00375       const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir =
00376         static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir);
00377       DIFFOP::ApplyIR (bfel, mir, x, flux, lh);
00378     }
00379 
00380     virtual void
00381     ApplyTrans (const FiniteElement & bfel,
00382                 const BaseMappedIntegrationPoint & bmip,
00383                 FlatVector<double> flux,
00384                 FlatVector<double> x, 
00385                 LocalHeap & lh) const 
00386     {
00387       const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip =
00388         static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip);
00389       DIFFOP::ApplyTrans (bfel, mip, flux, x, lh);
00390     }    
00391   };
00392 
00393 
00394 
00395 
00396 
00397 
00398 
00399 
00400   
00401   // new design, code is still experimental ...
00402   
00403   template <typename DOP, typename F>
00404   class T_FunctionDiffOp : public DifferentialOperator
00405   {
00406 
00407     // possible conversion from vector to scalar 
00408     class V2VS 
00409     {
00410       FlatVector<> v;
00411     public:
00412       V2VS (FlatVector<> av) : v(av) { ; }
00413       
00414       template <int D>
00415       V2VS (Vec<D> av) : v(av) { ; }
00416       
00417       operator double () { return v(0); }
00418       
00419       operator FlatVector<> () { return v; }
00420       
00421       template <int D>
00422       operator Vec<D> () { return v; }
00423     };
00424     
00425     
00426     int dim;
00427     const F & func;
00428   public:
00429     
00430     T_FunctionDiffOp (const F & afunc, int adim) : func(afunc), dim(adim) { ; }
00431     
00432     virtual int Dim() const { return dim; }
00433     
00434     virtual void Apply (const FiniteElement & fel,
00435                         const BaseMappedIntegrationPoint & mip,
00436                         FlatVector<double> x, 
00437                         FlatVector<double> flux,
00438                         LocalHeap & lh) const 
00439     {
00440       Vec<DOP::DIM_DMAT> u;
00441       DOP::Apply (fel, mip, x, u, lh);
00442       flux = func(V2VS(u));
00443     }
00444   };
00445   
00446   
00447   template <typename DOP, typename F>
00448   DifferentialOperator * CreateFunctionDiffOp (const DOP & dop, 
00449                                                const F & func, int dim = 1)
00450   {
00451     return new T_FunctionDiffOp<DOP, F> (func, dim);
00452   }
00453 
00454 
00455   /* examples:
00456      
00457   double myexp (double x)  { return exp(x); }
00458   Vec<1> myexpVec (FlatVector<> x)  { return Vec<1> (exp(x(0))); }
00459 
00460   CreateFunctionDiffOp(DiffOpId<2>(), myexpVec));
00461   */
00462   
00463 }
00464 
00465 
00466 #endif