Helpers.h 12.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_HELPERS_H
#define EXTENSIONS_HELPERS_H
21

22
// #include "AMDiS.h"
23

24
25
26
27
28
29
#include "FixVec.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "DOFIterator.h"

30
31
32
33
#if HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/StdMpi.h"
#endif

34
35
36
#ifdef _WIN32
#include <io.h>
#else
37
#include <unistd.h>
38
#endif
39
40
41
42
#include <ios>
#include <iostream>
#include <fstream>
#include <string>
43

44
45
46
47
#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"
48
#include <boost/math/special_functions/fpclassify.hpp>
49
#include <boost/math/special_functions/atanh.hpp>
50

51
52
53
54
#include <time.h>
#include <stdarg.h>

#include "VectorOperations.h"
55
#include "Views.h"
56
57
58
59
60
61

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

using namespace std;
using namespace AMDiS; 

62
static long random_seed_initial_value = 0;
63

64
namespace Helpers {
65

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
66
67
  // math routines
  // ==============
68
  
69
  inline double cint(double x) {
70
    double intpart; 
71
72
    if (modf(x, &intpart) >= 0.5)
      return x >=0 ? ceil(x) : floor(x);
73
    else
74
      return x < 0 ? ceil(x) : floor(x);
75
76
  }
  
77
78
  inline double round (double r, double places) {
    double off = pow(10.0, places);
79
80
81
    return cint(r*off)/off;
  }
  
82
83
84
85
  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; }
86

87
88
89
  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; }
90

91
92
93
  template<typename T>
  T atanh(T x)
  {
94
95
96
97
98
99
100
101
102
103
    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;
104
105
  }
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
106
107
  // routines for conversion to string
  // =================================
108
109
  
  template<typename T>
110
111
  inline std::string toString(const T &value, 
			      ios_base::fmtflags flag = ios_base::scientific)
112
113
114
115
116
117
118
  {
    std::stringstream ss;
    ss.setf(flag);
    ss.precision(6);
    if (!(ss<<value))
      throw(std::runtime_error("toString() not possible"));
    return ss.str();
119
  }
120
  
121
  
122
  inline std::string toString(const int &value)
123
124
  {
    return boost::lexical_cast<std::string>(value);
125
  }
126
  
127
  
128
  template<typename T>
129
130
131
  inline std::string toString(const WorldVector<T> &value, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
132
133
134
135
136
137
138
139
140
141
  {
    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();
142
  }
143
  
144
  
145
  template<typename T>
146
147
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
Praetorius, Simon's avatar
Praetorius, Simon committed
148
149
150
151
			      std::string separator = ",",
			      ios_base::fmtflags flag = ios_base::scientific,
			      size_t precision = 6
 			    )
152
153
154
155
  {
    std::stringstream ss;
    if (brackets) ss<<"[";
    ss.setf(flag);
Praetorius, Simon's avatar
Praetorius, Simon committed
156
    ss.precision(precision);
157
    for (size_t i = 0; i < vec.size(); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
158
      if (!(ss<<(i>0?(separator+" "):std::string(""))<<vec[i]))
159
160
161
162
	throw(std::runtime_error("toString() not possible"));
    }
    if (brackets) ss<<"]";
    return ss.str();
163
  }
164
  
165
  
166
  template<class T>
167
  inline T fromString(const std::string& s)
168
169
170
171
172
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
173
  }
174

175
  
176
  inline std::string fillString(int length, char c, int numValues, ...)
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
  {
    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);
196
  }
197

198
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
199
  /// produce printable string of memory
200
  inline std::string memoryToString(double mem, int startIdx = 0) {
201
202
203
204
205
206
207
208
    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;
209
  }
210
211
  
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
212
213
  // some mesh routines
  // ===========================
214
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
215
  /// rotate scale and shift the mesh coordinates
216
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate);
217

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
218
  /// scale a vector of meshes by different values in all directions
219
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
220
221
  
  /// scale and shift by different values in all directions
222
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
223
224
  
  /// scale by different values in all directions
225
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
226
227
  
  /// scale and shift by the same values in all directions
228
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
229
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
230
231
232
233
234
  /**
   * 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]
   **/
235
  WorldVector<double> getMeshDimension(Mesh *mesh);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
236
237
  
  /// calculate the dimension of a mesh, by mesh traversal.
238
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
Praetorius, Simon's avatar
Praetorius, Simon committed
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
  inline double calcMeshSizes(Mesh* mesh, double& minH, double& maxH, int& minLevel, int& maxLevel) 
  {
    FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    minH = 1e15; maxH = 0.0;
    int k = 0;
    minLevel = 100;
    maxLevel = 0;
    while (elInfo) {
      maxLevel = std::max(maxLevel,elInfo->getLevel());
      minLevel = std::min(minLevel,elInfo->getLevel());
      coords = elInfo->getCoords();
      double h = 0.0;
      for (int i = 0; i < coords.size(); i++) {
        for (int j = 0; j < coords.size(); j++) {
	  if (i != j)
            h = std::max(h,norm(coords[i]-coords[j]));
        }
      }
      minH = std::min(h, minH);
      maxH = std::max(h, maxH);
      elInfo = stack.traverseNext(elInfo);
      k++;
    }
    minLevel += mesh->getMacroElementLevel();
    maxLevel += mesh->getMacroElementLevel();
    
    return minH;
  }
  
272
  /// read DOFVector from AMDiS .dat-files
273
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
274

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
275
276
  // some linear algebra methods
  // ===========================
277

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
278
  /// calculate the determant of a 3x3 matrix of type dense2D
279
  inline double determinant(mtl::dense2D<double> &A) // only for 3x3 - matrices
280
281
282
283
284
285
  {
    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;
286
  }
287
  
288
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
289
  /// calculate the determant of a 3x3 matrix of type WorldMatrix
290
  inline double determinant(WorldMatrix<double> &A) // only for 3x3 - matrices
291
292
293
294
295
296
  {
    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;
297
  }
298
  
299
  
300
301
  /// inverse power method to find minimal eigenvalue of matrix A and appropriate eigenvector x, using m iterations
  template <class LinearSolver, class Matrix, class EigenVector>
302
  inline double inverse_iteration(LinearSolver& solver, const Matrix& A, EigenVector& x, int m)
Praetorius, Simon's avatar
Praetorius, Simon committed
303
  {
304
305
306
307
308
309
310
311
312
313

    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
314
315
316
317
318
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
319
320
321
322
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

323
324
    return lambda;
  }
325

326
327
328
329
  
  /// calculate approximation to condition number of matrix A using the LinearSolver solver
  template <class LinSolver, class Matrix>
  inline double condition (LinSolver& solver, const Matrix& A, int m=10)
Praetorius, Simon's avatar
Praetorius, Simon committed
330
  {
331
332
333
334
335
336
337
338
    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;
339
      lambda_old = lambda;
340
341
      lambda = mtl::dot(y,x);
      relErr = abs((lambda-lambda_old)/lambda);
342
343
      if (relErr < 1.e-5)
	break;
344
345
346
347
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
348
  }
349
  
350
  
351
352
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
353
  inline void interpolOverLine(const Container &vec, 
354
355
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
356
357
358
359
360
			       std::vector<WorldVector<double> > &x, std::vector<double> &y)
  { FUNCNAME("Helpers::interpolOverLine()");

    WorldVector<double> p;

Praetorius, Simon's avatar
Praetorius, Simon committed
361
    TEST_EXIT(nPoints > 1)("Zu wenig Punkte fuer eine Interpolation!\n");
362
363
364
365
366
367
368
369
370
371
372
    
    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);
    }
  }
373
  
374

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
375
  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
376
  void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
377
378
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
379
  void calcMaxOnYAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
380
381

  /// calc normal vectors of surface from element normals by averaging
382
  void getNormalsWeighted(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals);
383
384
  
  /// calc normal vectors of surface from element normals by local approximation by quartic
385
  void getNormals(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals, DOFVector<WorldMatrix<double> > *gradNormals = NULL);
386
387
  
  /// calc curvature from given normal vectors
388
  void getCurvature(DOFVector<WorldVector<double> >* normals, DOFVector<double>* curvature);
389
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
390
391
  // misc routines
  // =============
392
  
393
  void plot(std::vector<double> &values, int numRows = 7, int numCols = 20, std::string symbol = "*");
394
  
395
396
  inline long getMicroTimestamp() 
  {
397
398
399
400
401
    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();
402
  }
403
  
404
405
406
  
  inline long getRandomSeed(int param = 0) 
  {
407
408
409
410
411
412
413
414
    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++;

415
    return value0 + value1 + value2 + param;
416
  }
417

418
419
  
} // end namespace Helpers
420
421


422
#endif // EXTENSIONS_HELPERS_H