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