Helpers.h 12.4 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
#include "operations/norm.hpp"
#include "operations/product.hpp"
58
59
60

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

61
62
namespace AMDiS { namespace extensions {
  
63
64
using namespace std;

65
static long random_seed_initial_value = 0;
66

67
namespace Helpers {
68

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

90
91
92
  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; }
93

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

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

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

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
221
  /// scale a vector of meshes by different values in all directions
222
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
223
224
  
  /// scale and shift by different values in all directions
225
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
226
227
  
  /// scale by different values in all directions
228
  void scaleMesh(Mesh *mesh, WorldVector<double> scale);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
229
230
  
  /// scale and shift by the same values in all directions
231
  void scaleMesh(Mesh *mesh, double shift=0.0, double scale=1.0);
232
  
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
233
234
235
236
237
  /**
   * 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]
   **/
238
  WorldVector<double> getMeshDimension(Mesh *mesh);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
239
240
  
  /// calculate the dimension of a mesh, by mesh traversal.
241
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
Praetorius, Simon's avatar
Praetorius, Simon committed
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
  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;
  }
  
275
  /// read DOFVector from AMDiS .dat-files
276
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
277

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
278
279
  // some linear algebra methods
  // ===========================
280

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

    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
319
320
321
322
323
      lambda_old = lambda;
      lambda = dot(y,x) / dot(y,y);
      relErr = abs((lambda - lambda_old) / lambda);
      if (relErr < 1.e-5)
	break;
324
325
326
327
      x = y / two_norm(y);
    }
    solver.setMultipleRhs(false);

328
329
    return lambda;
  }
330

331
332
333
334
  
  /// 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
335
  {
336
337
    using mtl::two_norm; using AMDiS::two_norm;
    using mtl::dot; using AMDiS::dot;
338
339
340
341
342
343
344
345
    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;
346
      lambda_old = lambda;
347
348
      lambda = dot(y,x);
      relErr = std::abs((lambda-lambda_old)/lambda);
349
350
      if (relErr < 1.e-5)
	break;
351
352
353
354
      x = y / two_norm(y) ;
    }
    result1 = std::abs(lambda / inverse_iteration(solver, A, x, m));
    return result1;
355
  }
356
  
357
  
358
359
  /// interpolate values of DOFVector along line from (x1,y1) -> (x2,y1) with nPoints points
  template<typename Container>
360
  inline void interpolOverLine(const Container &vec, 
361
362
			       double x1, double y1, double x2, double y2, 
			       int nPoints,
363
364
365
366
367
			       std::vector<WorldVector<double> > &x, std::vector<double> &y)
  { FUNCNAME("Helpers::interpolOverLine()");

    WorldVector<double> p;

Praetorius, Simon's avatar
Praetorius, Simon committed
368
    TEST_EXIT(nPoints > 1)("Zu wenig Punkte fuer eine Interpolation!\n");
369
370
371
372
373
374
375
376
377
378
379
    
    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);
    }
  }
380
  
381

Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
382
  /// calculate maxima of DOFVector along x-axis, using interpolOverLine
383
  void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
Praetorius, Simon's avatar
Helpers    
Praetorius, Simon committed
384
385
  
  /// calculate maxima of DOFVector along y-axis, using interpolOverLine
386
  void calcMaxOnYAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima);
387
388

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

422
    return value0 + value1 + value2 + param;
423
  }
424

425
426
  
} // end namespace Helpers
427

428
429
430
} }

using namespace AMDiS::extensions;
431

432
#endif // EXTENSIONS_HELPERS_H