NGSolve
4.9
|
00001 #ifndef FILE_BDBINTEGRATOR 00002 #define FILE_BDBINTEGRATOR 00003 00004 /*********************************************************************/ 00005 /* File: bdbintegrator.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 25. Mar. 2000 */ 00008 /*********************************************************************/ 00009 00010 00011 namespace ngfem 00012 { 00013 00014 00015 00016 00017 00018 00024 template <class DMO> 00025 class DMatOp 00026 { 00027 public: 00028 // typedef double TSCAL; 00030 enum { SYMMETRIC = 1 }; 00031 00032 template <typename FEL, typename MIR, typename MAT> 00033 void GenerateMatrixIR (const FEL & fel, const MIR & mir, 00034 const FlatArray<MAT> & dmats, 00035 LocalHeap & lh) const 00036 { 00037 for (int i = 0; i < mir.IR().GetNIP(); i++) 00038 static_cast<const DMO*>(this) -> GenerateMatrix (fel, mir[i], dmats[i], lh); 00039 } 00040 00041 00043 template <typename FEL, typename MIP, typename VEC, typename MAT> 00044 void GenerateLinearizedMatrix (const FEL & fel, const MIP & mip, VEC & vec, 00045 MAT & mat, LocalHeap & lh) const 00046 { 00047 static_cast<const DMO*>(this) -> GenerateMatrix (fel, mip, mat, lh); 00048 } 00049 00051 template <typename FEL, typename MIP, class TVX, class TVY> 00052 void Apply (const FEL & fel, const MIP & mip, 00053 const TVX & x, TVY & y, 00054 LocalHeap & lh) const 00055 { 00056 Mat<DMO::DIM_DMAT, DMO::DIM_DMAT, double> mat; 00057 static_cast<const DMO*>(this) -> GenerateMatrix (fel, mip, mat, lh); 00058 y = mat * x; 00059 } 00060 00061 template <typename FEL, typename MIP, class TVX> 00062 void Apply1 (const FEL & fel, const MIP & mip, 00063 TVX & x, LocalHeap & lh) const 00064 { 00065 Vec<DMO::DIM_DMAT, typename TVX::TSCAL> y; 00066 static_cast<const DMO*>(this) -> Apply (fel, mip, x, y, lh); 00067 x = y; 00068 } 00069 00070 template <typename FEL, typename MIR, typename TVX> 00071 void ApplyIR (const FEL & fel, const MIR & mir, 00072 TVX & x, LocalHeap & lh) const 00073 { 00074 for (int i = 0; i < mir.Size(); i++) 00075 static_cast<const DMO*>(this) -> Apply1 (fel, mir[i], x.Row(i), lh); 00076 } 00077 00078 00079 template <typename FEL, typename MIP, class TVX, class TVY> 00080 void ApplyInv (const FEL & fel, const MIP & mip, 00081 const TVX & x, TVY & y, 00082 LocalHeap & lh) const 00083 { 00084 Mat<DMO::DIM_DMAT, DMO::DIM_DMAT, double> mat; 00085 Mat<DMO::DIM_DMAT, DMO::DIM_DMAT, double> inv; 00086 00087 static_cast<const DMO*>(this) -> GenerateMatrix (fel, mip, mat, lh); 00088 CalcInverse (mat, inv); 00089 y = inv * x; 00090 } 00091 00093 template <typename FEL, typename MIP, class TVX, class TVY> 00094 void ApplyTrans (const FEL & fel, const MIP & mip, 00095 const TVX & x, TVY & y, 00096 LocalHeap & lh) const 00097 { 00098 Mat<DMO::DIM_DMAT, DMO::DIM_DMAT, double> mat; 00099 static_cast<const DMO*>(this) -> GenerateMatrix (fel, mip, mat, lh); 00100 y = Trans (mat) * x; 00101 } 00102 00104 template <typename FEL, typename MIP, class TVX> 00105 double Energy (const FEL & fel, const MIP & mip, 00106 const TVX & x, LocalHeap & lh) const 00107 { 00108 TVX y; 00109 static_cast<const DMO*>(this) -> Apply (fel, mip, x, y, lh); 00110 return 0.5 * InnerProduct (x,y); 00111 } 00112 }; 00113 00114 00115 00116 00117 00118 00119 #ifdef WIN32 00120 #define __restrict__ __restrict 00121 #endif 00122 00123 template <int M> NGS_DLL_HEADER 00124 void FastMat (int n, Complex * ba, Complex * pb, Complex * pc); 00125 00126 template <int M> NGS_DLL_HEADER 00127 void FastMat (int n, double * __restrict__ ba, double * __restrict__ pb, double * __restrict__ pc); 00128 00129 00130 00131 00132 00133 00134 00144 template <class DIFFOP, class DMATOP, class FEL = FiniteElement> 00145 class T_BDBIntegrator : public virtual BilinearFormIntegrator 00146 { 00147 protected: 00148 DMATOP dmatop; 00149 public: 00150 // typedef double TSCAL; // necessary ? 00151 00152 enum { DIM_SPACE = DIFFOP::DIM_SPACE }; 00153 enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT }; 00154 enum { DIM_DMAT = DIFFOP::DIM_DMAT }; 00155 enum { DIM = DIFFOP::DIM }; 00156 // typedef typename DMATOP::TSCAL TSCAL; 00157 00159 T_BDBIntegrator (Array<CoefficientFunction*> & coeffs) 00160 : dmatop(coeffs) 00161 { ; } 00162 00164 T_BDBIntegrator (const DMATOP & admat) 00165 : dmatop(admat) 00166 { ; } 00167 00168 00170 virtual ~T_BDBIntegrator () 00171 { ; } 00172 00174 virtual bool BoundaryForm () const 00175 { return int(DIM_SPACE) > int(DIM_ELEMENT); } 00176 00177 virtual int DimElement () const 00178 { return DIM_ELEMENT; } 00179 00180 virtual int DimSpace () const 00181 { return DIM_SPACE; } 00182 00183 virtual int DimFlux () const 00184 { return DIM_DMAT; } 00185 00186 DMATOP & DMat () { return dmatop; } 00187 const DMATOP & DMat () const { return dmatop; } 00188 00189 virtual void CheckElement (const FiniteElement & el) const 00190 { 00191 if (!dynamic_cast<const FEL*> (&el) ) 00192 throw Exception (string ("Element does not match integrator\n") + 00193 string ("element type is ") + typeid(el).name() + 00194 string (" expected type is ") + typeid(FEL).name() + 00195 string (" integrator is ") + Name()); 00196 } 00197 00198 00199 virtual void ApplyDMat (const FiniteElement & bfel, 00200 const BaseMappedIntegrationPoint & bmip, 00201 const FlatVector<double> & elx, 00202 FlatVector<double> & eldx, 00203 LocalHeap & lh) const 00204 { 00205 dmatop.Apply(static_cast<const FEL&> (bfel), 00206 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> &>(bmip), 00207 elx, eldx ,lh); 00208 } 00209 00210 virtual void ApplyDMat (const FiniteElement & bfel, 00211 const BaseMappedIntegrationPoint & bmip, 00212 const FlatVector<Complex> & elx, 00213 FlatVector<Complex> & eldx, 00214 LocalHeap & lh) const 00215 { 00216 dmatop.Apply(static_cast<const FEL&> (bfel), 00217 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> &>(bmip), 00218 elx,eldx,lh); 00219 } 00220 00221 virtual void ApplyDMat (const FiniteElement & bfel, 00222 const BaseMappedIntegrationRule & bmir, 00223 const FlatMatrix<double> & elx, 00224 FlatMatrix<double> & eldx, 00225 LocalHeap & lh) const 00226 { 00227 const FEL & fel = static_cast<const FEL&> (bfel); 00228 const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir = 00229 static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir); 00230 00231 for (int i = 0; i < mir.Size(); i++) 00232 dmatop.Apply (fel, mir[i], elx.Row(i), eldx.Row(i), lh); 00233 } 00234 00235 00236 virtual void ApplyDMatInv (const FiniteElement & bfel, 00237 const BaseMappedIntegrationRule & bmir, 00238 const FlatMatrix<double> & elx, 00239 FlatMatrix<double> & eldx, 00240 LocalHeap & lh) const 00241 { 00242 const FEL & fel = static_cast<const FEL&> (bfel); 00243 const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir = 00244 static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir); 00245 00246 for (int i = 0; i < mir.Size(); i++) 00247 dmatop.ApplyInv (fel, mir[i], elx.Row(i), eldx.Row(i), lh); 00248 } 00249 00250 00251 00252 00253 00254 00255 virtual void ApplyDMat (const FiniteElement & bfel, 00256 const BaseMappedIntegrationRule & bmir, 00257 const FlatMatrix<Complex> & elx, 00258 FlatMatrix<Complex> & eldx, 00259 LocalHeap & lh) const 00260 { 00261 const FEL & fel = static_cast<const FEL&> (bfel); 00262 const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir = 00263 static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir); 00264 00265 for (int i = 0; i < mir.Size(); i++) 00266 dmatop.Apply (fel, mir[i], elx.Row(i), eldx.Row(i), lh); 00267 } 00268 00269 00270 00271 virtual void 00272 CalcElementMatrix (const FiniteElement & bfel, 00273 const ElementTransformation & eltrans, 00274 FlatMatrix<double> & elmat, 00275 LocalHeap & lh) const 00276 { 00277 T_CalcElementMatrix<double> (bfel, eltrans, elmat, lh); 00278 } 00279 00280 virtual void 00281 CalcElementMatrix (const FiniteElement & bfel, 00282 const ElementTransformation & eltrans, 00283 FlatMatrix<Complex> & elmat, 00284 LocalHeap & lh) const 00285 { 00286 T_CalcElementMatrix<Complex> (bfel, eltrans, elmat, lh); 00287 } 00288 00289 00290 #ifdef TEXT_BOOK_VERSION 00291 00292 template <typename TSCAL> 00293 void T_CalcElementMatrix (const FiniteElement & bfel, 00294 const ElementTransformation & eltrans, 00295 FlatMatrix<TSCAL> & elmat, 00296 LocalHeap & lh) const 00297 { 00298 00299 try 00300 { 00301 const FEL & fel = static_cast<const FEL&> (bfel); 00302 int ndof = fel.GetNDof(); 00303 00304 elmat = 0; 00305 00306 FlatMatrixFixHeight<DIM_DMAT, TSCAL> bmat (ndof * DIM, lh); 00307 FlatMatrixFixHeight<DIM_DMAT, TSCAL> dbmat (ndof * DIM, lh); 00308 Mat<DIM_DMAT,DIM_DMAT,TSCAL> dmat; 00309 00310 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00311 00312 for (int i = 0 ; i < ir.GetNIP(); i++) 00313 { 00314 HeapReset hr(lh); 00315 00316 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> 00317 mip(ir[i], eltrans, lh); 00318 00319 dmatop.GenerateMatrix (fel, mip, dmat, lh); 00320 00321 DIFFOP::GenerateMatrix (fel, mip, bmat, lh); 00322 double fac = mip.GetMeasure() * mip.IP().Weight(); 00323 00324 dbmat = fac * (dmat * bmat); 00325 elmat += Trans (bmat) * dbmat; 00326 } 00327 } 00328 00329 catch (Exception & e) 00330 { 00331 e.Append ("in CalcElementMatrix - textbook, type = "); 00332 e.Append (typeid(*this).name()); 00333 e.Append ("\n"); 00334 throw; 00335 } 00336 catch (exception & e) 00337 { 00338 Exception e2(e.what()); 00339 e2.Append ("\nin CalcElementMatrix - textbook, type = "); 00340 e2.Append (typeid(*this).name()); 00341 e2.Append ("\n"); 00342 throw e2; 00343 } 00344 } 00345 #endif 00346 00347 00348 00349 00350 #ifdef __SSE3__ 00351 #define BLOCK_VERSION 00352 #endif 00353 00354 #ifdef BLOCK_VERSION 00355 00356 template <typename TSCAL> 00357 void T_CalcElementMatrix (const FiniteElement & fel, 00358 const ElementTransformation & eltrans, 00359 FlatMatrix<TSCAL> & elmat, 00360 LocalHeap & lh) const 00361 { 00362 // static int timer = NgProfiler::CreateTimer (string ("Elementmatrix, ") + Name()); 00363 // NgProfiler::RegionTimer reg (timer); 00364 00365 try 00366 { 00367 // const FEL & fel = *static_cast<const FEL*> (&bfel); 00368 int ndof = fel.GetNDof(); 00369 00370 elmat = 0; 00371 00372 00373 enum { BLOCK = 2 * (12 / DIM_DMAT + 1) }; 00374 00375 HeapReset hr1(lh); 00376 00377 FlatMatrixFixHeight<DIM_DMAT, TSCAL> bmat (ndof * DIM, lh); 00378 FlatMatrixFixHeight<DIM_DMAT, TSCAL> dbmat (ndof * DIM, lh); 00379 00380 FlatMatrixFixHeight<2*DIM_DMAT, TSCAL> bmat2 (ndof * DIM, lh); 00381 FlatMatrixFixHeight<2*DIM_DMAT, TSCAL> dbmat2 (ndof * DIM, lh); 00382 00383 FlatMatrixFixHeight<DIM_DMAT*BLOCK, TSCAL> bbmat (ndof * DIM, lh); 00384 FlatMatrixFixHeight<DIM_DMAT*BLOCK, TSCAL> bdbmat (ndof * DIM, lh); 00385 Mat<DIM_DMAT,DIM_DMAT, TSCAL> dmat; 00386 00387 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00388 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 00389 00390 FlatArray<Mat<DIM_DMAT,DIM_DMAT,TSCAL> >dmats(ir.GetNIP(), lh); 00391 dmatop.GenerateMatrixIR (fel, mir, dmats, lh); 00392 00393 int i = 0; 00394 for (int i1 = 0; i1 < ir.GetNIP() / BLOCK; i1++) 00395 { 00396 for (int i2 = 0; i2 < BLOCK; i++, i2++) 00397 { 00398 HeapReset hr(lh); 00399 00400 DIFFOP::GenerateMatrix (fel, mir[i], bmat, lh); 00401 // dmatop.GenerateMatrix (fel, mir[i], dmat, lh); 00402 dmat = dmats[i]; 00403 dmat *= mir[i].GetWeight(); 00404 00405 bbmat.Rows(i2*DIM_DMAT, (i2+1)*DIM_DMAT) = bmat; 00406 bdbmat.Rows(i2*DIM_DMAT, (i2+1)*DIM_DMAT) = dmat * bmat; 00407 } 00408 00409 if (DMATOP::SYMMETRIC) 00410 FastMat<DIM_DMAT*BLOCK> (elmat.Height(), &bdbmat(0,0), &bbmat(0,0), &elmat(0,0)); 00411 else 00412 elmat += Trans (bbmat) * bdbmat; 00413 } 00414 00415 for ( ; i < ir.GetNIP()-1; i+=2) 00416 { 00417 HeapReset hr(lh); 00418 00419 DIFFOP::GenerateMatrix (fel, mir[i], bmat, lh); 00420 // dmatop.GenerateMatrix (fel, mir[i], dmat, lh); 00421 dmat = dmats[i]; 00422 dmat *= mir[i].GetWeight(); 00423 00424 bmat2.Rows(0,DIM_DMAT) = bmat; 00425 dbmat2.Rows(0,DIM_DMAT) = dmat * bmat; 00426 00427 DIFFOP::GenerateMatrix (fel, mir[i+1], bmat, lh); 00428 // dmatop.GenerateMatrix (fel, mir[i+1], dmat, lh); 00429 dmat = dmats[i+1]; 00430 dmat *= mir[i+1].GetWeight(); 00431 00432 bmat2.Rows(DIM_DMAT,2*DIM_DMAT) = bmat; 00433 dbmat2.Rows(DIM_DMAT,2*DIM_DMAT) = dmat * bmat; 00434 00435 if (DMATOP::SYMMETRIC) 00436 FastMat<2*DIM_DMAT> (elmat.Height(), &dbmat2(0,0), &bmat2(0,0), &elmat(0,0)); 00437 else 00438 elmat += Trans (bmat) * dbmat; 00439 } 00440 00441 for ( ; i < ir.GetNIP(); i++) 00442 { 00443 HeapReset hr(lh); 00444 00445 DIFFOP::GenerateMatrix (fel, mir[i], bmat, lh); 00446 // dmatop.GenerateMatrix (fel, mir[i], dmat, lh); 00447 dmat = dmats[i]; 00448 dmat *= mir[i].GetWeight(); 00449 00450 dbmat = dmat * bmat; 00451 00452 if (DMATOP::SYMMETRIC) 00453 // elmat += Symmetric (Trans (dbmat) * bmat); 00454 FastMat<DIM_DMAT> (elmat.Height(), &dbmat(0,0), &bmat(0,0), &elmat(0,0)); 00455 else 00456 elmat += Trans (bmat) * dbmat; 00457 } 00458 00459 if (DMATOP::SYMMETRIC) 00460 { 00461 for (int i = 0; i < elmat.Height(); i++) 00462 for (int j = 0; j < i; j++) 00463 elmat(j,i) = elmat(i,j); 00464 } 00465 } 00466 00467 catch (Exception & e) 00468 { 00469 e.Append ("in CalcElementMatrix - blockversion, type = "); 00470 e.Append (typeid(*this).name()); 00471 e.Append ("\n"); 00472 throw; 00473 } 00474 catch (exception & e) 00475 { 00476 Exception e2(e.what()); 00477 e2.Append ("\nin CalcElementMatrix - blockversion, type = "); 00478 e2.Append (typeid(*this).name()); 00479 e2.Append ("\n"); 00480 throw e2; 00481 } 00482 } 00483 00484 00485 00486 #else // blockversion 00487 // one matrix matrix multiplication 00489 template <typename TSCAL> 00490 void T_CalcElementMatrix (const FiniteElement & fel, 00491 const ElementTransformation & eltrans, 00492 FlatMatrix<TSCAL> & elmat, 00493 LocalHeap & lh) const 00494 { 00495 static Timer timer (string ("Elementmatrix, ") + Name(), 2); 00496 static Timer timer2 (string ("Elementmatrix, ") + Name() + ", Lapack", 3); 00497 RegionTimer reg (timer); 00498 00499 try 00500 { 00501 // const FEL & fel = static_cast<const FEL&> (bfel); 00502 int ndof = fel.GetNDof(); 00503 00504 HeapReset hr(lh); 00505 00506 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00507 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 00508 00509 FlatMatrixFixHeight<DIM_DMAT, TSCAL> bmat (ndof * DIM, lh); 00510 Mat<DIM_DMAT,DIM_DMAT,TSCAL> dmat; 00511 00512 FlatMatrix<TSCAL> bbmat (ndof * DIM, DIM_DMAT*ir.GetNIP(), lh); 00513 FlatMatrix<TSCAL> bdbmat (ndof * DIM, DIM_DMAT*ir.GetNIP(), lh); 00514 00515 for (int i = 0; i < ir.GetNIP(); i++) 00516 { 00517 HeapReset hr(lh); 00518 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = mir[i]; 00519 00520 DIFFOP::GenerateMatrix (fel, mip, bmat, lh); 00521 dmatop.GenerateMatrix (fel, mip, dmat, lh); 00522 dmat *= mip.GetWeight(); 00523 00524 bbmat.Cols(i*DIM_DMAT, (i+1)*DIM_DMAT) = Trans (bmat); 00525 bdbmat.Cols(i*DIM_DMAT, (i+1)*DIM_DMAT) = Trans (dmat * bmat); 00526 } 00527 00528 RegionTimer reg2 (timer2); 00529 00530 if (ndof < 20) 00531 { 00532 if (DMATOP::SYMMETRIC) 00533 elmat = Symmetric ( bdbmat * Trans (bbmat)); 00534 else 00535 elmat = bbmat * Trans (bdbmat); 00536 } 00537 else 00538 elmat = bbmat * Trans(bdbmat) | Lapack; 00539 00540 timer.AddFlops (long(elmat.Height())*long(elmat.Width())*bbmat.Width()); 00541 } 00542 00543 00544 catch (Exception & e) 00545 { 00546 e.Append (string ("in CalcElementMatrix - lapack, type = ") + typeid(*this).name() + "\n"); 00547 throw; 00548 } 00549 catch (exception & e) 00550 { 00551 Exception e2(e.what()); 00552 e2.Append (string ("\nin CalcElementMatrix - lapack, type = ") + typeid(*this).name() + "\n"); 00553 throw e2; 00554 } 00555 } 00556 00557 #endif 00558 00559 00560 00561 00562 00563 virtual void 00564 CalcElementMatrixDiag (const FiniteElement & bfel, 00565 const ElementTransformation & eltrans, 00566 FlatVector<double> & diag, 00567 LocalHeap & lh) const 00568 { 00569 try 00570 { 00571 const FEL & fel = dynamic_cast<const FEL&> (bfel); 00572 int ndof = fel.GetNDof(); 00573 00574 diag.AssignMemory (ndof*DIM, lh); 00575 diag = 0.0; 00576 00577 FlatMatrixFixHeight<DIM_DMAT, double> bmat (ndof * DIM, lh); 00578 Mat<DIM_DMAT,DIM_DMAT> dmat; 00579 00580 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00581 00582 void * heapp = lh.GetPointer(); 00583 for (int i = 0; i < ir.GetNIP(); i++) 00584 { 00585 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip(ir[i], eltrans); 00586 00587 DIFFOP::GenerateMatrix (fel, mip, bmat, lh); 00588 dmatop.GenerateMatrix (fel, mip, dmat, lh); 00589 00590 double fac = 00591 fabs (mip.GetJacobiDet()) * mip.IP().Weight(); 00592 00593 for (int j = 0; j < diag.Size(); j++) 00594 { 00595 Vec<DIM_DMAT> hv = dmat * bmat.Col(j); 00596 diag(j) += fac * InnerProduct (bmat.Col(j), hv); 00597 } 00598 00599 lh.CleanUp (heapp); 00600 } 00601 } 00602 00603 catch (Exception & e) 00604 { 00605 e.Append ("in CalcElementMatrixDiag, type = "); 00606 e.Append (typeid(*this).name()); 00607 e.Append ("\n"); 00608 throw; 00609 } 00610 catch (exception & e) 00611 { 00612 Exception e2(e.what()); 00613 e2.Append ("\nin CalcElementMatrixDiag, type = "); 00614 e2.Append (typeid(*this).name()); 00615 e2.Append ("\n"); 00616 throw e2; 00617 } 00618 } 00619 00620 00621 /* 00623 virtual void 00624 ApplyElementMatrix (const FiniteElement & bfel, 00625 const ElementTransformation & eltrans, 00626 const FlatVector<double> & elx, 00627 FlatVector<double> & ely, 00628 void * precomputed, 00629 LocalHeap & lh) const 00630 { 00631 static Timer timer (string ("BDBIntegrator::ApplyElementMatrix, ") + Name()); 00632 RegionTimer reg (timer); 00633 00634 00635 const FEL & fel = static_cast<const FEL&> (bfel); 00636 int ndof = fel.GetNDof (); 00637 00638 ely = 0; 00639 00640 Vec<DIM_DMAT,double> hv1; 00641 Vec<DIM_DMAT,double> hv2; 00642 00643 FlatVector<double> hely (ndof*DIM, lh); 00644 00645 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00646 00647 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 00648 00649 for (int i = 0; i < ir.GetNIP(); i++) 00650 { 00651 HeapReset hr (lh); 00652 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = mir[i]; 00653 00654 DIFFOP::Apply (fel, mip, elx, hv1, lh); 00655 dmatop.Apply (fel, mip, hv1, hv2, lh); 00656 DIFFOP::ApplyTrans (fel, mip, hv2, hely, lh); 00657 00658 double fac = fabs (mip.GetJacobiDet()) * ir[i].Weight(); 00659 ely += fac * hely; 00660 } 00661 } 00662 00663 00665 virtual void 00666 ApplyElementMatrix (const FiniteElement & bfel, 00667 const ElementTransformation & eltrans, 00668 const FlatVector<Complex> & elx, 00669 FlatVector<Complex> & ely, 00670 void * precomputed, 00671 LocalHeap & lh) const 00672 { 00673 static Timer timer (string ("BDBIntegrator::ApplyElementMatrix<Complex>, ") + Name(), 2); 00674 RegionTimer reg (timer); 00675 00676 const FEL & fel = static_cast<const FEL&> (bfel); 00677 int ndof = fel.GetNDof (); 00678 00679 ely = 0; 00680 00681 Vec<DIM_DMAT,Complex> hv1; 00682 Vec<DIM_DMAT,Complex> hv2; 00683 00684 FlatVector<Complex> hely (ndof*DIM, lh); 00685 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00686 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 00687 00688 for (int i = 0; i < ir.GetNIP(); i++) 00689 { 00690 HeapReset hr (lh); 00691 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = mir[i]; 00692 00693 DIFFOP::Apply (fel, mip, elx, hv1, lh); 00694 dmatop.Apply (fel, mip, hv1, hv2, lh); 00695 DIFFOP::ApplyTrans (fel, mip, hv2, hely, lh); 00696 00697 double fac = fabs (mip.GetJacobiDet()) * mip.IP().Weight(); 00698 ely += fac * hely; 00699 } 00700 } 00701 */ 00702 00703 00704 00705 00706 00707 virtual void 00708 ApplyElementMatrix (const FiniteElement & bfel, 00709 const ElementTransformation & eltrans, 00710 const FlatVector<double> & elx, 00711 FlatVector<double> & ely, 00712 void * precomputed, 00713 LocalHeap & lh) const 00714 { 00715 T_ApplyElementMatrix<double> (bfel, eltrans, elx, ely, precomputed, lh); 00716 } 00717 00718 virtual void 00719 ApplyElementMatrix (const FiniteElement & bfel, 00720 const ElementTransformation & eltrans, 00721 const FlatVector<Complex> & elx, 00722 FlatVector<Complex> & ely, 00723 void * precomputed, 00724 LocalHeap & lh) const 00725 { 00726 T_ApplyElementMatrix<Complex> (bfel, eltrans, elx, ely, precomputed, lh); 00727 } 00728 00729 00730 template <typename TSCAL> 00731 void T_ApplyElementMatrix (const FiniteElement & fel, 00732 const ElementTransformation & eltrans, 00733 const FlatVector<TSCAL> & elx, 00734 FlatVector<TSCAL> & ely, 00735 void * precomputed, 00736 LocalHeap & lh) const 00737 { 00738 static Timer timer (string ("BDBIntegrator::ApplyElementMatrix<") 00739 + typeid(TSCAL).name() + ">" + Name(), 2); 00740 RegionTimer reg (timer); 00741 00742 // const FEL & fel = static_cast<const FEL&> (bfel); 00743 int ndof = fel.GetNDof (); 00744 00745 ely = 0; 00746 00747 Vec<DIM_DMAT,TSCAL> hv1; 00748 Vec<DIM_DMAT,TSCAL> hv2; 00749 00750 FlatVector<TSCAL> hely (ndof*DIM, lh); 00751 const IntegrationRule & ir = GetIntegrationRule (fel,eltrans.HigherIntegrationOrderSet()); 00752 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 00753 00754 for (int i = 0; i < ir.GetNIP(); i++) 00755 { 00756 HeapReset hr (lh); 00757 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = mir[i]; 00758 00759 DIFFOP::Apply (fel, mip, elx, hv1, lh); 00760 dmatop.Apply (fel, mip, hv1, hv2, lh); 00761 DIFFOP::ApplyTrans (fel, mip, hv2, hely, lh); 00762 00763 ely += mip.GetWeight() * hely; 00764 } 00765 } 00766 00767 00768 00769 00770 00771 00772 00774 virtual void 00775 ApplyMixedElementMatrix (const FiniteElement & bfel1, 00776 const FiniteElement & bfel2, 00777 const ElementTransformation & eltrans, 00778 const FlatVector<double> & elx, 00779 FlatVector<double> & ely, 00780 LocalHeap & lh) const 00781 { 00782 const FEL & fel1 = dynamic_cast<const FEL&> (bfel1); 00783 const FEL & fel2 = dynamic_cast<const FEL&> (bfel2); 00784 // int ndof1 = fel1.GetNDof (); 00785 int ndof2 = fel2.GetNDof (); 00786 00787 ely = 0; 00788 00789 Vec<DIM_DMAT,double> hv1; 00790 Vec<DIM_DMAT,double> hv2; 00791 00792 FlatVector<double> hely (ndof2*DIM, lh); 00793 00794 const IntegrationRule & ir = GetIntegrationRule (fel2,eltrans.HigherIntegrationOrderSet()); 00795 00796 void * heapp = lh.GetPointer(); 00797 for (int i = 0; i < ir.GetNIP(); i++) 00798 { 00799 const IntegrationPoint & ip = ir[i]; 00800 00801 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip (ir[i], eltrans); 00802 00803 DIFFOP::Apply (fel1, mip, elx, hv1, lh); 00804 dmatop.Apply (fel1, mip, hv1, hv2, lh); 00805 DIFFOP::ApplyTrans (fel2, mip, hv2, hely, lh); 00806 00807 double fac = fabs (mip.GetJacobiDet()) * ip.Weight(); 00808 ely += fac * hely; 00809 lh.CleanUp (heapp); 00810 } 00811 } 00812 00813 00814 00815 IntegrationRule GetIntegrationRule (const FiniteElement & fel, 00816 const bool use_higher_integration_order = false) const 00817 { 00818 int order = 2 * fel.Order(); 00819 00820 ELEMENT_TYPE et = fel.ElementType(); 00821 00822 if (et == ET_TET || et == ET_TRIG || et == ET_SEGM) 00823 order -= 2 * DIFFOP::DIFFORDER; 00824 00825 if (common_integration_order >= 0) 00826 order = common_integration_order; 00827 00828 if (integration_order >= 0) 00829 order = integration_order; 00830 00831 if(use_higher_integration_order && higher_integration_order > order) 00832 { 00833 //(*testout) << "changing order " << order << " to " << higher_integration_order << endl; 00834 order = higher_integration_order; 00835 } 00836 00837 // return SelectIntegrationRule (et, order); 00838 return IntegrationRule (et, order); 00839 } 00840 00841 00842 00843 00844 00845 00846 00847 00848 00849 00850 virtual void 00851 CalcFlux (const FiniteElement & fel, 00852 const BaseMappedIntegrationPoint & bmip, 00853 const FlatVector<double> & elx, 00854 FlatVector<double> & flux, 00855 bool applyd, 00856 LocalHeap & lh) const 00857 { 00858 // const FEL & fel = static_cast<const FEL&> (bfel); 00859 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = 00860 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00861 00862 00863 DIFFOP::Apply (fel, mip, elx, flux, lh); 00864 if (applyd) 00865 dmatop.Apply1 (fel, mip, flux, lh); 00866 } 00867 00868 virtual void 00869 CalcFlux (const FiniteElement & fel, 00870 const BaseMappedIntegrationRule & bmir, 00871 const FlatVector<double> & elx, 00872 FlatMatrix<double> & flux, 00873 bool applyd, 00874 LocalHeap & lh) const 00875 { 00876 // const FEL & fel = static_cast<const FEL&> (bfel); 00877 const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir = 00878 static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir); 00879 00880 00881 DIFFOP::ApplyIR (fel, mir, elx, flux, lh); 00882 if (applyd) 00883 dmatop.ApplyIR (fel, mir, flux, lh); 00884 } 00885 00886 virtual void 00887 CalcFlux (const FiniteElement & fel, 00888 const BaseMappedIntegrationPoint & bmip, 00889 const FlatVector<Complex> & elx, 00890 FlatVector<Complex> & flux, 00891 bool applyd, 00892 LocalHeap & lh) const 00893 { 00894 // const FEL & fel = static_cast<const FEL&> (bfel); 00895 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = 00896 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00897 00898 if (applyd) 00899 { 00900 Vec<DIM_DMAT,Complex> hv1; 00901 DIFFOP::Apply (fel, mip, elx, hv1, lh); 00902 dmatop.Apply (fel, mip, hv1, flux, lh); 00903 } 00904 else 00905 { 00906 DIFFOP::Apply (fel, mip, elx, flux, lh); 00907 } 00908 } 00909 00910 00911 00912 00913 00914 00915 virtual void 00916 CalcFluxMulti (const FiniteElement & bfel, 00917 const BaseMappedIntegrationPoint & bmip, 00918 int m, 00919 const FlatVector<double> & elx, 00920 FlatVector<double> & flux, 00921 bool applyd, 00922 LocalHeap & lh) const 00923 { 00924 const FEL & fel = static_cast<const FEL&> (bfel); 00925 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = 00926 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00927 00928 int ndof = fel.GetNDof(); 00929 FlatMatrixFixHeight<DIM_DMAT> bmat (ndof * DIM, lh); 00930 00931 DIFFOP::GenerateMatrix (fel, mip, bmat, lh); 00932 00933 if (applyd) 00934 { 00935 Vec<DIM_DMAT> hv1; 00936 Mat<DIM_DMAT,DIM_DMAT> dmat; 00937 dmatop.GenerateMatrix (fel, mip, dmat, lh); 00938 00939 for (int i = 0; i < m; i++) 00940 { 00941 SliceVector<double> slice_x (ndof*DIM, m, &const_cast<double&> (elx(i))); 00942 SliceVector<double> slice_flux (DIM_DMAT, m, &flux(i)); 00943 hv1 = bmat * slice_x; 00944 slice_flux = dmat * hv1; 00945 } 00946 } 00947 else 00948 { 00949 for (int i = 0; i < m; i++) 00950 { 00951 // SliceVector<double> slice_x (ndof*DIM, m, &const_cast<double&> (elx(i))); 00952 SliceVector<double> slice_x (ndof*DIM, m, &elx(i)); 00953 SliceVector<double> slice_flux (DIM_DMAT, m, &flux(i)); 00954 slice_flux = bmat * slice_x; 00955 } 00956 } 00957 } 00958 00959 00960 00961 00962 00963 virtual void 00964 ApplyBTrans (const FiniteElement & fel, 00965 const BaseMappedIntegrationPoint & bmip, 00966 const FlatVector<double> & elx, 00967 FlatVector<double> & ely, 00968 LocalHeap & lh) const 00969 { 00970 // const FEL & fel = static_cast<const FEL&> (bfel); 00971 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip = 00972 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00973 00974 DIFFOP::ApplyTrans (fel, mip, elx, ely, lh); 00975 } 00976 00977 00978 virtual void 00979 ApplyBTrans (const FiniteElement & fel, 00980 const BaseMappedIntegrationPoint & bmip, 00981 const FlatVector<Complex> & elx, 00982 FlatVector<Complex> & ely, 00983 LocalHeap & lh) const 00984 { 00985 // const FEL & fel = static_cast<const FEL&> (bfel); 00986 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip = 00987 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00988 00989 DIFFOP::ApplyTrans (fel, mip, elx, ely, lh); 00990 } 00991 00992 00993 virtual void 00994 ApplyBTrans (const FiniteElement & fel, 00995 const BaseMappedIntegrationRule & bmir, 00996 const FlatMatrix<double> & elx, 00997 FlatVector<double> & ely, 00998 LocalHeap & lh) const 00999 { 01000 // const FEL & fel = static_cast<const FEL&> (bfel); 01001 const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir = 01002 static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir); 01003 01004 DIFFOP::ApplyTransIR (fel, mir, elx, ely, lh); 01005 } 01006 01008 virtual int GetDimension () const { return DIM; } 01009 }; 01010 01011 01012 01013 01014 01015 01016 01017 01018 01019 01020 01021 01022 01023 01024 01025 01026 01027 01028 01029 01030 01031 01032 01033 01034 template <class DIFFOP, class DMATOP, class FEL = FiniteElement> 01035 class T_NonlinearBDBIntegrator : public T_BDBIntegrator<DIFFOP, DMATOP, FEL> 01036 { 01037 protected: 01038 enum { DIM_SPACE = DIFFOP::DIM_SPACE }; 01039 enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT }; 01040 enum { DIM_DMAT = DIFFOP::DIM_DMAT }; 01041 enum { DIM = DIFFOP::DIM }; 01042 01043 using T_BDBIntegrator<DIFFOP,DMATOP,FEL>::GetIntegrationRule; 01044 01045 public: 01047 T_NonlinearBDBIntegrator (const DMATOP & admat) 01048 : T_BDBIntegrator<DIFFOP,DMATOP,FEL> (admat) 01049 { ; } 01050 01052 virtual ~T_NonlinearBDBIntegrator () 01053 { ; } 01054 01055 01056 01057 virtual void 01058 CalcLinearizedElementMatrix (const FiniteElement & fel, 01059 const ElementTransformation & eltrans, 01060 FlatVector<double> & elveclin, 01061 FlatMatrix<double> & elmat, 01062 LocalHeap & lh) const 01063 { 01064 static Timer maintimer (string ("NonlinearBDB, CalcLinearized, ") + this->Name(), 2); 01065 static Timer bdbtimer ("NonlinearBDB, bdb product", 3); 01066 RegionTimer reg(maintimer); 01067 01068 try 01069 { 01070 // const FEL & fel = dynamic_cast<const FEL&> (bfel); 01071 int ndof = fel.GetNDof(); 01072 01073 elmat = 0; 01074 01075 FlatMatrixFixHeight<DIM_DMAT, double> bmat (ndof * DIM, lh); 01076 FlatMatrixFixHeight<DIM_DMAT, double> dbmat (ndof * DIM, lh); 01077 Vec<DIM_DMAT,double> hvlin; 01078 01079 Mat<DIM_DMAT,DIM_DMAT> dmat; 01080 01081 const IntegrationRule & ir = GetIntegrationRule (fel); 01082 01083 void * heapp = lh.GetPointer(); 01084 for (int i = 0; i < ir.GetNIP(); i++) 01085 { 01086 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip(ir[i], eltrans); 01087 01088 DIFFOP::Apply (fel, mip, elveclin, hvlin, lh); 01089 01090 DIFFOP::GenerateMatrix (fel, mip, bmat, lh); 01091 01092 this->dmatop . GenerateLinearizedMatrix (fel, mip, hvlin, dmat, lh); 01093 01094 double fac = 01095 fabs (mip.GetJacobiDet()) * mip.IP().Weight(); 01096 01097 { 01098 NgProfiler::RegionTimer reg(bdbtimer); 01099 dbmat = fac * (dmat * bmat); 01100 elmat += Trans (bmat) * dbmat; 01101 } 01102 01103 lh.CleanUp (heapp); 01104 } 01105 } 01106 01107 catch (Exception & e) 01108 { 01109 e.Append ("in CalcLinearizedElementMatrix, type = "); 01110 e.Append (typeid(*this).name()); 01111 e.Append ("\n"); 01112 throw; 01113 } 01114 catch (exception & e) 01115 { 01116 Exception e2(e.what()); 01117 e2.Append ("\nin CalcLinearizedElementMatrix, type = "); 01118 e2.Append (typeid(*this).name()); 01119 e2.Append ("\n"); 01120 throw e2; 01121 } 01122 } 01123 01124 01125 01126 virtual void 01127 CalcLinearizedElementMatrix (const FiniteElement & bfel, 01128 const ElementTransformation & eltrans, 01129 FlatVector<Complex> & elveclin, 01130 FlatMatrix<Complex> & elmat, 01131 LocalHeap & lh) const 01132 { 01133 static Timer maintimer (string ("NonlinearBDB, CalcLinearized<Complex>, ") + this->Name(), 2); 01134 static Timer bdbtimer ("NonlinearBDB, bdb product", 3); 01135 RegionTimer reg(maintimer); 01136 01137 try 01138 { 01139 const FEL & fel = dynamic_cast<const FEL&> (bfel); 01140 int ndof = fel.GetNDof(); 01141 01142 elmat = 0; 01143 01144 FlatMatrixFixHeight<DIM_DMAT, Complex> bmat (ndof * DIM, lh); 01145 FlatMatrixFixHeight<DIM_DMAT, Complex> dbmat (ndof * DIM, lh); 01146 Vec<DIM_DMAT,Complex> hvlin; 01147 01148 Mat<DIM_DMAT,DIM_DMAT, Complex> dmat; 01149 01150 const IntegrationRule & ir = GetIntegrationRule (fel); 01151 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 01152 01153 01154 FlatMatrix<Complex> bbmat (ndof * DIM, DIM_DMAT*ir.GetNIP(), lh); 01155 FlatMatrix<Complex> bdbmat (ndof * DIM, DIM_DMAT*ir.GetNIP(), lh); 01156 01157 01158 for (int i = 0; i < ir.GetNIP(); i++) 01159 { 01160 HeapReset hr(lh); 01161 01162 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = mir[i]; 01163 01164 DIFFOP::Apply (fel, mip, elveclin, hvlin, lh); 01165 01166 DIFFOP::GenerateMatrix (fel, mip, bmat, lh); 01167 this->dmatop . GenerateLinearizedMatrix (fel, mip, hvlin, dmat, lh); 01168 01169 double fac = mip.GetWeight(); 01170 /* 01171 dbmat = dmat * bmat; 01172 elmat += fac * (Trans (bmat) * dbmat); 01173 */ 01174 01175 dmat *= fac; 01176 bbmat.Cols(i*DIM_DMAT, (i+1)*DIM_DMAT) = Trans (bmat); 01177 bdbmat.Cols(i*DIM_DMAT, (i+1)*DIM_DMAT) = Trans (dmat * bmat); 01178 } 01179 01180 01181 RegionTimer reg(bdbtimer); 01182 01183 if (ndof < 20) 01184 { 01185 if (DMATOP::SYMMETRIC) 01186 elmat = Symmetric ( bdbmat * Trans (bbmat)); 01187 else 01188 elmat = bbmat * Trans (bdbmat); 01189 } 01190 else 01191 // LapackMultABt (bbmat, bdbmat, elmat); 01192 LapackMult (bbmat, Trans(bdbmat), elmat); 01193 01194 01195 } 01196 01197 catch (Exception & e) 01198 { 01199 e.Append ("in CalcLinearizedElementMatrix, type = "); 01200 e.Append (typeid(*this).name()); 01201 e.Append ("\n"); 01202 throw; 01203 } 01204 catch (exception & e) 01205 { 01206 Exception e2(e.what()); 01207 e2.Append ("\nin CalcLinearizedElementMatrix, type = "); 01208 e2.Append (typeid(*this).name()); 01209 e2.Append ("\n"); 01210 throw e2; 01211 } 01212 } 01213 01214 01215 01216 01217 01218 01219 01221 virtual void 01222 ApplyLinearizedElementMatrix (const FiniteElement & fel, 01223 const ElementTransformation & eltrans, 01224 const FlatVector<double> & ellin, 01225 const FlatVector<double> & elx, 01226 FlatVector<double> & ely, 01227 LocalHeap & lh) const 01228 { 01229 // const FEL & fel = dynamic_cast<const FEL&> (bfel); 01230 int ndof = fel.GetNDof (); 01231 01232 ely = 0; 01233 01234 Vec<DIM_DMAT,double> hvlin; 01235 Vec<DIM_DMAT,double> hvx; 01236 Vec<DIM_DMAT,double> hvy; 01237 01238 FlatVector<double> hely (ndof*DIM, lh); 01239 01240 const IntegrationRule & ir = GetIntegrationRule (fel); 01241 01242 for (int i = 0; i < ir.GetNIP(); i++) 01243 { 01244 const IntegrationPoint & ip = ir[i]; 01245 01246 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip (ir[i], eltrans); 01247 01248 DIFFOP::Apply (fel, mip, ellin, hvlin, lh); 01249 DIFFOP::Apply (fel, mip, elx, hvx, lh); 01250 this->dmatop.ApplyLinearized (fel, mip, hvlin, hvx, hvy, lh); 01251 DIFFOP::ApplyTrans (fel, mip, hvy, hely, lh); 01252 01253 double fac = fabs (mip.GetJacobiDet()) * ip.Weight(); 01254 ely += fac * hely; 01255 } 01256 } 01257 01258 01259 01261 virtual void 01262 ApplyLinearizedElementMatrix (const FiniteElement & bfel, 01263 const ElementTransformation & eltrans, 01264 const FlatVector<Complex> & ellin, 01265 const FlatVector<Complex> & elx, 01266 FlatVector<Complex> & ely, 01267 LocalHeap & lh) const 01268 { 01269 const FEL & fel = dynamic_cast<const FEL&> (bfel); 01270 int ndof = fel.GetNDof (); 01271 01272 ely = 0; 01273 01274 Vec<DIM_DMAT,Complex> hvlin; 01275 Vec<DIM_DMAT,Complex> hvx; 01276 Vec<DIM_DMAT,Complex> hvy; 01277 01278 FlatVector<Complex> hely (ndof*DIM, lh); 01279 01280 /* 01281 int order = IntegrationOrder (fel); 01282 const IntegrationRule & ir = 01283 GetIntegrationRules().SelectIntegrationRule (fel.ElementType(), order); 01284 */ 01285 const IntegrationRule & ir = GetIntegrationRule (fel); 01286 01287 01288 01289 for (int i = 0; i < ir.GetNIP(); i++) 01290 { 01291 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip (ir[i], eltrans); 01292 01293 DIFFOP::Apply (fel, mip, ellin, hvlin, lh); 01294 DIFFOP::Apply (fel, mip, elx, hvx, lh); 01295 this->dmatop.ApplyLinearized (fel, mip, hvlin, hvx, hvy, lh); 01296 DIFFOP::ApplyTrans (fel, mip, hvy, hely, lh); 01297 01298 double fac = fabs (mip.GetJacobiDet()) * mip.IP().Weight(); 01299 ely += fac * hely; 01300 } 01301 } 01302 01303 01304 01305 virtual double Energy (const FiniteElement & bfel, 01306 const ElementTransformation & eltrans, 01307 const FlatVector<double> & elx, 01308 LocalHeap & lh) const 01309 { 01310 const FEL & fel = dynamic_cast<const FEL&> (bfel); 01311 // int ndof = fel.GetNDof (); 01312 01313 Vec<DIM_DMAT,double> hvx; 01314 const IntegrationRule & ir = GetIntegrationRule (fel); 01315 01316 double energy = 0; 01317 01318 for (int i = 0; i < ir.GetNIP(); i++) 01319 { 01320 const IntegrationPoint & ip = ir[i]; 01321 01322 MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> mip (ir[i], eltrans); 01323 DIFFOP::Apply (fel, mip, elx, hvx, lh); 01324 01325 double fac = fabs (mip.GetJacobiDet()) * ip.Weight(); 01326 energy += fac * this->dmatop.Energy (fel, mip, hvx, lh); 01327 } 01328 01329 return energy; 01330 } 01331 01332 }; 01333 01334 01335 01336 01337 01338 01339 01340 01341 01342 01343 01344 01345 01346 01347 01348 01349 01350 01351 01352 01353 01354 01355 01356 01357 01358 01363 template <class DIFFOP, class DVecOp, class FEL = FiniteElement> 01364 class T_BIntegrator : public LinearFormIntegrator 01365 { 01366 protected: 01367 DVecOp dvecop; 01368 01369 public: 01370 enum { DIM_SPACE = DIFFOP::DIM_SPACE }; 01371 enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT }; 01372 enum { DIM_DMAT = DIFFOP::DIM_DMAT }; 01373 enum { DIM = DIFFOP::DIM }; 01374 // typedef typename DVecOp::TSCAL TSCAL; 01375 01377 T_BIntegrator (const DVecOp & advec) 01378 : dvecop(advec) 01379 { ; } 01380 01382 virtual ~T_BIntegrator () 01383 { ; } 01384 01386 virtual void CheckElement (const FiniteElement & el) const 01387 { 01388 if (!dynamic_cast<const FEL*> (&el) ) 01389 throw Exception (string ("Element does not match integrator\n") + 01390 string ("element type is ") + typeid(el).name() + 01391 string (" expected type is ") + typeid(FEL).name() + 01392 string ("integrator is ") + Name()); 01393 } 01394 01396 virtual bool BoundaryForm () const 01397 { return int(DIM_SPACE) > int(DIM_ELEMENT); } 01398 01399 virtual int DimElement () const 01400 { return DIM_ELEMENT; } 01401 01402 virtual int DimSpace () const 01403 { return DIM_SPACE; } 01404 01405 01406 01407 virtual void 01408 CalcElementVector (const FiniteElement & bfel, 01409 const ElementTransformation & eltrans, 01410 FlatVector<double> & elvec, 01411 LocalHeap & lh) const 01412 { 01413 T_CalcElementVector<double> (bfel, eltrans, elvec, lh); 01414 } 01415 01416 virtual void 01417 CalcElementVector (const FiniteElement & bfel, 01418 const ElementTransformation & eltrans, 01419 FlatVector<Complex> & elvec, 01420 LocalHeap & lh) const 01421 { 01422 T_CalcElementVector<Complex> (bfel, eltrans, elvec, lh); 01423 } 01424 01425 01426 01428 01429 template <typename TSCAL> 01430 void T_CalcElementVector (const FiniteElement & fel, 01431 const ElementTransformation & eltrans, 01432 FlatVector<TSCAL> & elvec, 01433 LocalHeap & lh) const 01434 { 01435 try 01436 { 01437 // const FEL & fel = static_cast<const FEL&> (bfel); 01438 int ndof = fel.GetNDof(); 01439 01440 elvec = 0; 01441 FlatVector<TSCAL> hv(ndof * DIM, lh); 01442 Vec<DIM_DMAT, TSCAL> dvec; 01443 01444 IntegrationRule ir(fel.ElementType(), IntegrationOrder(fel)); 01445 MappedIntegrationRule<DIM_ELEMENT, DIM_SPACE> mir(ir, eltrans, lh); 01446 01447 // FlatMatrix<TSCAL> dvecs(ir.GetNIP(), DIM_DMAT, lh); 01448 FlatMatrixFixWidth<DIM_DMAT, TSCAL> dvecs(ir.GetNIP(), lh); 01449 dvecop.GenerateVectorIR (fel, mir, dvecs, lh); 01450 01451 for (int i = 0; i < ir.GetNIP(); i++) 01452 { 01453 HeapReset hr(lh); 01454 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = mir[i]; 01455 01456 // dvecop.GenerateVector (fel, mip, dvec, lh); 01457 DIFFOP::ApplyTrans (fel, mip, dvecs.Row(i), hv, lh); 01458 01459 elvec += mip.GetWeight() * hv; 01460 } 01461 } 01462 01463 catch (Exception & e) 01464 { 01465 e.Append ("in CalcElementVector, type = "); 01466 e.Append (typeid(*this).name()); 01467 e.Append ("\n"); 01468 throw; 01469 } 01470 catch (exception & e) 01471 { 01472 Exception e2(e.what()); 01473 e2.Append ("\nin CalcElementVector, type = "); 01474 e2.Append (typeid(*this).name()); 01475 e2.Append ("\n"); 01476 throw e2; 01477 } 01478 } 01479 01480 01481 01482 01483 01484 virtual void 01485 CalcElementVectorIndependent (const FiniteElement & gfel, 01486 const BaseMappedIntegrationPoint & s_mip, 01487 const BaseMappedIntegrationPoint & g_mip, 01488 FlatVector<double> & elvec, 01489 LocalHeap & lh, 01490 const bool curveint = false) const 01491 { 01492 T_CalcElementVectorIndependent (gfel, s_mip, g_mip, elvec, lh, curveint); 01493 } 01494 01495 virtual void 01496 CalcElementVectorIndependent (const FiniteElement & gfel, 01497 const BaseMappedIntegrationPoint & s_mip, 01498 const BaseMappedIntegrationPoint & g_mip, 01499 FlatVector<Complex> & elvec, 01500 LocalHeap & lh, 01501 const bool curveint = false) const 01502 { 01503 T_CalcElementVectorIndependent (gfel, s_mip, g_mip, elvec, lh, curveint); 01504 } 01505 01506 01507 template <typename TSCAL> 01508 void T_CalcElementVectorIndependent (const FiniteElement & gfel, 01509 const BaseMappedIntegrationPoint & s_mip, 01510 const BaseMappedIntegrationPoint & g_mip, 01511 FlatVector<TSCAL> & elvec, 01512 LocalHeap & lh, 01513 const bool curveint = false) const 01514 { 01515 const FEL & fel = dynamic_cast<const FEL&> (gfel); 01516 int ndof = fel.GetNDof(); 01517 01518 elvec.AssignMemory (ndof * DIM, lh); 01519 //elvec = 0; 01520 01521 Vec<DIM_DMAT, TSCAL> dvec; 01522 01523 const MappedIntegrationPoint< DIM_SPACE, DIM_SPACE > & d_g_mip 01524 (static_cast<const MappedIntegrationPoint< DIM_SPACE, DIM_SPACE > &>(g_mip)); 01525 01526 if(curveint) 01527 { 01528 const MappedIntegrationPoint< 1, DIM_SPACE > & d_s_mip 01529 (static_cast<const MappedIntegrationPoint< 1, DIM_SPACE > &>(s_mip)); 01530 01531 dvecop.GenerateVector (fel, d_s_mip, dvec, lh); 01532 } 01533 else 01534 { 01535 enum { HDIM = (DIM_SPACE > 1) ? DIM_SPACE-1 : 1 }; 01536 01537 const MappedIntegrationPoint< HDIM, DIM_SPACE > & d_s_mip 01538 (static_cast<const MappedIntegrationPoint< HDIM, DIM_SPACE > &>(s_mip)); 01539 01540 dvecop.GenerateVector (fel, d_s_mip, dvec, lh); 01541 } 01542 01543 01544 DIFFOP::ApplyTrans (fel, d_g_mip, dvec, elvec, lh); 01545 01546 //(*testout) << "dvec " << dvec << " elvec " << elvec << endl; 01547 01548 } 01549 01550 01551 01552 01553 int IntegrationOrder (const FiniteElement & fel) const 01554 { 01555 // int order = fel.Order()+2; // low order case 01556 int order = 2*fel.Order()+1; // high order case 01557 01558 ELEMENT_TYPE et = fel.ElementType(); 01559 01560 if (et == ET_TET || et == ET_TRIG || et == ET_SEGM) 01561 order -= DIFFOP::DIFFORDER; 01562 01563 /* 01564 if (common_integration_order >= 0) 01565 order = common_integration_order; 01566 */ 01567 if (integration_order >= 0) 01568 order = integration_order; 01569 01570 return order; 01571 } 01572 01573 01575 virtual int GetDimension () const { return DIM; } 01576 01578 virtual string Name () const { return "BDB integrator"; } 01579 }; 01580 01581 01582 01583 } 01584 01585 01586 #endif