Global.h 15.1 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 defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI)
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 defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_MPI)
135
#define PRINT_LINE(stream, line) stream << "[" << MPI::COMM_WORLD.Get_rank() << "] " << line
136
137
#elif defined(HAVE_OPENMP)
#define PRINT_LINE(stream, line) stream << "[" << omp_get_thread_num() << "] " << line
138
#else
139
#define PRINT_LINE(stream, line) stream << line
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
140
#endif
141
142
#endif

143
  /// Calculates factorial of i
144
145
  int fac(int i);

146
147
  void waitSec(int seconds);

148
  void processMemUsage(double& vm_usage, double& resident_set, bool inMegaByte = true);
149

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

161

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

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

178

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

187
  // ===== some simple template functions ======================================
188

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

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

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

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

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
257
258
  // ===== 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:
259
    /// Prints a formated message to the message stream
260
261
    static void print(const char *format, ...);

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

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

268
269
    ///
    static void catch_error_exit(const char *format, ...) {}
270

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

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

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

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

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

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

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

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

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

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

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

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

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

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

335
  public:
336
#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_OPENMP) || defined(HAVE_MPI)
337
338
339
340
341
342
    /// 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

343
  protected:
344
    /// Message stram
345
    static std::ostream *out;
346

347
    /// Error stream
348
    static std::ostream *error;
349

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

354
    /// Global info level
355
356
    static int msgInfo;

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

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
361
362
363
  // ===========================================================================
  // ===== message macros ======================================================
  // ===========================================================================
364

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

509
510
    static void clear();

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

  private:
517
    /// Dimension of the simulated world
518
519
    static int dimOfWorld;

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

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

#define COMMA ,

  const int RescheduleErrorCode = 23;

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

545
546
547
548
}

#endif // AMDIS_GLOBAL_H