NGSolve
4.9
|
00001 #ifndef FILE_TSCALARFE 00002 #define FILE_TSCALARFE 00003 00004 /*********************************************************************/ 00005 /* File: tscalarfe.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 25. Mar. 2000 */ 00008 /*********************************************************************/ 00009 00010 namespace ngfem 00011 { 00012 00013 00017 template <int DIM> 00018 class DShapeElement 00019 { 00020 FlatVec<DIM> data; 00021 public: 00023 DShapeElement (FlatVec<DIM> adata) : data(adata) { ; } 00024 00026 void operator= (AutoDiff<DIM> ad) 00027 { 00028 for (int i = 0; i < DIM; i++) 00029 data(i) = ad.DValue(i); 00030 } 00031 }; 00032 00033 00037 template <int DIM> 00038 class DShapeAssign 00039 { 00040 FlatMatrixFixWidth<DIM> dshape; 00041 public: 00043 DShapeAssign (FlatMatrixFixWidth<DIM> mat) : dshape(mat) { ; } 00044 00046 DShapeElement<DIM> operator[] (int i) const 00047 { return DShapeElement<DIM> (dshape.Row(i)); } 00048 00050 const DShapeAssign Addr (int i) const 00051 { return DShapeAssign (dshape.Rows(i, dshape.Height())); } 00052 }; 00053 00054 00055 00056 00057 00061 class EvaluateShapeElement 00062 { 00063 double coef; 00064 double * sum; 00065 public: 00067 EvaluateShapeElement (double acoef, double * asum) 00068 : coef(acoef), sum(asum) { ; } 00069 00071 void operator= (double val) 00072 { 00073 *sum += coef * val; 00074 } 00075 }; 00076 00077 00081 class EvaluateShape 00082 { 00083 const double * coefs; 00084 double * sum; 00085 public: 00087 EvaluateShape (FlatVector<> acoefs, double * asum) 00088 : coefs(&acoefs(0)), sum(asum) { ; } 00089 00091 EvaluateShape (const double * acoefs, double * asum) 00092 : coefs(acoefs), sum(asum) { ; } 00093 00095 EvaluateShapeElement operator[] (int i) const 00096 { return EvaluateShapeElement (coefs[i], sum); } 00097 00099 const EvaluateShape Addr (int i) const 00100 { return EvaluateShape (coefs+i, sum); } 00101 }; 00102 00103 00104 00105 00107 class EvaluateShapeTransElement 00108 { 00109 double & data; 00110 const double & fac; 00111 public: 00113 EvaluateShapeTransElement (double & adata, const double & afac) 00114 : data(adata), fac(afac) { ; } 00115 00117 void operator= (double ad) 00118 { 00119 data += ad * fac; 00120 } 00121 }; 00122 00123 class EvaluateShapeTrans 00124 { 00125 double * coefs; 00126 const double & fac; 00127 public: 00128 EvaluateShapeTrans (FlatVector<> acoefs, const double & afac) 00129 : coefs(&acoefs(0)), fac(afac) { ; } 00130 00131 EvaluateShapeTrans (double * acoefs, const double & afac) 00132 : coefs(acoefs), fac(afac) { ; } 00133 00134 EvaluateShapeTransElement operator[] (int i) const 00135 { return EvaluateShapeTransElement (coefs[i], fac); } 00136 00137 const EvaluateShapeTrans Addr (int i) const 00138 { return EvaluateShapeTrans (coefs+i, fac); } 00139 }; 00140 00141 00142 00143 00144 00145 00146 00147 00148 00149 template <int DIM> 00150 class EvaluateDShapeElement 00151 { 00152 double data; 00153 Vec<DIM> & sum; 00154 public: 00155 EvaluateDShapeElement (double adata, Vec<DIM> & asum) : data(adata), sum(asum) { ; } 00156 void operator= (AutoDiff<DIM> ad) 00157 { 00158 for (int i = 0; i < DIM; i++) 00159 sum(i) += ad.DValue(i) * data; 00160 } 00161 }; 00162 00163 template <int DIM> 00164 class EvaluateDShape 00165 { 00166 double * coefs; 00167 Vec<DIM> & sum; 00168 public: 00169 EvaluateDShape (FlatVector<> acoefs, Vec<DIM> & asum) 00170 : coefs(&acoefs(0)), sum(asum) { ; } 00171 00172 EvaluateDShape (double * acoefs, Vec<DIM> & asum) 00173 : coefs(acoefs), sum(asum) { ; } 00174 00175 EvaluateDShapeElement<DIM> operator[] (int i) const 00176 { return EvaluateDShapeElement<DIM> (coefs[i], sum); } 00177 00178 const EvaluateDShape Addr (int i) const 00179 { return EvaluateDShape (coefs+i, sum); } 00180 }; 00181 00182 00183 00184 00185 00186 00187 00188 template <int DIM> 00189 class EvaluateDShapeTransElement 00190 { 00191 double & data; 00192 const Vec<DIM> & fac; 00193 public: 00194 EvaluateDShapeTransElement (double & adata, const Vec<DIM> & afac) : data(adata), fac(afac) { ; } 00195 void operator= (AutoDiff<DIM> ad) 00196 { 00197 for (int i = 0; i < DIM; i++) 00198 data += ad.DValue(i) * fac(i); 00199 } 00200 }; 00201 00203 template <int DIM> 00204 class EvaluateDShapeTrans 00205 { 00206 double * coefs; 00207 const Vec<DIM> & fac; 00208 public: 00209 EvaluateDShapeTrans (FlatVector<> acoefs, const Vec<DIM> & afac) 00210 : coefs(&acoefs(0)), fac(afac) { ; } 00211 00212 EvaluateDShapeTrans (double * acoefs, const Vec<DIM> & afac) 00213 : coefs(acoefs), fac(afac) { ; } 00214 00215 EvaluateDShapeTransElement<DIM> operator[] (int i) const 00216 { return EvaluateDShapeTransElement<DIM> (coefs[i], fac); } 00217 00218 const EvaluateDShapeTrans Addr (int i) const 00219 { return EvaluateDShapeTrans (coefs+i, fac); } 00220 }; 00221 00222 00223 00224 00225 00226 00227 00232 template <class FEL, ELEMENT_TYPE ET, int NDOF, int ORDER> 00233 class T_ScalarFiniteElement : public ScalarFiniteElement<ET_trait<ET>::DIM> 00234 { 00235 00236 public: 00237 00238 protected: 00239 enum { DIM = ET_trait<ET>::DIM }; 00240 00241 T_ScalarFiniteElement () 00242 : ScalarFiniteElement<DIM> (ET, NDOF, ORDER) { ; } 00243 00244 virtual ~T_ScalarFiniteElement() { ; } 00245 00246 public: 00247 00248 virtual void CalcShape (const IntegrationPoint & ip, 00249 FlatVector<> shape) const 00250 { 00251 // double pt[DIM]; 00252 Vec<DIM> pt; 00253 for (int i = 0; i < DIM; i++) pt[i] = ip(i); 00254 FEL::T_CalcShape (&pt(0), shape); 00255 } 00256 00257 virtual double 00258 Evaluate (const IntegrationPoint & ip, FlatVector<double> x) const 00259 { 00260 // double pt[DIM]; 00261 Vec<DIM> pt; 00262 for (int i = 0; i < DIM; i++) pt[i] = ip(i); 00263 00264 double sum = 0.0; 00265 EvaluateShape eval(x, &sum); 00266 00267 FEL::T_CalcShape (&pt(0), eval); 00268 return sum; 00269 } 00270 00271 virtual void 00272 Evaluate (const IntegrationRule & ir, FlatVector<double> coefs, FlatVector<double> vals) const 00273 { 00274 // double pt[DIM]; 00275 Vec<DIM> pt; 00276 for (int i = 0; i < ir.GetNIP(); i++) 00277 { 00278 for (int j = 0; j < DIM; j++) pt[j] = ir[i](j); 00279 00280 vals(i) = 0.0; 00281 EvaluateShape eval(coefs, &vals(i)); 00282 FEL::T_CalcShape (&pt(0), eval); 00283 } 00284 } 00285 00286 00287 00288 static void CalcShapeStat (const IntegrationPoint & ip, 00289 FlatVector<> shape) 00290 { 00291 // double pt[DIM]; 00292 Vec<DIM> pt; 00293 for (int i = 0; i < DIM; i++) pt[i] = ip(i); 00294 FEL::T_CalcShape (&pt(0), shape); 00295 } 00296 00297 virtual void CalcDShape (const IntegrationPoint & ip, 00298 FlatMatrixFixWidth<DIM> dshape) const 00299 { 00300 // AutoDiff<DIM> adp[DIM]; 00301 Vec<DIM, AutoDiff<DIM> > adp; 00302 for (int i = 0; i < DIM; i++) 00303 adp[i] = AutoDiff<DIM> (ip(i), i); 00304 00305 DShapeAssign<DIM> ds(dshape); 00306 FEL::T_CalcShape (&adp(0), ds); 00307 } 00308 00309 static void CalcDShapeStat (const IntegrationPoint & ip, 00310 FlatMatrixFixWidth<DIM> dshape) 00311 { 00312 // AutoDiff<DIM> adp[DIM]; 00313 Vec<DIM, AutoDiff<DIM> > adp; 00314 for (int i = 0; i < DIM; i++) 00315 adp[i] = AutoDiff<DIM> (ip(i), i); 00316 00317 DShapeAssign<DIM> ds(dshape); 00318 FEL::T_CalcShape (&adp(0), ds); 00319 } 00320 00321 virtual void 00322 CalcMappedDShape (const MappedIntegrationPoint<DIM,DIM> & sip, 00323 FlatMatrixFixWidth<DIM> dshape) const 00324 { 00325 // AutoDiff<DIM> adp[DIM]; 00326 Vec<DIM, AutoDiff<DIM> > adp; 00327 for (int i = 0; i < DIM; i++) 00328 adp[i].Value() = sip.IP()(i); 00329 00330 for (int i = 0; i < DIM; i++) 00331 for (int j = 0; j < DIM; j++) 00332 adp[i].DValue(j) = sip.GetJacobianInverse()(i,j); 00333 00334 DShapeAssign<DIM> ds(dshape); 00335 FEL::T_CalcShape (&adp(0), ds); 00336 } 00337 }; 00338 00339 00340 00341 00342 00343 00344 00345 00346 00347 00348 00349 00350 00351 00352 00353 00354 00355 00356 00357 00358 00359 00360 00361 00362 00363 00364 00365 00371 template <class FEL, ELEMENT_TYPE ET> 00372 class NGS_DLL_HEADER T_ScalarFiniteElement2 : virtual public ScalarFiniteElement<ET_trait<ET>::DIM> 00373 { 00374 public: 00375 enum { DIM = ET_trait<ET>::DIM }; 00376 using ScalarFiniteElement<DIM>::eltype; 00377 00378 T_ScalarFiniteElement2 () { eltype = ET; } 00379 00380 virtual void CalcShape (const IntegrationPoint & ip, 00381 FlatVector<> shape) const; 00382 00383 virtual double Evaluate (const IntegrationPoint & ip, FlatVector<double> x) const; 00384 00385 virtual void Evaluate (const IntegrationRule & ir, FlatVector<double> coefs, FlatVector<double> vals) const; 00386 00387 virtual void EvaluateTrans (const IntegrationRule & ir, FlatVector<> vals, FlatVector<double> coefs) const; 00388 00389 virtual void EvaluateGrad (const IntegrationRule & ir, FlatVector<double> coefs, FlatMatrixFixWidth<DIM> vals) const; 00390 00391 virtual void EvaluateGradTrans (const IntegrationRule & ir, FlatMatrixFixWidth<DIM> vals, FlatVector<double> coefs) const; 00392 00393 virtual void CalcDShape (const IntegrationPoint & ip, 00394 FlatMatrixFixWidth<DIM> dshape) const; 00395 00396 virtual void 00397 CalcMappedDShape (const MappedIntegrationPoint<DIM,DIM> & sip, 00398 FlatMatrixFixWidth<DIM> dshape) const; 00399 00400 00401 00402 }; 00403 00404 template <typename TFA> 00405 inline void SetZero (TFA & shape, int first, int next) 00406 { 00407 for (int i = first; i < next; i++) 00408 shape[i] = 0.0; 00409 } 00410 00411 inline void SetZero (EvaluateShape & shape, int first, int next) 00412 { 00413 ; 00414 } 00415 00416 inline void SetZero (EvaluateShapeTrans & shape, int first, int next) 00417 { 00418 ; 00419 } 00420 00421 00422 } 00423 00424 00425 #endif