dune-localfunctions
2.2.0
|
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