Global.h 15.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * 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 Vey, Thomas Witkowski, Andreas Naumann, 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.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20

21
22
23
24
25
26
27
28
29
30
31
32
33
34

/** \file Global.h */

/** \mainpage AMDiS
 * @{ <img src="vis.png"> @}
 */

/** \defgroup Common Common
 */

#ifndef AMDIS_GLOBAL_H
#define AMDIS_GLOBAL_H

#if (__GNUC__) && (__GNUC__ > 2) 
35
#define OPENMODE std::ios::openmode
36
#else
37
#define OPENMODE std::ios::open_mode
38
39
#endif

40
41
42
43
44
45
46
#ifndef HAS_Cpp11 //__cplusplus < 199711L
  #define override 
#endif
#ifndef HAS_Cpp11 //__cplusplus < 199711L
  #define nullptr NULL 
#endif

47
/** \brief current AMDiS version */
48
49
50
#ifndef AMDIS_VERSION
#define AMDIS_VERSION  "AMDiS: Version 0.9.1"
#endif
51
52
53

#include <string>
#include <vector>
54
#include <set>
55
56
57
58
59
60
#include <fstream>
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <functional>
#include <float.h>
61
#include <time.h>
62

63
64
65
66
67
68
#ifdef _WIN32
#include <io.h>
#else
#include <unistd.h>
#endif

69
#if HAVE_PARALLEL_DOMAIN_AMDIS
70
#include <mpi.h>
71
72
#endif

73
#include <boost/algorithm/string.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
74
#include <boost/algorithm/string/trim.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
75
#include <boost/math/special_functions/fpclassify.hpp>
76
#include "boost/tuple/tuple.hpp"
77
#include "AMDiS_fwd.h"
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
78
#include "OpenMP.h"
79

80
namespace AMDiS {
81

82
83
  using namespace std;

84
85
  extern const char *funcName;

86
  const int amdisRevisionNumber = 1700;
87

88
  /// Used by matrix vector multiplication
89
90
91
92
  typedef enum { NoTranspose,
		 Transpose,
		 ConjugateTranspose } MatrixTranspose;

93
  /// Speciefies the norm used by Estimator.
94
95
96
  typedef enum { NO_NORM = 0, 
		 H1_NORM = 1, 
		 L2_NORM = 2 } Norm;
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
  
  // may be used in the future
  struct DofIndex 
  {
    typedef signed int size_type;
    
    DofIndex() : idx(0) {}
    DofIndex(size_type i) : idx(i) {}
    
    DofIndex& operator=(const size_type &rhs) 
    {
      idx = rhs;
      return *this;
    }
    
    operator size_type() const
    {
      return idx;
    }
    
    size_type idx;
  };
  
  std::ostream& operator<<(std::ostream& os, const DofIndex& di);  
  std::istream& operator>>(std::istream& is, DofIndex& di);
122

123
  /// Datatype for degrees of freedom 
124
125
//   typedef signed int DegreeOfFreedom;
  typedef DofIndex::size_type DegreeOfFreedom;
126

127
  /// Defines type for a vector of DOF pointers.
128
  typedef vector<const DegreeOfFreedom*> DofContainer;
129

130
131
  typedef std::set<const DegreeOfFreedom*> DofContainerSet;

132
  /// Defines a type for global edge identification via its DOFs.
133
  typedef pair<DegreeOfFreedom, DegreeOfFreedom> DofEdge;
134
135
136

  /// Defines a tzpe for global face identiication via its DOFs.
  typedef boost::tuple<DegreeOfFreedom, DegreeOfFreedom, DegreeOfFreedom> DofFace;
137

138

139
  /// Returns the GeoIndex of d for dimension dim.
140
141
#define INDEX_OF_DIM(d, dim) (static_cast<GeoIndex>((d == dim) ? CENTER : d + 1))

142
  /// Returns the dimension of GeoIndex ind for dimension dim
143
144
#define DIM_OF_INDEX(ind, dim) ((static_cast<int>(ind) == 0) ? dim : static_cast<int>(ind) - 1)

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
145
146
147
148

#if SUPPRESS_OUTPUT
#define PRINT_LINE(stream, line)
#else
149
#if HAVE_PARALLEL_DOMAIN_AMDIS
150
#define PRINT_LINE(stream, line) stream << "[" << MPI::COMM_WORLD.Get_rank() << "] " << line
151
#else
152
#define PRINT_LINE(stream, line) stream << line
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
153
#endif
154
155
#endif

156
  /// Calculates factorial of i
157
158
  int fac(int i);

159
160
  void waitSec(int seconds);

161
  void processMemUsage(double& vm_usage, double& resident_set, bool inMegaByte = true);
162

163
  /// Content comparision of two pointers. Used e.g. for find_if
164
  template<typename T>
165
  struct comparePtrContents : public binary_function<T*, T*, bool>
166
  {
167
    /// Overrides binary_function::operator()
168
169
    bool operator()(T* a, T* b) const 
    {
170
      return (*a == *b);
171
    }
172
173
  };

174
  
Praetorius, Simon's avatar
Praetorius, Simon committed
175
176
177
  /// check for file existence
  inline bool file_exists(const std::string filename)
  {
178
179
180
#ifdef _WIN32
    return _access(filename.c_str(), 0) == 0;
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
181
    return access(filename.c_str(), F_OK) == 0;
182
#endif
183
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
184
185
186
187
188
189
  
  /// check for inf and nan values
  inline bool isNumber(double val)
  {
    return !boost::math::isnan(val) && !boost::math::isinf(val);
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
190

191
  
Praetorius, Simon's avatar
Praetorius, Simon committed
192
193
194
195
196
197
  /// trim std::string
  inline std::string trim(const std::string& oldStr)
  {
    std::string swap(oldStr);
    boost::algorithm::trim(swap);
    return swap;
198
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
199

200
201
202
203
204
205
206
207
208
209
210
211
  // ===== some simple template functions ====================================== 

  template<typename T> inline T abs(T a) 
  {
    return  ((a) >= 0 ? (a) : -(a));
  }

  template<typename T> inline T sqr(T a) 
  {
    return  ((a)*(a));
  }

212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
  template<typename T> inline void nullify(T &a)
  {
    a = 0.0;
  }

  template<typename T> inline void nullify(std::vector<T> &a)
  {
    typename std::vector<T>::iterator it;
    for (it = a.begin(); it != a.end(); it++)
      nullify(*it);
  }

  template<typename T> inline void nullify(mtl::dense_vector<T> &a)
  {
    T null;
    nullify(null);
    a = null;
  }

  template<typename T> inline void nullify(WorldVector<T> &a)
  {
233
234
235
236
237
238
239
    T null; nullify(null);
    a = null;
  }

  template<typename T> inline void nullify(WorldMatrix<T> &a)
  {
    T null; nullify(null);
240
241
242
    a = null;
  }

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
  // ===== some predefined values ===============================================
  const double m_e = 		2.7182818284590452354;
  const double m_log2e = 		1.4426950408889634074;
  const double m_log10e = 	0.43429448190325182765;
  const double m_ln2 = 		0.69314718055994530942;
  const double m_ln10 = 		2.30258509299404568402;
  const double m_pi = 		3.14159265358979323846;
  const double m_pi_2 = 		1.57079632679489661923;
  const double m_pi_4 = 		0.78539816339744830962;
  const double m_1_pi = 		0.31830988618379067154;
  const double m_2_pi = 		0.63661977236758134308;
  const double m_2_sqrtpi = 	1.12837916709551257390;
  const double m_sqrt2 = 		1.41421356237309504880;
  const double m_sqrt1_2 = 	0.70710678118654752440;

  // ===== tolerance for floating point comparison ==============================
#define DBL_TOL DBL_EPSILON
#define FLT_TOL FLT_EPSILON


  /** \brief
   * Manages the output of messages, warnings, errors, ...
   * Used by the macros FUNCNAME, ERROR, ERROR_EXIT, WARNING, TEST, MSG, INFO,
   * PRINT_INFO, WAIT, WAIT_REALLY.
   * Don't use this class directly but only via these macros!
   */
  class Msg
  {
  public:
272
    /// Prints a formated message to the message stream
273
274
    static void print(const char *format, ...);

275
    /// Prints a formated message to the error stream
276
277
    static void print_error(const char *format, ...);

278
    /// Prints a formated message to the error stream and exits 
279
280
    static void print_error_exit(const char *format, ...);

281
282
    ///
    static void catch_error_exit(const char *format, ...) {}
283

284
    /// Prints an error message with funcname, file, and line to the error stream
285
286
287
288
    static void print_error_funcname(const char *funcname,
				     const char *file, 
				     int line);

289
    /// Prints a warning to the message stream
290
291
    static void print_warn(const char *format, ...);

292
    /// Prints a warning with funcname, file, and line to the message stream
293
294
295
296
    static void print_warn_funcname(const char *funcname,
				    const char *file, 
				    int line);

297
    /// Prints the funcname to the message stream
298
299
    static void print_funcname(const char *funcname);

300
    /// Changes the message stream
301
    static void change_out(ostream*);
302

303
    /// Changes the error stream 
304
    static void change_error_out(ofstream *fp);
305

306
    /// Creates a filestream and sets the error stream to this filestream
307
308
    static void open_error_file(const char *filename, OPENMODE);

309
    /// Sets \ref msgInfo
310
311
312
313
    static void setMsgInfo(int info) 
    { 
      msgInfo = info; 
    }
314

315
    /// Returns \ref msgInfo
316
317
318
319
    static int  getMsgInfo() 
    { 
      return msgInfo; 
    }
320

321
    /// Sets \ref msgWait
322
323
324
325
    static void setMsgWait(bool wait) 
    { 
      msgWait = wait; 
    }
326

327
    /// Returns \ref msgWait
328
329
330
331
    static bool getMsgWait() 
    { 
      return msgWait; 
    }
332

333
    /// Waits for enter if w is true
334
335
    static void wait(bool w);

336
    /// Returns \ref out
337
    static ostream *getOutStream() 
338
339
340
    { 
      return out; 
    }
341

342
    /// Returns \ref error
343
    static ostream *getErrorStream() 
344
345
346
    { 
      return error; 
    }
347

348
349
350
351
352
353
354
355
  public:
#if HAVE_PARALLEL_DOMAIN_AMDIS
    /// In parallel computations, when this variable is true, only the 0 rank will
    /// print messages to the output stream. Error messages and warnings are always
    /// printed from all ranks.
    static bool outputMainRank;
#endif

356
  protected:
357
    /// Message stram
358
    static ostream *out;
359

360
    /// Error stream
361
    static ostream *error;
362

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
363
364
    /// Remember funcName to avoid multiple output of funcName within the same
    /// function call
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
365
    static ThreadPrivate<const char*> oldFuncName;
366

367
    /// Global info level
368
369
    static int msgInfo;

370
    /// Spezifies whether to wait when WAIT is called
371
372
373
    static bool msgWait;
  };

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
374
375
376
  // ===========================================================================
  // ===== message macros ======================================================
  // ===========================================================================
377

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
378
379
  /// Should be the first call in every functions. It defines the current 
  /// function name nn for message output via MSG, WARNING, ...
380
#define FUNCNAME(nn) const char *funcName; funcName = nn;
381
382
383
384
385
386
  
#if (DEBUG == 0) 
  #define FUNCNAME_DBG(nn)
#else
  #define FUNCNAME_DBG(nn) const char *funcName; funcName = nn;
#endif
387

388
  /// prints an error message 
389
390
391
#define ERROR Msg::print_error_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_error

392
  /// prints an error message and exits 
393
394
395
#define ERROR_EXIT Msg::print_error_funcname(funcName,__FILE__, __LINE__), \
    Msg::print_error_exit

396
  /// prints a warning
397
398
399
#define WARNING Msg::print_warn_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_warn

400
  /// if test is false, an error message is printed
401
402
#define TEST(test) if ((test));else ERROR

403
  /// if test is false, an error message is printed and the program exits
Thomas Witkowski's avatar
Thomas Witkowski committed
404
#define TEST_EXIT(test) if ((test));else ERROR_EXIT
405

406
  /// In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop.
407
408
#if (DEBUG == 0) 
  #define TEST_EXIT_DBG(test) if (false) Msg::catch_error_exit
409
  #define DBG_VAR(var)
410
411
#else
#define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT
412
  #define DBG_VAR(var) var
413
#endif
414

415
  /// prints a message
416
417
#define MSG Msg::print_funcname(funcName), Msg::print

418
419
420
421
422
423
#if (DEBUG == 0) 
  #define MSG_DBG
#else
  #define MSG_DBG Msg::print_funcname(funcName), Msg::print
#endif

424
  /// prints a message, if min(Msg::msgInfo, info) >= noinfo
425
#define INFO(info,noinfo)						\
426
  if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) MSG
427

428
  /// prints a message, if min(Msg::msgInfo, info) >= noinfo
429
#define PRINT_INFO(info,noinfo)						\
430
  if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) Msg::print
431
432


Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
433
434
435
  /// If the value of Msg::wait is not zero the macro will produce the message 
  /// 'wait for <enter> ...' and will continue after pressing the return or enter
  /// key. Otherwise the program continues without a message.
436
437
#define WAIT Msg::wait(Msg::getMsgWait())

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
438
439
  /// produces the message 'wait for <enter> ...' and will continue after
  /// pressing the return or enter key.
440
441
442
443
#define WAIT_REALLY Msg::wait(true)

#define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC)

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
444
445
  /// internal used indices to represent the different geometrical objects.
  /// Used as parameter for getGeo() and as template parameter for FixVec. 
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
  typedef enum
    {
      CENTER   = 0, /**< in 1d the center is at the edge, in 2d at the face, in 3d 
		     * at the interior of an element. So a Line has no edge but
		     * only a center, a Triangle has no face but only a center.
		     */
      VERTEX   = 1, /**< index for element vertices.
		     * number of vertices is equal to number of parts and 
		     * neighbours.
		     */
      EDGE     = 2, /**< index for element edges */
      FACE     = 3, /**< index for element faces */
      DIMEN    =-1, /**< index for problem dimension */
      PARTS    =-2, /**< index for parts of an element (vertices in 1d, edges in 2d
		     * , faces in 3d). Number of element parts is equal to number
		     * of vertices and neighbours. 
		     */
      NEIGH    =-3, /**< index for neighbours of an element.
		     * Number of element neighbours is equal to number of 
		     * vertices and parts.
		     */
      WORLD    =-4, /**< index for world dimension */
      BOUNDARY =-5, /**< index for boundary nodes of an element. This could be
		     * vertices, edges or faces.
		     */
      PROJECTION=-6 /**< index for element and boundary projections */
    } GeoIndex;

#define MAXPART FACE
#define MINPART PROJECTION


  /** \ingroup Common
   * \brief
   * Static class wich holds global information about the world and all types of
   * elements.
   */
  class Global
  {
  public:
Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
486
487
    /// returns a pointer to \ref referenceElement [dim]. With this pointer you
    /// can get information about the element via Element's getGeo method.
488
489
    static const Element *getReferenceElement(int dim) 
    {
490
      FUNCNAME("Global::getReferenceElement()");
491
      TEST_EXIT(dim > 0 && dim < 4)("invalid dim: %d\n", dim);
492
      return referenceElement[dim];
493
    }
494

495
    /// returns geometrical information. Currently this is only dimOfWorld.
496
497
    static inline int getGeo(GeoIndex p) 
    {
498
499
500
      if (WORLD == p) 
	return dimOfWorld; 

501
      ERROR_EXIT("Illegal request for geometry data: part = %d!\n", p);
502
      return 0;
503
    }
504

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
505
506
    /// Returns geometrical information about elements of the dimension dim.
    /// getGeo(VERTEX, 3) returns 4 because a Tetrahedron has 4 vertices.
507
508
    static inline int getGeo(GeoIndex p, int dim) 
    {
509
      TEST_EXIT_DBG(p >= MINPART && p <= MAXPART)
510
	("Calling for invalid geometry value %d\n",p);
511
      TEST_EXIT_DBG(dim >= 0 && dim < 4)
512
513
514
515
516
517
	("invalid dim: %d\n", dim);
      TEST_EXIT_DBG((dim != 0) || (p == PARTS || p == VERTEX || p == EDGE || p == FACE))
	("dim = 0\n");

      return geoIndexTable[dim][p - MINPART];
    }
518

519
    /// Inits the Global class with the help of Parameters.
520
521
    static void init();

522
523
    static void clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
524
  private:
Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
525
526
    /// Global is a pure static class. So the constructor is private to avoid
    /// instantiation.
Thomas Witkowski's avatar
Thomas Witkowski committed
527
    Global();
528
529

  private:
530
    /// Dimension of the simulated world
531
532
    static int dimOfWorld;

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
533
534
535
536
537
    /// contains a pointer to a Line, a Triangle, and a Tetrahedron.
    /// This allows the access to information of the concrete elements via
    /// the dimension index.
    /// referenceElement[3]->getGeo(VERTEX) gives the number of vertices of a
    /// Tetrahedron wich is 4 => no switch statement is necessary.
538
    static Element *referenceElement[4];
539

540
    /// Stores the precalculated results that should be returned by Global::getGeo.
541
    static vector<vector<int> > geoIndexTable;
542
543
544
545
546
547
  };

#define COMMA ,

  const int RescheduleErrorCode = 23;

548
549
550
551
552
553
554
555
556
  /**
   * \ingroup Assembler
   * \brief
   * Specifies the type of a FirstOrderTerm 
   */
  enum FirstOrderType {
    GRD_PSI,
    GRD_PHI
  };
557

558
559
560
561
}

#endif // AMDIS_GLOBAL_H