NGSolve
4.9
|
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