NGSolve  4.9
fem/bdbequations.hpp
00001 #ifndef FILE_BDBEQUATIONS
00002 #define FILE_BDBEQUATIONS
00003 
00004 /*********************************************************************/
00005 /* File:   bdbequations.hpp                                          */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   25. Mar. 2000                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngfem
00011 {
00012 
00013 /* 
00014    realizations of bdb integrators for many equations.
00015    The differential operators provide the B-matrix,
00016    the DMatOps provide the coefficient tensors
00017 */
00018 
00019 
00021 template <int D, typename FEL = ScalarFiniteElement<D> >
00022 class DiffOpGradient : public DiffOp<DiffOpGradient<D, FEL> >
00023 {
00024 public:
00025   enum { DIM = 1 };
00026   enum { DIM_SPACE = D };
00027   enum { DIM_ELEMENT = D };
00028   enum { DIM_DMAT = D };
00029   enum { DIFFORDER = 1 };
00030   
00031   static const FEL & Cast (const FiniteElement & fel) 
00032   { return static_cast<const FEL&> (fel); }
00033 
00035   template <typename MIP, typename MAT>
00036   static void GenerateMatrix (const FiniteElement & fel, 
00037                               const MIP & mip,
00038                               MAT & mat, LocalHeap & lh)
00039   {
00040     mat = Trans (mip.GetJacobianInverse ()) * 
00041       Trans (Cast(fel).GetDShape(mip.IP(),lh));
00042   }
00043 
00044   static void GenerateMatrix (const FiniteElement & fel, 
00045                               const MappedIntegrationPoint<D,D> & mip,
00046                               FlatMatrixFixHeight<D> & mat, LocalHeap & lh)
00047   {
00048     Cast(fel).CalcMappedDShape (mip, mat.Trans());
00049   }
00050 
00052   template <typename MIP, class TVX, class TVY>
00053   static void Apply (const FiniteElement & fel, const MIP & mip,
00054                      const TVX & x, TVY & y,
00055                      LocalHeap & lh) 
00056   {
00057     typedef typename TVX::TSCAL TSCAL;
00058 
00059     Vec<D,TSCAL> hv = Trans (Cast(fel).GetDShape(mip.IP(), lh)) * x;
00060     y = Trans (mip.GetJacobianInverse()) * hv;
00061   }
00062 
00063   template <class MIR>
00064   static void ApplyIR (const FiniteElement & fel, const MIR & mir,
00065                        const FlatVector<double> & x, FlatMatrix<double> & y,
00066                        LocalHeap & lh)
00067   {
00068     FlatMatrixFixWidth<D> grad(mir.Size(), &y(0));
00069     Cast(fel).EvaluateGrad (mir.IR(), x, grad);
00070     for  (int i = 0; i < mir.Size(); i++)
00071       {
00072         Vec<D> hv = grad.Row(i);
00073         grad.Row(i) = Trans (mir[i].GetJacobianInverse()) * hv;
00074       }
00075 
00076     // for  (int i = 0; i < mir.Size(); i++)
00077     //   grad.Row(i) = Trans (mir[i].GetJacobianInverse()) * Vec<D> (grad.Row(i));
00078   }
00079 
00080 
00082   template <typename MIP, class TVX, class TVY>
00083   static void ApplyTrans (const FiniteElement & fel, const MIP & mip,
00084                           const TVX & x, TVY & y,
00085                           LocalHeap & lh) 
00086   {
00087     typedef typename TVX::TSCAL TSCAL;
00088 
00089     Vec<D,TSCAL> hv = mip.GetJacobianInverse() * x;
00090     y = Cast(fel).GetDShape(mip.IP(),lh) * hv;
00091   }
00092 
00093 
00095   template <typename MIP, class TVX>
00096   static void Transform (const MIP & mip, TVX & x)
00097   {
00098     Vec<D> hx = Trans (mip.GetJacobianInverse()) * x;
00099     x = hx; 
00100   }
00101 
00103   template <typename MIP, class TVX>
00104   static void TransformT (const MIP & mip, TVX & x)
00105   {
00106     Vec<D> hx = mip.GetJacobianInverse() * x;
00107     x = hx; 
00108   }
00109   
00110   /*
00112   template <class TVD, class TVY>
00113   static void ApplyGrid (const FiniteElement & fel, 
00114                          const IntegrationRuleTP<D> & ir,
00115                          const TVY & hv, 
00116                          TVD & dvecs, 
00117                          LocalHeap & lh)
00118   {
00119     FlatMatrix<double> dmat(dvecs.Size(), D, const_cast<double*> (&dvecs[0](0)));
00120     Cast(fel).EvaluateDShapeGrid (ir, hv, dmat, lh);
00121   }
00122 
00124   template <class TVD, class TVY>
00125   static void ApplyTransGrid (const FiniteElement & fel, 
00126                               const IntegrationRuleTP<D> & ir,
00127                               const TVD & dvecs, 
00128                               TVY & hv, LocalHeap & lh)
00129   {
00130     FlatMatrix<double> dmat(dvecs.Size(), D, const_cast<double*> (&dvecs[0](0)));
00131     Cast(fel).EvaluateDShapeGridTrans (ir, dmat, hv, lh);
00132   }
00133   */
00134 };
00135 
00136 
00137 
00138 
00140 template <int D, typename FEL = ScalarFiniteElement<D-1> >
00141 class DiffOpGradientBoundary : public DiffOp<DiffOpGradientBoundary<D, FEL> >
00142 {
00143 public:
00144   enum { DIM = 1 };
00145   enum { DIM_SPACE = D };
00146   enum { DIM_ELEMENT = D-1 };
00147   enum { DIM_DMAT = D };
00148   enum { DIFFORDER = 1 };
00149 
00151   template <typename AFEL, typename MIP, typename MAT>
00152   static void GenerateMatrix (const AFEL & fel, const MIP & mip,
00153                               MAT & mat, LocalHeap & lh)
00154   {
00155     mat = Trans (mip.GetJacobianInverse ()) * 
00156       Trans (static_cast<const FEL&>(fel).GetDShape(mip.IP(),lh));
00157   }
00158 };
00159 
00160 
00161 
00162 
00163 
00165 template <int D>
00166 class DiffOpGradientRotSym : 
00167   public DiffOp<DiffOpGradientRotSym<D> >
00168 {
00169 public:
00170   enum { DIM = 1 };
00171   enum { DIM_SPACE = D };
00172   enum { DIM_ELEMENT = D };
00173   enum { DIM_DMAT = D };
00174   enum { DIFFORDER = 1 };
00175 
00177   template <typename FEL, typename MIP, typename MAT>
00178   static void GenerateMatrix (const FEL & fel, const MIP & mip,
00179                               MAT & mat, LocalHeap & lh)
00180   {
00181     typedef typename MAT::TSCAL TSCAL;
00182 
00183     mat =  Trans (mip.GetJacobianInverse ()) * 
00184       Trans (fel.GetDShape(mip.IP(),lh));
00185 
00186     int i;
00187     double cx = mip.GetPoint()(0);
00188     if (cx == 0) cx = 1e-10;
00189     for (int i = 0; i < mat.Width(); i++)
00190       mat(0,i) += fel.GetShape(mip.IP(), lh)(i) / cx;
00191 
00192     // do the rot
00193     for (int i = 0; i < mat.Width(); i++)
00194       {
00195         TSCAL hv = mat(0,i);
00196         mat(0,i) = mat(1,i);
00197         mat(1,i) = -hv;
00198       }
00199   }
00200 };
00201 
00202 
00203 
00205 template <int D, typename FEL = ScalarFiniteElement<D> >
00206 class DiffOpId : public DiffOp<DiffOpId<D, FEL> >
00207 {
00208 public:
00209   enum { DIM = 1 };
00210   enum { DIM_SPACE = D };
00211   enum { DIM_ELEMENT = D };
00212   enum { DIM_DMAT = 1 };
00213   enum { DIFFORDER = 0 };
00214 
00215   static const FEL & Cast (const FiniteElement & fel) 
00216   { return static_cast<const FEL&> (fel); }
00217 
00218   template <typename MIP, typename MAT>
00219   static void GenerateMatrix (const FiniteElement & fel, const MIP & mip,
00220                               MAT & mat, LocalHeap & lh)
00221   {
00222     mat.Row(0) = Cast(fel).GetShape(mip.IP(), lh);
00223   }
00224 
00225   static void GenerateMatrix (const FiniteElement & fel, 
00226                               const MappedIntegrationPoint<D,D> & mip,
00227                               FlatMatrixFixHeight<1> & mat, LocalHeap & lh)
00228   {
00229     Cast(fel).CalcShape (mip.IP(), FlatVector<> (fel.GetNDof(), &mat(0,0)));
00230   }
00231 
00232 
00233   template <typename MIP, class TVX, class TVY>
00234   static void Apply (const FiniteElement & fel, const MIP & mip,
00235                      const TVX & x, TVY & y,
00236                      LocalHeap & lh) 
00237   {
00238     y = Trans (Cast(fel).GetShape (mip.IP(), lh)) * x;
00239   }
00240 
00241   static void Apply (const FiniteElement & fel, const MappedIntegrationPoint<D,D> & mip,
00242                      const FlatVector<double> & x, FlatVector<double> & y,
00243                      LocalHeap & lh) 
00244   {
00245     y(0) = Cast(fel).Evaluate(mip.IP(), x);
00246   }
00247 
00248   template <class MIR>
00249   static void ApplyIR (const FiniteElement & fel, const MIR & mir,
00250                        const FlatVector<double> & x, FlatMatrix<double> & y,
00251                        LocalHeap & lh)
00252   {
00253     Cast(fel).Evaluate (mir.IR(), x, FlatVector<> (mir.Size(), &y(0,0)));
00254   }
00255 
00256 
00257 
00258   template <typename MIP, class TVX, class TVY>
00259   static void ApplyTrans (const FiniteElement & fel, const MIP & mip,
00260                           const TVX & x, TVY & y,
00261                           LocalHeap & lh) 
00262   {
00263     y = Cast(fel).GetShape (mip.IP(), lh) * x;
00264   }
00265 
00266 
00267 
00268   template <typename MIP, class TVX>
00269   static void Transform (const MIP & mip, TVX & x)
00270   {
00271     // do nothing
00272     ; 
00273   }
00274 
00275   template <typename MIP, class TVX>
00276   static void TransformT (const MIP & mip, TVX & x)
00277   {
00278     // do nothing
00279     ; 
00280   }
00281 
00282   /*
00283   template <class TVD, class TVY>
00284   static void ApplyGrid (const FiniteElement & fel, 
00285                          const IntegrationRuleTP<D> & ir,
00286                          const TVY & hv, 
00287                          TVD & dvecs, 
00288                          LocalHeap & lh)
00289   {
00290     FlatVector<double> dvec(dvecs.Size(), const_cast<double*> (&dvecs[0](0)));
00291     Cast(fel).EvaluateShapeGrid (ir, hv, dvec, lh);
00292   }
00293 
00294 
00295   template <class TVD, class TVY>
00296   static void ApplyTransGrid (const FiniteElement & fel, 
00297                               const IntegrationRuleTP<D> & ir,
00298                               const TVD & dvecs, 
00299                               TVY & hv, LocalHeap & lh)
00300   {
00301     FlatVector<double> dvec(dvecs.Size(), const_cast<double*> (&dvecs[0](0)));
00302     Cast(fel).EvaluateShapeGridTrans (ir, dvec, hv, lh);
00303   }
00304   */
00305 };
00306 
00307 
00308 
00310 template <int D, int SYSDIM>
00311 class DiffOpIdSys : public DiffOp<DiffOpIdSys<D,SYSDIM> >
00312 {
00313 public:
00314   enum { DIM = SYSDIM };
00315   enum { DIM_SPACE = D };
00316   enum { DIM_ELEMENT = D };
00317   enum { DIM_DMAT = SYSDIM };
00318   enum { DIFFORDER = 0 };
00319 
00320   template <typename FEL, typename MIP, typename MAT>
00321   static void GenerateMatrix (const FEL & fel, const MIP & mip,
00322                               MAT & mat, LocalHeap & lh)
00323   {
00324     const FlatVector<> shape = fel.GetShape (mip.IP(), lh);
00325   
00326 
00327     typedef typename MAT::TSCAL TSCAL;
00328     mat = TSCAL(0.); 
00329     for (int j = 0; j < shape.Height(); j++)
00330       for (int i = 0; i < SYSDIM; i++)
00331         mat(i, j*SYSDIM+i) = shape(j);
00332   }
00333 };
00334 
00335 
00336 
00337 
00338 
00340 template <int D, typename FEL = ScalarFiniteElement<D-1> >
00341 class DiffOpIdBoundary : public DiffOp<DiffOpIdBoundary<D, FEL> >
00342 {
00343 public:
00344   enum { DIM = 1 };
00345   enum { DIM_SPACE = D };
00346   enum { DIM_ELEMENT = D-1 };
00347   enum { DIM_DMAT = 1 };
00348   enum { DIFFORDER = 0 };
00349 
00350   template <typename AFEL, typename MIP, typename MAT>
00351   static void GenerateMatrix (const AFEL & fel, const MIP & mip,
00352                               MAT & mat, LocalHeap & lh)
00353   {
00354     mat.Row(0) = static_cast<const FEL&>(fel).GetShape(mip.IP(), lh);
00355   }
00356 
00357   template <typename AFEL, typename MIP, class TVX, class TVY>
00358   static void Apply (const AFEL & fel, const MIP & mip,
00359                      const TVX & x, TVY & y,
00360                      LocalHeap & lh) 
00361   {
00362     y = Trans (static_cast<const FEL&>(fel).GetShape (mip.IP(), lh)) * x;
00363   }
00364 
00365   static void Apply (const ScalarFiniteElement<D-1> & fel, const MappedIntegrationPoint<D-1,D> & mip,
00366                      const FlatVector<double> & x, FlatVector<double> & y,
00367                      LocalHeap & lh) 
00368   {
00369     y(0) = static_cast<const FEL&>(fel).Evaluate(mip.IP(), x);
00370   }
00371 
00372 
00373   template <typename AFEL, typename MIP, class TVX, class TVY>
00374   static void ApplyTrans (const AFEL & fel, const MIP & mip,
00375                           const TVX & x, TVY & y,
00376                           LocalHeap & lh) 
00377   {
00378     y = static_cast<const FEL&>(fel).GetShape (mip.IP(), lh) * x;
00379   }
00380 
00381 
00382 
00383   template <typename MIP, class TVX>
00384   static void Transform (const MIP & mip, TVX & x)
00385   {
00386     // do nothing
00387     ; 
00388   }
00389 
00390   template <typename MIP, class TVX>
00391   static void TransformT (const MIP & mip, TVX & x)
00392   {
00393     // do nothing
00394     ; 
00395   }
00396 
00397   /*
00398   template <typename AFEL, class TVD, class TVY>
00399   static void ApplyGrid (const AFEL & fel, 
00400                          const IntegrationRuleTP<D-1> & ir,
00401                          const TVY & hv, 
00402                          TVD & dvecs, 
00403                          LocalHeap & lh)
00404   {
00405     FlatVector<double> dvec(dvecs.Size(), const_cast<double*> (&dvecs[0](0)));
00406     static_cast<const FEL&>(fel).EvaluateShapeGrid (ir, hv, dvec, lh);
00407   }
00408 
00409 
00410   template <typename AFEL, class TVD, class TVY>
00411   static void ApplyTransGrid (const AFEL & fel, 
00412                               const IntegrationRuleTP<D-1> & ir,
00413                               const TVD & dvecs, 
00414                               TVY & hv, LocalHeap & lh)
00415   {
00416     FlatVector<double> dvec(dvecs.Size(), const_cast<double*> (&dvecs[0](0)));
00417     static_cast<const FEL&>(fel).EvaluateShapeGridTrans (ir, dvec, hv, lh);
00418   }
00419   */
00420 };
00421 
00422 
00423 
00425 template <int D, int SYSDIM>
00426 class DiffOpIdBoundarySys : public DiffOp<DiffOpIdBoundarySys<D,SYSDIM> >
00427 {
00428 public:
00429   enum { DIM = SYSDIM };
00430   enum { DIM_SPACE = D };
00431   enum { DIM_ELEMENT = D-1 };
00432   enum { DIM_DMAT = SYSDIM };
00433   enum { DIFFORDER = 0 };
00434 
00435   template <typename FEL, typename MIP, typename MAT>
00436   static void GenerateMatrix (const FEL & fel, const MIP & mip,
00437                               MAT & mat, LocalHeap & lh)
00438   {
00439     mat = 0.; 
00440     const FlatVector<> shape = fel.GetShape (mip.IP(), lh);
00441     for (int j = 0; j < shape.Height(); j++)
00442       for (int i = 0; i < SYSDIM; i++)
00443         mat(i, j*SYSDIM+i) = shape(j);
00444   }
00445 };
00446 
00447 
00448 
00449 
00450 
00451 
00452 
00453 
00454 
00455 
00456 
00457 
00458 
00459 
00460 
00461 
00462 
00463 
00464 
00465 
00466 
00467 
00468 
00469 
00471 template <int DIM>
00472 class DiagDMat : public DMatOp<DiagDMat<DIM> >
00473 {
00474   CoefficientFunction * coef;
00475 public:
00476   // typedef SCAL TSCAL;
00477   enum { DIM_DMAT = DIM };
00478   DiagDMat (CoefficientFunction * acoef) : coef(acoef) { ; }
00479 
00480   DiagDMat (Array<CoefficientFunction*> & acoefs) : coef(acoefs[0]) { ; }
00481 
00482   template <typename FEL, typename MIP, typename MAT>
00483   void GenerateMatrix (const FEL & fel, const MIP & mip,
00484                        MAT & mat, LocalHeap & lh) const
00485   {
00486     typedef typename MAT::TSCAL TRESULT;
00487     TRESULT val = coef -> T_Evaluate<TRESULT> (mip);
00488     /*
00489     mat = TRESULT(0);
00490     for (int i = 0; i < DIM; i++)
00491       mat(i, i) = val;
00492     */
00493     mat = val * Id<DIM>();
00494   }  
00495 
00496   template <typename FEL, typename MIR, typename MAT>
00497   void GenerateMatrixIR (const FEL & fel, const MIR & mir,
00498                          const FlatArray<MAT> & mats, LocalHeap & lh) const
00499   {
00500     typedef typename MAT::TSCAL TRESULT;
00501     FlatMatrix<TRESULT> vals(mir.IR().GetNIP(), 1, lh);
00502     coef -> Evaluate (mir, vals);
00503     
00504     for (int j = 0; j < mir.IR().GetNIP(); j++)
00505       {
00506         mats[j] = vals(j,0) * Id<DIM>();
00507         /*
00508         mats[j] = TRESULT(0.0);
00509         for (int i = 0; i < DIM; i++)
00510           mats[j](i, i) = vals(j,0);
00511         */
00512       }
00513   }  
00514 
00515   template <typename FEL, class VECX, class VECY>
00516   void Apply (const FEL & fel, const BaseMappedIntegrationPoint & mip,
00517               const VECX & x, VECY & y, LocalHeap & lh) const
00518   {
00519     typedef typename VECY::TSCAL TRESULT;
00520     TRESULT val = coef -> T_Evaluate<TRESULT> (mip);
00521     for (int i = 0; i < DIM; i++)
00522       y(i) = val * x(i);
00523   }
00524 
00525   template <typename FEL, class VECX>
00526   void Apply1 (const FEL & fel, const BaseMappedIntegrationPoint & mip,
00527                const VECX & x, LocalHeap & lh) const
00528   {
00529     typedef typename VECX::TSCAL TSCAL;
00530     x *= coef -> T_Evaluate<TSCAL> (mip);
00531   }
00532 
00533   template <typename FEL, typename MIR, typename TVX>
00534   void ApplyIR (const FEL & fel, const MIR & mir,
00535                 TVX & x, LocalHeap & lh) const
00536   {
00537     typedef typename TVX::TSCAL TSCAL;
00538     FlatMatrix<TSCAL> values(mir.Size(), 1, lh);
00539     coef -> Evaluate (mir, values);
00540     for (int i = 0; i < mir.Size(); i++)
00541       x.Row(i) *=  values(i, 0);
00542   }
00543 };
00544 
00545 
00546 
00547 
00548 
00549 
00550 
00551 
00553 template <int N> 
00554 class OrthoDMat
00555 {
00556 };
00557 
00558 template <> class OrthoDMat<1> : public DMatOp<OrthoDMat<1> >
00559 {
00560   CoefficientFunction * coef;
00561 public:
00562   enum { DIM_DMAT = 1 };
00563   OrthoDMat (CoefficientFunction * acoef) : coef(acoef) { ; }
00564 
00565   template <typename FEL, typename MIP, typename MAT>
00566   void GenerateMatrix (const FEL & fel, const MIP & mip,
00567                        MAT & mat, LocalHeap & lh) const
00568   {
00569     mat(0,0) = Evaluate (*coef, mip);
00570   }  
00571 
00572   template <typename FEL, typename MIP, class VECX, class VECY>
00573   void Apply (const FEL & fel, const MIP & mip,
00574               const VECX & x, VECY & y, LocalHeap & lh) const
00575   {
00576     y(0) = Evaluate (*coef, mip) * x(0);
00577   }
00578 
00579   template <typename FEL, typename MIP, class VECX, class VECY>
00580   void ApplyTrans (const FEL & fel, const MIP & mip,
00581                    const VECX & x, VECY & y, LocalHeap & lh) const
00582   {
00583     y(0) = Evaluate (*coef, mip) * x(0);
00584   }
00585 
00586 };
00587 
00588 
00589 template <> class OrthoDMat<2>: public DMatOp<OrthoDMat<2> >
00590 {
00591   CoefficientFunction * coef1;
00592   CoefficientFunction * coef2;
00593 public:
00594   enum { DIM_DMAT = 2 };
00595   OrthoDMat (CoefficientFunction * acoef1,
00596         CoefficientFunction * acoef2)
00597     : coef1(acoef1), coef2(acoef2) { ; }
00598   OrthoDMat (CoefficientFunction * acoef1,
00599         CoefficientFunction * acoef2,
00600         CoefficientFunction * acoef3)
00601     : coef1(acoef1), coef2(acoef2) { ; }
00602   
00603   template <typename FEL, typename MIP, typename MAT>
00604   void GenerateMatrix (const FEL & fel, const MIP & mip,
00605                        MAT & mat, LocalHeap & lh) const
00606   {
00607     mat = 0;
00608     mat(0,0) = Evaluate (*coef1, mip);
00609     mat(1,1) = Evaluate (*coef2, mip);
00610   }  
00611 
00612   template <typename FEL, typename MIP>
00613   void GetEigensystem (const FEL & fel, const MIP & mip, 
00614                   Array<double> & eigenvalues,
00615                   Array<double> & eigenvectors,
00616                   LocalHeap & lh) const
00617   {
00618     eigenvalues[0] = Evaluate (*coef1, mip);
00619     eigenvalues[1] = Evaluate (*coef2, mip);
00620     eigenvectors[0] = eigenvectors[3] = 1.;
00621     eigenvectors[1] = eigenvectors[2] = 0.;
00622   }
00623 
00624   template <typename FEL, typename MIP, class VECX, class VECY>
00625   void Apply (const FEL & fel, const MIP & mip,
00626               const VECX & x, VECY & y, LocalHeap & lh) const
00627   {
00628     y(0) = Evaluate (*coef1, mip) * x(0);
00629     y(1) = Evaluate (*coef2, mip) * x(1);
00630   }
00631 
00632   template <typename FEL, typename MIP, class VECX, class VECY>
00633   void ApplyTrans (const FEL & fel, const MIP & mip,
00634                    const VECX & x, VECY & y, LocalHeap & lh) const
00635   {
00636     y(0) = Evaluate (*coef1, mip) * x(0);
00637     y(1) = Evaluate (*coef2, mip) * x(1);
00638   }
00639 };
00640 
00641 template <> class OrthoDMat<3> : public DMatOp<OrthoDMat<3> >
00642 {
00643   CoefficientFunction * coef1;
00644   CoefficientFunction * coef2;
00645   CoefficientFunction * coef3;
00646 public:
00647   enum { DIM_DMAT = 3 };
00648   OrthoDMat (CoefficientFunction * acoef1,
00649              CoefficientFunction * acoef2)
00650     : coef1(acoef1), coef2(acoef2), coef3(acoef2) { ; }
00651   OrthoDMat (CoefficientFunction * acoef1,
00652              CoefficientFunction * acoef2,
00653              CoefficientFunction * acoef3)
00654     : coef1(acoef1), coef2(acoef2), coef3(acoef3) { ; }
00655   
00656   template <typename FEL, typename MIP, typename MAT>
00657   void GenerateMatrix (const FEL & fel, const MIP & mip,
00658                        MAT & mat, LocalHeap & lh) const
00659   {
00660     mat = 0;
00661     mat(0,0) = Evaluate (*coef1, mip);
00662     mat(1,1) = Evaluate (*coef2, mip);
00663     mat(2,2) = Evaluate (*coef3, mip);
00664   }  
00665 
00666   
00667   template <typename FEL, typename MIP>
00668   void GetEigensystem (const FEL & fel, const MIP & mip, 
00669                   Array<double> & eigenvalues,
00670                   Array<double> & eigenvectors,
00671                   LocalHeap & lh) const
00672   {
00673     
00674     eigenvalues[0] = Evaluate(*coef1,mip);
00675     eigenvalues[1] = Evaluate(*coef2,mip);
00676     eigenvalues[2] = Evaluate(*coef3,mip);
00677 
00678     eigenvectors = 0.;
00679     eigenvectors[0] = eigenvectors[4] = eigenvectors[8] = 1.;
00680   }
00681 
00682 
00683   template <typename FEL, typename MIP, class VECX, class VECY>
00684   void Apply (const FEL & fel, const MIP & mip,
00685               const VECX & x, VECY & y, LocalHeap & lh) const
00686   {
00687     y(0) = Evaluate (*coef1, mip) * x(0);
00688     y(1) = Evaluate (*coef2, mip) * x(1);
00689     y(2) = Evaluate (*coef3, mip) * x(2);
00690   }
00691 
00692   template <typename FEL, typename MIP, class VECX, class VECY>
00693   void ApplyTrans (const FEL & fel, const MIP & mip,
00694                    const VECX & x, VECY & y, LocalHeap & lh) const
00695   {
00696     y(0) = Evaluate (*coef1, mip) * x(0);
00697     y(1) = Evaluate (*coef2, mip) * x(1);
00698     y(2) = Evaluate (*coef3, mip) * x(2);
00699   }
00700 
00701   void SetCoefficientFunctions( CoefficientFunction * acoef1,
00702                                 CoefficientFunction * acoef2,
00703                                 CoefficientFunction * acoef3 )
00704   {
00705     // NOTE: alte coefficient-functions werden nicht geloescht!
00706     coef1 = acoef1;
00707     coef2 = acoef2;
00708     coef3 = acoef3;
00709   }
00710 };
00711 
00712 
00713 
00714 
00715 
00716 
00717 
00718 
00719 
00721 template <int N> 
00722 class SymDMat : public DMatOp<SymDMat<N> >
00723 {
00724 };
00725 
00726 template <> class SymDMat<1> : public DMatOp<SymDMat<1> >
00727 {
00728   CoefficientFunction * coef;
00729 public:
00730   enum { DIM_DMAT = 1 };
00731 
00732   SymDMat (CoefficientFunction * acoef) : coef(acoef) { ; }
00733 
00734   template <typename FEL, typename MIP, typename MAT>
00735   void GenerateMatrix (const FEL & fel, const MIP & mip,
00736                        MAT & mat, LocalHeap & lh) const
00737   {
00738     mat(0,0) = Evaluate (*coef, mip);
00739   }  
00740 };
00741 
00742 
00743 template <> class SymDMat<2> : public DMatOp<SymDMat<2> >
00744 {
00745   CoefficientFunction * coef00;
00746   CoefficientFunction * coef01;
00747   CoefficientFunction * coef11;
00748 public:
00749   enum { DIM_DMAT = 2 };
00750 
00751   SymDMat (CoefficientFunction * acoef00,
00752            CoefficientFunction * acoef01,
00753            CoefficientFunction * acoef11)
00754     : coef00(acoef00), coef01(acoef01), coef11(acoef11) { ; }
00755   
00756   template <typename FEL, typename MIP, typename MAT>
00757   void GenerateMatrix (const FEL & fel, const MIP & mip,
00758                        MAT & mat, LocalHeap & lh) const
00759   {
00760     mat = 0;
00761     mat(0,0) = Evaluate (*coef00, mip);
00762     mat(0,1) = mat(1,0) = Evaluate (*coef01, mip);
00763     mat(1,1) = Evaluate (*coef11, mip);
00764   }  
00765 };
00766 
00767 template <> class SymDMat<3> : public DMatOp<SymDMat<3> >
00768 {
00769   CoefficientFunction * coef00;
00770   CoefficientFunction * coef10;
00771   CoefficientFunction * coef11;
00772   CoefficientFunction * coef20;
00773   CoefficientFunction * coef21;
00774   CoefficientFunction * coef22;
00775 public:
00776   enum { DIM_DMAT = 3 };
00777 
00778   SymDMat (CoefficientFunction * acoef00,
00779            CoefficientFunction * acoef10,
00780            CoefficientFunction * acoef11,
00781            CoefficientFunction * acoef20,
00782            CoefficientFunction * acoef21,
00783            CoefficientFunction * acoef22)
00784     : coef00(acoef00), coef10(acoef10), coef11(acoef11),
00785       coef20(acoef20), coef21(acoef21), coef22(acoef22) { ; }
00786   
00787   template <typename FEL, typename MIP, typename MAT>
00788   void GenerateMatrix (const FEL & fel, const MIP & mip,
00789                        MAT & mat, LocalHeap & lh) const
00790   {
00791     mat = 0;
00792     mat(0,0) = Evaluate (*coef00, mip);
00793     mat(1,0) = mat(0,1) = Evaluate (*coef10, mip);
00794     mat(1,1) = Evaluate (*coef11, mip);
00795     mat(2,0) = mat(0,2) = Evaluate (*coef20, mip);
00796     mat(2,1) = mat(1,2) = Evaluate (*coef21, mip);
00797     mat(2,2) = Evaluate (*coef22, mip);
00798   }  
00799 };
00800 
00801 
00802 
00803 
00804 
00805 
00806 
00807 
00808 
00810 template <int DIM>
00811 class NormalDMat : public DMatOp<NormalDMat<DIM> >
00812 {
00813   CoefficientFunction * coef;
00814 public:
00815   enum { DIM_DMAT = DIM };
00816   NormalDMat (CoefficientFunction * acoef) : coef(acoef) { ; }
00817 
00818   template <typename FEL, typename MIP, typename MAT>
00819   void GenerateMatrix (const FEL & fel, const MIP & mip,
00820                        MAT & mat, LocalHeap & lh) const
00821   {
00822     mat = 0;
00823     double val = Evaluate (*coef, mip);
00824     Vec<DIM> nv = mip.GetNV();
00825     for (int i = 0; i < DIM; i++)
00826       for (int j = 0; j < DIM; j++)
00827         mat(i, j) = val * nv(i) * nv(j);
00828   }  
00829 };
00830 
00831 
00832 
00833 
00834 
00835 
00836 
00837 template <int N>
00838 class NGS_DLL_HEADER DVec  // : public DVecBase<N,T>
00839 {
00840   CoefficientFunction * coefs[N];
00841   bool vectorial;
00842 public:
00843   DVec (Array<CoefficientFunction*> & acoeffs)
00844   {
00845     vectorial = (N > 1) && (N == acoeffs[0]->Dimension());
00846 
00847     if (vectorial)
00848       coefs[0] = acoeffs[0];
00849     else
00850       for (int i = 0; i < N; i++)
00851         coefs[i] = acoeffs[i];
00852   }
00853   
00854   DVec (CoefficientFunction * acoef)
00855   { 
00856     vectorial = (N > 1) && (N == acoef->Dimension());
00857     coefs[0] = acoef;
00858   }
00859 
00860   DVec (CoefficientFunction * acoef1,
00861         CoefficientFunction * acoef2,
00862         CoefficientFunction * acoef3 = NULL,
00863         CoefficientFunction * acoef4 = NULL,
00864         CoefficientFunction * acoef5 = NULL,
00865         CoefficientFunction * acoef6 = NULL)
00866   { 
00867     vectorial = false;
00868 
00869     coefs[0] = acoef1;
00870     coefs[1] = acoef2;
00871     if (N >= 3) coefs[2] = acoef3;
00872     if (N >= 4) coefs[3] = acoef4;
00873     if (N >= 5) coefs[4] = acoef5;
00874     if (N >= 6) coefs[5] = acoef6;
00875   }
00876     
00877 
00878 
00879 
00880   template <typename FEL, typename MIP, typename VEC>
00881   void GenerateVector (const FEL & fel, const MIP & mip,
00882                        VEC & vec, LocalHeap & lh) const
00883   {
00884     typedef typename VEC::TSCAL TSCAL;
00885 
00886     if (vectorial)
00887       {
00888         coefs[0] -> Evaluate (mip, vec);
00889       }
00890     else
00891       for (int i = 0; i < N; i++)
00892         {
00893           CoefficientFunction * hp = coefs[i];
00894           vec(i) = hp -> T_Evaluate<TSCAL> (mip);
00895 
00896           // gcc 4.6 complains, why ???? (JS)
00897           // vec(i) = coefs[i] -> T_Evaluate<TSCAL> (mip);  
00898         }
00899   }  
00900 
00901 
00902 
00903   template <typename FEL, typename MIR, typename VEC>
00904   void GenerateVectorIR (const FEL & fel, const MIR & mir,
00905                          VEC & vecs, LocalHeap & lh) const
00906   {
00907     typedef typename VEC::TSCAL TSCAL;
00908 
00909     if (vectorial || N == 1)
00910       {
00911         coefs[0] -> Evaluate (mir, vecs);
00912       }
00913     else
00914       for (int j = 0; j < mir.Size(); j++)
00915         for (int i = 0; i < N; i++)
00916           {
00917             CoefficientFunction * hp = coefs[i];
00918             vecs(j,i) = hp -> T_Evaluate<TSCAL> (mir[j]);
00919             
00920             // gcc 4.6 complains, why ???? (JS)
00921             // vec(i) = (coefs[i]) -> T_Evaluate<TSCAL> (mip);  
00922           }
00923   }  
00924 };
00925 
00926 
00927 
00928 
00929 
00930 template <int N, typename T = double>  
00931 class DVecN
00932 {
00933   CoefficientFunction * coef;
00934 public:
00935   typedef T TSCAL;
00936   DVecN (CoefficientFunction * acoef)
00937     : coef(acoef) { ; }
00938   
00939   template <typename FEL, typename MIP, typename VEC>
00940   void GenerateVector (const FEL & fel, const MIP & mip,
00941                        VEC & vec, LocalHeap & lh) const
00942   {
00943     Vec<N> hv;
00944     coef -> Evaluate (mip, hv);
00945     for (int i = 0; i < N; i++)
00946       vec(i) = hv(i);
00947   }  
00948 
00949   template <typename FEL, typename MIR, typename VEC>
00950   void GenerateVectorIR (const FEL & fel, const MIR & mir,
00951                          VEC & vecs, LocalHeap & lh) const
00952   {
00953     for (int i = 0; i < mir.Size(); i++)
00954       GenerateVector (fel, mir[i], vecs.Row(i), lh);
00955   }
00956 };
00957 
00958 
00959 
00960 template <int N>
00961 class TVec
00962 {
00963   CoefficientFunction * coef;
00964 
00965 public:
00966   typedef double TSCAL;
00967   TVec (CoefficientFunction * acoef) : coef(acoef) {;}
00968 
00969   
00970 
00971   template <typename FEL, typename MIP, typename VEC>
00972   void GenerateVector (const FEL & fel, const MIP & mip,
00973                        VEC & vec, LocalHeap & lh) const
00974   {
00975     vec = 0.0;
00976 
00977     typedef typename VEC::TSCAL TSCAL;
00978     
00979     TSCAL length = 0.;
00980     for(int i=0; i<N; i++)
00981       {
00982         //vec(i) = mip.GetJacobian()(i,0);
00983         vec(i) = mip.GetTV()(i);
00984         length += vec(i)*vec(i);
00985       }
00986     //(*testout) << "point " << mip.GetPoint() << " tv " << vec;
00987     vec *= Evaluate (*coef, mip)/sqrt(length);
00988     //(*testout) << " retval " << vec << endl;
00989   }
00990 
00991   template <typename FEL, typename MIR, typename VEC>
00992   void GenerateVectorIR (const FEL & fel, const MIR & mir,
00993                          VEC & vecs, LocalHeap & lh) const
00994   {
00995     for (int i = 0; i < mir.Size(); i++)
00996       GenerateVector (fel, mir[i], vecs.Row(i), lh);
00997   }
00998 
00999 };
01000 
01001 
01002 
01003 
01004 
01005 
01006 
01007 
01008 
01009 
01010 
01012 template <int DIM>
01013 class RotSymLaplaceDMat : public DMatOp<RotSymLaplaceDMat<DIM> >
01014 {
01015   CoefficientFunction * coef;
01016 public:
01017   enum { DIM_DMAT = DIM };
01018 
01019   RotSymLaplaceDMat (CoefficientFunction * acoef) : coef(acoef) { ; }
01020 
01021   template <typename FEL, typename MIP, typename MAT>
01022   void GenerateMatrix (const FEL & fel, const MIP & mip,
01023                        MAT & mat, LocalHeap & lh) const
01024   {
01025     mat = 0;
01026     const double r = mip.GetPoint()(0);
01027     double val = r*Evaluate (*coef, mip);
01028     for (int i = 0; i < DIM; i++)
01029       mat(i, i) = val;
01030   }  
01031   
01032 
01033   template <typename FEL, typename MIP, class VECX, class VECY>
01034   void Apply (const FEL & fel, const MIP & mip,
01035               const VECX & x, VECY & y, LocalHeap & lh) const
01036   {
01037     const double r = mip.GetPoint()(0);
01038     double val = r*Evaluate (*coef, mip);
01039     y = val * x;
01040   }
01041 };
01042 
01043 
01044 
01045 
01046 
01047 
01049 template <int D, typename FEL = ScalarFiniteElement<D> >
01050 class DiffOpNormal : public DiffOp<DiffOpNormal<D, FEL> >
01051 {
01052 public:
01053   enum { DIM = D };
01054   enum { DIM_SPACE = D };
01055   enum { DIM_ELEMENT = D-1 };
01056   enum { DIM_DMAT = 1 };
01057   enum { DIFFORDER = 0 };
01058 
01059   template <typename MIP, typename MAT>
01060   static void GenerateMatrix (const FiniteElement & fel, const MIP & mip,
01061                               MAT & mat, LocalHeap & lh)
01062   {
01063     FlatVector<> shape = static_cast<const FEL&> (fel).GetShape (mip.IP(), lh);
01064     Vec<D> nv = mip.GetNV();
01065     //Vec<D> p = mip.GetPoint();
01066     for (int j = 0; j < shape.Size(); j++)
01067       for (int i = 0; i < D; i++)
01068         mat(0, j*D+i) =  shape(j) * nv(i);
01069 
01070     //(*testout) << "mip = " << p << ", nv = " << nv << endl;
01071     //p(0) = 0.0;
01072     //p /= L2Norm(p);
01073     //nv /= L2Norm(nv);
01074     //(*testout) << "normalized, mip = " << p << ", nv = " << nv << endl;
01075     //(*testout) << "mat = " << mat << endl;
01076   }
01077 };
01078 
01079 
01080 
01081 
01082 
01084 template <int D>
01085 class NormalRobinIntegrator 
01086   : public T_BDBIntegrator<DiffOpNormal<D>, DiagDMat<1>, ScalarFiniteElement<D-1> >
01087 {
01088 public:
01089   NormalRobinIntegrator (CoefficientFunction * coeff)
01090     : T_BDBIntegrator<DiffOpNormal<D>, DiagDMat<1>, ScalarFiniteElement<D-1> > (DiagDMat<1> (coeff))
01091   { ; }
01092 
01093 
01094   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01095   {
01096     return new NormalRobinIntegrator (coeffs[0]);
01097   }
01098 
01099   virtual bool BoundaryForm () const { return 1; }
01100   virtual string Name () const { return "NormalRobin"; }
01101 };
01102 
01103 
01104 
01105 
01106 
01107 
01108 // ********************************* Scalar integrators: ********************
01109 
01111 template <int D, typename FEL = ScalarFiniteElement<D> >
01112 class NGS_DLL_HEADER LaplaceIntegrator 
01113   : public T_BDBIntegrator<DiffOpGradient<D>, DiagDMat<D>, FEL>
01114 {
01115 public:
01117   LaplaceIntegrator (CoefficientFunction * coeff);
01118   LaplaceIntegrator (Array<CoefficientFunction*> & coeffs);
01119   virtual ~LaplaceIntegrator ();
01120   virtual string Name () const { return "Laplace"; }
01121 };
01122 
01123 
01124 
01126 template <int D, typename FEL = ScalarFiniteElement<D-1> >
01127 class LaplaceBoundaryIntegrator 
01128   : public T_BDBIntegrator<DiffOpGradientBoundary<D>, DiagDMat<D>, FEL>
01129 {
01130 public:
01132   LaplaceBoundaryIntegrator (CoefficientFunction * coeff);
01133 
01134   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01135   {
01136     return new LaplaceBoundaryIntegrator (coeffs[0]);
01137   }
01138 
01140   virtual string Name () const { return "Laplace-Boundary"; }
01141 };
01142 
01143 
01144 
01146 template <int D, typename FEL = ScalarFiniteElement<D> >
01147 class RotSymLaplaceIntegrator 
01148   : public T_BDBIntegrator<DiffOpGradient<D>, RotSymLaplaceDMat<D>, FEL>
01149 {
01150 public:
01152   RotSymLaplaceIntegrator (CoefficientFunction * coeff);
01153 
01154   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01155   {
01156     return new RotSymLaplaceIntegrator (coeffs[0]);
01157   }
01158   
01160   virtual string Name () const { return "RotSymLaplace"; }
01161 };
01162 
01163 
01165 template <int D, typename FEL = ScalarFiniteElement<D> >
01166 class OrthoLaplaceIntegrator 
01167   : public T_BDBIntegrator<DiffOpGradient<D>, OrthoDMat<D>, FEL>
01168 {
01169 public:
01171   OrthoLaplaceIntegrator (CoefficientFunction * coeff1, CoefficientFunction * coeff2)
01172     : T_BDBIntegrator<DiffOpGradient<D>, OrthoDMat<D>, FEL> (OrthoDMat<D> (coeff1, coeff2))
01173   { ; }
01174   OrthoLaplaceIntegrator (CoefficientFunction * coeff1, CoefficientFunction * coeff2, CoefficientFunction * coeff3)
01175     : T_BDBIntegrator<DiffOpGradient<D>, OrthoDMat<D>, FEL> (OrthoDMat<D> (coeff1, coeff2, coeff3))
01176   { ; }
01177   
01178   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01179   {
01180     if(coeffs.Size() == 2)
01181       return new OrthoLaplaceIntegrator (coeffs[0], coeffs[1]);
01182     else
01183       return new OrthoLaplaceIntegrator (coeffs[0], coeffs[1], coeffs[2]);
01184   }
01185 
01186 
01188   virtual string Name () const { return "OrthoLaplace"; }
01189 };
01190 
01191 
01192 
01194 template <int D, typename FEL = ScalarFiniteElement<D> >
01195 class NGS_DLL_HEADER MassIntegrator 
01196   : public T_BDBIntegrator<DiffOpId<D>, DiagDMat<1>, FEL >
01197 {
01198 public:
01200   MassIntegrator (CoefficientFunction * coeff);
01202   MassIntegrator (Array<CoefficientFunction*> & coeffs);
01204   virtual ~MassIntegrator ();
01206   virtual string Name () const { return "Mass"; }
01207 };
01208 
01209 
01210 
01211 
01212 
01213 
01215   template <int D, typename FEL = ScalarFiniteElement<D-1> >
01216 class NGS_DLL_HEADER RobinIntegrator 
01217   : public T_BDBIntegrator<DiffOpIdBoundary<D>, DiagDMat<1>, FEL>
01218 {
01219   typedef T_BDBIntegrator<DiffOpIdBoundary<D>, DiagDMat<1>, FEL> BASE;
01220 public:
01222   RobinIntegrator (CoefficientFunction * coeff);
01224   RobinIntegrator (Array<CoefficientFunction*> & coeffs);
01226   virtual ~RobinIntegrator ();
01228   virtual bool BoundaryForm () const { return 1; }
01230   virtual string Name () const { return "Robin"; }
01231 };
01232 
01233 
01234 
01235 
01236 
01237 /*
01238 template <int D>
01239 class NormalRobinIntegrator 
01240   : public T_BDBIntegrator<DiffOpIdBoundary<D,D>, NormalDMat<D>, ScalarFiniteElement<D> >
01241 {
01242 public:
01243   NormalRobinIntegrator (CoefficientFunction * coeff)
01244     : T_BDBIntegrator<DiffOpIdBoundary<D,D>, NormalDMat<D>, ScalarFiniteElement<D> > (NormalDMat<D> (coeff))
01245   { ; }
01246 
01247   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01248   {
01249     return new NormalRobinIntegrator (coeffs[0]);
01250   }
01251 
01252   virtual bool BoundaryForm () const { return 1; }
01253   virtual string Name () const { return "NormalRobin"; }
01254 };
01255 */
01256 
01257 
01259 template <int D> 
01260 class DiffOpDiv : public DiffOp<DiffOpDiv<D> >
01261 {
01262 public:
01263   enum { DIM = D };
01264   enum { DIM_SPACE = D };
01265   enum { DIM_ELEMENT = D };
01266   enum { DIM_DMAT = 1 };
01267   enum { DIFFORDER = 1 };
01268 
01269   template <typename FEL, typename MIP, typename MAT>
01270   static void GenerateMatrix (const FEL & fel, const MIP & mip,
01271                               MAT & mat, LocalHeap & lh)
01272   {
01273     int nd = fel.GetNDof();
01274 
01275     FlatMatrix<> grad (D, nd, lh);
01276     grad = Trans (mip.GetJacobianInverse ()) * 
01277       Trans (fel.GetDShape(mip.IP(), lh));
01278     
01279     mat = 0;
01280     for (int i = 0; i < nd; i++)
01281       for (int j = 0; j < DIM; j++)
01282         mat(0, DIM*i+j) = grad(j, i);
01283   }
01284 };
01285 
01286 
01287 template <int D, typename FEL = ScalarFiniteElement<D> >
01288 class DivDivIntegrator 
01289   : public T_BDBIntegrator<DiffOpDiv<D>, DiagDMat<1>, FEL>
01290 {
01291 public:
01293   DivDivIntegrator (CoefficientFunction * coeff)
01294     : T_BDBIntegrator<DiffOpDiv<D>, DiagDMat<1>, FEL> (DiagDMat<1> (coeff))
01295   {
01296     ;
01297   }
01298   
01299   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01300   {
01301     return new DivDivIntegrator (coeffs[0]);
01302   }
01303 
01305   virtual string Name () const { return "DivDiv"; }
01306 };
01307 
01308 
01309 
01310 
01311 
01312 
01313 
01314 
01315 
01316 
01317 
01318 
01320 class DiffOpCurl : public DiffOp<DiffOpCurl>
01321 {
01322 public:
01323   enum { DIM = 2 };
01324   enum { DIM_SPACE = 2 };
01325   enum { DIM_ELEMENT = 2 };
01326   enum { DIM_DMAT = 1 };
01327   enum { DIFFORDER = 1 };
01328 
01329   template <typename FEL, typename MIP, typename MAT>
01330   static void GenerateMatrix (const FEL & fel, const MIP & mip,
01331                               MAT & mat, LocalHeap & lh)
01332   {
01333     int nd = fel.GetNDof();
01334 
01335     FlatMatrix<> grad (2, nd, lh);
01336     grad = Trans (mip.GetJacobianInverse ()) * 
01337       Trans (fel.GetDShape(mip.IP(), lh));
01338     
01339     mat = 0;
01340     for (int i = 0; i < nd; i++)
01341       {
01342         mat(0, DIM*i  ) = grad(1, i);
01343         mat(0, DIM*i+1) = -grad(0, i);
01344       }
01345   }
01346 };
01347 
01348 
01349 template <typename FEL = ScalarFiniteElement<2> >
01350 class CurlCurlIntegrator 
01351   : public T_BDBIntegrator<DiffOpCurl, DiagDMat<1>, FEL>
01352 {
01353 public:
01355   CurlCurlIntegrator (CoefficientFunction * coeff)
01356     : T_BDBIntegrator<DiffOpCurl, DiagDMat<1>, FEL> (DiagDMat<1> (coeff))
01357   {
01358     ;
01359   }
01360   
01361   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01362   {
01363     return new CurlCurlIntegrator (coeffs[0]);
01364   }
01365 
01367   virtual string Name () const { return "CurlCurl"; }
01368 };
01369 
01370 
01371 
01372 
01374 class DiffOpCurl3d : public DiffOp<DiffOpCurl3d>
01375 {
01376 public:
01377   enum { DIM = 3 };
01378   enum { DIM_SPACE = 3 };
01379   enum { DIM_ELEMENT = 3 };
01380   enum { DIM_DMAT = 3 };
01381   enum { DIFFORDER = 1 };
01382 
01383   template <typename FEL, typename MIP, typename MAT>
01384   static void GenerateMatrix (const FEL & fel, const MIP & mip,
01385                               MAT & mat, LocalHeap & lh)
01386   {
01387     int nd = fel.GetNDof();
01388 
01389     FlatMatrix<> grad (3, nd, lh);
01390     grad = Trans (mip.GetJacobianInverse ()) * 
01391       Trans (fel.GetDShape(mip.IP(), lh));
01392     
01393     mat = 0;
01394     for (int i = 0; i < nd; i++)
01395       {
01396         mat(0, DIM*i+2) =  grad(1, i);
01397         mat(0, DIM*i+1) = -grad(2, i);
01398         mat(1, DIM*i+0) =  grad(2, i);
01399         mat(1, DIM*i+2) = -grad(0, i);
01400         mat(2, DIM*i+1) =  grad(0, i);
01401         mat(2, DIM*i+0) = -grad(1, i);
01402       }
01403   }
01404 };
01405 
01406 
01407 template <typename FEL = ScalarFiniteElement<3> >
01408 class CurlCurl3dIntegrator 
01409   : public T_BDBIntegrator<DiffOpCurl3d, DiagDMat<3>, FEL>
01410 {
01411 public:
01413   CurlCurl3dIntegrator (CoefficientFunction * coeff)
01414     : T_BDBIntegrator<DiffOpCurl3d, DiagDMat<3>, FEL> (DiagDMat<3> (coeff))
01415   {
01416     ;
01417   }
01418   
01419   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01420   {
01421     return new CurlCurl3dIntegrator (coeffs[0]);
01422   }
01423 
01425   virtual string Name () const { return "CurlCurl3d"; }
01426 };
01427 
01428 
01429 
01430 
01431 
01432 
01433 
01434 
01435 
01436 
01437 
01438 
01439 
01440 
01441 
01442 
01443 
01444 
01445 
01446 
01447 
01448 
01449 
01450 
01451 
01452 /* ************************** Linearform Integrators ************************* */
01453 
01455 template <int D, typename FEL = ScalarFiniteElement<D>  >
01456 class NGS_DLL_HEADER SourceIntegrator 
01457   : public T_BIntegrator<DiffOpId<D>, DVec<1>, FEL>
01458 {
01459 public:
01460   SourceIntegrator (CoefficientFunction * coeff);
01461   SourceIntegrator (Array<CoefficientFunction*> & coeffs);
01462   virtual ~SourceIntegrator ();
01463   virtual string Name () const { return "Source"; }
01464 };
01465 
01466 
01468 template <int D, typename FEL = ScalarFiniteElement<D-1> >
01469 class NGS_DLL_HEADER NeumannIntegrator 
01470   : public T_BIntegrator<DiffOpIdBoundary<D>, DVec<1>, FEL>
01471 {
01472 public:
01474   NeumannIntegrator (CoefficientFunction * coeff);
01475   NeumannIntegrator (Array<CoefficientFunction*> & coeffs);
01476   virtual ~NeumannIntegrator ();
01478   virtual bool BoundaryForm () const { return 1; }
01480   virtual string Name () const { return "Neumann"; }
01481 };
01482 
01483 
01484 
01485 
01487 template <int D, typename FEL = ScalarFiniteElement<D-1> >
01488 class NormalNeumannIntegrator 
01489   : public T_BIntegrator<DiffOpNormal<D>, DVec<1>, FEL>
01490 {
01491 public:
01493   NormalNeumannIntegrator (CoefficientFunction * coeff)
01494     : T_BIntegrator<DiffOpNormal<D>, DVec<1>, FEL> (DVec<1> (coeff))
01495   { ; }
01496 
01497   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01498   {
01499     return new NormalNeumannIntegrator (coeffs[0]);
01500   }
01501 
01503   virtual bool BoundaryForm () const { return 1; }
01505   virtual string Name () const { return "NormalNeumann"; }
01506 };
01507 
01508 
01509 
01510 
01511 
01512 
01514 template <int D, typename FEL = ScalarFiniteElement<D>  >
01515 class GradSourceIntegrator 
01516   : public T_BIntegrator<DiffOpGradient<D>, DVec<D>, FEL>
01517 {
01518 public:
01520   GradSourceIntegrator (CoefficientFunction * coeff1, 
01521                         CoefficientFunction * coeff2, 
01522                         CoefficientFunction * coeff3)
01523     : T_BIntegrator<DiffOpGradient<D>, DVec<D>, FEL> (DVec<D> (coeff1, coeff2, coeff3))
01524   { ; }
01525   
01526   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
01527   {
01528     return new GradSourceIntegrator (coeffs[0],coeffs[1],coeffs[2]);
01529   }
01530 
01532   virtual string Name () const { return "GradSource"; }
01533 };
01534 
01535 
01536 
01537 
01538 
01539 
01540 
01541 
01542 
01543 }
01544 
01545 #endif