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