Utility.h 15.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
//======================================================================================================================
//
//  This file is part of waLBerla. waLBerla is free software: you can 
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of 
//  the License, or (at your option) any later version.
//  
//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
//  for more details.
//  
//  You should have received a copy of the GNU General Public License along
//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file Utility.h
//! \ingroup core
//! \author Klaus Iglberger
//! \brief Header file for mathematical functions and constants
//
//======================================================================================================================

#pragma once

#include "MathTrait.h"
#include "core/DataTypes.h"

#include <boost/math/constants/constants.hpp>
#include <cmath>
#include <cstddef>
#include <limits>


namespace walberla {
namespace math {

//======================================================================================================================
//
//  MATHEMATICAL CONSTANTS
//
//======================================================================================================================

//! Definition of the mathematical value \f$ \pi \f$.
const real_t PI = boost::math::constants::pi<real_t>();



//======================================================================================================================
//
//  MATHEMATICAL UTILITY FUNCTIONS
//
//======================================================================================================================

//**********************************************************************************************************************
/*!\name Mathematical utility functions */
//@{
template< typename T >
inline const T sign( T a );

template< typename T >
inline const typename boost::enable_if_c<  boost::is_unsigned<T>::value, T >::type abs( T a );
template< typename T >
inline const typename boost::enable_if_c< !boost::is_unsigned<T>::value, T >::type abs( T a );


template< typename T1, typename T2 >
inline const typename MathTrait<T1,T2>::High min( const T1& a, const T2& b );

template< typename T1, typename T2, typename T3 >
inline const typename MathTrait< typename MathTrait<T1,T2>::High, T3 >::High min( const T1& a, const T2& b, const T3& c );

template< typename T1, typename T2 >
inline const typename MathTrait<T1,T2>::High max( const T1& a, const T2& b );

template< typename T1, typename T2, typename T3 >
inline const typename MathTrait< typename MathTrait<T1,T2>::High, T3 >::High max( const T1& a, const T2& b, const T3& c );

template< typename T >
inline const T sqr( const T& a );

template< typename T1, typename T2 >
inline bool equal( T1 a, T2 b );

inline real_t round( real_t a );
//@}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Sign function.
// \ingroup math
//
// \param a The signed value.
// \return 1 if the value is greater than or equal to zero, -1 if the value is smaller than zero.
//
// The sign function only works for signed built-in data types. The attempt to use unsigned data
// types or user-defined class types will result in a compile time error.
 */
template< typename T >
inline const T sign( T a )
{
102
   WALBERLA_STATIC_ASSERT( std::numeric_limits<T>::is_signed );
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
   return ( a < T(0) )?( T(-1) ):( T(1) );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Absolute value function.
// \ingroup math
//
// \param a The value.
// \return The value if it is greater than or equal to zero, -1 times the value if the value is smaller than zero.
 */
template< typename T >
inline const typename boost::enable_if_c<  boost::is_unsigned<T>::value, T >::type abs( T a )
{
   return a;
}

template< typename T >
inline const typename boost::enable_if_c< !boost::is_unsigned<T>::value, T >::type abs( T a )
{
   return std::abs( a );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Minimum function for two arguments.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return The minimum of the two values.
//
// This function returns the minimum of the two given data values. The return type of the
// function is determined by the data types of the given arguments (for further detail see
// the MathTrait class description).
 */
template< typename T1, typename T2 >
inline const typename MathTrait<T1,T2>::High min( const T1& a, const T2& b )
{
   // The return type of the function is only a copy of the one of the arguments for two reasons:
   //  - in case the data types T1 and T2 are equal, a reference return type could result in a
   //    bug if combined with literals
   //  - in case the two data types are unequal, the result of the comparison could be converted
   //    to the more significant data type, which results in a local temporary value
   // The copy operation might cause a performance decrease for class types, which is probably
   // avoided if the function is inlined.
   return ( a < b )?( a ):( b );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Minimum function for three arguments.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \param c Third value
// \return The minimum of the three values.
//
// This function returns the minimum of the three given data values. The return type of the
// function is determined by the data types of the given arguments (for further detail see
// the MathTrait class description).
 */
template< typename T1, typename T2, typename T3 >
inline const typename MathTrait< typename MathTrait<T1,T2>::High, T3 >::High min( const T1& a, const T2& b, const T3& c )
{
   // The return type of the function is only a copy of the one of the arguments for two reasons:
   //  - in case the data types T1, T2, and T3 are equal, a reference return type could result in
   //    a bug if combined with literals
   //  - in case the three data types are unequal, the result of the comparison could be converted
   //    to the more significant data type, which results in a local temporary value
   // The copy operation might cause a performance decrease for class types, which is probably
   // avoided if the function is inlined.
   return ( a < b )?( ( a < c )?( a ):( c ) ):( ( b < c )?( b ):( c ) );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Maximum function for two arguments.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return The maximum of the two values.
//
// This function returns the maximum of the two given data values. The return type of the
// function is determined by the data types of the given arguments (for further detail see
// the MathTrait class description).
 */
template< typename T1, typename T2 >
inline const typename MathTrait<T1,T2>::High max( const T1& a, const T2& b )
{
   // The return type of the function is only a copy of the one of the arguments for two reasons:
   //  - in case the data types T1 and T2 are equal, a reference return type could result in a
   //    bug if combined with literals
   //  - in case the two data types are unequal, the result of the comparison could be converted
   //    to the more significant data type, which results in a local temporary value
   // The copy operation might cause a performance decrease for class types, which is probably
   // avoided if the function is inlined.
   return ( a > b )?( a ):( b );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Maximum function for three arguments.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \param c Third value.
// \return The maximum of the three values.
//
// This function returns the maximum of the three given data values. The return type of the
// function is determined by the data types of the given arguments (for further detail see
// the MathTrait class description).
 */
template< typename T1, typename T2, typename T3 >
inline const typename MathTrait< typename MathTrait<T1,T2>::High, T3 >::High max( const T1& a, const T2& b, const T3& c )
{
   // The return type of the function is only a copy of the one of the arguments for two reasons:
   //  - in case the data types T1, T2, and T3 are equal, a reference return type could result in
   //    a bug if combined with literals
   //  - in case the three data types are unequal, the result of the comparison could be converted
   //    to the more significant data type, which results in a local temporary value
   // The copy operation might cause a performance decrease for class types, which is probably
   // avoided if the function is inlined.
   return ( a > b )?( ( a > c )?( a ):( c ) ):( ( b > c )?( b ):( c ) );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Square function.
// \ingroup math
//
// \param a Value to be squared.
// \return The square of the input value.
 */
template< typename T >
inline const T sqr( const T& a )
{
   return a*a;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*! \cond internal */
/*!\brief Equality check for integral data types.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return \a true if the two values are equal, \a false if not.
//
// Equal function for the comparison of two integral values.
 */
template< typename T >
inline bool equal_backend( T a, T b )
{
   return a == b;
}
/*! \endcond */
//**********************************************************************************************************************


//**********************************************************************************************************************
/*! \cond internal */
/*!\brief Equality check for single precision floating point values.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return \a true if the two values are equal, \a false if not.
//
// Equal function for the comparison of two single precision floating point numbers. Due to the
// limited machine accuracy, a direct comparison of two floating point numbers should be avoided.
// This functions offers the possibility to compare two floating-point values with a certain
// accuracy margin.
 */
template<>
inline bool equal_backend<float>( float a, float b )
{
   return std::fabs( a - b ) < 1E-8;
}
/*! \endcond */
//**********************************************************************************************************************


//**********************************************************************************************************************
/*! \cond internal */
/*!\brief Equality check for double precision floating point values.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return \a true if the two values are equal, \a false if not.
//
// Equal function for the comparison of two double precision floating point numbers. Due to the
// limited machine accuracy, a direct comparison of two floating point numbers should be avoided.
// This functions offers the possibility to compare two floating-point values with a certain
// accuracy margin.
 */
template<>
inline bool equal_backend<double>( double a, double b )
{
   return std::fabs( a - b ) < 1E-8;
}
/*! \endcond */
//**********************************************************************************************************************


//**********************************************************************************************************************
/*! \cond internal */
/*!\brief Equality check for long double precision floating point values.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return \a true if the two values are equal, \a false if not.
//
// Equal function for the comparison of two long double precision floating point numbers. Due
// to the limited machine accuracy, a direct comparison of two floating point numbers should be
// avoided. This functions offers the possibility to compare two floating-point values with a
// certain accuracy margin.
 */
template<>
inline bool equal_backend<long double>( long double a, long double b )
{
   return std::fabs( a - b ) < 1E-10;
}
/*! \endcond */
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Generic equality check.
// \ingroup math
//
// \param a First value.
// \param b Second value.
// \return \a true if the two values are equal, \a false if not.
//
// Generic equal function for the comparison of two numeric values. Depending on the types of
// the two arguments, a special comparison for floating point values is selected that takes
// the limited machine accuracy into account.
 */
template< typename T1, typename T2 >
inline bool equal( T1 a, T2 b )
{
   typedef typename MathTrait<T1,T2>::High High;
   return equal_backend<High>( a, b );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\brief Rounding function.
// \ingroup math
//
// \param  a Floating point value.
// \return The floating point value rounded towards the next whole number.

 */
inline real_t round( real_t a )
{
   return std::floor( a + real_t(0.5) );
}
//**********************************************************************************************************************



} // namespace math
} // namespace walberla