Helpers.h 12.3 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
/** \file Helpers.h */

#ifndef MY_HELPERS_H
#define MY_HELPERS_H

#include "AMDiS.h"

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

#include <unistd.h>
#include <ios>
#include <iostream>
#include <fstream>
#include <string>
#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"
#include <time.h>
#include <stdarg.h>

#include "VectorOperations.h"
25
#include "Views.h"
26
27
28
29
30
31

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

using namespace std;
using namespace AMDiS; 

32
static long random_seed_initial_value = 0;
33

34
namespace Helpers {
35

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

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
61
62
  // routines for conversion to string
  // =================================
63
64
  
  template<typename T>
65
66
  inline std::string toString(const T &value, 
			      ios_base::fmtflags flag = ios_base::scientific)
67
68
69
70
71
72
73
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
74
  }
75
  
76
  inline std::string toString(const int &value)
77
78
  {
    return boost::lexical_cast<std::string>(value);
79
  }
80
81
  
  template<typename T>
82
83
84
  inline std::string toString(const WorldVector<T> &value, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
85
86
87
88
89
90
91
92
93
94
  {
    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();
95
  }
96
97
  
  template<typename T>
98
99
100
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
101
102
103
104
105
106
107
108
109
110
111
  {
    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();
112
  }
113
114
  
  template<class T>
115
  inline T fromString(const std::string& s)
116
117
118
119
120
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
121
  }
122

123
  inline std::string fillString(int length, char c, int numValues, ...)
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
  {
    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);
143
  }
144

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
145
  /// produce printable string of memory
146
  inline std::string memoryToString(double mem, int startIdx = 0) {
147
148
149
150
151
152
153
154
    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;
155
  }
156
157
  
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
158
159
  // some mesh routines
  // ===========================
160
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
161
  /// rotate scale and shift the mesh coordinates
162
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
163

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
164
  /// scale a vector of meshes by different values in all directions
165
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
166
167
  
  /// scale and shift by different values in all directions
168
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
169
170
  
  /// scale by different values in all directions
171
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
172
173
  
  /// scale and shift by the same values in all directions
174
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
175
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
176
177
178
179
180
  /**
   * 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]
   **/
181
  WorldVector<double> getMeshDimension(Mesh *mesh);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
182
183
  
  /// calculate the dimension of a mesh, by mesh traversal.
184
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
185
186

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

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
189
190
  // some linear algebra methods
  // ===========================
191

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
192
  /// calculate the determant of a 3x3 matrix of type dense2D
193
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
194
195
196
197
198
199
  {
    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;
200
  }
201
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
202
  /// calculate the determant of a 3x3 matrix of type WorldMatrix
203
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
204
205
206
207
208
209
  {
    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;
210
  }
211
212
213
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
214
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
215
216
217
218
219
220
221
222
223
224
225
  { 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
226
227
228
229
230
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
231
232
233
234
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

235
236
    return lambda;
  }
237
238
239

  /// calculate approximation to condition number of matrix A using the OEMSolver solver
  template <class LinearSolver, class Matrix>
240
  inline double condition (LinearSolver& solver, const Matrix& A, int m=10)
241
242
243
244
245
246
247
248
249
250
  { 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;
251
      lambda_old = lambda;
252
253
      lambda = mtl::dot(y,x);
      relErr = abs((lambda-lambda_old)/lambda);
254
255
      if (relErr < 1.e-5)
	break;
256
257
258
259
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
260
  }
261
262
263
  
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
264
  inline void interpolOverLine(const Container &vec, 
265
266
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
			       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);
    }
  }
285

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
286
  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
287
  static void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
288
289
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
290
291
292
293
294
295
296
297
298
299
300
  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
301
302
  // misc routines
  // =============
303
  
304
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
305
306
307
308
309
310

  // process_mem_usage(double &, double &) - takes two doubles by reference,
  // attempts to read the system-dependent data for a process' virtual memory
  // size and resident set size, and return the results in KB.
  //
  // On failure, returns 0.0, 0.0
311
  void process_mem_usage(double& vm_usage, double& resident_set);
312
  
313
  inline long getMicroTimestamp() {
314
315
316
317
318
319
320
    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();
  };
  
321
  inline long getRandomSeed(int param = 0) {
322
323
324
325
326
327
328
329
    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++;

330
    return value0 + value1 + value2 + param;
331
332
  };

333
334
  
} // end namespace Helpers
335

336
337
338
339
340
341
342
343
namespace AMDiS {

  /** \ingroup DOFAdministration
    * \brief
    * Implements a DOFIterator for a DOFIndexed<T> object
    */
  template<typename T>
  class DOFConstIterator : public DOFIteratorBase
344
  {
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
  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;
    }
365

366
367
368
369
370
    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }
371

372
373
374
375
    inline bool operator!=(const DOFIterator<T>& rhs)
    {
      if (this->iteratedObject != rhs.iteratedObject)
	return true;
376

377
378
      if (this->it != rhs.it)
	return true;
379

380
381
      return false;
    }
382

383
384
385
386
    inline bool operator==(const DOFIterator<T>& rhs)
    {
      return !(this->operator==(rhs));
    }
387

388
389
390
391
392
393
  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }
394

395
396
397
398
399
    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }
400

401
402
403
404
405
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }
406

407
408
409
410
411
412
413
414
415
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

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

417
418
419
    /// Iterator for \ref iteratedObject
    typename std::vector<T>::const_iterator it;
  };
420

421
} // end namespace AMDiS
422
423

#endif