dune-localfunctions  2.2.1
field.hh
Go to the documentation of this file.
1 #ifndef DUNE_LOCALFUNCTIONS_UTILITY_FIELD_HH
2 #define DUNE_LOCALFUNCTIONS_UTILITY_FIELD_HH
3 
4 #if HAVE_ALGLIB
5 #include <alglib/amp.h>
6 #warning ALGLIB support is deprecated and will be dropped after DUNE 2.2 (cf. FS#931)
7 #endif
8 
9 #include <dune/common/gmpfield.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
12 
13 namespace Dune
14 {
15 
16  // Unity
17  // -----
18 
29  template< class Field >
30  struct Unity
31  {
32  operator Field () const
33  {
34  return Field( 1 );
35  }
36  };
37 
38  template< class Field >
39  Field operator+ ( const Unity< Field > &u, const Field &f )
40  {
41  return (Field)u + f;
42  }
43 
44  template< class Field >
45  Field operator- ( const Unity< Field > &u, const Field &f )
46  {
47  return (Field)u - f;
48  }
49 
50  template< class Field >
51  Field operator* ( const Unity< Field > &u, const Field &f )
52  {
53  return f;
54  }
55 
56  template< class Field >
57  Field operator/ ( const Unity< Field > &u, const Field &f )
58  {
59  return (Field)u / f;
60  }
61 
62 
63 
64  // Zero
65  // ----
66 
78  template< class Field >
79  struct Zero
80  {
81  operator Field () const
82  {
83  return Field( 0 );
84  }
85  static const Field epsilon()
86  {
87  return Field(1e-12);
88  }
89  };
90 #if HAVE_ALGLIB
91  template< unsigned int precision >
92  struct Zero< amp::ampf< precision > >
93  {
94  typedef amp::ampf< precision > Field;
95  operator Field () const
96  {
97  return Field( 0 );
98  }
99  static const Field epsilon()
100  {
101  return Field(1e-20);
102  }
103  };
104 #endif
105 #if HAVE_GMP
106  template< unsigned int precision >
107  struct Zero< GMPField< precision > >
108  {
109  typedef GMPField< precision > Field;
110  operator Field () const
111  {
112  return Field( 0 );
113  }
114  static const Field epsilon()
115  {
116  return Field(1e-20);
117  }
118  };
119 #endif
120 
121  template< class Field >
122  inline bool operator == ( const Zero< Field > &, const Field &f )
123  {
124  return ( f < Zero<Field>::epsilon() && f > -Zero<Field>::epsilon() );
125  }
126 
127  template< class Field >
128  inline bool operator == ( const Field &f, const Zero< Field > &z)
129  {
130  return ( z == f );
131  }
132 
133  template< class Field >
134  inline bool operator< ( const Zero< Field > &, const Field &f )
135  {
136  return f > Zero<Field>::epsilon();
137  }
138 
139  template< class Field >
140  inline bool operator< ( const Field &f, const Zero< Field > & )
141  {
142  return f < -Zero<Field>::epsilon();
143  }
144 
145  template< class Field >
146  inline bool operator> ( const Zero< Field > &z, const Field &f )
147  {
148  return f < z;
149  }
150 
151  template< class Field >
152  inline bool operator> ( const Field &f, const Zero< Field > &z )
153  {
154  return z < f;
155  }
156 
157 
158  // field_cast
159  // ----------
160 
173  template< class F2, class F1 >
174  inline void field_cast ( const F1 &f1, F2 &f2 )
175  {
176  f2 = f1;
177  }
178 
179 #if HAVE_ALGLIB
180  template< unsigned int precision >
181  inline void field_cast ( const amp::ampf< precision > &f1, double &f2 )
182  {
183  f2 = f1.toDouble();
184  }
185 
186  template< unsigned int precision >
187  inline void field_cast ( const amp::ampf< precision > &f1, long double &f2 )
188  {
189  f2 = f1.toDouble();
190  }
191 #endif
192 
193 #if HAVE_GMP
194  template< unsigned int precision >
195  inline void field_cast ( const Dune::GMPField< precision > &f1, double &f2 )
196  {
197  f2 = f1.get_d();
198  }
199 
200  template< unsigned int precision >
201  inline void field_cast ( const Dune::GMPField< precision > &f1, long double &f2 )
202  {
203  f2 = f1.get_d();
204  }
205 #endif
206 
207  template< class F2, class F1, int dim >
208  inline void field_cast ( const Dune::FieldVector< F1, dim > &f1, Dune::FieldVector< F2, dim > &f2 )
209  {
210  for( int d = 0; d < dim; ++d )
211  field_cast( f1[ d ], f2[ d ] );
212  }
213  template< class F2, class F1 >
214  inline void field_cast ( const Dune::FieldVector< F1, 1 > &f1, F2 &f2 )
215  {
216  field_cast( f1[ 0 ], f2 );
217  }
218  template< class F2, class F1 >
219  inline void field_cast ( const F1 &f1, Dune::FieldVector< F2, 1 > &f2 )
220  {
221  field_cast( f1, f2[ 0 ] );
222  }
223 
224  template< class F2, class F1, int rdim, int cdim >
225  inline void field_cast ( const Dune::FieldMatrix< F1, rdim, cdim > &f1, Dune::FieldMatrix< F2, rdim, cdim > &f2 )
226  {
227  for( int r = 0; r < rdim; ++r )
228  field_cast( f1[ r ], f2[ r ] );
229  }
230  template< class F2, class F1 >
231  inline void field_cast ( const Dune::FieldMatrix<F1,1,1> &f1, Dune::FieldMatrix< F2, 1,1 > &f2 )
232  {
233  field_cast( f1[ 0 ][ 0 ], f2[ 0 ][ 0 ] );
234  }
235  template< class F2, class F1 >
236  inline void field_cast ( const Dune::FieldMatrix< F1, 1,1 > &f1, F2 &f2 )
237  {
238  field_cast( f1[ 0 ][ 0 ], f2 );
239  }
240  template< class F2, class F1 >
241  inline void field_cast ( const F1 &f1, Dune::FieldMatrix< F2, 1,1 > &f2 )
242  {
243  field_cast( f1, f2[ 0 ][ 0 ] );
244  }
245  template< class F2, class F1 >
246  inline void field_cast ( const Dune::FieldVector<F1,1> &f1, Dune::FieldMatrix< F2, 1,1 > &f2 )
247  {
248  field_cast( f1[ 0 ], f2[ 0 ][ 0 ] );
249  }
250  template< class F2, class F1 >
251  inline void field_cast ( const Dune::FieldMatrix<F1,1,1> &f1, Dune::FieldVector< F2, 1 > &f2 )
252  {
253  field_cast( f1[ 0 ][ 0 ], f2[ 0 ] );
254  }
255 
256  template< class F2, class F1 >
257  inline void field_cast ( const Dune::FieldVector< F1, 1 > &f1, Dune::FieldVector<F2, 1> &f2 )
258  {
259  field_cast( f1[ 0 ], f2[ 0 ] );
260  }
261 
262  template< class F2,class V >
263  struct FieldCast
264  {
265  typedef F2 type;
266  };
267  template< class F2,class F1,int dim >
268  struct FieldCast< F2, Dune::FieldVector<F1,dim> >
269  {
270  typedef Dune::FieldVector<F2,dim> type;
271  };
272  template< class F2,class F1,int dim1, int dim2>
273  struct FieldCast< F2, Dune::FieldMatrix<F1,dim1,dim2> >
274  {
275  typedef Dune::FieldMatrix<F2,dim1,dim2> type;
276  };
277  template< class F2,class V >
278  inline typename FieldCast<F2,V>::type field_cast ( const V &f1 )
279  {
280  typename FieldCast<F2,V>::type f2;
281  field_cast( f1, f2 );
282  return f2;
283  }
284 
285 
286  // Precision
287  // this is not a perfect solution to obtain the
288  // precision of a field - definition is not clear
289  // to be removed
290  // ---------
291 
292  template <class Field>
293  struct Precision;
294 
295  template<>
296  struct Precision< double >
297  {
298  static const unsigned int value = 64;
299  };
300 
301  template<>
302  struct Precision< long double >
303  {
304  static const unsigned int value = 80;
305  };
306 
307  template<>
308  struct Precision< float >
309  {
310  static const unsigned int value = 32;
311  };
312 
313 #if HAVE_ALGLIB
314  template< unsigned int precision >
315  struct Precision< amp::ampf< precision > >
316  {
317  static const unsigned int value = precision;
318  };
319 #endif
320 #if HAVE_GMP
321  template< unsigned int precision >
322  struct Precision< GMPField< precision > >
323  {
324  static const unsigned int value = precision;
325  };
326 #endif
327 
328  // ComputeField
329  // ------------
330 
331  template <class Field,unsigned int sum>
332  struct ComputeField
333  {
334  typedef Field Type;
335  };
336 
337 #if HAVE_ALGLIB
338  template< unsigned int precision, unsigned int sum >
339  struct ComputeField< amp::ampf< precision >, sum >
340  {
341  typedef amp::ampf<precision+sum> Type;
342  };
343 #endif
344 #if HAVE_GMP
345  template< unsigned int precision, unsigned int sum >
346  struct ComputeField< GMPField< precision >, sum >
347  {
348  typedef GMPField<precision+sum> Type;
349  };
350 #endif
351 } // namespace Dune
352 
353 // to be moved to different location...
354 namespace std
355 {
356 
357 #if HAVE_ALGLIB
358  template< unsigned int precision >
359  inline ostream &
360  operator<< ( ostream &out,
361  const amp::ampf< precision > &value )
362  {
363  return out << value.toDec();
364  }
365 
366  template< unsigned int precision >
367  inline amp::ampf< precision > sqrt ( const amp::ampf< precision > &a )
368  {
369  return amp::sqrt( a );
370  }
371 
372  template< unsigned int precision >
373  inline amp::ampf< precision > abs ( const amp::ampf< precision > &a )
374  {
375  return amp::abs( a );
376  }
377 #endif // #if HAVE_ALGLIB
378 
379 } // namespace std
380 
381 namespace Dune
382 {
383 
384  namespace GenericGeometry
385  {
386 
387  // FieldHelper
388  // -----------
389 
390  template< class Field >
391  struct FieldHelper;
392 
393 #if HAVE_ALGLIB
394  template< unsigned int precision >
395  struct FieldHelper< amp::ampf< precision > >
396  {
397  typedef amp::ampf< precision > Field;
398 
399  static Field abs ( const Field &x ) { return amp::abs( x ); }
400  };
401 #endif // #if HAVE_ALGLIB
402 
403  } // namespace GenericGeometry
404 
405 } // namespace Dune
406 
407 #endif // #ifndef DUNE_LOCALFUNCTIONS_UTILITY_FIELD_HH