NGSolve  4.9
fem/bdbintegrator.hpp
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