NGSolve  4.9
fem/coefficient.hpp
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