NGSolve
4.9
|
00001 #ifndef FILE_DIFFOP 00002 #define FILE_DIFFOP 00003 00004 /*********************************************************************/ 00005 /* File: diffop.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 24. Nov. 2009 */ 00008 /*********************************************************************/ 00009 00010 namespace ngfem 00011 { 00012 00013 00014 00021 template<class DOP> 00022 class NGS_DLL_HEADER DiffOp 00023 { 00024 public: 00025 // enum { DIM_ELEMENT = TDOP::DIM_ELEMENT }; 00026 // enum { DIM_SPACE = TDOP::DIM_SPACE }; 00027 00035 template <typename FEL, typename MIP, typename MAT> 00036 static void GenerateMatrix (const FEL & fel, const MIP & mip, 00037 MAT & mat, LocalHeap & lh) 00038 { 00039 ; 00040 } 00041 00046 template <typename FEL, typename MIP, class TVX, class TVY> 00047 static void Apply (const FEL & fel, const MIP & mip, 00048 const TVX & x, TVY & y, 00049 LocalHeap & lh) 00050 { 00051 typedef typename TVY::TSCAL TSCAL; 00052 00053 HeapReset hr(lh); 00054 00055 FlatMatrixFixHeight<DOP::DIM_DMAT, TSCAL> mat(DOP::DIM*fel.GetNDof(), lh); 00056 DOP::GenerateMatrix (fel, mip, mat, lh); 00057 y = mat * x; 00058 } 00059 00061 template <typename FEL, class MIR, class TVX, class TVY> 00062 static void ApplyIR (const FEL & fel, const MIR & mir, 00063 const TVX & x, TVY & y, 00064 LocalHeap & lh) 00065 { 00066 for (int i = 0; i < mir.Size(); i++) 00067 Apply (fel, mir[i], x, y.Row(i), lh); 00068 } 00069 00070 00072 template <typename FEL, typename MIP, class TVX, class TVY> 00073 static void ApplyTrans (const FEL & fel, const MIP & mip, 00074 const TVX & x, TVY & y, 00075 LocalHeap & lh) 00076 { 00077 typedef typename TVY::TSCAL TSCAL; 00078 00079 HeapReset hr(lh); 00080 00081 FlatMatrixFixHeight<DOP::DIM_DMAT, TSCAL> mat(DOP::DIM*fel.GetNDof(), lh); 00082 DOP::GenerateMatrix (fel, mip, mat, lh); 00083 y = Trans (mat) * x; 00084 } 00085 00086 00088 template <typename FEL, typename MIR, class TVX, class TVY> 00089 static void ApplyTransIR (const FEL & fel, const MIR & mir, 00090 const TVX & x, TVY & y, 00091 LocalHeap & lh) 00092 { 00093 typedef typename TVY::TSCAL TSCAL; 00094 00095 HeapReset hr(lh); 00096 FlatVector<TSCAL> hy(y.Size(), lh); 00097 00098 y = 0.0; 00099 for (int i = 0; i < mir.Size(); i++) 00100 { 00101 ApplyTrans (fel, mir[i], x.Row(i), hy, lh); 00102 y += hy; 00103 } 00104 } 00105 00106 00107 00108 00109 00111 template <typename MIP, class TVX> 00112 static void Transform (const MIP & mip, TVX & x) 00113 { 00114 ; 00115 } 00116 00118 template <typename MIP, class TVX> 00119 static void TransformT (const MIP & mip, TVX & x) 00120 { 00121 ; 00122 } 00123 00124 /* 00126 template <typename FEL, class TVD, class TVY, int D> 00127 static void ApplyGrid (const FEL & fel, 00128 const IntegrationRuleTP<D> & ir, 00129 const TVY & hv, 00130 TVD & dvecs, 00131 LocalHeap & lh) 00132 { 00133 ; 00134 } 00135 00137 template <typename FEL, class TVD, class TVY, int D> 00138 static void ApplyTransGrid (const FEL & fel, 00139 const IntegrationRuleTP<D> & ir, 00140 const TVD & dvecs, 00141 TVY & hv, LocalHeap & lh) 00142 { 00143 ; 00144 } 00145 */ 00146 }; 00147 00148 00149 00150 00151 00152 00153 00159 class DifferentialOperator 00160 { 00161 public: 00163 NGS_DLL_HEADER virtual ~DifferentialOperator (); 00165 virtual int Dim() const = 0; 00167 virtual bool Boundary() const { return false; } 00169 virtual void 00170 CalcMatrix (const FiniteElement & fel, 00171 const BaseMappedIntegrationPoint & mip, 00172 FlatMatrix<double> mat, 00173 LocalHeap & lh) const 00174 { 00175 cerr << "DifferentialOperator::CalcMatrix called for base class" << endl; 00176 } 00177 00178 virtual void 00179 Apply (const FiniteElement & fel, 00180 const BaseMappedIntegrationPoint & mip, 00181 FlatVector<double> x, 00182 FlatVector<double> flux, 00183 LocalHeap & lh) const 00184 { 00185 FlatMatrix<> mat(Dim(), fel.GetNDof(), lh); 00186 CalcMatrix (fel, mip, mat, lh); 00187 flux = mat * x; 00188 } 00189 00190 virtual void 00191 Apply (const FiniteElement & fel, 00192 const BaseMappedIntegrationPoint & mip, 00193 FlatVector<Complex> x, 00194 FlatVector<Complex> flux, 00195 LocalHeap & lh) const 00196 { 00197 FlatMatrix<> mat(Dim(), fel.GetNDof(), lh); 00198 CalcMatrix (fel, mip, mat, lh); 00199 flux = mat * x; 00200 } 00201 00202 virtual void 00203 Apply (const FiniteElement & fel, 00204 const BaseMappedIntegrationRule & mir, 00205 FlatVector<double> x, 00206 FlatMatrix<double> flux, 00207 LocalHeap & lh) const 00208 { 00209 for (int i = 0; i < mir.Size(); i++) 00210 Apply (fel, mir[i], x, flux.Row(i), lh); 00211 } 00212 00213 virtual void 00214 Apply (const FiniteElement & fel, 00215 const BaseMappedIntegrationRule & mir, 00216 FlatVector<Complex> x, 00217 FlatMatrix<Complex> flux, 00218 LocalHeap & lh) const 00219 { 00220 for (int i = 0; i < mir.Size(); i++) 00221 Apply (fel, mir[i], x, flux.Row(i), lh); 00222 } 00223 00224 00225 00226 00227 virtual void 00228 ApplyTrans (const FiniteElement & fel, 00229 const BaseMappedIntegrationPoint & mip, 00230 FlatVector<double> flux, 00231 FlatVector<double> x, 00232 LocalHeap & lh) const 00233 { 00234 FlatMatrix<> mat(Dim(), fel.GetNDof(), lh); 00235 CalcMatrix (fel, mip, mat, lh); 00236 flux = mat * x; 00237 } 00238 00239 00240 virtual void 00241 ApplyTrans (const FiniteElement & fel, 00242 const BaseMappedIntegrationPoint & mip, 00243 FlatVector<Complex> flux, 00244 FlatVector<Complex> x, 00245 LocalHeap & lh) const 00246 { 00247 FlatMatrix<> mat(Dim(), fel.GetNDof(), lh); 00248 CalcMatrix (fel, mip, mat, lh); 00249 flux = mat * x; 00250 } 00251 00252 00253 virtual void 00254 ApplyTrans (const FiniteElement & fel, 00255 const BaseMappedIntegrationRule & mir, 00256 FlatMatrix<double> flux, 00257 FlatVector<double> x, 00258 LocalHeap & lh) const 00259 { 00260 FlatVector<double> hx(x.Size(), lh); 00261 x = 0.0; 00262 for (int i = 0; i < mir.Size(); i++) 00263 { 00264 ApplyTrans (fel, mir[i], flux.Row(i), hx, lh); 00265 x += hx; 00266 } 00267 } 00268 00269 virtual void 00270 ApplyTrans (const FiniteElement & fel, 00271 const BaseMappedIntegrationRule & mir, 00272 FlatMatrix<Complex> flux, 00273 FlatVector<Complex> x, 00274 LocalHeap & lh) const 00275 { 00276 FlatVector<Complex> hx(x.Size(), lh); 00277 x = 0.0; 00278 for (int i = 0; i < mir.Size(); i++) 00279 { 00280 ApplyTrans (fel, mir[i], flux.Row(i), hx, lh); 00281 x += hx; 00282 } 00283 } 00284 }; 00285 00286 00287 class BlockDifferentialOperator : public DifferentialOperator 00288 { 00289 const DifferentialOperator & diffop; 00290 int dim; 00291 int comp; 00292 public: 00293 BlockDifferentialOperator (const DifferentialOperator & adiffop, 00294 int adim, int acomp = -1) 00295 : diffop(adiffop), dim(adim), comp(acomp) { ; } 00296 00298 virtual int Dim() const { return dim*diffop.Dim(); } 00299 virtual bool Boundary() const { return diffop.Boundary(); } 00300 00301 00302 00303 virtual void 00304 CalcMatrix (const FiniteElement & fel, 00305 const BaseMappedIntegrationPoint & mip, 00306 FlatMatrix<double> mat, 00307 LocalHeap & lh) const 00308 { 00309 FlatMatrix<double> mat1(diffop.Dim(), fel.GetNDof(), lh); 00310 diffop.CalcMatrix (fel, mip, mat1, lh); 00311 mat = 0; 00312 00313 if (comp == -1) 00314 for (int i = 0; i < mat1.Height(); i++) 00315 for (int j = 0; j < mat1.Width(); j++) 00316 for (int k = 0; k < dim; k++) 00317 mat(dim*i+k, dim*j+k) = mat1(i,j); 00318 else 00319 for (int i = 0; i < mat1.Height(); i++) 00320 for (int j = 0; j < mat1.Width(); j++) 00321 mat(dim*i+comp, dim*j+comp) = mat1(i,j); 00322 } 00323 00324 }; 00325 00326 00327 00331 template <typename DIFFOP> 00332 class T_DifferentialOperator : public DifferentialOperator 00333 { 00334 protected: 00335 enum { DIM_SPACE = DIFFOP::DIM_SPACE }; 00336 enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT }; 00337 enum { DIM_DMAT = DIFFOP::DIM_DMAT }; 00338 enum { DIM = DIFFOP::DIM }; 00339 00340 public: 00341 virtual int Dim() const { return DIFFOP::DIM_DMAT; } 00342 virtual bool Boundary() const { return int(DIM_SPACE) > int(DIM_ELEMENT); } 00343 00344 virtual void 00345 CalcMatrix (const FiniteElement & bfel, 00346 const BaseMappedIntegrationPoint & bmip, 00347 FlatMatrix<double> mat, 00348 LocalHeap & lh) const 00349 { 00350 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = 00351 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00352 DIFFOP::GenerateMatrix (bfel, mip, mat, lh); 00353 } 00354 00355 virtual void 00356 Apply (const FiniteElement & bfel, 00357 const BaseMappedIntegrationPoint & bmip, 00358 FlatVector<double> x, 00359 FlatVector<double> flux, 00360 LocalHeap & lh) const 00361 { 00362 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = 00363 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00364 DIFFOP::Apply (bfel, mip, x, flux, lh); 00365 } 00366 00367 00368 virtual void 00369 Apply (const FiniteElement & bfel, 00370 const BaseMappedIntegrationRule & bmir, 00371 FlatVector<double> x, 00372 FlatMatrix<double> flux, 00373 LocalHeap & lh) const 00374 { 00375 const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> & mir = 00376 static_cast<const MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE>&> (bmir); 00377 DIFFOP::ApplyIR (bfel, mir, x, flux, lh); 00378 } 00379 00380 virtual void 00381 ApplyTrans (const FiniteElement & bfel, 00382 const BaseMappedIntegrationPoint & bmip, 00383 FlatVector<double> flux, 00384 FlatVector<double> x, 00385 LocalHeap & lh) const 00386 { 00387 const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE> & mip = 00388 static_cast<const MappedIntegrationPoint<DIM_ELEMENT,DIM_SPACE>&> (bmip); 00389 DIFFOP::ApplyTrans (bfel, mip, flux, x, lh); 00390 } 00391 }; 00392 00393 00394 00395 00396 00397 00398 00399 00400 00401 // new design, code is still experimental ... 00402 00403 template <typename DOP, typename F> 00404 class T_FunctionDiffOp : public DifferentialOperator 00405 { 00406 00407 // possible conversion from vector to scalar 00408 class V2VS 00409 { 00410 FlatVector<> v; 00411 public: 00412 V2VS (FlatVector<> av) : v(av) { ; } 00413 00414 template <int D> 00415 V2VS (Vec<D> av) : v(av) { ; } 00416 00417 operator double () { return v(0); } 00418 00419 operator FlatVector<> () { return v; } 00420 00421 template <int D> 00422 operator Vec<D> () { return v; } 00423 }; 00424 00425 00426 int dim; 00427 const F & func; 00428 public: 00429 00430 T_FunctionDiffOp (const F & afunc, int adim) : func(afunc), dim(adim) { ; } 00431 00432 virtual int Dim() const { return dim; } 00433 00434 virtual void Apply (const FiniteElement & fel, 00435 const BaseMappedIntegrationPoint & mip, 00436 FlatVector<double> x, 00437 FlatVector<double> flux, 00438 LocalHeap & lh) const 00439 { 00440 Vec<DOP::DIM_DMAT> u; 00441 DOP::Apply (fel, mip, x, u, lh); 00442 flux = func(V2VS(u)); 00443 } 00444 }; 00445 00446 00447 template <typename DOP, typename F> 00448 DifferentialOperator * CreateFunctionDiffOp (const DOP & dop, 00449 const F & func, int dim = 1) 00450 { 00451 return new T_FunctionDiffOp<DOP, F> (func, dim); 00452 } 00453 00454 00455 /* examples: 00456 00457 double myexp (double x) { return exp(x); } 00458 Vec<1> myexpVec (FlatVector<> x) { return Vec<1> (exp(x(0))); } 00459 00460 CreateFunctionDiffOp(DiffOpId<2>(), myexpVec)); 00461 */ 00462 00463 } 00464 00465 00466 #endif