Helpers.h 11.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
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
102
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
382
383
384
385
386
387
388
389
390
/** \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"

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

using namespace std;
using namespace AMDiS; 

static long random_seed_initial_value=0;

class Helpers {
public:

  /// math routines
  /// ==============
  
  static double cint(double x){
    double intpart; 
    if (modf(x,&intpart)>=.5)
      return x>=0?ceil(x):floor(x);
    else
      return x<0?ceil(x):floor(x);
  }
  
  static double round(double r, double places){
    double off= pow(10.0, places);
    return cint(r*off)/off;
  }
  
  static bool isNumber(double val) { return !isnan(val) && !isinf(val); };
  static double signum(double val, double tol=1.e-10) { return (val>tol?1.0:(val<-tol?-1.0:0.0)); }
  static int signum(int val) { return (val>0?1:(val<0?-1:0)); }
  static int signum(bool val) { return (val?1:-1); }
  static int toInt(bool val) { return (val?1:0); }

  /// routines for conversion to string
  /// =================================
  
  template<typename T>
  static std::string toString(const T &value, ios_base::fmtflags flag= ios_base::scientific)
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
  };
  
  static std::string toString(const int &value)
  {
    return boost::lexical_cast<std::string>(value);
  };
  
  template<typename T>
  static std::string toString(const WorldVector<T> &value, bool brackets= true, ios_base::fmtflags flag= ios_base::scientific)
  {
    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();
  };
  
  template<typename T>
  static std::string toString(const std::vector<T> &vec, bool brackets= true, ios_base::fmtflags flag= ios_base::scientific)
  {
    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();
  };
  
  template<class T>
  static T fromString(const std::string& s)
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
  };

  static std::string fillString(int length, char c, int numValues, ...)
  {
    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);
  };

  // process printable string of memory
  static std::string memoryToString(double mem, int startIdx=0) {
    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;
  };  
  
  
  /// some mesh routines
  /// ===========================
  
  static void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);

  static void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
  // scale and shift by different values in all directions
  static void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
  // scale by different values in all directions
  static void scaleMesh(Mesh *mesh, WorldVector<double> scale);
  // scale and shift by the same values in all directions
  static void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
  
  /// calculate the dimension of a mesh
  static WorldVector<double> getMeshDimension(Mesh *mesh);
  static void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);

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

  /// copy DOFVectors that live on different meshes
  static void copyDOF(std::vector<DOFVector<double>*> dof1, std::vector<DOFVector<double>*> dof2, WorldVector<double> shift, double fillValue = 0.0);
  
  static void copyDOF(DOFVector<double>* dof1, DOFVector<double>* dof2) {
    WorldVector<double> shift; shift.set(0.0);
    std::vector<DOFVector<double>*> dof1Vec; dof1Vec.push_back(dof1);
    std::vector<DOFVector<double>*> dof2Vec; dof2Vec.push_back(dof2);
    copyDOF(dof1Vec, dof2Vec, shift);
  };

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

  static double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
  {
    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;
  };
  
  static double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
  {
    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;
  };
  
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
  static double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
  { 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
      lambda_old=lambda;
      lambda = dot(y,x)/dot(y,y);
      relErr = abs((lambda-lambda_old)/lambda);
      if(relErr<1.e-5) break;
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

    return lambda ;
  };

  /// calculate approximation to condition number of matrix A using the OEMSolver solver
  template <class LinearSolver, class Matrix>
  static double condition(LinearSolver& solver, const Matrix& A, int m=10)
  { 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;
      lambda_old=lambda;
      lambda = mtl::dot(y,x);
      relErr = abs((lambda-lambda_old)/lambda);
      if(relErr<1.e-5) break;
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
  };
  
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
  static void interpolOverLine(const Container &vec, 
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
			       std::vector<WorldVector<double> > &x, std::vector<double> &y);

  /// 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
  /// =============
  
  static void plot(std::vector<double> &values, int numRows=7, int numCols=20, std::string symbol="*");

  // 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
  static void process_mem_usage(double& vm_usage, double& resident_set);
  
  static long getMicroTimestamp() {
    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();
  };
  
  static long getRandomSeed(int param=0) {
    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++;

    return value0+value1+value2+param;
  };
};



/** \ingroup DOFAdministration
  * \brief
  * Implements a DOFIterator for a DOFIndexed<T> object
  */
template<typename T>
class DOFConstIterator : public DOFIteratorBase
{
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;
  }

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

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

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

    return false;
  }

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

protected:
  /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
  inline void goToBeginOfIteratedObject()
  {
    it = iteratedObject->begin();
  }

  /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
  inline void goToEndOfIteratedObject()
  {
    it = iteratedObject->end();
  }

  /// Implementation of DOFIteratorBase::incObjectIterator()
  inline void incObjectIterator()
  {
    ++it;
  }

  /// Implementation of DOFIteratorBase::incObjectIterator()
  inline void decObjectIterator()
  {
    --it;
  }

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

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

#include "Helpers.hh"


#endif