Helpers.h 12.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
/** \file Helpers.h */

#ifndef MY_HELPERS_H
#define MY_HELPERS_H

#include "AMDiS.h"

#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/StdMpi.h"
#endif

12
13
14
#ifdef _WIN32
#include <io.h>
#else
15
#include <unistd.h>
16
#endif
17
18
19
20
#include <ios>
#include <iostream>
#include <fstream>
#include <string>
21

22
23
24
25
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
#include "boost/lexical_cast.hpp"
#include "boost/date_time/posix_time/posix_time.hpp"
26
27
#include <boost/math/special_functions/fpclassify.hpp>

28
29
30
31
#include <time.h>
#include <stdarg.h>

#include "VectorOperations.h"
32
#include "Views.h"
33
34
35
36
37
38

#define TEST_WARNING(test) if ((test));else WARNING

using namespace std;
using namespace AMDiS; 

39
static long random_seed_initial_value = 0;
40

41
namespace Helpers {
42

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
43
44
  // math routines
  // ==============
45
  
46
  inline double cint(double x) {
47
    double intpart; 
48
49
    if (modf(x, &intpart) >= 0.5)
      return x >=0 ? ceil(x) : floor(x);
50
    else
51
      return x < 0 ? ceil(x) : floor(x);
52
53
  }
  
54
55
  inline double round (double r, double places) {
    double off = pow(10.0, places);
56
57
58
    return cint(r*off)/off;
  }
  
59
60
61
62
  inline double signum(double val, double tol = DBL_TOL) { return val > tol ? 1.0 : (val < -tol ? -1.0 : 0.0); }
  inline int signum(int val) { return val > 0 ? 1 : (val < 0 ? -1 : 0); }
  inline int signum(bool val) { return val ? 1 : -1; }
  inline int toInt(bool val) { return val ? 1 : 0; }
63
  inline bool isNumber(double val) { return !boost::math::isnan(val) && !boost::math::isinf(val); }
64
65
66
  inline bool toBool(double val, double tol = DBL_TOL) { return val > tol || val < -tol ? true : false; }
  inline bool isPositiv(double val, double tol = DBL_TOL) { return val > -tol; }
  inline bool isStrictPositiv(double val, double tol = DBL_TOL) { return val > tol; }
67

68
69
70
71
72
73
74
75
76
77
  template<typename T>
  T atanh(T x)
  {
#ifndef _WIN32
    return atanh(x);
#else
    return 0.5*log((1.0+x)/(1.0-x));
#endif
  }
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
78
79
  // routines for conversion to string
  // =================================
80
81
  
  template<typename T>
82
83
  inline std::string toString(const T &value, 
			      ios_base::fmtflags flag = ios_base::scientific)
84
85
86
87
88
89
90
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
91
  }
92
  
93
  inline std::string toString(const int &value)
94
95
  {
    return boost::lexical_cast<std::string>(value);
96
  }
97
98
  
  template<typename T>
99
100
101
  inline std::string toString(const WorldVector<T> &value, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
102
103
104
105
106
107
108
109
110
111
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    
    if (brackets) ss<<"]";
    return ss.str();
112
  }
113
114
  
  template<typename T>
115
116
117
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
118
119
120
121
122
123
124
125
126
127
128
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
    ss.precision(6);
    for (size_t i = 0; i < vec.size(); ++i) {
      if (!(ss<<(i>0?", ":"")<<vec[i]))
	throw(std::runtime_error("toString() not possible"));
    }
    if (brackets) ss<<"]";
    return ss.str();
129
  }
130
131
  
  template<class T>
132
  inline T fromString(const std::string& s)
133
134
135
136
137
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
138
  }
139

140
  inline std::string fillString(int length, char c, int numValues, ...)
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
  {
    va_list values;
    int value;
    va_start(values, numValues);

    int len = length;
    for(int i=0; i<numValues; i++) {
      value = va_arg(values, int);
      double tmp = static_cast<double>(value);
      while (true) {
	len --;
	tmp /= 10.0;
	if (tmp < 1.0 - DBL_TOL)
	  break;
      }
    }
    va_end(values);

    return std::string(len, c);
160
  }
161

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
162
  /// produce printable string of memory
163
  inline std::string memoryToString(double mem, int startIdx = 0) {
164
165
166
167
168
169
170
171
    int idx = startIdx;
    double mem_ = mem;
    while (mem_/1024.0 > 1.0) {
      mem_ /= 1024.0;
      idx++;
    }
    std::string result = toString(mem_)+" "+(idx==0?"B":(idx==1?"KB":(idx==2?"MB":"GB")));
    return result;
172
  }
173
174
  
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
175
176
  // some mesh routines
  // ===========================
177
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
178
  /// rotate scale and shift the mesh coordinates
179
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
180

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
181
  /// scale a vector of meshes by different values in all directions
182
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
183
184
  
  /// scale and shift by different values in all directions
185
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
186
187
  
  /// scale by different values in all directions
188
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
189
190
  
  /// scale and shift by the same values in all directions
191
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
192
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
193
194
195
196
197
  /**
   * calculate the dimension of a mesh, by mesh traversal.
   * The mthod determines the maximal mesh coordinates and assumes symmetry of the
   * mesh around the origin, i.e. mesh = [-rsult, result]
   **/
198
  WorldVector<double> getMeshDimension(Mesh *mesh);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
199
200
  
  /// calculate the dimension of a mesh, by mesh traversal.
201
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
202
203

  /// read DOFVector from AMDiS .dat-files
204
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
205

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
206
207
  // some linear algebra methods
  // ===========================
208

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
209
  /// calculate the determant of a 3x3 matrix of type dense2D
210
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
211
212
213
214
215
216
  {
    if(num_rows(A)==3 && num_cols(A)==3) {
      double det = A(0,0)*A(1,1)*A(2,2)+A(0,1)*A(1,2)*A(2,0)+A(0,2)*A(1,0)*A(2,1);
      return det-(A(0,2)*A(1,1)*A(2,0)+A(0,1)*A(1,0)*A(2,2)+A(0,0)*A(1,2)*A(2,1));
    }
    return 1.0;
217
  }
218
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
219
  /// calculate the determant of a 3x3 matrix of type WorldMatrix
220
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
221
222
223
224
225
226
  {
    if(A.getNumRows()==3 && A.getNumCols()==3) {
      double det = A[0][0]*A[1][1]*A[2][2]+A[0][1]*A[1][2]*A[2][0]+A[0][2]*A[1][0]*A[2][1];
      return det-(A[0][2]*A[1][1]*A[2][0]+A[0][1]*A[1][0]*A[2][2]+A[0][0]*A[1][2]*A[2][1]);
    }
    return 1.0;
227
  }
228
229
230
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
231
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
232
233
234
235
236
237
238
239
240
241
242
  { FUNCNAME("Helpers::inverse_iteration()");

    EigenVector y( size(x) ), res( size(x) );
    double lambda = 0.0, lambda_old = 0.0, relErr;
    mtl::seed<double> seed;
    random(x, seed);
    x /= two_norm(x);

    solver.setMultipleRhs(true);
    for (int i = 0; i < m; ++i) {
      solver.solveSystem(A, y, x ); // solve Ay=x --> y
243
244
245
246
247
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
248
249
250
251
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

252
253
    return lambda;
  }
254
255
256

  /// calculate approximation to condition number of matrix A using the OEMSolver solver
  template <class LinearSolver, class Matrix>
257
  inline double condition (LinearSolver& solver, const Matrix& A, int m=10)
258
259
260
261
262
263
264
265
266
267
  { FUNCNAME("Helpers::condition()");

    mtl::dense_vector<double>   x(num_rows(A)), y(num_rows(A));
    mtl::seed<double>           seed;
    random(x, seed);

    double lambda = 0.0, lambda_old = 0.0, relErr, result1;
    x /= two_norm(x);
    for (int i = 0; i < m; ++i) {
      y = A * x;
268
      lambda_old = lambda;
269
270
      lambda = mtl::dot(y,x);
      relErr = abs((lambda-lambda_old)/lambda);
271
272
      if (relErr < 1.e-5)
	break;
273
274
275
276
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
277
  }
278
279
280
  
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
281
  inline void interpolOverLine(const Container &vec, 
282
283
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
			       std::vector<WorldVector<double> > &x, std::vector<double> &y)
  { FUNCNAME("Helpers::interpolOverLine()");

    WorldVector<double> p;

    if (nPoints <= 1)
      throw(std::runtime_error("Zu wenig Punkte fuer eine Interpolation!"));
    
    for (int i = 0; i < nPoints; ++i) {
      double lambda = static_cast<double>(i)/static_cast<double>(nPoints - 1.0);
      p[0] = lambda*x2 + (1.0-lambda)*x1;
      p[1] = lambda*y2 + (1.0-lambda)*y1;
      
      double value = evalAtPoint(vec, p);
      x.push_back(p);
      y.push_back(value);
    }
  }
302

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
303
  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
304
  static void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
305
306
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
307
308
309
310
311
312
313
314
315
316
317
  static void calcMaxOnYAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);

  /// calc normal vectors of surface from element normals by averaging
  static void getNormalsWeighted(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals);
  
  /// calc normal vectors of surface from element normals by local approximation by quartic
  static void getNormals(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals, DOFVector<WorldMatrix<double> > *gradNormals = NULL);
  
  /// calc curvature from given normal vectors
  static void getCurvature(DOFVector<WorldVector<double> >* normals, DOFVector<double>* curvature);
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
318
319
  // misc routines
  // =============
320
  
321
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
322
  
323
  inline long getMicroTimestamp() {
324
325
326
327
328
    using namespace boost::posix_time;
    ptime t0(min_date_time);
    ptime now = microsec_clock::local_time();
    time_duration td = now-t0;
    return td.total_microseconds();
329
  }
330
  
331
  inline long getRandomSeed(int param = 0) {
332
333
334
335
336
337
338
339
    using namespace boost::posix_time;
    ptime t0(min_date_time);
    ptime now = microsec_clock::local_time();
    time_duration td = now-t0;
    long value0 = td.total_microseconds();
    long value1 = clock();
    long value2 = random_seed_initial_value++;

340
    return value0 + value1 + value2 + param;
341
  }
342

343
344
  
} // end namespace Helpers
345

346
347
348
349
350
351
352
353
namespace AMDiS {

  /** \ingroup DOFAdministration
    * \brief
    * Implements a DOFIterator for a DOFIndexed<T> object
    */
  template<typename T>
  class DOFConstIterator : public DOFIteratorBase
354
  {
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
  public:
    /// Constructs a DOFIterator for cont of type t
    DOFConstIterator(const DOFIndexed<T> *obj, DOFIteratorType t)
      : DOFIteratorBase(dynamic_cast<DOFAdmin*>(obj->getFeSpace()->getAdmin()), t),
	iteratedObject(obj)
    {}

    /// Constructs a DOFIterator for cont of type t
    DOFConstIterator(DOFAdmin *admin,
		const DOFIndexed<T> *obj,
		DOFIteratorType t)
      : DOFIteratorBase(admin, t),
	iteratedObject(obj)
    {}

    /// Dereference operator
    inline const T& operator*()
    {
      return *it;
    }
375

376
377
378
379
380
    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }
381

382
383
384
385
    inline bool operator!=(const DOFIterator<T>& rhs)
    {
      if (this->iteratedObject != rhs.iteratedObject)
	return true;
386

387
388
      if (this->it != rhs.it)
	return true;
389

390
391
      return false;
    }
392

393
394
395
396
    inline bool operator==(const DOFIterator<T>& rhs)
    {
      return !(this->operator==(rhs));
    }
397

398
399
400
401
402
403
  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }
404

405
406
407
408
409
    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }
410

411
412
413
414
415
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }
416

417
418
419
420
421
422
423
424
425
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

  protected:
    /// Object that is iterated
    const DOFIndexed<T> *iteratedObject;
426

427
428
429
    /// Iterator for \ref iteratedObject
    typename std::vector<T>::const_iterator it;
  };
430

431
} // end namespace AMDiS
432
433

#endif