QuOc

 

discreteFunction.h

Go to the documentation of this file.
00001 /* Copyright (c) 2001-2010 AG Rumpf, INS, Universitaet Bonn                      */
00002 /*                                                                               */
00003 /* The contents of this file are subject to the terms of the Common Development  */
00004 /* and Distribution License Version 1.0 (the "License"); you may not use         */
00005 /* this file except in compliance with the License. You may obtain a copy of     */
00006 /* the License at http://www.sun.com/cddl/                                       */
00007 /*                                                                               */
00008 /* Software distributed under the License is distributed on an "AS IS" basis,    */
00009 /* WITHOUT WARRANTY OF ANY KIND, either express or implied.                      */
00010 /*                                                                               */
00011 #ifndef __DISCRETEFUNCTION_H
00012 #define __DISCRETEFUNCTION_H
00013 
00014 #include <vec.h>
00015 #include <multiVector.h>
00016 #include <smallMat.h>
00017 #include <pointerClasses.h>
00018 
00019 namespace aol {
00020 
00021 /** Piecewise linear interpolation of discrete data given at not necessarily equidistant sample points
00022  *  \author Schwen
00023  */
00024 template < typename DomType, typename DiscValType, typename ContValType = DiscValType >
00025 class DiscreteValueInterpolator {
00026 protected:
00027   typedef std::map < DomType, DiscValType > MapType;
00028   MapType _vals;
00029 
00030 public:
00031   // destructor, default copy constructor and assignment operators are OK
00032 
00033   DiscreteValueInterpolator ( ) {
00034   }
00035 
00036   DiscreteValueInterpolator ( const aol::Vector<DomType> &domPoints, const aol::Vector<DiscValType> &valPoints ) {
00037     if ( domPoints.size() != valPoints.size() ) {
00038       throw aol::Exception ( "vector sizes do not match", __FILE__, __LINE__ );
00039     }
00040 
00041     for ( int i = 0; i < domPoints.size(); ++i ) {
00042       _vals[ domPoints[i] ] = valPoints[i];
00043     }
00044   }
00045 
00046   ContValType evaluate ( const DomType &point ) const {
00047     if ( point < _vals.begin()->first ) {
00048       throw aol::Exception ( "Value too small, not in interpolatable range.", __FILE__, __LINE__ );
00049     } else if ( point > _vals.rbegin()->first ) {
00050       throw aol::Exception ( "Value too large, not in interpolatable range.", __FILE__, __LINE__ );
00051     }
00052 
00053     const typename MapType::const_iterator found = _vals.find ( point );
00054     if ( found != _vals.end() ) {
00055       return ( static_cast< ContValType > ( found->second ) );
00056     }
00057 
00058     typename MapType::const_iterator ubd1 = _vals.upper_bound ( point ), ubd0 = ubd1; --ubd0;
00059     const double v = static_cast<double> ( point - ubd1->first ) / ( ubd0->first - ubd1->first );
00060 
00061     if ( ( v < 0 ) || ( v > 1 ) ) {
00062       cerr << point << " " << ubd0->first << " " << ubd0->second << " " << ubd1->first << " " << ubd1->second << " " << v << endl;
00063       throw aol::Exception ( "barycentric coordinate outside [0,1]", __FILE__, __LINE__ );
00064     }
00065 
00066     return ( v * ubd0->second + ( 1.0 - v ) * ubd1->second );
00067   }
00068 
00069   void setVals ( const MapType &Vals ) {
00070     _vals = Vals;
00071   }
00072 
00073   void invertIfMonotonic ( ) {
00074     aol::Vector<DiscValType> differences;
00075     for ( typename MapType::const_iterator it = _vals.begin(); it != _vals.end(); ++it ) {
00076       typename MapType::const_iterator itp = it; ++itp;
00077       if ( itp != _vals.end() ) {
00078         differences.pushBack ( itp->second - it->second );
00079       }
00080     }
00081     for ( int i = 0; i < differences.size() - 1; ++i ) {
00082       if ( aol::signum ( differences[i+1] ) * aol::signum ( differences[i] ) != aol::NumberTrait<DiscValType>::one ) {
00083         throw aol::Exception ( "cannot invert DiscreteValueInterpolator that is not strictly monotonic", __FILE__, __LINE__ );
00084       }
00085     }
00086 
00087     MapType inverted;
00088     for ( typename MapType::const_iterator it = _vals.begin(); it != _vals.end(); ++it ) {
00089       inverted[ it->second ] = it->first;
00090     }
00091 
00092     _vals = inverted;
00093   }
00094 
00095   void dump ( ostream &out = cout ) const {
00096     for ( typename MapType::const_iterator it = _vals.begin(); it != _vals.end(); ++it ) {
00097       out << it->first << " " << it->second << endl;
00098     }
00099   }
00100 
00101 };
00102 
00103 
00104 /**
00105  * interface class for discrete functions (FE-functions)
00106  * \author md
00107  */
00108 template <typename ConfiguratorType, typename Imp>
00109 class DiscreteFunctionInterface {
00110 
00111 public:
00112   typedef typename ConfiguratorType::RealType RealType;
00113 
00114 protected:
00115   DeleteFlagPointer<const ConfiguratorType> _configPtr;
00116   const typename ConfiguratorType::InitType &_initializer;
00117 
00118   // barton-nackman
00119   typedef typename ConfiguratorType::ElementType ElementType;
00120   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
00121   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
00122 
00123 public:
00124   //! constructor from grid ( = initializer)
00125   DiscreteFunctionInterface ( const typename ConfiguratorType::InitType &Initializer ) :
00126       _configPtr ( new ConfiguratorType(Initializer), /* deleteFlag = */ true ),
00127       _initializer ( Initializer ) {}
00128 
00129   //! constructor from configurator
00130   DiscreteFunctionInterface ( const ConfiguratorType & config ) :
00131       _configPtr ( &config, /* deleteFlag = */ false ),
00132       _initializer ( config.getInitializer() ) {}
00133 
00134   //! copy constructor
00135   //! We have to make sure that the copy survives, after the copied object is destroyed.
00136   DiscreteFunctionInterface ( const DiscreteFunctionInterface<ConfiguratorType, Imp> & other ) :
00137       _configPtr ( other._configPtr.getDeleteFlag() ? new ConfiguratorType( *(other._configPtr.get()) ) : other._configPtr.get(),
00138                    /* deleteFlag = */ other._configPtr.getDeleteFlag() ),
00139       _initializer ( other._initializer ) {}
00140 
00141   typedef typename ConfiguratorType::VecType  VecType;
00142   typedef typename ConfiguratorType::DomVecType  DomVecType;
00143 
00144   const ConfiguratorType& getConfigurator( ) const { return *_configPtr; }
00145 
00146   /*************************************************************************
00147    ***** interface begins here >>> *****************************************
00148    *************************************************************************/
00149   RealType evaluateAtQuadPoint ( const ElementType &El, int QuadPoint ) const {
00150     return asImp().evaluateAtQuadPoint ( El, QuadPoint );
00151   }
00152 
00153   RealType evaluate ( const ElementType &El, const DomVecType& RefCoord ) const {
00154     return asImp().evaluate ( El, RefCoord );
00155   }
00156 
00157   void evaluateGradientAtQuadPoint ( const ElementType &El, int QuadPoint, VecType& Grad ) const {
00158     asImp().evaluateGradientAtQuadPoint ( El, QuadPoint, Grad );
00159   }
00160 
00161   void evaluateGradient ( const ElementType &El, const DomVecType& RefCoord, VecType& Grad ) const {
00162     asImp().evaluateGradient ( El, RefCoord, Grad );
00163   }
00164 
00165   /*************************************************************************
00166    ***** <<< interface ends here *******************************************
00167    *************************************************************************/
00168 
00169   RealType evaluateGlobal ( const DomVecType& Coord ) const {
00170     DomVecType localCoord;
00171     ElementType el = getConfigurator().getEmptyElement();
00172     getConfigurator().getLocalCoords ( Coord, el, localCoord );
00173 
00174     return evaluate ( el, localCoord );
00175   }
00176 
00177   RealType evaluateGradientGlobal ( const DomVecType& Coord, VecType& Grad ) const {
00178     DomVecType localCoord;
00179     ElementType el = getConfigurator().getEmptyElement();
00180     getConfigurator().getLocalCoords ( Coord, el, localCoord );
00181 
00182     return evaluateGradient ( el, localCoord, Grad );
00183   }
00184 
00185 
00186 };
00187 
00188 /**
00189  * special implementation for discrete functions, that make use of the basefunctionset
00190  * \author md
00191  */
00192 template <typename ConfiguratorType>
00193 class DiscreteFunctionDefault : public DiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType> > {
00194 public:
00195   typedef typename ConfiguratorType::RealType RealType;
00196   typedef typename ConfiguratorType::VecType  VecType;
00197   typedef typename ConfiguratorType::DomVecType  DomVecType;
00198   typedef typename ConfiguratorType::ElementType ElementType;
00199 
00200   const aol::Vector<RealType> & _dofs;
00201 
00202   typedef ConfiguratorType ConfType;
00203 
00204   //! constructor from grid ( = initializer)
00205   DiscreteFunctionDefault ( const typename ConfiguratorType::InitType &Initializer, const aol::Vector<RealType> &Dofs )
00206       : DiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType> > ( Initializer ), _dofs ( Dofs ) {}
00207 
00208   //! constructor from configurator
00209   DiscreteFunctionDefault ( const ConfiguratorType & config, const aol::Vector<RealType> & Dofs )
00210     : DiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType> > ( config ), _dofs ( Dofs ) {}
00211 
00212   // implicit copy constructor does the right thing.
00213   // DiscreteFunctionDefault ( const DiscreteFunctionDefault<ConfiguratorType> &DiscrFunc );
00214 
00215   RealType evaluate ( const ElementType &El, const DomVecType& RefCoord ) const {
00216     const typename ConfiguratorType::BaseFuncSetType &bfs = this->getConfigurator().getBaseFunctionSet ( El );
00217     RealType w = 0.;
00218     for ( int b = 0; b < this->getConfigurator().getNumLocalDofs ( El ); b++ ) {
00219       w += this->_dofs[ this->getConfigurator().localToGlobal ( El, b ) ] * bfs.evaluate ( b, RefCoord );
00220     }
00221     return w;
00222   }
00223 
00224   RealType evaluateAtQuadPoint ( const ElementType &El, int QuadPoint ) const {
00225     const typename ConfiguratorType::BaseFuncSetType &bfs = this->getConfigurator().getBaseFunctionSet ( El );
00226     RealType w = 0.;
00227     for ( int b = 0; b < static_cast<int> ( this->getConfigurator().getNumLocalDofs ( El ) ); b++ ) {
00228       w += this->_dofs[ this->getConfigurator().localToGlobal ( El, b ) ] * bfs.evaluate ( b, QuadPoint );
00229     }
00230     return w;
00231   }
00232 
00233   void evaluateGradientAtQuadPoint ( const ElementType &El, int QuadPoint, aol::Vec<ConfiguratorType::DomDim, RealType>& Grad ) const {
00234     const typename ConfiguratorType::BaseFuncSetType &bfs = this->getConfigurator().getBaseFunctionSet ( El );
00235     VecType v;
00236     Grad.setZero();
00237     for ( int b = 0; b < this->getConfigurator().getNumLocalDofs ( El ); b++ ) {
00238       v = bfs.evaluateGradient ( b, QuadPoint );
00239       v *= this->_dofs[ this->getConfigurator().localToGlobal ( El, b ) ];
00240       Grad += v;
00241     }
00242   }
00243 
00244   void evaluateGradient ( const ElementType &El, const DomVecType& RefCoord, aol::Vec<ConfiguratorType::DomDim, RealType>& Grad ) const {
00245     const typename ConfiguratorType::BaseFuncSetType &bfs = this->getConfigurator().getBaseFunctionSet ( El );
00246     VecType v;
00247     Grad.setZero();
00248 
00249     for ( int b = 0; b < this->getConfigurator().getNumLocalDofs ( El ); b++ ) {
00250       bfs.evaluateGradient ( b, RefCoord, v );
00251       v *= this->_dofs[ this->getConfigurator().localToGlobal ( El, b ) ];
00252       Grad += v;
00253     }
00254   }
00255   const aol::Vector<RealType> & getDofs() const {
00256     return _dofs;
00257   }
00258 };
00259 
00260 /**
00261  * special implementation for discrete functions, that make use of the basefunctionset
00262  * \author md
00263  */
00264 template <typename ConfiguratorType, typename DiscreteFunctionType, typename Imp>
00265 class ComposedDiscreteFunctionInterface : public DiscreteFunctionInterface < ConfiguratorType,
00266       ComposedDiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionType, Imp> > {
00267 public:
00268   typedef typename ConfiguratorType::RealType RealType;
00269   typedef typename ConfiguratorType::VecType  VecType;
00270   typedef typename ConfiguratorType::DomVecType  DomVecType;
00271   typedef typename ConfiguratorType::ElementType ElementType;
00272   typedef ConfiguratorType ConfType;
00273 
00274   inline Imp& asImp() { return static_cast<Imp&> ( *this ); }
00275   inline const Imp& asImp() const { return static_cast<const Imp&> ( *this ); }
00276 
00277   const DiscreteFunctionType &_discreteFunction;
00278 
00279   ComposedDiscreteFunctionInterface ( const typename ConfiguratorType::InitType &Initializer, const DiscreteFunctionType &DiscreteFunction )
00280       : DiscreteFunctionInterface<ConfiguratorType, ComposedDiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionType, Imp> > ( Initializer ),
00281       _discreteFunction ( DiscreteFunction ) {}
00282 
00283   //! copy constructor
00284   ComposedDiscreteFunctionInterface ( const ComposedDiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionType, Imp> &DiscrFunc )
00285       : DiscreteFunctionInterface<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType> > ( DiscrFunc._initializer ),
00286       _discreteFunction ( DiscrFunc._discreteFunction ) {}
00287 
00288   RealType evaluateCompositionFunction ( RealType s ) const {
00289     return asImp().evaluateCompositionFunction ( s );
00290   }
00291 
00292   RealType evaluateCompositionFunctionDerivative ( RealType s ) const {
00293     return asImp().evaluateCompositionFunctionDerivative ( s );
00294   }
00295 
00296   RealType evaluateGlobal ( const DomVecType& Coord ) const {
00297     return evaluateCompositionFunction ( _discreteFunction.evaluateGlobal ( Coord ) );
00298   }
00299 
00300   RealType evaluate ( const ElementType &El, const DomVecType& RefCoord ) const {
00301     return evaluateCompositionFunction ( _discreteFunction.evaluate ( El, RefCoord ) );
00302   }
00303 
00304   RealType evaluateAtQuadPoint ( const ElementType &El, int QuadPoint ) const {
00305     return evaluateCompositionFunction ( _discreteFunction.evaluate ( El, QuadPoint ) );
00306   }
00307 
00308   void evaluateGradientAtQuadPoint ( const ElementType &El, int QuadPoint, aol::Vec<ConfiguratorType::DomDim, RealType>& Grad ) const {
00309     _discreteFunction.evaluateGradientAtQuadPointWeight ( El, QuadPoint, Grad );
00310     Grad *= evaluateCompositionFunctionDerivative ( _discreteFunction.evaluateAtQuadPoint ( El, QuadPoint ) );
00311   }
00312 
00313   void evaluateGradient ( const ElementType &El, const DomVecType& RefCoord, aol::Vec<ConfiguratorType::DomDim, RealType>& Grad ) const {
00314     _discreteFunction.evaluateGradient ( El, RefCoord, Grad );
00315     Grad *= evaluateCompositionFunctionDerivative ( _discreteFunction.evaluate ( El, RefCoord ) );
00316   }
00317 
00318 
00319 };
00320 
00321 
00322 
00323 
00324 
00325 template <typename ConfiguratorType, typename DiscFuncType, int DimRange>
00326 class DiscreteVectorFunction {
00327 public:
00328   typedef typename ConfiguratorType::RealType RealType;
00329   typedef typename ConfiguratorType::VecType  VecType;
00330   typedef typename ConfiguratorType::DomVecType  DomVecType;
00331   typedef typename ConfiguratorType::ElementType ElementType;
00332 
00333   mutable aol::auto_container<DimRange, DiscFuncType> _discrFuncs;
00334   const aol::MultiVector<RealType> & _dofsReference;
00335 protected:
00336 public:
00337   //! constructor from grid
00338   DiscreteVectorFunction ( const typename ConfiguratorType::InitType &Initializer,
00339                            const aol::MultiVector<RealType> &Dofs,
00340                            bool allowMoreComponentsThanNeeded ) :
00341       _dofsReference ( Dofs ) {
00342     if ( !allowMoreComponentsThanNeeded ) {
00343       if ( Dofs.numComponents() != DimRange )
00344         throw aol::Exception ( "you have to pass a multivector with the correct amount of components.", __FILE__, __LINE__ );
00345     } else {
00346       if ( Dofs.numComponents() < DimRange )
00347         throw aol::Exception ( "you have to pass a multivector with at least as much components as the discrete function has.", __FILE__, __LINE__ );
00348     }
00349     for ( int c = 0; c < DimRange; c++ ) {
00350       _discrFuncs.take_over ( c, new DiscFuncType ( Initializer, Dofs[c] ) );
00351     }
00352   }
00353 
00354   //! constructor from configurator
00355   DiscreteVectorFunction ( const ConfiguratorType &Configurator,
00356                            const aol::MultiVector<RealType> &Dofs,
00357                            bool allowMoreComponentsThanNeeded ) :
00358       _dofsReference ( Dofs ) {
00359     if ( !allowMoreComponentsThanNeeded ) {
00360       if ( Dofs.numComponents() != DimRange )
00361         throw aol::Exception ( "you have to pass a multivector with the correct amount of components.", __FILE__, __LINE__ );
00362     } else {
00363       if ( Dofs.numComponents() < DimRange )
00364         throw aol::Exception ( "you have to pass a multivector with at least as much components as the discrete function has.", __FILE__, __LINE__ );
00365     }
00366     for ( int c = 0; c < DimRange; c++ ) {
00367       _discrFuncs.take_over ( c, new DiscFuncType ( Configurator, Dofs[c] ) );
00368     }
00369   }
00370 
00371 
00372   void evaluate ( const ElementType &El, const DomVecType& RefCoord, aol::Vec<DimRange, RealType> &Value ) const {
00373     for ( int c = 0; c < DimRange; c++ ) {
00374       Value[c] = _discrFuncs[c].evaluate ( El, RefCoord );
00375     }
00376   }
00377 
00378   void evaluateAtQuadPoint ( const ElementType &El, int QuadPoint, aol::Vec<DimRange, RealType> &Value ) const {
00379     for ( int c = 0; c < DimRange; c++ ) {
00380       Value[c] = _discrFuncs[c].evaluateAtQuadPoint ( El, QuadPoint );
00381     }
00382   }
00383 
00384   void evaluateGradient ( const ElementType &El, const DomVecType& RefCoord, Mat<DimRange, ConfiguratorType::DomDim, RealType> &Mat ) const {
00385     VecType v;
00386     for ( int c = 0; c < DimRange; c++ ) {
00387       _discrFuncs[c].evaluateGradient ( El, RefCoord, v );
00388       Mat.setRow ( c, v );
00389     }
00390   }
00391 
00392   void evaluateGradientAtQuadPoint ( const ElementType &El, int QuadPoint, Mat<DimRange, ConfiguratorType::DomDim, RealType> &Mat ) const {
00393     VecType v;
00394     for ( int c = 0; c < DimRange; c++ ) {
00395       _discrFuncs[c].evaluateGradientAtQuadPoint ( El, QuadPoint, v );
00396       Mat.setRow ( c, v );
00397     }
00398   }
00399 
00400 
00401   const DiscFuncType& operator[] ( int i ) const {
00402     return _discrFuncs[i];
00403   }
00404 
00405   DiscFuncType& operator[] ( int i ) {
00406     return _discrFuncs[i];
00407   }
00408 
00409 
00410   const aol::MultiVector<RealType> & getDofs() const {
00411     return _dofsReference;
00412   }
00413 };
00414 
00415 
00416 template <typename ConfiguratorType, int DimRange>
00417 class DiscreteVectorFunctionDefault : public DiscreteVectorFunction<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType>, DimRange > {
00418 
00419 public:
00420   typedef typename ConfiguratorType::RealType RealType;
00421 
00422   //! constructor from grid
00423   DiscreteVectorFunctionDefault ( const typename ConfiguratorType::InitType &Initializer,
00424                                   const aol::MultiVector<RealType> &Dofs,
00425                                   bool allowMoreDofsThanNeeded = false )
00426       : DiscreteVectorFunction<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType>, DimRange > ( Initializer, Dofs, allowMoreDofsThanNeeded ) { }
00427 
00428   //! constructor from configurator
00429   DiscreteVectorFunctionDefault ( const ConfiguratorType &Configurator,
00430                                   const aol::MultiVector<RealType> &Dofs,
00431                                   bool allowMoreDofsThanNeeded = false )
00432       : DiscreteVectorFunction<ConfiguratorType, DiscreteFunctionDefault<ConfiguratorType>, DimRange > ( Configurator, Dofs, allowMoreDofsThanNeeded ) { }
00433 };
00434 
00435 } // end namespace aol
00436 
00437 #endif

Generated on Fri Sep 9 2011 21:07:39 for QuocMesh by doxygen 1.7.1