NGSolve
4.9
|
00001 #ifndef FILE_HCURLFE 00002 #define FILE_HCURLFE 00003 00004 /*********************************************************************/ 00005 /* File: hcurlfe.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 16. Apr. 2000 */ 00008 /*********************************************************************/ 00009 00010 00011 00012 namespace ngfem 00013 { 00014 00015 00016 /* 00017 HCurl Finite Element Definitions 00018 */ 00019 00020 00021 template <int D> class HDivFiniteElement; 00022 00023 00024 template <int D> 00025 class DIM_CURL_TRAIT 00026 { 00027 public: 00028 enum { DIM = (D*(D-1))/2 }; 00029 }; 00030 00031 00032 template <> class DIM_CURL_TRAIT<1> 00033 { 00034 public: 00035 enum { DIM = 0 }; // should be 0; set to 1 to make gcc34 feel ok 00036 }; 00037 00038 00042 template <int D> 00043 class NGS_DLL_HEADER HCurlFiniteElement : public FiniteElement 00044 { 00045 00046 public: 00047 enum { DIM = D }; 00048 enum { DIM_CURL = DIM_CURL_TRAIT<D>::DIM }; 00049 00050 00051 public: 00053 HCurlFiniteElement () { ; } 00054 00056 HCurlFiniteElement (ELEMENT_TYPE aeltype, int andof, int aorder) 00057 : FiniteElement (aeltype, andof, aorder) { ; } 00058 00059 virtual ~HCurlFiniteElement () { ; } 00060 00061 virtual string ClassName() const; 00062 00064 virtual void CalcShape (const IntegrationPoint & ip, 00065 FlatMatrixFixWidth<DIM> shape) const = 0; 00066 00068 virtual void CalcCurlShape (const IntegrationPoint & ip, 00069 FlatMatrixFixWidth<DIM_CURL> curlshape) const; 00070 00072 virtual void CalcMappedShape (const MappedIntegrationPoint<DIM,DIM> & mip, 00073 FlatMatrixFixWidth<DIM> shape) const; 00074 00076 virtual void CalcMappedCurlShape (const MappedIntegrationPoint<DIM,DIM> & mip, 00077 FlatMatrixFixWidth<DIM_CURL> curlshape) const; 00078 00079 00081 const FlatMatrixFixWidth<DIM> GetShape (const IntegrationPoint & ip, 00082 LocalHeap & lh) const 00083 { 00084 FlatMatrixFixWidth<DIM> shape(ndof, lh); 00085 CalcShape (ip, shape); 00086 return shape; 00087 } 00088 00090 const FlatMatrixFixWidth<DIM_CURL> GetCurlShape (const IntegrationPoint & ip, 00091 LocalHeap & lh) const 00092 { 00093 FlatMatrixFixWidth<DIM_CURL> curlshape(ndof, lh); 00094 CalcCurlShape (ip, curlshape); 00095 return curlshape; 00096 } 00097 00098 template <typename TVX> 00099 Vec<DIM_CURL, typename TVX::TSCAL> 00100 EvaluateCurlShape (const IntegrationPoint & ip, 00101 const TVX & x, LocalHeap & lh) const 00102 { 00103 return Trans (GetCurlShape(ip, lh)) * x; 00104 } 00105 00106 virtual Vec<DIM_CURL> 00107 EvaluateCurlShape (const IntegrationPoint & ip, 00108 FlatVector<double> x, LocalHeap & lh) const 00109 { 00110 return Trans (GetCurlShape(ip, lh)) * x; 00111 } 00112 00113 00114 protected: 00116 virtual void CalcShape1 (const IntegrationPoint & ip, 00117 FlatMatrixFixWidth<D> shape) const 00118 { ; } 00120 virtual void CalcShape2 (const IntegrationPoint & ip, 00121 FlatMatrixFixWidth<D> shape) const 00122 { ; } 00123 00124 virtual void CalcShape3 (const IntegrationPoint & ip, 00125 FlatMatrixFixWidth<D> shape) const 00126 { ; } 00127 00128 virtual void CalcShape4 (const IntegrationPoint & ip, 00129 FlatMatrixFixWidth<D> shape) const 00130 { ; } 00131 00133 void ComputeEdgeMoments (int enr, ScalarFiniteElement<1> & testfe, 00134 FlatMatrix<> moments, int order, int shape = 1) const; 00136 void ComputeFaceMoments (int fnr, HDivFiniteElement<2> & testfe, 00137 FlatMatrix<> moments, int order, int shape = 1) const; 00139 void ComputeVolMoments (HDivFiniteElement<3> & testfe, 00140 FlatMatrix<> moments, int order, int shape = 1) const; 00141 }; 00142 00143 00144 00145 00146 00147 00148 00149 00150 00151 00152 00153 00154 00155 00156 // hv.DValue() = (grad u) x (grad v) 00157 inline AutoDiff<3> Cross (const AutoDiff<3> & u, 00158 const AutoDiff<3> & v) 00159 { 00160 AutoDiff<3> hv; 00161 hv.Value() = 0.0; 00162 hv.DValue(0) = u.DValue(1)*v.DValue(2)-u.DValue(2)*v.DValue(1); 00163 hv.DValue(1) = u.DValue(2)*v.DValue(0)-u.DValue(0)*v.DValue(2); 00164 hv.DValue(2) = u.DValue(0)*v.DValue(1)-u.DValue(1)*v.DValue(0); 00165 return hv; 00166 } 00167 00168 inline AutoDiff<1> Cross (const AutoDiff<2> & u, 00169 const AutoDiff<2> & v) 00170 { 00171 AutoDiff<1> hv; 00172 hv.Value() = 0.0; 00173 hv.DValue(0) = u.DValue(0)*v.DValue(1)-u.DValue(1)*v.DValue(0); 00174 return hv; 00175 } 00176 00177 00178 00179 00180 00181 00182 00183 template <int DIM> 00184 class Du 00185 { 00186 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00187 00188 public: 00189 const AutoDiff<DIM> & u; 00190 00191 Du (const AutoDiff<DIM> & au) 00192 : u(au) { ; } 00193 00194 Vec<DIM> Value () const 00195 { 00196 Vec<DIM> val; 00197 for (int j = 0; j < DIM; j++) 00198 val(j) = u.DValue(j); 00199 return val; 00200 } 00201 00202 Vec<DIM_CURL> CurlValue () const 00203 { 00204 return Vec<DIM> (0.0); 00205 } 00206 }; 00207 00208 00209 00210 00211 template <int DIM> 00212 class uDv 00213 { 00214 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00215 00216 public: 00217 const AutoDiff<DIM> & u, v; 00218 00219 uDv (const AutoDiff<DIM> & au, 00220 const AutoDiff<DIM> & av) 00221 : u(au), v(av) { ; } 00222 00223 Vec<DIM> Value () const 00224 { 00225 Vec<DIM> val; 00226 for (int j = 0; j < DIM; j++) 00227 val(j) = u.Value() * v.DValue(j); 00228 return val; 00229 } 00230 00231 Vec<DIM_CURL> CurlValue () const 00232 { 00233 AutoDiff<DIM_CURL> hd = Cross (u, v); 00234 Vec<DIM_CURL> val; 00235 for (int i = 0; i < DIM_CURL; i++) 00236 val(i) = hd.DValue(i); 00237 return val; 00238 } 00239 }; 00240 00241 00242 00243 template <int DIM> 00244 class uDv_minus_vDu 00245 { 00246 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00247 00248 public: 00249 const AutoDiff<DIM> & u, v; 00250 00251 uDv_minus_vDu (const AutoDiff<DIM> & au, 00252 const AutoDiff<DIM> & av) 00253 : u(au), v(av) { ; } 00254 00255 Vec<DIM> Value () const 00256 { 00257 Vec<DIM> val; 00258 for (int j = 0; j < DIM; j++) 00259 val(j) = u.Value() * v.DValue(j) - v.Value() * u.DValue(j); 00260 return val; 00261 } 00262 00263 Vec<DIM_CURL> CurlValue () const 00264 { 00265 AutoDiff<DIM_CURL> hd = Cross (u, v); 00266 Vec<DIM_CURL> val; 00267 for (int i = 0; i < DIM_CURL; i++) 00268 val(i) = 2 * hd.DValue(i); 00269 return val; 00270 } 00271 }; 00272 00273 00274 template <int DIM> 00275 class wuDv_minus_wvDu 00276 { 00277 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00278 00279 public: 00280 const AutoDiff<DIM> & u, v, w; 00281 00282 wuDv_minus_wvDu (const AutoDiff<DIM> & au, 00283 const AutoDiff<DIM> & av, 00284 const AutoDiff<DIM> & aw) 00285 : u(au), v(av), w(aw) { ; } 00286 00287 Vec<DIM> Value () const 00288 { 00289 Vec<DIM> val; 00290 for (int j = 0; j < DIM; j++) 00291 val(j) = w.Value() * (u.Value() * v.DValue(j) - v.Value() * u.DValue(j)); 00292 return val; 00293 } 00294 00295 Vec<DIM_CURL> CurlValue () const 00296 { 00297 AutoDiff<DIM_CURL> hd = Cross (u*w, v) + Cross(u, v*w); 00298 Vec<DIM_CURL> val; 00299 for (int i = 0; i < DIM_CURL; i++) 00300 val(i) = hd.DValue(i); 00301 return val; 00302 } 00303 }; 00304 00305 00306 00307 00308 00309 00310 template <int DIM> 00311 class HCurlShapeElement 00312 { 00313 FlatVec<DIM> data; 00314 public: 00315 HCurlShapeElement (FlatVec<DIM> adata) : data(adata) { ; } 00316 00317 template <typename F1> 00318 void operator= (F1 form1) { data = form1.Value(); } 00319 }; 00320 00321 template <int DIM> 00322 class HCurlCurlShapeElement 00323 { 00324 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00325 FlatVec<DIM_CURL> data; 00326 public: 00327 HCurlCurlShapeElement (FlatVec<DIM_CURL> adata) : data(adata) { ; } 00328 00329 template <typename F1> 00330 void operator= (F1 form1) { data = form1.CurlValue(); } 00331 }; 00332 00333 00334 template <int DIM> 00335 class HCurlEvaluateCurlElement 00336 { 00337 const double & coef; 00338 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00339 Vec<DIM_CURL> & sum; 00340 public: 00341 HCurlEvaluateCurlElement (const double & acoef, Vec<DIM_CURL> & asum) 00342 : coef(acoef), sum(asum) { ; } 00343 00344 template <typename F1> 00345 void operator= (F1 form1) { sum += coef * form1.CurlValue(); } 00346 }; 00347 00348 00349 00350 template <int DIM> 00351 class HCurlShapeAssign 00352 { 00353 FlatMatrixFixWidth<DIM> shape; 00354 public: 00355 HCurlShapeAssign (FlatMatrixFixWidth<DIM> mat) : shape(mat) { ; } 00356 00357 HCurlShapeElement<DIM> operator[] (int i) const 00358 { return HCurlShapeElement<DIM> (shape.Row(i)); } 00359 }; 00360 00361 template <int DIM> 00362 class HCurlCurlShapeAssign 00363 { 00364 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00365 FlatMatrixFixWidth<DIM_CURL> cshape; 00366 public: 00367 HCurlCurlShapeAssign (FlatMatrixFixWidth<DIM_CURL> mat) : cshape(mat) { ; } 00368 00369 HCurlCurlShapeElement<DIM> operator[] (int i) const 00370 { return HCurlCurlShapeElement<DIM> (cshape.Row(i)); } 00371 }; 00372 00373 template <int DIM> 00374 class HCurlEvaluateCurl 00375 { 00376 FlatVector<> coefs; 00377 enum { DIM_CURL = (DIM * (DIM-1))/2 }; 00378 Vec<DIM_CURL> sum; 00379 public: 00380 HCurlEvaluateCurl (FlatVector<> acoefs) : coefs(acoefs), sum(0.0) { ; } 00381 00382 HCurlEvaluateCurlElement<DIM> operator[] (int i) 00383 { return HCurlEvaluateCurlElement<DIM> (coefs(i), sum); } 00384 00385 Vec<DIM_CURL> Sum() { return sum; } 00386 }; 00387 00388 00389 00390 00391 00392 00393 00399 template <class FEL, ELEMENT_TYPE ET, int NDOF, int ORDER> 00400 class T_HCurlFiniteElement : public HCurlFiniteElement<ET_trait<ET>::DIM> 00401 { 00402 00403 public: 00404 00405 protected: 00406 enum { DIM = ET_trait<ET>::DIM }; 00407 using HCurlFiniteElement<DIM>::DIM_CURL; 00408 00409 T_HCurlFiniteElement () 00410 : HCurlFiniteElement<DIM> (ET, NDOF, ORDER) { ; } 00411 00412 virtual ~T_HCurlFiniteElement() { ; } 00413 00414 public: 00415 00416 virtual void CalcShape (const IntegrationPoint & ip, 00417 FlatMatrixFixWidth<DIM> shape) const 00418 { 00419 AutoDiff<DIM> adp[DIM]; 00420 for (int i = 0; i < DIM; i++) 00421 adp[i] = AutoDiff<DIM> (ip(i), i); 00422 00423 HCurlShapeAssign<DIM> ds(shape); 00424 FEL::T_CalcShape (adp, ds); 00425 } 00426 00427 virtual void 00428 CalcMappedShape (const MappedIntegrationPoint<DIM,DIM> & mip, 00429 FlatMatrixFixWidth<DIM> shape) const 00430 { 00431 AutoDiff<DIM> adp[DIM]; 00432 00433 for (int i = 0; i < DIM; i++) 00434 adp[i].Value() = mip.IP()(i); 00435 00436 for (int i = 0; i < DIM; i++) 00437 for (int j = 0; j < DIM; j++) 00438 adp[i].DValue(j) = mip.GetJacobianInverse()(i,j); 00439 00440 HCurlShapeAssign<DIM> ds(shape); 00441 FEL::T_CalcShape (adp, ds); 00442 } 00443 00444 00445 virtual void CalcCurlShape (const IntegrationPoint & ip, 00446 FlatMatrixFixWidth<DIM_CURL> curlshape) const 00447 { 00448 AutoDiff<DIM> adp[DIM]; 00449 for (int i = 0; i < DIM; i++) 00450 adp[i] = AutoDiff<DIM> (ip(i), i); 00451 00452 HCurlCurlShapeAssign<DIM> ds(curlshape); 00453 FEL::T_CalcShape (adp, ds); 00454 } 00455 00456 virtual void 00457 CalcMappedCurlShape (const MappedIntegrationPoint<DIM,DIM> & mip, 00458 FlatMatrixFixWidth<DIM_CURL> curlshape) const 00459 { 00460 AutoDiff<DIM> adp[DIM]; 00461 00462 for (int i = 0; i < DIM; i++) 00463 adp[i].Value() = mip.IP()(i); 00464 00465 for (int i = 0; i < DIM; i++) 00466 for (int j = 0; j < DIM; j++) 00467 adp[i].DValue(j) = mip.GetJacobianInverse()(i,j); 00468 00469 HCurlCurlShapeAssign<DIM> ds(curlshape); 00470 FEL::T_CalcShape (adp, ds); 00471 } 00472 00473 virtual Vec <DIM_CURL> 00474 EvaluateCurlShape (const IntegrationPoint & ip, 00475 FlatVector<double> x, 00476 LocalHeap & lh) const 00477 { 00478 AutoDiff<DIM> adp[DIM]; 00479 for (int i = 0; i < DIM; i++) 00480 adp[i] = AutoDiff<DIM> (ip(i), i); 00481 00482 HCurlEvaluateCurl<DIM> ds(x); 00483 FEL::T_CalcShape (adp, ds); 00484 return ds.Sum(); 00485 } 00486 00487 }; 00488 00489 00490 00491 00492 template <int D> 00493 extern void ComputeGradientMatrix (const ScalarFiniteElement<D> & h1fe, 00494 const HCurlFiniteElement<D> & hcurlfe, 00495 FlatMatrix<> gradient); 00496 00497 } 00498 00499 00500 00501 00502 #endif