NGSolve  4.9
comp/bilinearform.hpp
00001 #ifndef FILE_BILINEARFORM
00002 #define FILE_BILINEARFORM
00003 
00004 /*********************************************************************/
00005 /* File:   bilinearform.hpp                                          */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   25. Mar. 2000                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngcomp
00011 {
00012 
00013   class LinearForm;
00014   class Preconditioner;
00015 
00023   class NGS_DLL_HEADER BilinearForm : public NGS_Object
00024   {
00025   protected:
00027     const FESpace & fespace;
00029     const FESpace * fespace2;
00030 
00032     bool nonassemble;
00034     bool diagonal;
00036     bool multilevel;
00038     bool galerkin;
00040     bool hermitean;
00042     bool symmetric;
00044     double eps_regularization; 
00046     double unuseddiag;
00047 
00049     BilinearForm * low_order_bilinear_form;
00050 
00052     LinearForm * linearform;
00053 
00055     Array<Preconditioner*> preconditioners; 
00056 
00058     Array<BaseMatrix*> mats;
00059 
00061     Array<BilinearFormIntegrator*> parts;
00062 
00064     Array<bool> parts_deletable;
00065 
00066     /*
00067     Array<BilinearFormIntegrator*> independent_parts;
00068     Array<bool> independent_parts_deletable;
00069     Array< Vec<2,int> > independent_meshindex;
00070     */
00071 
00073     bool timing;
00075     bool print;
00077     bool printelmat;
00079     bool elmat_ev;
00081     bool eliminate_internal;
00083     bool keep_internal;
00085     bool store_inner; 
00086     
00088     bool precompute;
00090     Array<void*> precomputed_data;
00091 
00092 
00093   public:
00095     BilinearForm (const FESpace & afespace,
00096                   const string & aname, 
00097                   const Flags & flags);
00098 
00100     BilinearForm (const FESpace & afespace, 
00101                   const FESpace & afespace2, 
00102                   const string & aname,
00103                   const Flags & flags);
00105     virtual ~BilinearForm ();
00106   
00107 
00109     void AddIntegrator (BilinearFormIntegrator * bfi, bool deletable = true);
00110 
00111     /*
00112     void AddIndependentIntegrator (BilinearFormIntegrator * bfi,
00113                                    const int master_surface,
00114                                    const int slave,
00115                                    const bool deletable = true);
00116     */
00117   
00119     int NumIntegrators () const 
00120     {
00121       return parts.Size(); 
00122     }
00123 
00124     /*
00125     int NumIndependentIntegrators(void) const {return independent_parts.Size();}
00126     */
00127 
00129     BilinearFormIntegrator * GetIntegrator (int i) const { return parts[i]; }
00130 
00131     // const BilinearFormIntegrator * GetIntegrator (int i) const { return parts[i]; }
00132 
00133 
00134     /*
00135     BilinearFormIntegrator * GetIndependentIntegrator (const int i) 
00136     { return independent_parts[i]; }
00137 
00138     const BilinearFormIntegrator * GetIndependentIntegrator (const int i) const 
00139     { return independent_parts[i]; }
00140 
00141     int GetIndependentMasterSurfaceIndex(const int i) const
00142     {
00143       return independent_meshindex[i](0);
00144     }
00145     int GetIndependentSlaveIndex(const int i) const
00146     {
00147       return independent_meshindex[i](1);
00148     }
00149     */
00150 
00151 
00153     void SetLinearForm (LinearForm * alf) { linearform = alf; }
00154     
00156     void SetPreconditioner (Preconditioner * pre) 
00157     { preconditioners.Append (pre); }
00158 
00160     virtual MatrixGraph * GetGraph (int level, bool symmetric);
00161 
00163     void Assemble (LocalHeap & lh);
00164 
00167     void ReAssemble (LocalHeap & lh, bool reallocate = 0);
00168 
00171     virtual void AssembleLinearization (const BaseVector & lin,
00172                                         LocalHeap & lh, 
00173                                         bool reallocate = 0) = 0;
00174 
00176     void ApplyMatrix (const BaseVector & x,
00177                       BaseVector & y) const
00178     {
00179       y = 0;
00180       AddMatrix (1, x, y);
00181     }
00182 
00184     virtual void AddMatrix (double val, const BaseVector & x,
00185                             BaseVector & y) const = 0;
00186   
00188     virtual void AddMatrix (Complex val, const BaseVector & x,
00189                             BaseVector & y) const = 0;
00190   
00192     virtual void ApplyLinearizedMatrixAdd (double val,
00193                                            const BaseVector & lin,
00194                                            const BaseVector & x,
00195                                            BaseVector & y) const = 0;
00197     virtual void ApplyLinearizedMatrixAdd (Complex val,
00198                                            const BaseVector & lin,
00199                                            const BaseVector & x,
00200                                            BaseVector & y) const = 0;
00201 
00203     virtual double Energy (const BaseVector & x) const = 0;
00204 
00206     BaseMatrix & GetMatrix () const { return *mats.Last(); }
00207 
00208     // operator const BaseMatrix& () const { return GetMatrix(); }
00209 
00211     BaseMatrix & GetMatrix (int level) const { return *mats[level]; }
00212 
00213 
00214     // BaseMatrix & GetMatrix ()  { return *mats.Last(); }
00215     // BaseMatrix & GetMatrix (int level)  { return *mats[level]; }
00216 
00219     virtual BaseMatrix & GetHarmonicExtension () const = 0;
00220 
00224     virtual BaseMatrix & GetHarmonicExtensionTrans () const = 0;
00225 
00227     virtual BaseMatrix & GetInnerSolve () const = 0;
00228 
00230     virtual BaseMatrix & GetInnerMatrix () const = 0;
00231 
00233     bool HasLowOrderBilinearForm () const { return low_order_bilinear_form != NULL; }
00234 
00236     BilinearForm & GetLowOrderBilinearForm() const
00237     {
00238       return *low_order_bilinear_form;
00239     }
00240 
00241 
00243     bool UsesEliminateInternal () const { return eliminate_internal; }
00244 
00246     bool UsesKeepInternal () const { return keep_internal; }
00247 
00249     bool UsesStoreInner () const { return store_inner; }
00250 
00251 
00253     const FESpace & GetFESpace() const { return fespace; }
00254 
00255 
00257     bool MixedSpaces () const { return fespace2 != NULL; }
00258 
00260     const FESpace & GetFESpace2() const { return *fespace2; }
00261 
00263     int GetNLevels() const { return mats.Size(); }
00264 
00266     bool IsSymmetric(void) const { return symmetric; }
00267 
00269     void SetNonAssemble (bool na = true) { nonassemble = na; }
00270 
00272     void SetGalerkin (bool agalerkin = true) { galerkin = agalerkin; }
00273 
00275     void SetDiagonal (bool adiagonal = true) { diagonal = adiagonal; }
00276 
00278     void SetSymmetric (bool asymmetric = true) { symmetric = asymmetric; }
00279 
00281     void SetHermitean (bool ahermitean = true) { hermitean = ahermitean; }
00282 
00284     void SetMultiLevel (bool amultilevel = 1) { multilevel = amultilevel; }
00285 
00287     void SetTiming (bool at) { timing = at; }
00288 
00289     void SetEliminateInternal (bool eliminate) 
00290     { eliminate_internal = eliminate; }
00291 
00292     void SetKeepInternal (bool keep)
00293     { keep_internal = keep; }
00294 
00295     void SetStoreInner (bool storei) 
00296     { store_inner = storei; }
00297 
00298     void SetPrint (bool ap);
00299     void SetPrintElmat (bool ap);
00300     void SetElmatEigenValues (bool ee);
00301 
00303     void GalerkinProjection ();
00304 
00306     virtual void ComputeInternal (BaseVector & u, const BaseVector & f, LocalHeap & lh) const = 0;
00307   
00309     void SetEpsRegularization(double val) { eps_regularization = val; }
00310 
00312     void SetUnusedDiag (double val) { unuseddiag = val; }
00313 
00315     bool UseGalerkin () const { return galerkin; }
00316 
00318     virtual string GetClassName () const
00319     {
00320       return "BilinearForm";
00321     }
00322 
00324     virtual void PrintReport (ostream & ost);
00325 
00327     virtual void MemoryUsage (Array<MemoryUsageStruct*> & mu) const;
00328 
00330     virtual BaseVector * CreateVector() const = 0;
00331 
00333     virtual void CleanUpLevel() { ; }
00334 
00335   private:
00337     virtual void DoAssemble (LocalHeap & lh) = 0;
00338 
00340     virtual void AllocateMatrix () = 0;
00341   };
00342 
00343 
00344 
00345 
00346 
00347 
00351   template <class SCAL>
00352   class NGS_DLL_HEADER S_BilinearForm : public BilinearForm
00353   {
00354   protected:
00355 
00356     ElementByElementMatrix<SCAL> * harmonicext;
00357     BaseMatrix * harmonicexttrans;
00358     ElementByElementMatrix<SCAL> * innersolve;
00359     ElementByElementMatrix<SCAL> * innermatrix;
00360 
00361         
00362   public:
00364     S_BilinearForm (const FESpace & afespace, const string & aname,
00365                     const Flags & flags)
00366       : BilinearForm (afespace, aname, flags) 
00367     { 
00368       harmonicext = NULL;
00369       harmonicexttrans = NULL;
00370       innersolve = NULL;
00371       innermatrix = NULL;
00372     }
00373 
00375     S_BilinearForm (const FESpace & afespace, 
00376                     const FESpace & afespace2,
00377                     const string & aname, const Flags & flags)
00378       : BilinearForm (afespace, afespace2, aname, flags) 
00379     {
00380       harmonicext = NULL;
00381       harmonicexttrans = NULL;
00382       innersolve = NULL;
00383       innermatrix = NULL;
00384     }
00385 
00386 
00387     ~S_BilinearForm();
00388 
00390     void AddMatrix1 (SCAL val, const BaseVector & x,
00391                      BaseVector & y) const;
00392 
00393     virtual void AddMatrix (double val, const BaseVector & x,
00394                             BaseVector & y) const
00395     {
00396       AddMatrix1 (val, x, y);
00397     }
00398 
00399 
00400     virtual void AddMatrix (Complex val, const BaseVector & x,
00401                             BaseVector & y) const
00402     {
00403       AddMatrix1 (ConvertTo<SCAL> (val), x, y);
00404     }
00405 
00406 
00407     virtual void LapackEigenSystem(FlatMatrix<SCAL> & elmat, LocalHeap & lh) const 
00408     {
00409       ;
00410     }
00411   
00412     void ApplyLinearizedMatrixAdd1 (SCAL val,
00413                                     const BaseVector & lin,
00414                                     const BaseVector & x,
00415                                     BaseVector & y) const;
00416   
00417     virtual void ApplyLinearizedMatrixAdd (double val,
00418                                            const BaseVector & lin,
00419                                            const BaseVector & x,
00420                                            BaseVector & y) const
00421     {
00422       ApplyLinearizedMatrixAdd1 (val, lin, x, y);
00423     }
00424   
00425     virtual void ApplyLinearizedMatrixAdd (Complex val,
00426                                            const BaseVector & lin,
00427                                            const BaseVector & x,
00428                                            BaseVector & y) const
00429     {
00430       ApplyLinearizedMatrixAdd1 (ConvertTo<SCAL> (val), lin, x, y);
00431     }
00432   
00433 
00434     virtual double Energy (const BaseVector & x) const;
00435 
00436     virtual void ComputeInternal (BaseVector & u, const BaseVector & f, LocalHeap & lh) const;
00437 
00439     virtual void DoAssemble (LocalHeap & lh);
00441     // virtual void DoAssembleIndependent (BitArray & useddof, LocalHeap & lh);
00443     virtual void AssembleLinearization (const BaseVector & lin,
00444                                         LocalHeap & lh, 
00445                                         bool reallocate = 0);
00447     virtual void AddElementMatrix (const Array<int> & dnums1,
00448                                    const Array<int> & dnums2,
00449                                    const FlatMatrix<SCAL> & elmat,
00450                                    bool inner_element, int elnr,
00451                                    LocalHeap & lh) = 0;
00452 
00453     virtual void ApplyElementMatrix(const BaseVector & x,
00454                                     BaseVector & y,
00455                                     const SCAL & val,
00456                                     const Array<int> & dnums,
00457                                     const ElementTransformation & eltrans,
00458                                     const int elnum,
00459                                     const int type,
00460                                     int & cnt,
00461                                     LocalHeap & lh,
00462                                     const FiniteElement * fel,
00463                                     const SpecialElement * sel = NULL) const
00464     { cerr << "ApplyElementMatrix called for baseclass" << endl;}
00465 
00466     virtual void AddDiagElementMatrix (const Array<int> & dnums1,
00467                                        const FlatVector<SCAL> & diag,
00468                                        bool inner_element, int elnr,
00469                                        LocalHeap & lh);
00470 
00471 
00472     BaseMatrix & GetHarmonicExtension () const 
00473     { 
00474       return *harmonicext; 
00475     }
00477     BaseMatrix & GetHarmonicExtensionTrans () const
00478     { 
00479       return *harmonicexttrans; 
00480     }
00482     BaseMatrix & GetInnerSolve () const
00483     { 
00484       return *innersolve; 
00485     }
00487     BaseMatrix & GetInnerMatrix () const
00488     { 
00489       return *innermatrix; 
00490     }
00491 
00492 
00493   };
00494 
00495 
00496 
00497   template <class TM, class TV = typename mat_traits<TM>::TV_COL>
00498   class NGS_DLL_HEADER T_BilinearForm : public S_BilinearForm<typename mat_traits<TM>::TSCAL>
00499   {
00500   public:
00501     typedef typename mat_traits<TM>::TSCAL TSCAL;
00502     typedef TV TV_COL;
00503     typedef SparseMatrix<TM,TV,TV> TMATRIX;
00504     
00505   protected:
00506 
00507   public:
00509     T_BilinearForm (const FESpace & afespace, const string & aname, const Flags & flags);
00511     T_BilinearForm (const FESpace & afespace, 
00512                     const FESpace & afespace2,
00513                     const string & aname,
00514                     const Flags & flags);
00516     virtual ~T_BilinearForm ();
00517 
00519     virtual void AllocateMatrix ();
00520 
00522     virtual BaseVector * CreateVector() const;
00523 
00525     virtual void CleanUpLevel();
00526 
00528     virtual void AddElementMatrix (const Array<int> & dnums1,
00529                                    const Array<int> & dnums2,
00530                                    const FlatMatrix<TSCAL> & elmat,
00531                                    bool inner_element, int elnr,
00532                                    LocalHeap & lh);
00533 
00534     virtual void ApplyElementMatrix(const BaseVector & x,
00535                                     BaseVector & y,
00536                                     const TSCAL & val,
00537                                     const Array<int> & dnums,
00538                                     const ElementTransformation & eltrans,
00539                                     const int elnum,
00540                                     const int type,
00541                                     int & cnt,
00542                                     LocalHeap & lh,
00543                                     const FiniteElement * fel,
00544                                     const SpecialElement * sel = NULL) const;
00545 
00546     virtual void LapackEigenSystem(FlatMatrix<TSCAL> & elmat, LocalHeap & lh) const;
00547   };
00548 
00549 
00550 
00551 
00552 
00553 
00554   template <class TM, class TV = typename mat_traits<TM>::TV_COL>
00555   class NGS_DLL_HEADER T_BilinearFormSymmetric : public S_BilinearForm<typename mat_traits<TM>::TSCAL>
00556   {
00557 
00558   public:
00559     typedef typename mat_traits<TM>::TSCAL TSCAL;
00560     typedef TV TV_COL;
00561     typedef SparseMatrixSymmetric<TM,TV> TMATRIX;
00562     
00563   protected:
00564     
00565 
00566   public:
00567     T_BilinearFormSymmetric (const FESpace & afespace, const string & aname,
00568                              const Flags & flags);
00569     virtual ~T_BilinearFormSymmetric ();
00570 
00571     virtual void AllocateMatrix ();
00572     virtual void CleanUpLevel();
00573 
00574     virtual BaseVector * CreateVector() const;
00575 
00576     virtual void AddElementMatrix (const Array<int> & dnums1,
00577                                    const Array<int> & dnums2,
00578                                    const FlatMatrix<TSCAL> & elmat,
00579                                    bool inner_element, int elnr,
00580                                    LocalHeap & lh);
00581     virtual void ApplyElementMatrix(const BaseVector & x,
00582                                     BaseVector & y,
00583                                     const TSCAL & val,
00584                                     const Array<int> & dnums,
00585                                     const ElementTransformation & eltrans,
00586                                     const int elnum,
00587                                     const int type,
00588                                     int & cnt,
00589                                     LocalHeap & lh,
00590                                     const FiniteElement * fel,
00591                                     const SpecialElement * sel = NULL) const;
00592 
00593     virtual void LapackEigenSystem(FlatMatrix<TSCAL> & elmat, LocalHeap & lh) const;
00594   };
00595 
00596 
00597 
00598 
00599 
00600 
00601 
00602   template <class TM>
00603   class NGS_DLL_HEADER T_BilinearFormDiagonal : public S_BilinearForm<typename mat_traits<TM>::TSCAL>
00604   {
00605 
00606   public:
00607     typedef typename mat_traits<TM>::TSCAL TSCAL;
00608     typedef typename mat_traits<TM>::TV_COL TV_COL;
00609     typedef SparseMatrixSymmetric<TM> TMATRIX;
00610 
00611   protected:
00612 
00613   public:
00614     T_BilinearFormDiagonal (const FESpace & afespace, const string & aname,
00615                             const Flags & flags);
00616     virtual ~T_BilinearFormDiagonal ();
00617 
00618     virtual void AllocateMatrix ();
00619     virtual BaseVector * CreateVector() const;
00620 
00621     virtual void AddElementMatrix (const Array<int> & dnums1,
00622                                    const Array<int> & dnums2,
00623                                    const FlatMatrix<TSCAL> & elmat,
00624                                    bool inner_element, int elnr,
00625                                    LocalHeap & lh);
00626 
00627     virtual void AddDiagElementMatrix (const Array<int> & dnums1,
00628                                        const FlatVector<TSCAL> & diag,
00629                                        bool inner_element, int elnr,
00630                                        LocalHeap & lh);
00631     virtual void ApplyElementMatrix(const BaseVector & x,
00632                                     BaseVector & y,
00633                                     const TSCAL & val,
00634                                     const Array<int> & dnums,
00635                                     const ElementTransformation & eltrans,
00636                                     const int elnum,
00637                                     const int type,
00638                                     int & cnt,
00639                                     LocalHeap & lh,
00640                                     const FiniteElement * fel,
00641                                     const SpecialElement * sel = NULL) const;
00642   };
00643 
00644 
00645 
00646 
00652   extern NGS_DLL_HEADER BilinearForm * CreateBilinearForm (const FESpace * space,
00653                                                            const string & name,
00654                                                            const Flags & flags);
00655 
00656 
00657 
00663   class NGS_DLL_HEADER BilinearFormApplication : public BaseMatrix
00664   {
00665   protected:
00667     const BilinearForm * bf;
00668 
00669   public:
00671     BilinearFormApplication (const BilinearForm * abf);
00673     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00675     virtual void MultAdd (double val, const BaseVector & v, BaseVector & prod) const;
00677     virtual void MultAdd (Complex val, const BaseVector & v, BaseVector & prod) const;
00679     virtual BaseVector * CreateVector () const;
00681     virtual int VHeight() const
00682     {
00683       return bf->GetFESpace().GetNDof(); 
00684     }
00686     virtual int VWidth() const
00687     {
00688       return bf->GetFESpace().GetNDof(); 
00689     }
00690   };
00691 
00696   class NGS_DLL_HEADER LinearizedBilinearFormApplication : public BilinearFormApplication 
00697   {
00698   protected:
00699     const BaseVector * veclin;
00700 
00701   public:
00703     LinearizedBilinearFormApplication (const BilinearForm * abf,
00704                                        const BaseVector * aveclin);
00705 
00707     virtual void Mult (const BaseVector & v, BaseVector & prod) const;
00709     virtual void MultAdd (double val, const BaseVector & v, BaseVector & prod) const;
00711     virtual void MultAdd (Complex val, const BaseVector & v, BaseVector & prod) const;
00712   };
00713 
00714 
00718   template<class SCAL>
00719   class ElementByElement_BilinearForm : public S_BilinearForm<SCAL>
00720   {
00721   public:
00722     ElementByElement_BilinearForm (const FESpace & afespace, const string & aname,
00723                     const Flags & flags);
00724     virtual ~ElementByElement_BilinearForm ();
00725     
00726     virtual void AllocateMatrix ();
00727     virtual BaseVector * CreateVector() const;
00728     
00729     virtual void AddElementMatrix (const Array<int> & dnums1,
00730                                   const Array<int> & dnums2,
00731                                   const FlatMatrix<SCAL> & elmat,
00732                                   bool inner_element, int elnr,
00733                                   LocalHeap & lh);    
00734   };
00735   
00736   
00737   
00738 }
00739 
00740 
00741 #endif