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

Praetorius, Simon's avatar
Helpers  
Praetorius, Simon committed
272
273
  // some linear algebra methods
  // ===========================
274

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

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

320
321
    return lambda;
  }
322

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

    WorldVector<double> p;

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

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

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

412
    return value0 + value1 + value2 + param;
413
  }
414

415
416
  
} // end namespace Helpers
417

418
419
420
421
422
423
424
425
namespace AMDiS {

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

448
449
450
451
452
    /// Dereference operator
    inline const T* operator->()
    {
      return &(*it);
    }
453

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

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

462
463
      return false;
    }
464

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

470
471
472
473
474
475
  protected:
    /// Implementation of DOFIteratorBase::goToBeginOfIteratedObject()
    inline void goToBeginOfIteratedObject()
    {
      it = iteratedObject->begin();
    }
476

477
478
479
480
481
    /// Implementation of DOFIteratorBase::goToEndOfIteratedObject()
    inline void goToEndOfIteratedObject()
    {
      it = iteratedObject->end();
    }
482

483
484
485
486
487
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void incObjectIterator()
    {
      ++it;
    }
488

489
490
491
492
493
494
495
496
497
    /// Implementation of DOFIteratorBase::incObjectIterator()
    inline void decObjectIterator()
    {
      --it;
    }

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

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

503
} // end namespace AMDiS
504

505
#endif // EXTENSIONS_HELPERS_H