Helpers.h 14 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
23

#include "AMDiS.h"

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
148
  inline std::string toString(const std::vector<T> &vec, 
			      bool brackets = true, 
			      ios_base::fmtflags flag = ios_base::scientific)
149
150
151
152
153
154
155
156
157
158
159
  {
    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();
160
  }
161
  
162
  
163
  template<class T>
164
  inline T fromString(const std::string& s)
165
166
167
168
169
  {
    std::istringstream stream(s);
    T t;
    stream >> t;
    return t;
170
  }
171

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

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

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

    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;
  }
  
271
  /// read DOFVector from AMDiS .dat-files
272
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
273

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

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

322
323
    return lambda;
  }
324

325
326
327
328
  
  /// 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)
329
330
331
332
333
334
335
336
337
338
  { 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;
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
361
362
363
364
365
366
367
368
369
370
371
372
373
			       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);
    }
  }
374
  
375

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

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

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

419
420
  
} // end namespace Helpers
421

422
423
424
425
426
427
428
429
namespace AMDiS {

  /** \ingroup DOFAdministration
    * \brief
    * Implements a DOFIterator for a DOFIndexed<T> object
    */
  template<typename T>
  class DOFConstIterator : public DOFIteratorBase
430
  {
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
  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;
    }
451

452
453
454
455
456
    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }
457

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

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

466
467
      return false;
    }
468

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

474
475
476
477
478
479
  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }
480

481
482
483
484
485
    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }
486

487
488
489
490
491
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }
492

493
494
495
496
497
498
499
500
501
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

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

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

507
} // end namespace AMDiS
508

509
#endif // EXTENSIONS_HELPERS_H