Helpers.h 12.4 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
#include <boost/math/special_functions/fpclassify.hpp>
27
#include <boost/math/special_functions/atanh.hpp>
28

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

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

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

using namespace std;
using namespace AMDiS; 

40
static long random_seed_initial_value = 0;
41

42
namespace Helpers {
43

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
44
45
  // math routines
  // ==============
46
  
47
  inline double cint(double x) {
48
    double intpart; 
49
50
    if (modf(x, &intpart) >= 0.5)
      return x >=0 ? ceil(x) : floor(x);
51
    else
52
      return x < 0 ? ceil(x) : floor(x);
53
54
  }
  
55
56
  inline double round (double r, double places) {
    double off = pow(10.0, places);
57
58
59
    return cint(r*off)/off;
  }
  
60
61
62
63
  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; }
64
  inline bool isNumber(double val) { return !boost::math::isnan(val) && !boost::math::isinf(val); }
65
66
67
  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; }
68

69
70
71
  template<typename T>
  T atanh(T x)
  {
72
73
74
75
76
77
78
79
80
81
    T result;
    try {
      result = boost::math::atanh(x);
    } catch(...) { // domain_error, or overflow_error
      if (x > 0.0)
	result = 10.0;
      else
	result = -10.0;
    }
    return result;
82
83
  }
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
84
85
  // routines for conversion to string
  // =================================
86
87
  
  template<typename T>
88
89
  inline std::string toString(const T &value, 
			      ios_base::fmtflags flag = ios_base::scientific)
90
91
92
93
94
95
96
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
97
  }
98
  
99
  inline std::string toString(const int &value)
100
101
  {
    return boost::lexical_cast<std::string>(value);
102
  }
103
104
  
  template<typename T>
105
106
107
  inline std::string toString(const WorldVector<T> &value, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
108
109
110
111
112
113
114
115
116
117
  {
    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();
118
  }
119
120
  
  template<typename T>
121
122
123
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
124
125
126
127
128
129
130
131
132
133
134
  {
    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();
135
  }
136
137
  
  template<class T>
138
  inline T fromString(const std::string& s)
139
140
141
142
143
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
144
  }
145

146
  inline std::string fillString(int length, char c, int numValues, ...)
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
  {
    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);
166
  }
167

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
168
  /// produce printable string of memory
169
  inline std::string memoryToString(double mem, int startIdx = 0) {
170
171
172
173
174
175
176
177
    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;
178
  }
179
180
  
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
181
182
  // some mesh routines
  // ===========================
183
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
184
  /// rotate scale and shift the mesh coordinates
185
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
186

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
187
  /// scale a vector of meshes by different values in all directions
188
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
189
190
  
  /// scale and shift by different values in all directions
191
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
192
193
  
  /// scale by different values in all directions
194
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
195
196
  
  /// scale and shift by the same values in all directions
197
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
198
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
199
200
201
202
203
  /**
   * 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]
   **/
204
  WorldVector<double> getMeshDimension(Mesh *mesh);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
205
206
  
  /// calculate the dimension of a mesh, by mesh traversal.
207
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
Praetorius, Simon's avatar
Praetorius, Simon committed
208
  
209
  /// read DOFVector from AMDiS .dat-files
210
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
211

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
212
213
  // some linear algebra methods
  // ===========================
214

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
215
  /// calculate the determant of a 3x3 matrix of type dense2D
216
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
217
218
219
220
221
222
  {
    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;
223
  }
224
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
225
  /// calculate the determant of a 3x3 matrix of type WorldMatrix
226
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
227
228
229
230
231
232
  {
    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;
233
  }
234
235
236
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
237
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
238
239
240
241
242
243
244
245
246
247
248
  { 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
249
250
251
252
253
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
254
255
256
257
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

258
259
    return lambda;
  }
260
261
262

  /// calculate approximation to condition number of matrix A using the OEMSolver solver
  template <class LinearSolver, class Matrix>
263
  inline double condition (LinearSolver& solver, const Matrix& A, int m=10)
264
265
266
267
268
269
270
271
272
273
  { 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;
274
      lambda_old = lambda;
275
276
      lambda = mtl::dot(y,x);
      relErr = abs((lambda-lambda_old)/lambda);
277
278
      if (relErr < 1.e-5)
	break;
279
280
281
282
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
283
  }
284
285
286
  
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
287
  inline void interpolOverLine(const Container &vec, 
288
289
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
			       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);
    }
  }
308

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
309
  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
310
  static void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
311
312
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
313
314
315
316
317
318
319
320
321
322
323
  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
324
325
  // misc routines
  // =============
326
  
327
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
328
  
329
  inline long getMicroTimestamp() {
330
331
332
333
334
    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();
335
  }
336
  
337
  inline long getRandomSeed(int param = 0) {
338
339
340
341
342
343
344
345
    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++;

346
    return value0 + value1 + value2 + param;
347
  }
348

349
350
  
} // end namespace Helpers
351

352
353
354
355
356
357
358
359
namespace AMDiS {

  /** \ingroup DOFAdministration
    * \brief
    * Implements a DOFIterator for a DOFIndexed<T> object
    */
  template<typename T>
  class DOFConstIterator : public DOFIteratorBase
360
  {
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
  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;
    }
381

382
383
384
385
386
    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }
387

388
389
390
391
    inline bool operator!=(const DOFIterator<T>& rhs)
    {
      if (this->iteratedObject != rhs.iteratedObject)
	return true;
392

393
394
      if (this->it != rhs.it)
	return true;
395

396
397
      return false;
    }
398

399
400
401
402
    inline bool operator==(const DOFIterator<T>& rhs)
    {
      return !(this->operator==(rhs));
    }
403

404
405
406
407
408
409
  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }
410

411
412
413
414
415
    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }
416

417
418
419
420
421
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }
422

423
424
425
426
427
428
429
430
431
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

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

433
434
435
    /// Iterator for \ref iteratedObject
    typename std::vector<T>::const_iterator it;
  };
436

437
} // end namespace AMDiS
438
439

#endif