NGSolve  4.9
fem/hdiv_equations.hpp
00001 #ifndef FILE_HDIV_EQUATIONS
00002 #define FILE_HDIV_EQUATIONS
00003 
00004 /*********************************************************************/
00005 /* File:   hdiv_equations.hpp                                        */
00006 /* Author: Joachim Schoeberl, Almedin Becirovic                      */
00007 /* Date:   10. Feb. 2002                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngfem
00011 {
00012 
00013 
00014 /*
00015   Finite Element Integrators for H(div)
00016 
00017   Mapping with Piola transformation
00018 
00019   Requires H(div) finite elements
00020 */
00021 
00022 
00023 
00025 template <int D, typename FEL = HDivFiniteElement<D> >
00026 class DiffOpIdHDiv : public DiffOp<DiffOpIdHDiv<D, FEL> >
00027 {
00028 public:
00029   enum { DIM = 1 };
00030   enum { DIM_SPACE = D };
00031   enum { DIM_ELEMENT = D };
00032   enum { DIM_DMAT = D };
00033   enum { DIFFORDER = 0 };
00034     
00035   template <typename AFEL, typename SIP, typename MAT>
00036   static void GenerateMatrix (const AFEL & fel, const SIP & sip,
00037                               MAT & mat, LocalHeap & lh)
00038   {
00039     mat = (1.0/sip.GetJacobiDet()) * 
00040       (sip.GetJacobian() * Trans (static_cast<const FEL&> (fel).GetShape(sip.IP(), lh)));
00041   }
00042     
00043   template <typename AFEL, typename SIP, class TVX, class TVY>
00044   static void Apply (const AFEL & fel, const SIP & sip,
00045                      const TVX & x, TVY & y,
00046                      LocalHeap & lh) 
00047   {
00048     typedef typename TVX::TSCAL TSCAL;
00049       
00050     Vec<D,TSCAL> hv = Trans (static_cast<const FEL&>(fel).GetShape(sip.IP(), lh)) * x;
00051     hv *= (1.0/sip.GetJacobiDet());
00052     y = sip.GetJacobian() * hv;
00053   }
00054 
00055 
00056   template <typename AFEL, class MIR>
00057   static void ApplyIR (const AFEL & fel, const MIR & mir,
00058                        const FlatVector<double> & x, FlatMatrix<double> & y,
00059                        LocalHeap & lh)
00060   {
00061     static_cast<const FEL&>(fel).Evaluate (mir.IR(), x, FlatMatrixFixWidth<D> (y.Height(), &y(0,0)));
00062     for (int i = 0; i < mir.Size(); i++)
00063       {
00064         Vec<D> hy = (1.0/mir[i].GetJacobiDet()) * mir[i].GetJacobian() * y.Row(i);
00065         y.Row(i) = hy;
00066       }
00067   }
00068 
00069 
00070 
00071   template <typename AFEL, typename SIP, class TVX, class TVY>
00072   static void ApplyTrans (const AFEL & fel, const SIP & sip,
00073                           const TVX & x, TVY & y,
00074                           LocalHeap & lh) 
00075   {
00076     typedef typename TVX::TSCAL TSCAL;
00077 
00078     Vec<D,TSCAL> hv = Trans (sip.GetJacobian()) * x;
00079     hv *= (1.0/sip.GetJacobiDet());
00080     y = static_cast<const FEL&>(fel).GetShape(sip.IP(),lh) * hv;
00081   }
00082 };
00083 
00084 
00085 
00086 
00088 template <int D, typename FEL = HDivFiniteElement<D> >
00089 class DiffOpDivHDiv : public DiffOp<DiffOpDivHDiv<D, FEL> >
00090 {
00091 public:
00092   enum { DIM = 1 };
00093   enum { DIM_SPACE = D };
00094   enum { DIM_ELEMENT = D };
00095   enum { DIM_DMAT = 1 };
00096   enum { DIFFORDER = 1 };
00097 
00098   template <typename AFEL, typename SIP, typename MAT>
00099   static void GenerateMatrix (const AFEL & fel, const SIP & sip,
00100                               MAT & mat, LocalHeap & lh)
00101   {
00102     
00103     mat = 1.0/sip.GetJacobiDet() *
00104       Trans (static_cast<const FEL&>(fel).GetDivShape(sip.IP(), lh));
00105   }
00106 
00107   template <typename AFEL, typename SIP>
00108   static void GenerateMatrix (const AFEL & fel, const SIP & sip,
00109                               FlatVector<double> & mat, LocalHeap & lh)
00110   {
00111     mat = 1.0/sip.GetJacobiDet() * 
00112       (static_cast<const FEL&>(fel).GetDivShape(sip.IP(), lh));
00113   }
00114 
00115 
00116   template <typename AFEL, typename SIP, class TVX, class TVY>
00117   static void Apply (const AFEL & fel, const SIP & sip,
00118                      const TVX & x, TVY & y,
00119                      LocalHeap & lh) 
00120   {
00121     typedef typename TVX::TSCAL TSCAL;
00122       
00123     Vec<DIM,TSCAL> hv = Trans (static_cast<const FEL&>(fel).GetDivShape(sip.IP(), lh)) * x;
00124     y = (1.0/sip.GetJacobiDet()) * hv;
00125   }
00126 
00127   template <typename AFEL, typename SIP, class TVX, class TVY>
00128   static void ApplyTrans (const AFEL & fel, const SIP & sip,
00129                           const TVX & x, TVY & y,
00130                           LocalHeap & lh) 
00131   {
00132     typedef typename TVX::TSCAL TSCAL;
00133     Vec<DIM,TSCAL> hv = x;
00134     hv *= (1.0/sip.GetJacobiDet());
00135     y = static_cast<const FEL&>(fel).GetDivShape(sip.IP(),lh) * hv;
00136   }
00137 };
00138 
00139 
00140 
00141 
00143 template <int D, typename FEL = HDivNormalFiniteElement<D-1> >
00144 class DiffOpIdHDivBoundary : public DiffOp<DiffOpIdHDivBoundary<D, FEL> >
00145 {
00146 public:
00147   enum { DIM = 1 };
00148   enum { DIM_SPACE = D };
00149   enum { DIM_ELEMENT = D-1 };
00150   enum { DIM_DMAT = 1 };
00151   enum { DIFFORDER = 0 };
00152 
00153   template <typename AFEL, typename SIP, typename MAT>
00154   static void GenerateMatrix (const AFEL & fel, const SIP & sip,
00155                               MAT & mat, LocalHeap & lh)
00156   {
00157     mat =  (1.0/sip.GetJacobiDet())*
00158       Trans(static_cast<const FEL&> (fel).GetShape (sip.IP(), lh));
00159   }
00160 
00161   template <typename AFEL, typename SIP, class TVX, class TVY>
00162   static void Apply (const AFEL & fel, const SIP & sip,
00163                      const TVX & x, TVY & y,
00164                      LocalHeap & lh)
00165   {
00166     y = (1.0/sip.GetJacobiDet())*
00167       (Trans (static_cast<const FEL&> (fel).GetShape (sip.IP(), lh)) * x);
00168   }
00169 
00170   template <typename AFEL, typename SIP, class TVX, class TVY>
00171   static void ApplyTrans (const AFEL & fel, const SIP & sip,
00172                           const TVX & x, TVY & y,
00173                           LocalHeap & lh)
00174   {
00175     y = static_cast<const FEL&> (fel).GetShape (sip.IP(), lh)*((1.0/sip.GetJacobiDet())* x);
00176   }
00177 };
00178 
00179 
00181 template <int D, typename FEL = HDivNormalFiniteElement<D-1> >
00182 class DiffOpIdVecHDivBoundary : public DiffOp<DiffOpIdVecHDivBoundary<D,FEL> >
00183 {
00184 public:
00185   enum { DIM = 1 };
00186   enum { DIM_SPACE = D };
00187   enum { DIM_ELEMENT = D-1 };
00188   enum { DIM_DMAT = D };
00189   enum { DIFFORDER = 0 };
00190 
00191   template <typename AFEL, typename SIP, typename MAT>
00192   static void GenerateMatrix (const AFEL & fel, const SIP & sip,
00193                               MAT & mat, LocalHeap & lh)
00194   {
00195     mat =  (1.0/sip.GetJacobiDet())*
00196       Trans(static_cast<const FEL&> (fel).GetShape (sip.IP(), lh)) 
00197       * Trans (sip.GetNV());;
00198   }
00199 
00200   template <typename AFEL, typename SIP, class TVX, class TVY> 
00201   static void Apply (const AFEL & fel, const SIP & sip,
00202                      const TVX & x, TVY & y,
00203                      LocalHeap & lh)
00204   {
00205     y = ( (1.0/sip.GetJacobiDet())*(InnerProduct (static_cast<const FEL&> (fel).GetShape (sip.IP(), lh), x) )) * sip.GetNV();
00206   }
00207 
00208   template <typename AFEL, typename SIP, class TVX, class TVY>
00209   static void ApplyTrans (const AFEL & fel, const SIP & sip,
00210                           const TVX & x, TVY & y,
00211                           LocalHeap & lh)
00212   {
00213     y = ((1.0/sip.GetJacobiDet())* InnerProduct (x, sip.GetNV()) ) * static_cast<const FEL&> (fel).GetShape (sip.IP(), lh);
00214   }
00215 };
00216 
00217 
00218 
00219 
00221 template <int D>
00222 class MassHDivIntegrator
00223   : public T_BDBIntegrator<DiffOpIdHDiv<D>, DiagDMat<D>, HDivFiniteElement<D> >
00224 {
00225 public:
00226   NGS_DLL_HEADER MassHDivIntegrator (CoefficientFunction * coeff);
00227 
00228   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
00229   {
00230     return new MassHDivIntegrator (coeffs[0]);
00231   }
00232 
00233   virtual string Name () const { return "MassHDiv"; }
00234 };
00235 
00236   
00237 
00239 template <int D>
00240 class DivDivHDivIntegrator
00241   : public T_BDBIntegrator<DiffOpDivHDiv<D>, DiagDMat<1>, HDivFiniteElement<D> >
00242 {
00243 public:
00244   NGS_DLL_HEADER DivDivHDivIntegrator (CoefficientFunction * coeff);
00245 
00246   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
00247   {
00248     return new DivDivHDivIntegrator (coeffs[0]);
00249   }
00250     
00251   virtual string Name () const { return "DivDivHDiv"; }
00252 };
00253 
00254 
00255 
00256 
00258 template <int D>
00259 class NGS_DLL_HEADER DivSourceHDivIntegrator 
00260   : public T_BIntegrator<DiffOpDivHDiv<D>, DVec<1>, HDivFiniteElement<D> >
00261 {
00262 public:
00263   DivSourceHDivIntegrator (CoefficientFunction * coeff);
00264 
00265   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
00266   {
00267     return new DivSourceHDivIntegrator (coeffs[0]);
00268   }
00269 
00270   virtual string Name () const { return "DivSourceHDiv"; }
00271 };
00272 
00273 /*
00275 template <int D>
00276 class SourceHDivIntegrator 
00277   : public T_BIntegrator<DiffOpIdHDiv<D>, DVec<D>, HDivFiniteElement<D> >
00278 {
00279 public:
00280   SourceHDivIntegrator (CoefficientFunction * coeff1,
00281                         CoefficientFunction * coeff2,
00282                         CoefficientFunction * coeff3);
00283 
00284   SourceHDivIntegrator (CoefficientFunction * coeff1,
00285                         CoefficientFunction * coeff2);
00286 
00287   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
00288   {
00289     if(D == 2)
00290       return new SourceHDivIntegrator (coeffs[0], coeffs[1]);
00291     else if (D == 3)
00292       return new SourceHDivIntegrator (coeffs[0], coeffs[1], coeffs[2]);
00293   }
00294     
00296   virtual string Name () const { return "SourceHDiv"; }
00297   };
00298 */
00299 
00300 
00301   template <int D> class SourceHDivIntegrator;
00302 
00303 
00304 template <int D>
00305 class NGS_DLL_HEADER BaseSourceHDivIntegrator 
00306   : public T_BIntegrator<DiffOpIdHDiv<D>, DVec<D>, HDivFiniteElement<D> >
00307 {
00308 public:
00309   BaseSourceHDivIntegrator(const DVec<D> & dvec)
00310     : T_BIntegrator<DiffOpIdHDiv<D>, DVec<D>, HDivFiniteElement<D> > (dvec) { ; }
00311 
00312   static Integrator * Create (Array<CoefficientFunction*> & coeffs);
00313   /*
00314   {
00315     if(D == 2)
00316       return new SourceHDivIntegrator<2> (coeffs[0], coeffs[1]);
00317     else if (D == 3)
00318       return new SourceHDivIntegrator<3> (coeffs[0], coeffs[1], coeffs[2]);
00319   }
00320   */
00322   virtual string Name () const { return "SourceHDiv"; }
00323 };
00324 
00325 
00326 
00327 template <>
00328 class NGS_DLL_HEADER SourceHDivIntegrator<2>
00329   : public BaseSourceHDivIntegrator<2> 
00330 {
00331 public:
00332   SourceHDivIntegrator (CoefficientFunction * coeff1,
00333                         CoefficientFunction * coeff2);
00334 };
00335 
00336 template <>
00337 class NGS_DLL_HEADER SourceHDivIntegrator<3>
00338   : public BaseSourceHDivIntegrator<3> 
00339 {
00340 public:
00341   SourceHDivIntegrator (CoefficientFunction * coeff1,
00342                         CoefficientFunction * coeff2,
00343                         CoefficientFunction * coeff3);
00344 };
00345 
00346 
00347 
00348 template <int D>
00349 class NGS_DLL_HEADER SourceHDivIntegratorN 
00350   : public T_BIntegrator<DiffOpIdHDiv<D>, DVecN<D>, HDivFiniteElement<D> >
00351 {
00352 public:
00353   SourceHDivIntegratorN(CoefficientFunction * coeff)
00354     : T_BIntegrator<DiffOpIdHDiv<D>, DVecN<D>, HDivFiniteElement<D> > (DVecN<D> (coeff)) { ; }
00355   virtual string Name () const { return "VecSourceHDiv"; }
00356 };
00357 
00358 
00359 
00360 
00361 
00362 
00363 
00365 template <int D, typename FEL = HDivNormalFiniteElement<D-1> >
00366 class NGS_DLL_HEADER NeumannHDivIntegrator
00367   : public T_BIntegrator<DiffOpIdHDivBoundary<D>, DVec<1>, FEL>
00368 {
00369 public:
00371   NeumannHDivIntegrator (CoefficientFunction * coeff)
00372     : T_BIntegrator<DiffOpIdHDivBoundary<D>, DVec<1>, FEL> (DVec<1> (coeff))
00373   { ; }
00374 
00375   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
00376   {
00377     return new NeumannHDivIntegrator (coeffs[0]);
00378   }
00379 
00381   virtual bool BoundaryForm () const { return 1; }
00383   virtual string Name () const { return "NeumannHDiv"; }
00384 };
00385 
00386 
00388 template <int D>
00389 class RobinHDivIntegrator
00390   : public T_BDBIntegrator<DiffOpIdHDivBoundary<D>, DiagDMat<1>, HDivNormalFiniteElement<D-1> >
00391 {
00392 public:
00393   NGS_DLL_HEADER RobinHDivIntegrator (CoefficientFunction * coeff)
00394     : T_BDBIntegrator<DiffOpIdHDivBoundary<D>, DiagDMat<1>, HDivNormalFiniteElement<D-1> > (DiagDMat<1> (coeff))
00395   { ; }
00396 
00397 
00398   static Integrator * Create (Array<CoefficientFunction*> & coeffs)
00399   {
00400     return new RobinHDivIntegrator (coeffs[0]);
00401 
00402   }
00403 
00404   virtual bool BoundaryForm () const { return 1; }
00405   virtual string Name () const { return "RobinHDiv"; }
00406 };
00407 
00408 
00409 }
00410 
00411 
00412 #endif
00413 
00414 
00415