Global.h 14 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// 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

/** \file Global.h */

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

/** \defgroup Common Common
 */

#ifndef AMDIS_GLOBAL_H
#define AMDIS_GLOBAL_H

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

/** \brief current AMDiS version */
40
41
42
#ifndef AMDIS_VERSION
#define AMDIS_VERSION  "AMDiS: Version 0.9.1"
#endif
43
44
45

#include <string>
#include <vector>
46
#include <set>
47
48
49
50
51
52
#include <fstream>
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <functional>
#include <float.h>
53
#include <time.h>
54

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

59
#include <boost/algorithm/string.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
60
#include <boost/algorithm/string/trim.hpp>
61
#include "boost/tuple/tuple.hpp"
62
#include "AMDiS_fwd.h"
63

64
namespace AMDiS {
65

66
67
  using namespace std;

68
69
  extern const char *funcName;

70
  const int amdisRevisionNumber = 1700;
71

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

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

82
  /// Datatype for degrees of freedom 
83
84
  typedef signed int DegreeOfFreedom;

85
  /// Defines type for a vector of DOF pointers.
86
  typedef vector<const DegreeOfFreedom*> DofContainer;
87

88
89
  typedef std::set<const DegreeOfFreedom*> DofContainerSet;

90
  /// Defines a type for global edge identification via its DOFs.
91
  typedef pair<DegreeOfFreedom, DegreeOfFreedom> DofEdge;
92
93
94

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

96

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

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

103
#if HAVE_PARALLEL_DOMAIN_AMDIS
104
#define PRINT_LINE(stream, line) stream << "[" << MPI::COMM_WORLD.Get_rank() << "] " << line
105
#else
106
#define PRINT_LINE(stream, line) stream << line
107
108
#endif

109
  /// Calculates factorial of i
110
111
  int fac(int i);

112
113
  void waitSec(int seconds);

114
  void processMemUsage(double& vm_usage, double& resident_set, bool inMegaByte = true);
115

116
  /// Content comparision of two pointers. Used e.g. for find_if
117
  template<typename T>
118
  struct comparePtrContents : public binary_function<T*, T*, bool>
119
  {
120
    /// Overrides binary_function::operator()
121
122
    bool operator()(T* a, T* b) const 
    {
123
      return (*a == *b);
124
    }
125
126
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
  /// check for file existence
  inline bool file_exists(const std::string filename)
  {
    return access(filename.c_str(), F_OK) == 0;
  };

  /// trim std::string
  inline std::string trim(const std::string& oldStr)
  {
    std::string swap(oldStr);
    boost::algorithm::trim(swap);
    return swap;
  };

141
142
143
144
145
146
147
148
149
150
151
152
  // ===== 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));
  }

153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
  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)
  {
    T null;
    nullify(null);
    a = null;
  }

179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
  // ===== 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:
208
    /// Prints a formated message to the message stream
209
210
    static void print(const char *format, ...);

211
    /// Prints a formated message to the error stream
212
213
    static void print_error(const char *format, ...);

214
    /// Prints a formated message to the error stream and exits 
215
216
    static void print_error_exit(const char *format, ...);

217
218
    ///
    static void catch_error_exit(const char *format, ...) {}
219

220
    /// Prints an error message with funcname, file, and line to the error stream
221
222
223
224
    static void print_error_funcname(const char *funcname,
				     const char *file, 
				     int line);

225
    /// Prints a warning to the message stream
226
227
    static void print_warn(const char *format, ...);

228
    /// Prints a warning with funcname, file, and line to the message stream
229
230
231
232
    static void print_warn_funcname(const char *funcname,
				    const char *file, 
				    int line);

233
    /// Prints the funcname to the message stream
234
235
    static void print_funcname(const char *funcname);

236
    /// Changes the message stream
237
    static void change_out(ostream*);
238

239
    /// Changes the error stream 
240
    static void change_error_out(ofstream *fp);
241

242
    /// Creates a filestream and sets the error stream to this filestream
243
244
    static void open_error_file(const char *filename, OPENMODE);

245
    /// Sets \ref msgInfo
246
247
248
249
    static void setMsgInfo(int info) 
    { 
      msgInfo = info; 
    }
250

251
    /// Returns \ref msgInfo
252
253
254
255
    static int  getMsgInfo() 
    { 
      return msgInfo; 
    }
256

257
    /// Sets \ref msgWait
258
259
260
261
    static void setMsgWait(bool wait) 
    { 
      msgWait = wait; 
    }
262

263
    /// Returns \ref msgWait
264
265
266
267
    static bool getMsgWait() 
    { 
      return msgWait; 
    }
268

269
    /// Waits for enter if w is true
270
271
    static void wait(bool w);

272
    /// Returns \ref out
273
    static ostream *getOutStream() 
274
275
276
    { 
      return out; 
    }
277

278
    /// Returns \ref error
279
    static ostream *getErrorStream() 
280
281
282
    { 
      return error; 
    }
283

284
285
286
287
288
289
290
291
  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

292
  protected:
293
    /// Message stram
294
    static ostream *out;
295

296
    /// Error stream
297
    static ostream *error;
298

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
299
300
    /// Remember funcName to avoid multiple output of funcName within the same
    /// function call
301
302
    static const char *oldFuncName;

303
    /// Global info level
304
305
    static int msgInfo;

306
    /// Spezifies whether to wait when WAIT is called
307
308
309
    static bool msgWait;
  };

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
310
311
312
  // ===========================================================================
  // ===== message macros ======================================================
  // ===========================================================================
313

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

318
  /// prints an error message 
319
320
321
#define ERROR Msg::print_error_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_error

322
  /// prints an error message and exits 
323
324
325
#define ERROR_EXIT Msg::print_error_funcname(funcName,__FILE__, __LINE__), \
    Msg::print_error_exit

326
  /// prints a warning
327
328
329
#define WARNING Msg::print_warn_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_warn

330
  /// if test is false, an error message is printed
331
332
#define TEST(test) if ((test));else ERROR

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

336
  /// In debug mode, it corresponds to ERROR_EXIT, otherwise it is noop.
337
338
#if (DEBUG == 0) 
  #define TEST_EXIT_DBG(test) if (false) Msg::catch_error_exit
339
  #define DBG_VAR(var)
340
341
#else
#define TEST_EXIT_DBG(test) if ((test));else ERROR_EXIT
342
  #define DBG_VAR(var) var
343
#endif
344

345
  /// prints a message
346
347
#define MSG Msg::print_funcname(funcName), Msg::print

348
349
350
351
352
353
#if (DEBUG == 0) 
  #define MSG_DBG
#else
  #define MSG_DBG Msg::print_funcname(funcName), Msg::print
#endif

354
  /// prints a message, if min(Msg::msgInfo, info) >= noinfo
355
#define INFO(info,noinfo)						\
356
  if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) MSG
357

358
  /// prints a message, if min(Msg::msgInfo, info) >= noinfo
359
#define PRINT_INFO(info,noinfo)						\
360
  if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) Msg::print
361
362


Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
363
364
365
  /// 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.
366
367
#define WAIT Msg::wait(Msg::getMsgWait())

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
368
369
  /// produces the message 'wait for <enter> ...' and will continue after
  /// pressing the return or enter key.
370
371
372
373
#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
374
375
  /// internal used indices to represent the different geometrical objects.
  /// Used as parameter for getGeo() and as template parameter for FixVec. 
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
  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
416
417
    /// returns a pointer to \ref referenceElement [dim]. With this pointer you
    /// can get information about the element via Element's getGeo method.
418
419
    static const Element *getReferenceElement(int dim) 
    {
420
      FUNCNAME("Global::getReferenceElement()");
421
      TEST_EXIT(dim > 0 && dim < 4)("invalid dim: %d\n", dim);
422
      return referenceElement[dim];
423
    }
424

425
    /// returns geometrical information. Currently this is only dimOfWorld.
426
427
    static inline int getGeo(GeoIndex p) 
    {
428
429
430
      if (WORLD == p) 
	return dimOfWorld; 

431
      ERROR_EXIT("Illegal request for geometry data: part = %d!\n", p);
432
      return 0;
433
    }
434

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
435
436
    /// Returns geometrical information about elements of the dimension dim.
    /// getGeo(VERTEX, 3) returns 4 because a Tetrahedron has 4 vertices.
437
438
    static inline int getGeo(GeoIndex p, int dim) 
    {
439
      TEST_EXIT_DBG(p >= MINPART && p <= MAXPART)
440
	("Calling for invalid geometry value %d\n",p);
441
      TEST_EXIT_DBG(dim >= 0 && dim < 4)
442
443
444
445
446
447
	("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];
    }
448

449
    /// Inits the Global class with the help of Parameters.
450
451
    static void init();

452
453
    static void clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
454
  private:
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
455
456
    /// Global is a pure static class. So the constructor is private to avoid
    /// instantiation.
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    Global();
458
459

  private:
460
    /// Dimension of the simulated world
461
462
    static int dimOfWorld;

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
463
464
465
466
467
    /// 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.
468
    static Element *referenceElement[4];
469

470
    /// Stores the precalculated results that should be returned by Global::getGeo.
471
    static vector<vector<int> > geoIndexTable;
472
473
474
475
476
477
  };

#define COMMA ,

  const int RescheduleErrorCode = 23;

478
479
480
481
482
483
484
485
486
  /**
   * \ingroup Assembler
   * \brief
   * Specifies the type of a FirstOrderTerm 
   */
  enum FirstOrderType {
    GRD_PSI,
    GRD_PHI
  };
487

488
489
490
491
}

#endif // AMDIS_GLOBAL_H