NGSolve
4.9
|
00001 #ifndef FILE_COEFFICIENT 00002 #define FILE_COEFFICIENT 00003 00004 /*********************************************************************/ 00005 /* File: coefficient.hh */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 25. Mar. 2000 */ 00008 /*********************************************************************/ 00009 00010 namespace ngfem 00011 { 00012 00018 class NGS_DLL_HEADER CoefficientFunction 00019 { 00020 public: 00022 CoefficientFunction (); 00024 virtual ~CoefficientFunction (); 00025 00027 virtual int NumRegions () { return INT_MAX; } 00029 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const = 0; 00030 00032 virtual void Evaluate (const BaseMappedIntegrationRule & ir, FlatMatrix<double> values) const; 00033 00034 virtual void Evaluate (const BaseMappedIntegrationRule & ir, FlatMatrix<Complex> values) const; 00035 00036 00038 virtual Complex EvaluateComplex (const BaseMappedIntegrationPoint & ip) const 00039 { 00040 return Evaluate (ip); 00041 } 00042 00043 00044 template <typename SCAL> 00045 inline SCAL T_Evaluate (const BaseMappedIntegrationPoint & ip) const 00046 { 00047 return SCAL (Evaluate (ip)); 00048 } 00049 00050 00051 virtual double Evaluate (const BaseMappedIntegrationPoint & ip, const double & t) const 00052 { 00053 return Evaluate(ip); 00054 } 00055 virtual double EvaluateDeri (const BaseMappedIntegrationPoint & ip, const double & t) const 00056 { 00057 return 0; 00058 } 00059 00060 // to be changed 00061 virtual double Evaluate (const BaseMappedIntegrationPoint & ip, 00062 const complex<double> & t) const 00063 { return Evaluate(ip,t.real()); } 00064 // to be changed 00065 virtual double EvaluateDeri (const BaseMappedIntegrationPoint & ip, 00066 const complex<double> & t) const 00067 { return EvaluateDeri(ip,t.real()); } 00068 00069 00070 virtual double EvaluateConst () const 00071 { 00072 throw Exception (string ("EvaluateConst called for non-const coefficient function ")+ 00073 typeid(*this).name()); 00074 } 00075 00076 virtual bool IsComplex() const { return false; } 00077 virtual int Dimension() const { return 1; } 00078 00079 virtual void Evaluate(const BaseMappedIntegrationPoint & ip, 00080 FlatVector<> result) const 00081 { 00082 double f = Evaluate (ip); 00083 result(0) = f; 00084 } 00085 00086 virtual void Evaluate(const BaseMappedIntegrationPoint & ip, 00087 FlatVector<Complex> result) const 00088 { 00089 Complex f = EvaluateComplex (ip); 00090 result(0) = f; 00091 } 00092 00093 00094 virtual void PrintReport (ostream & ost) 00095 { 00096 ost << "Base-Class CoefficientFunction" << endl; 00097 } 00098 }; 00099 00100 00101 template <> 00102 inline double CoefficientFunction :: 00103 T_Evaluate<double> (const BaseMappedIntegrationPoint & ip) const 00104 { 00105 return Evaluate (ip); 00106 } 00107 00108 template <> 00109 inline Complex CoefficientFunction :: 00110 T_Evaluate<Complex> (const BaseMappedIntegrationPoint & ip) const 00111 { 00112 return EvaluateComplex (ip); 00113 } 00114 00115 00116 00117 /* 00118 template <int S, int R> 00119 inline double Evaluate (const CoefficientFunction & fun, 00120 const MappedIntegrationPoint<S,R> & ip) 00121 { 00122 return fun.Evaluate(ip); 00123 } 00124 */ 00125 00126 inline double Evaluate (const CoefficientFunction & fun, 00127 const BaseMappedIntegrationPoint & ip) 00128 { 00129 return fun.Evaluate(ip); 00130 } 00131 00132 00133 00134 00136 class NGS_DLL_HEADER ConstantCoefficientFunction : public CoefficientFunction 00137 { 00139 double val; 00140 public: 00142 ConstantCoefficientFunction (double aval); 00144 virtual ~ConstantCoefficientFunction (); 00146 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const 00147 { 00148 return val; 00149 } 00150 00151 virtual double EvaluateConst () const 00152 { 00153 return val; 00154 } 00155 }; 00156 00157 00158 00160 class NGS_DLL_HEADER DomainConstantCoefficientFunction : public CoefficientFunction 00161 { 00163 Array<double> val; 00164 public: 00166 DomainConstantCoefficientFunction (const Array<double> & aval); 00168 virtual int NumRegions () { return val.Size(); } 00170 virtual ~DomainConstantCoefficientFunction (); 00172 00173 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const; 00174 00175 virtual void Evaluate (const BaseMappedIntegrationRule & ir, FlatMatrix<double> values) const; 00176 00177 virtual void Evaluate (const BaseMappedIntegrationRule & ir, FlatMatrix<Complex> values) const; 00178 00179 00180 virtual double EvaluateConst () const 00181 { 00182 return val[0]; 00183 } 00184 00185 double operator[] (int i) const { return val[i]; } 00186 }; 00187 00188 00189 00190 00192 template <int DIM> 00193 class NGS_DLL_HEADER DomainVariableCoefficientFunction : public CoefficientFunction 00194 { 00195 Array<EvalFunction*> fun; 00196 Array<CoefficientFunction*> depends_on; 00197 int numarg; 00198 public: 00200 DomainVariableCoefficientFunction (const Array<EvalFunction*> & afun); 00201 DomainVariableCoefficientFunction (const Array<EvalFunction*> & afun, 00202 const Array<CoefficientFunction*> & adepends_on); 00203 00205 virtual ~DomainVariableCoefficientFunction (); 00207 virtual int NumRegions () { return (fun.Size() == 1) ? INT_MAX : fun.Size(); } 00209 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const; 00210 00211 virtual Complex EvaluateComplex (const BaseMappedIntegrationPoint & ip) const; 00212 00213 EvalFunction & GetEvalFunction(const int index) 00214 { 00215 return *(fun[index]); 00216 } 00217 00218 virtual bool IsComplex() const 00219 { 00220 for (int i = 0; i < fun.Size(); i++) 00221 if (fun[i]->IsResultComplex()) return true; 00222 return false; 00223 } 00224 virtual int Dimension() const { return fun[0]->Dimension(); } 00225 00226 virtual void Evaluate(const BaseMappedIntegrationPoint & ip, 00227 FlatVector<> result) const; 00228 00229 virtual void Evaluate(const BaseMappedIntegrationPoint & ip, 00230 FlatVector<Complex> result) const; 00231 00232 virtual void Evaluate (const BaseMappedIntegrationRule & ir, 00233 FlatMatrix<double> values) const; 00234 }; 00235 00237 template <int DIM> 00238 class NGS_DLL_HEADER DomainInternalCoefficientFunction : public CoefficientFunction 00239 { 00241 int matnr; 00243 double (*f)(const double*); 00244 public: 00246 DomainInternalCoefficientFunction (int amatnr, double (*af)(const double*)) 00247 : matnr(amatnr), f(af) { ; } 00249 virtual ~DomainInternalCoefficientFunction () { ; } 00251 /* 00252 template <int S, int R> 00253 double Evaluate (const MappedIntegrationPoint<S,R> & ip) 00254 { 00255 int elind = ip.GetTransformation().GetElementIndex(); 00256 if (elind != matnr && matnr > 0) return 0; 00257 00258 return f(&ip.GetPoint()(0)); 00259 } 00260 */ 00261 00262 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const 00263 { 00264 int elind = ip.GetTransformation().GetElementIndex(); 00265 if (elind != matnr && matnr > 0) return 0; 00266 00267 return f(&static_cast<const DimMappedIntegrationPoint<DIM>&>(ip).GetPoint()(0)); 00268 // return f(&ip.GetPoint().REval(0)); 00269 } 00270 00271 }; 00272 00273 00274 00275 00276 00282 class IntegrationPointCoefficientFunction : public CoefficientFunction 00283 { 00284 int elems, ips_per_elem; 00286 Array<double> values; 00287 public: 00289 IntegrationPointCoefficientFunction (int aelems, int size) 00290 : elems(aelems), ips_per_elem(size), values(aelems*size) { ; } 00292 IntegrationPointCoefficientFunction (int aelems, int size, double val) 00293 : elems(aelems), ips_per_elem(size), values(aelems*size) 00294 { 00295 values = val; 00296 } 00298 IntegrationPointCoefficientFunction (int aelems, int size, Array<double> & avalues) 00299 : elems(aelems), ips_per_elem(size), values(avalues) 00300 { 00301 if ( avalues.Size() < aelems * size ) 00302 { 00303 cout << "Warning: IntegrationPointCoefficientFunction, constructor: sizes don't match!" << endl; 00304 values.SetSize(aelems*size); 00305 } 00306 } 00307 00309 virtual ~IntegrationPointCoefficientFunction () { ; } 00311 /* 00312 template <int S, int R> 00313 double Evaluate (const MappedIntegrationPoint<S,R> & ip) 00314 { 00315 int ipnr = ip.GetIPNr(); 00316 int elnr = ip.GetTransformation().GetElementNr(); 00317 00318 if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 00319 { 00320 ostringstream ost; 00321 ost << "IntegrationPointCoefficientFunction: ip = " 00322 << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 00323 << ips_per_elem << "/ 0 - " << elems << "!" << endl; 00324 throw Exception (ost.str()); 00325 } 00326 00327 return values[elnr*ips_per_elem+ipnr]; 00328 } 00329 */ 00330 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const 00331 { 00332 int ipnr = ip.GetIPNr(); 00333 int elnr = ip.GetTransformation().GetElementNr(); 00334 00335 if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 00336 { 00337 ostringstream ost; 00338 ost << "IntegrationPointCoefficientFunction: ip = " 00339 << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 00340 << ips_per_elem << "/ 0 - " << elems << "!" << endl; 00341 throw Exception (ost.str()); 00342 } 00343 00344 return values[elnr*ips_per_elem+ipnr]; 00345 } 00346 00347 00348 // direct access to the values at the integration points 00349 double & operator() (int elnr, int ipnr) 00350 { 00351 if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 00352 { 00353 ostringstream ost; 00354 ost << "IntegrationPointCoefficientFunction: ip = " 00355 << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 00356 << ips_per_elem-1 << " / 0 - " << elems-1 << "!" << endl; 00357 throw Exception (ost.str()); 00358 } 00359 00360 return values[elnr*ips_per_elem+ipnr]; 00361 } 00362 00363 double operator() (int elnr, int ipnr) const 00364 { 00365 if ( ipnr < 0 || ipnr >= ips_per_elem || elnr < 0 || elnr >= elems ) 00366 { 00367 ostringstream ost; 00368 ost << "IntegrationPointCoefficientFunction: ip = " 00369 << ipnr << " / elem = " << elnr << ". Ranges: 0 - " 00370 << ips_per_elem-1 << " / 0 - " << elems-1 << "!" << endl; 00371 throw Exception (ost.str()); 00372 } 00373 00374 return values[elnr*ips_per_elem+ipnr]; 00375 } 00376 00377 00378 00379 int GetNumIPs() const { return ips_per_elem; } 00380 int GetNumElems() const { return elems; } 00381 00382 00383 00384 void ReSetValues( Array<double> & avalues ) 00385 { 00386 if ( avalues.Size() < values.Size() ) 00387 { 00388 throw Exception("IntegrationPointCoefficientFunction::ReSetValues - sizes don't match!"); 00389 } 00390 values = avalues; 00391 } 00392 00393 }; 00394 00395 00397 class PolynomialCoefficientFunction : public CoefficientFunction 00398 { 00399 private: 00400 Array < Array< Array<double>* >* > polycoeffs; 00401 Array < Array<double>* > polybounds; 00402 00403 private: 00404 double EvalPoly(const double t, const Array<double> & coeffs) const; 00405 double EvalPolyDeri(const double t, const Array<double> & coeffs) const; 00406 00407 public: 00408 PolynomialCoefficientFunction(const Array < Array<double>* > & polycoeffs_in); 00409 PolynomialCoefficientFunction(const Array < Array< Array<double>* >* > & polycoeffs_in, const Array < Array<double>* > & polybounds_in); 00410 00411 virtual ~PolynomialCoefficientFunction(); 00412 00413 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const; 00414 00415 virtual double Evaluate (const BaseMappedIntegrationPoint & ip, const double & t) const; 00416 00417 virtual double EvaluateDeri (const BaseMappedIntegrationPoint & ip, const double & t) const; 00418 00419 virtual double EvaluateConst () const; 00420 }; 00421 00422 00423 00425 00426 class FileCoefficientFunction : public CoefficientFunction 00427 { 00428 private: 00429 Array < Array < double > * > ValuesAtIps; 00430 00431 ofstream outfile; 00432 00433 string valuesfilename; 00434 string infofilename; 00435 string ipfilename; 00436 00437 int maxelnum, maxipnum, totalipnum; 00438 00439 bool writeips; 00440 00441 private: 00442 void EmptyValues(void); 00443 00444 public: 00445 FileCoefficientFunction(); 00446 00447 FileCoefficientFunction(const string & filename); 00448 00449 FileCoefficientFunction(const string & aipfilename, 00450 const string & ainfofilename, 00451 const string & avaluesfilename, 00452 const bool loadvalues = false); 00453 00454 virtual ~FileCoefficientFunction(); 00455 00456 virtual double Evaluate (const BaseMappedIntegrationPoint & ip) const; 00457 00458 void LoadValues(const string & filename); 00459 inline void LoadValues(void){ LoadValues(valuesfilename); } 00460 00461 void StartWriteIps(const string & filename); 00462 inline void StartWriteIps(void){ StartWriteIps(ipfilename); } 00463 00464 void StopWriteIps(const string & infofilename); 00465 inline void StopWriteIps(void){ StopWriteIps(infofilename); } 00466 00467 void Reset(void); 00468 00469 00470 00471 00472 }; 00473 00474 } 00475 00476 00477 #endif