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
  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;
258
259
      for (int i = 0; i < coords.getSize(); i++) {
        for (int j = 0; j < coords.getSize(); j++) {
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
	  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