dune-localfunctions  2.2.0
field.hh
Go to the documentation of this file.
00001 #ifndef DUNE_LOCALFUNCTIONS_UTILITY_FIELD_HH
00002 #define DUNE_LOCALFUNCTIONS_UTILITY_FIELD_HH
00003 
00004 #if HAVE_ALGLIB
00005 #include <alglib/amp.h>
00006 #warning ALGLIB support is deprecated and will be dropped after DUNE 2.2 (cf. FS#931)
00007 #endif
00008 
00009 #include <dune/common/gmpfield.hh>
00010 #include <dune/common/fvector.hh>
00011 #include <dune/common/fmatrix.hh>
00012 
00013 namespace Dune
00014 {
00015 
00016   // Unity
00017   // -----
00018 
00029   template< class Field >
00030   struct Unity
00031   {
00032     operator Field () const
00033     {
00034       return Field( 1 );
00035     }
00036   };
00037 
00038   template< class Field >
00039   Field operator+ ( const Unity< Field > &u, const Field &f )
00040   {
00041     return (Field)u + f;
00042   }
00043 
00044   template< class Field >
00045   Field operator- ( const Unity< Field > &u, const Field &f )
00046   {
00047     return (Field)u - f;
00048   }
00049 
00050   template< class Field >
00051   Field operator* ( const Unity< Field > &u, const Field &f )
00052   {
00053     return f;
00054   }
00055 
00056   template< class Field >
00057   Field operator/ ( const Unity< Field > &u, const Field &f )
00058   {
00059     return (Field)u / f;
00060   }
00061 
00062 
00063 
00064   // Zero
00065   // ----
00066 
00078   template< class Field >
00079   struct Zero
00080   {
00081     operator Field () const
00082     {
00083       return Field( 0 );
00084     }
00085     static const Field epsilon()
00086     {
00087       return Field(1e-12);
00088     }
00089   };
00090 #if HAVE_ALGLIB
00091   template< unsigned int precision >
00092   struct Zero< amp::ampf< precision > >
00093   {
00094     typedef amp::ampf< precision > Field;
00095     operator Field () const
00096     {
00097       return Field( 0 );
00098     }
00099     static const Field epsilon()
00100     {
00101       return Field(1e-20);
00102     }
00103   };
00104 #endif
00105 #if HAVE_GMP
00106   template< unsigned int precision >
00107   struct Zero< GMPField< precision > >
00108   {
00109     typedef GMPField< precision > Field;
00110     operator Field () const
00111     {
00112       return Field( 0 );
00113     }
00114     static const Field epsilon()
00115     {
00116       return Field(1e-20);
00117     }
00118   };
00119 #endif
00120 
00121   template< class Field >
00122   inline bool operator == ( const Zero< Field > &, const Field &f )
00123   {
00124     return ( f < Zero<Field>::epsilon() && f > -Zero<Field>::epsilon() );
00125   }
00126 
00127   template< class Field >
00128   inline bool operator == ( const Field &f, const Zero< Field > &z)
00129   {
00130     return ( z == f );
00131   }
00132 
00133   template< class Field >
00134   inline bool operator< ( const Zero< Field > &, const Field &f )
00135   {
00136     return f > Zero<Field>::epsilon();
00137   }
00138 
00139   template< class Field >
00140   inline bool operator< ( const Field &f, const Zero< Field > & )
00141   {
00142     return f < -Zero<Field>::epsilon();
00143   }
00144 
00145   template< class Field >
00146   inline bool operator> ( const Zero< Field > &z, const Field &f )
00147   {
00148     return f < z;
00149   }
00150 
00151   template< class Field >
00152   inline bool operator> ( const Field &f, const Zero< Field > &z )
00153   {
00154     return z < f;
00155   }
00156 
00157 
00158   // field_cast
00159   // ----------
00160 
00173   template< class F2, class F1 >
00174   inline void field_cast ( const F1 &f1, F2 &f2 )
00175   {
00176     f2 = f1;
00177   }
00178 
00179 #if HAVE_ALGLIB
00180   template< unsigned int precision >
00181   inline void field_cast ( const amp::ampf< precision > &f1, double &f2 )
00182   {
00183     f2 = f1.toDouble();
00184   }
00185 
00186   template< unsigned int precision >
00187   inline void field_cast ( const amp::ampf< precision > &f1, long double &f2 )
00188   {
00189     f2 = f1.toDouble();
00190   }
00191 #endif
00192 
00193 #if HAVE_GMP
00194   template< unsigned int precision >
00195   inline void field_cast ( const Dune::GMPField< precision > &f1, double &f2 )
00196   {
00197     f2 = f1.get_d();
00198   }
00199 
00200   template< unsigned int precision >
00201   inline void field_cast ( const Dune::GMPField< precision > &f1, long double &f2 )
00202   {
00203     f2 = f1.get_d();
00204   }
00205 #endif
00206 
00207   template< class F2, class F1, int dim >
00208   inline void field_cast ( const Dune::FieldVector< F1, dim > &f1, Dune::FieldVector< F2, dim > &f2 )
00209   {
00210     for( int d = 0; d < dim; ++d )
00211       field_cast( f1[ d ], f2[ d ] );
00212   }
00213   template< class F2, class F1 >
00214   inline void field_cast ( const Dune::FieldVector< F1, 1 > &f1, F2 &f2 )
00215   {
00216     field_cast( f1[ 0 ], f2 );
00217   }
00218   template< class F2, class F1 >
00219   inline void field_cast ( const F1 &f1, Dune::FieldVector< F2, 1 > &f2 )
00220   {
00221     field_cast( f1, f2[ 0 ] );
00222   }
00223 
00224   template< class F2, class F1, int rdim, int cdim >
00225   inline void field_cast ( const Dune::FieldMatrix< F1, rdim, cdim > &f1, Dune::FieldMatrix< F2, rdim, cdim > &f2 )
00226   {
00227     for( int r = 0; r < rdim; ++r )
00228       field_cast( f1[ r ], f2[ r ] );
00229   }
00230   template< class F2, class F1 >
00231   inline void field_cast ( const Dune::FieldMatrix<F1,1,1> &f1, Dune::FieldMatrix< F2, 1,1 > &f2 )
00232   {
00233     field_cast( f1[ 0 ][ 0 ], f2[ 0 ][ 0 ] );
00234   }
00235   template< class F2, class F1 >
00236   inline void field_cast ( const Dune::FieldMatrix< F1, 1,1 > &f1, F2 &f2 )
00237   {
00238     field_cast( f1[ 0 ][ 0 ], f2 );
00239   }
00240   template< class F2, class F1 >
00241   inline void field_cast ( const F1 &f1, Dune::FieldMatrix< F2, 1,1 > &f2 )
00242   {
00243     field_cast( f1, f2[ 0 ][ 0 ] );
00244   }
00245   template< class F2, class F1 >
00246   inline void field_cast ( const Dune::FieldVector<F1,1> &f1, Dune::FieldMatrix< F2, 1,1 > &f2 )
00247   {
00248     field_cast( f1[ 0 ], f2[ 0 ][ 0 ] );
00249   }
00250   template< class F2, class F1 >
00251   inline void field_cast ( const Dune::FieldMatrix<F1,1,1> &f1, Dune::FieldVector< F2, 1 > &f2 )
00252   {
00253     field_cast( f1[ 0 ][ 0 ], f2[ 0 ] );
00254   }
00255 
00256   template< class F2, class F1 >
00257   inline void field_cast ( const Dune::FieldVector< F1, 1 > &f1, Dune::FieldVector<F2, 1> &f2 )
00258   {
00259     field_cast( f1[ 0 ], f2[ 0 ] );
00260   }
00261 
00262   template< class F2,class V > 
00263   struct FieldCast
00264   {
00265     typedef F2 type;
00266   };
00267   template< class F2,class F1,int dim > 
00268   struct FieldCast< F2, Dune::FieldVector<F1,dim> >
00269   {
00270     typedef Dune::FieldVector<F2,dim> type;
00271   };
00272   template< class F2,class F1,int dim1, int dim2> 
00273   struct FieldCast< F2, Dune::FieldMatrix<F1,dim1,dim2> >
00274   {
00275     typedef Dune::FieldMatrix<F2,dim1,dim2> type;
00276   };
00277   template< class F2,class V > 
00278   inline typename FieldCast<F2,V>::type field_cast ( const V &f1 )
00279   {
00280     typename FieldCast<F2,V>::type f2;
00281     field_cast( f1, f2 );
00282     return f2;
00283   }
00284 
00285 
00286   // Precision
00287   // this is not a perfect solution to obtain the
00288   // precision of a field - definition is not clear
00289   // to be removed
00290   // ---------
00291 
00292   template <class Field>
00293   struct Precision;
00294 
00295   template<>
00296   struct Precision< double >
00297   {
00298     static const unsigned int value = 64;
00299   };
00300 
00301   template<>
00302   struct Precision< long double >
00303   {
00304     static const unsigned int value = 80;
00305   };
00306 
00307   template<>
00308   struct Precision< float >
00309   {
00310     static const unsigned int value = 32;
00311   };
00312 
00313 #if HAVE_ALGLIB
00314   template< unsigned int precision >
00315   struct Precision< amp::ampf< precision > >
00316   {
00317     static const unsigned int value = precision;
00318   };
00319 #endif
00320 #if HAVE_GMP
00321   template< unsigned int precision >
00322   struct Precision< GMPField< precision > >
00323   {
00324     static const unsigned int value = precision;
00325   };
00326 #endif
00327 
00328   // ComputeField
00329   // ------------
00330 
00331   template <class Field,unsigned int sum>
00332   struct ComputeField 
00333   {
00334     typedef Field Type;
00335   };
00336 
00337 #if HAVE_ALGLIB
00338   template< unsigned int precision, unsigned int sum >
00339   struct ComputeField< amp::ampf< precision >, sum >
00340   {
00341     typedef amp::ampf<precision+sum> Type;
00342   };
00343 #endif
00344 #if HAVE_GMP
00345   template< unsigned int precision, unsigned int sum >
00346   struct ComputeField< GMPField< precision >, sum >
00347   {
00348     typedef GMPField<precision+sum> Type;
00349   };
00350 #endif
00351 } // namespace Dune
00352 
00353 // to be moved to different location...
00354 namespace std
00355 {
00356 
00357 #if HAVE_ALGLIB
00358   template< unsigned int precision >
00359   inline ostream &
00360   operator<< ( ostream &out, 
00361                const amp::ampf< precision > &value )
00362   {
00363     return out << value.toDec();
00364   }
00365 
00366   template< unsigned int precision >
00367   inline amp::ampf< precision > sqrt ( const amp::ampf< precision > &a )
00368   {
00369     return amp::sqrt( a );
00370   }
00371 
00372   template< unsigned int precision >
00373   inline amp::ampf< precision > abs ( const amp::ampf< precision > &a )
00374   {
00375     return amp::abs( a );
00376   }
00377 #endif // #if HAVE_ALGLIB
00378 
00379 } // namespace std
00380 
00381 namespace Dune
00382 {
00383 
00384   namespace GenericGeometry
00385   {
00386 
00387     // FieldHelper
00388     // -----------
00389 
00390     template< class Field >
00391     struct FieldHelper;
00392 
00393 #if HAVE_ALGLIB
00394     template< unsigned int precision >
00395     struct FieldHelper< amp::ampf< precision > >
00396     {
00397       typedef amp::ampf< precision > Field;
00398 
00399       static Field abs ( const Field &x ) { return amp::abs( x ); }
00400     };
00401 #endif // #if HAVE_ALGLIB
00402 
00403   } // namespace GenericGeometry
00404 
00405 } // namespace Dune
00406 
00407 #endif // #ifndef DUNE_LOCALFUNCTIONS_UTILITY_FIELD_HH