Global.h 14.9 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
 ******************************************************************************/
20

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

/** \file Global.h */

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

/** \defgroup Common Common
 */

#ifndef AMDIS_GLOBAL_H
#define AMDIS_GLOBAL_H

34
#include "Config.h"
35

36
37
38
39
#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif
#include <cmath>
40
41
#include <string>
#include <vector>
42
#include <set>
43
44
#include <fstream>
#include <iostream>
45
#include <cstdio>
46
#include <functional>
47
48
#include <cfloat>
#include <ctime>
49

50
#ifdef _MSC_VER
51
52
53
54
55
#include <io.h>
#else
#include <unistd.h>
#endif

56
#if HAVE_PARALLEL_DOMAIN_AMDIS
57
#include <mpi.h>
58
59
#endif

60
#include <boost/algorithm/string.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
61
#include <boost/algorithm/string/trim.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
62
#include <boost/math/special_functions/fpclassify.hpp>
63
#include "boost/tuple/tuple.hpp"
64
#include "AMDiS_fwd.h"
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
65
#include "OpenMP.h"
66

67
namespace AMDiS {
68
69
70

  extern const char *funcName;

71
  const int amdisRevisionNumber = 1700; // TODO: how to update this value
72

73
  /// Used by matrix vector multiplication
74
75
76
77
  typedef enum { NoTranspose,
		 Transpose,
		 ConjugateTranspose } MatrixTranspose;

78
  /// Speciefies the norm used by Estimator.
79
80
  typedef enum { NO_NORM = 0,
		 H1_NORM = 1,
81
		 L2_NORM = 2 } Norm;
82

83
  // may be used in the future
84
  struct DofIndex
85
86
  {
    typedef signed int size_type;
87

88
89
    DofIndex() : idx(0) {}
    DofIndex(size_type i) : idx(i) {}
90
91

    DofIndex& operator=(const size_type &rhs)
92
93
94
95
    {
      idx = rhs;
      return *this;
    }
96

97
98
99
100
    operator size_type() const
    {
      return idx;
    }
101

102
103
    size_type idx;
  };
104
105

  std::ostream& operator<<(std::ostream& os, const DofIndex& di);
106
  std::istream& operator>>(std::istream& is, DofIndex& di);
107

108
  /// Datatype for degrees of freedom
109
110
//   typedef signed int DegreeOfFreedom;
  typedef DofIndex::size_type DegreeOfFreedom;
111

112
  /// Defines type for a vector of DOF pointers.
113
  typedef std::vector<const DegreeOfFreedom*> DofContainer;
114

115
116
  typedef std::set<const DegreeOfFreedom*> DofContainerSet;

117
  /// Defines a type for global edge identification via its DOFs.
118
  typedef std::pair<DegreeOfFreedom, DegreeOfFreedom> DofEdge;
119
120
121

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

123

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

127
  /// Returns the dimension of GeoIndex ind for dimension dim
128
129
#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
130
131
132
133

#if SUPPRESS_OUTPUT
#define PRINT_LINE(stream, line)
#else
134
#if HAVE_PARALLEL_DOMAIN_AMDIS
135
#define PRINT_LINE(stream, line) stream << "[" << MPI::COMM_WORLD.Get_rank() << "] " << line
136
#else
137
#define PRINT_LINE(stream, line) stream << line
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
138
#endif
139
140
#endif

141
  /// Calculates factorial of i
142
143
  int fac(int i);

144
145
  void waitSec(int seconds);

146
  void processMemUsage(double& vm_usage, double& resident_set, bool inMegaByte = true);
147

148
  /// Content comparision of two pointers. Used e.g. for find_if
149
  template<typename T>
150
  struct comparePtrContents : public std::binary_function<T*, T*, bool>
151
  {
152
    /// Overrides binary_function::operator()
153
    bool operator()(T* a, T* b) const
154
    {
155
      return (*a == *b);
156
    }
157
158
  };

159

Praetorius, Simon's avatar
Praetorius, Simon committed
160
161
162
  /// check for file existence
  inline bool file_exists(const std::string filename)
  {
163
#ifdef _MSC_VER
164
165
    return _access(filename.c_str(), 0) == 0;
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
166
    return access(filename.c_str(), F_OK) == 0;
167
#endif
168
  }
169

Praetorius, Simon's avatar
Praetorius, Simon committed
170
171
172
173
174
  /// 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
175

176

Praetorius, Simon's avatar
Praetorius, Simon committed
177
178
179
180
181
182
  /// trim std::string
  inline std::string trim(const std::string& oldStr)
  {
    std::string swap(oldStr);
    boost::algorithm::trim(swap);
    return swap;
183
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
184

185
  // ===== some simple template functions ======================================
186

187
  template<typename T> inline T abs(T a)
188
189
190
191
  {
    return  ((a) >= 0 ? (a) : -(a));
  }

192
  template<typename T> inline T sqr(T a)
193
194
195
196
  {
    return  ((a)*(a));
  }

197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
  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)
  {
218
219
220
221
222
223
224
    T null; nullify(null);
    a = null;
  }

  template<typename T> inline void nullify(WorldMatrix<T> &a)
  {
    T null; nullify(null);
225
226
227
    a = null;
  }

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
  // ===== 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:
257
    /// Prints a formated message to the message stream
258
259
    static void print(const char *format, ...);

260
    /// Prints a formated message to the error stream
261
262
    static void print_error(const char *format, ...);

263
    /// Prints a formated message to the error stream and exits
264
265
    static void print_error_exit(const char *format, ...);

266
267
    ///
    static void catch_error_exit(const char *format, ...) {}
268

269
    /// Prints an error message with funcname, file, and line to the error stream
270
    static void print_error_funcname(const char *funcname,
271
				     const char *file,
272
273
				     int line);

274
    /// Prints a warning to the message stream
275
276
    static void print_warn(const char *format, ...);

277
    /// Prints a warning with funcname, file, and line to the message stream
278
    static void print_warn_funcname(const char *funcname,
279
				    const char *file,
280
281
				    int line);

282
    /// Prints the funcname to the message stream
283
284
    static void print_funcname(const char *funcname);

285
    /// Changes the message stream
286
    static void change_out(std::ostream*);
287

288
    /// Changes the error stream
289
    static void change_error_out(std::ofstream *fp);
290

291
    /// Creates a filestream and sets the error stream to this filestream
292
293
    static void open_error_file(const char *filename, OPENMODE);

294
    /// Sets \ref msgInfo
295
296
297
    static void setMsgInfo(int info)
    {
      msgInfo = info;
298
    }
299

300
    /// Returns \ref msgInfo
301
302
303
    static int  getMsgInfo()
    {
      return msgInfo;
304
    }
305

306
    /// Sets \ref msgWait
307
308
309
    static void setMsgWait(bool wait)
    {
      msgWait = wait;
310
    }
311

312
    /// Returns \ref msgWait
313
314
315
    static bool getMsgWait()
    {
      return msgWait;
316
    }
317

318
    /// Waits for enter if w is true
319
320
    static void wait(bool w);

321
    /// Returns \ref out
322
323
324
    static std::ostream *getOutStream()
    {
      return out;
325
    }
326

327
    /// Returns \ref error
328
329
330
    static std::ostream *getErrorStream()
    {
      return error;
331
    }
332

333
334
335
336
337
338
339
340
  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

341
  protected:
342
    /// Message stram
343
    static std::ostream *out;
344

345
    /// Error stream
346
    static std::ostream *error;
347

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
348
349
    /// Remember funcName to avoid multiple output of funcName within the same
    /// function call
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
350
    static ThreadPrivate<const char*> oldFuncName;
351

352
    /// Global info level
353
354
    static int msgInfo;

355
    /// Spezifies whether to wait when WAIT is called
356
357
358
    static bool msgWait;
  };

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
359
360
361
  // ===========================================================================
  // ===== message macros ======================================================
  // ===========================================================================
362

363
  /// Should be the first call in every functions. It defines the current
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
364
  /// function name nn for message output via MSG, WARNING, ...
365
#define FUNCNAME(nn) const char *funcName; funcName = nn;
366

Praetorius, Simon's avatar
Praetorius, Simon committed
367
#ifdef NDEBUG
368
369
370
371
  #define FUNCNAME_DBG(nn)
#else
  #define FUNCNAME_DBG(nn) const char *funcName; funcName = nn;
#endif
372

373
  /// prints an error message
374
375
376
#define ERROR Msg::print_error_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_error

377
  /// prints an error message and exits
378
379
380
#define ERROR_EXIT Msg::print_error_funcname(funcName,__FILE__, __LINE__), \
    Msg::print_error_exit

381
  /// prints a warning
382
383
384
#define WARNING Msg::print_warn_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_warn

385
  /// if test is false, an error message is printed
386
387
#define TEST(test) if ((test));else ERROR

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

391
  /// In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop.
Praetorius, Simon's avatar
Praetorius, Simon committed
392
#ifdef NDEBUG
393
  #define TEST_EXIT_DBG(test) if (false) Msg::catch_error_exit
394
  #define DBG_VAR(var)
395
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
396
  #define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT
397
  #define DBG_VAR(var) var
398
#endif
399

400
  /// prints a message
401
402
#define MSG Msg::print_funcname(funcName), Msg::print

Praetorius, Simon's avatar
Praetorius, Simon committed
403
#ifdef NDEBUG
404
405
406
407
408
  #define MSG_DBG
#else
  #define MSG_DBG Msg::print_funcname(funcName), Msg::print
#endif

409
  /// prints a message, if min(Msg::msgInfo, info) >= noinfo
410
#define INFO(info,noinfo)						\
411
  if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) MSG
412

413
  /// prints a message, if min(Msg::msgInfo, info) >= noinfo
414
#define PRINT_INFO(info,noinfo)						\
415
  if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) Msg::print
416
417


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

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
423
424
  /// produces the message 'wait for <enter> ...' and will continue after
  /// pressing the return or enter key.
425
426
#define WAIT_REALLY Msg::wait(true)

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
427
  /// internal used indices to represent the different geometrical objects.
428
  /// Used as parameter for getGeo() and as template parameter for FixVec.
429
430
  typedef enum
    {
431
      CENTER   = 0, /**< in 1d the center is at the edge, in 2d at the face, in 3d
432
433
434
435
		     * 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.
436
		     * number of vertices is equal to number of parts and
437
438
439
440
441
442
443
		     * 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
444
		     * of vertices and neighbours.
445
446
		     */
      NEIGH    =-3, /**< index for neighbours of an element.
447
		     * Number of element neighbours is equal to number of
448
449
450
451
452
453
		     * 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.
		     */
454
      PROJECTION=-6, /**< index for element and boundary projections */
455

456
      NO_INDEX =-127
457
458
459
460
461
462
463
464
465
466
467
468
469
470
    } 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
471
472
    /// returns a pointer to \ref referenceElement [dim]. With this pointer you
    /// can get information about the element via Element's getGeo method.
473
    static const Element *getReferenceElement(int dim)
474
    {
475
      FUNCNAME("Global::getReferenceElement()");
476
      TEST_EXIT(dim > 0 && dim < 4)("invalid dim: %d\n", dim);
477
      return referenceElement[dim];
478
    }
479

480
    /// returns geometrical information. Currently this is only dimOfWorld.
481
    static inline int getGeo(GeoIndex p)
482
    {
483
484
      if (WORLD == p)
	return dimOfWorld;
485

486
      ERROR_EXIT("Illegal request for geometry data: part = %d!\n", p);
487
      return 0;
488
    }
489

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
490
491
    /// Returns geometrical information about elements of the dimension dim.
    /// getGeo(VERTEX, 3) returns 4 because a Tetrahedron has 4 vertices.
492
    static inline int getGeo(GeoIndex p, int dim)
493
    {
494
      TEST_EXIT_DBG(p >= MINPART && p <= MAXPART)
495
	("Calling for invalid geometry value %d\n",p);
496
      TEST_EXIT_DBG(dim >= 0 && dim < 4)
497
498
499
500
501
502
	("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];
    }
503

504
    /// Inits the Global class with the help of Parameters.
505
506
    static void init();

507
508
    static void clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
509
  private:
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
510
511
    /// Global is a pure static class. So the constructor is private to avoid
    /// instantiation.
Thomas Witkowski's avatar
Thomas Witkowski committed
512
    Global();
513
514

  private:
515
    /// Dimension of the simulated world
516
517
    static int dimOfWorld;

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
518
519
520
521
522
    /// 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.
523
    static Element *referenceElement[4];
524

525
    /// Stores the precalculated results that should be returned by Global::getGeo.
526
    static std::vector<std::vector<int> > geoIndexTable;
527
528
529
530
531
532
  };

#define COMMA ,

  const int RescheduleErrorCode = 23;

533
534
535
  /**
   * \ingroup Assembler
   * \brief
536
   * Specifies the type of a FirstOrderTerm
537
538
539
540
541
   */
  enum FirstOrderType {
    GRD_PSI,
    GRD_PHI
  };
542

543
544
545
546
}

#endif // AMDIS_GLOBAL_H