Helpers.h 11.7 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
36
37
38

  /// math routines
  /// ==============
  
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
61
62
63
64

  /// routines for conversion to string
  /// =================================
  
  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
145

  // process 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
158
159
160
  
  
  /// some mesh routines
  /// ===========================
  
161
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
162

163
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
164
  // scale and shift by different values in all directions
165
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
166
  // scale by different values in all directions
167
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
168
  // scale and shift by the same values in all directions
169
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
170
171
  
  /// calculate the dimension of a mesh
172
173
  WorldVector<double> getMeshDimension(Mesh *mesh);
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
174
175

  /// read DOFVector from AMDiS .dat-files
176
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
177
178
179
180

  /// some linear algebra methods
  /// ===========================

181
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
182
183
184
185
186
187
  {
    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;
188
  }
189
  
190
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
191
192
193
194
195
196
  {
    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;
197
  }
198
199
200
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
201
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
202
203
204
205
206
207
208
209
210
211
212
  { 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
213
214
215
216
217
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
218
219
220
221
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

222
223
    return lambda;
  }
224
225
226

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

  /// calculate maxima of DOFVector along line, using interpolOverLine
  static void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
  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);
  
  /// misc routines
  /// =============
  
289
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
290
291
292
293
294
295

  // 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
296
  void process_mem_usage(double& vm_usage, double& resident_set);
297
  
298
  inline long getMicroTimestamp() {
299
300
301
302
303
304
305
    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();
  };
  
306
  inline long getRandomSeed(int param = 0) {
307
308
309
310
311
312
313
314
    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++;

315
    return value0 + value1 + value2 + param;
316
317
  };

318
319
  
} // end namespace Helpers
320

321
322
323
324
325
326
327
328
namespace AMDiS {

  /** \ingroup DOFAdministration
    * \brief
    * Implements a DOFIterator for a DOFIndexed<T> object
    */
  template<typename T>
  class DOFConstIterator : public DOFIteratorBase
329
  {
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
  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;
    }
350

351
352
353
354
355
    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }
356

357
358
359
360
    inline bool operator!=(const DOFIterator<T>& rhs)
    {
      if (this->iteratedObject != rhs.iteratedObject)
	return true;
361

362
363
      if (this->it != rhs.it)
	return true;
364

365
366
      return false;
    }
367

368
369
370
371
    inline bool operator==(const DOFIterator<T>& rhs)
    {
      return !(this->operator==(rhs));
    }
372

373
374
375
376
377
378
  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }
379

380
381
382
383
384
    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }
385

386
387
388
389
390
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }
391

392
393
394
395
396
397
398
399
400
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

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

402
403
404
    /// Iterator for \ref iteratedObject
    typename std::vector<T>::const_iterator it;
  };
405

406
} // end namespace AMDiS
407
408

#endif