Global.h 14.9 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

/** \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

#include <string>
#include <vector>
38
#include <set>
39 40 41 42 43 44
#include <fstream>
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <functional>
#include <float.h>
45
#include <time.h>
46

47
#ifdef _MSC_VER
48 49 50 51 52
#include <io.h>
#else
#include <unistd.h>
#endif

53
#if HAVE_PARALLEL_DOMAIN_AMDIS
54
#include <mpi.h>
55 56
#endif

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

64
namespace AMDiS {
65 66 67

  extern const char *funcName;

68
  const int amdisRevisionNumber = 1700; // TODO: how to update this value
69

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

75
  /// Speciefies the norm used by Estimator.
76 77 78
  typedef enum { NO_NORM = 0, 
		 H1_NORM = 1, 
		 L2_NORM = 2 } Norm;
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
  
  // 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);
104

105
  /// Datatype for degrees of freedom 
106 107
//   typedef signed int DegreeOfFreedom;
  typedef DofIndex::size_type DegreeOfFreedom;
108

109
  /// Defines type for a vector of DOF pointers.
110
  typedef std::vector<const DegreeOfFreedom*> DofContainer;
111

112 113
  typedef std::set<const DegreeOfFreedom*> DofContainerSet;

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

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

120

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

124
  /// Returns the dimension of GeoIndex ind for dimension dim
125 126
#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
127 128 129 130

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

138
  /// Calculates factorial of i
139 140
  int fac(int i);

141 142
  void waitSec(int seconds);

143
  void processMemUsage(double& vm_usage, double& resident_set, bool inMegaByte = true);
144

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

156
  
Praetorius, Simon's avatar
Praetorius, Simon committed
157 158 159
  /// check for file existence
  inline bool file_exists(const std::string filename)
  {
160
#ifdef _MSC_VER
161 162
    return _access(filename.c_str(), 0) == 0;
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
163
    return access(filename.c_str(), F_OK) == 0;
164
#endif
165
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
166 167 168 169 170 171
  
  /// 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
172

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

182 183 184 185 186 187 188 189 190 191 192 193
  // ===== 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));
  }

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

  template<typename T> inline void nullify(WorldMatrix<T> &a)
  {
    T null; nullify(null);
222 223 224
    a = null;
  }

225 226 227 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
  // ===== 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:
254
    /// Prints a formated message to the message stream
255 256
    static void print(const char *format, ...);

257
    /// Prints a formated message to the error stream
258 259
    static void print_error(const char *format, ...);

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

263 264
    ///
    static void catch_error_exit(const char *format, ...) {}
265

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

271
    /// Prints a warning to the message stream
272 273
    static void print_warn(const char *format, ...);

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

279
    /// Prints the funcname to the message stream
280 281
    static void print_funcname(const char *funcname);

282
    /// Changes the message stream
283
    static void change_out(std::ostream*);
284

285
    /// Changes the error stream 
286
    static void change_error_out(std::ofstream *fp);
287

288
    /// Creates a filestream and sets the error stream to this filestream
289 290
    static void open_error_file(const char *filename, OPENMODE);

291
    /// Sets \ref msgInfo
292 293 294 295
    static void setMsgInfo(int info) 
    { 
      msgInfo = info; 
    }
296

297
    /// Returns \ref msgInfo
298 299 300 301
    static int  getMsgInfo() 
    { 
      return msgInfo; 
    }
302

303
    /// Sets \ref msgWait
304 305 306 307
    static void setMsgWait(bool wait) 
    { 
      msgWait = wait; 
    }
308

309
    /// Returns \ref msgWait
310 311 312 313
    static bool getMsgWait() 
    { 
      return msgWait; 
    }
314

315
    /// Waits for enter if w is true
316 317
    static void wait(bool w);

318
    /// Returns \ref out
319
    static std::ostream *getOutStream() 
320 321 322
    { 
      return out; 
    }
323

324
    /// Returns \ref error
325
    static std::ostream *getErrorStream() 
326 327 328
    { 
      return error; 
    }
329

330 331 332 333 334 335 336 337
  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

338
  protected:
339
    /// Message stram
340
    static std::ostream *out;
341

342
    /// Error stream
343
    static std::ostream *error;
344

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

349
    /// Global info level
350 351
    static int msgInfo;

352
    /// Spezifies whether to wait when WAIT is called
353 354 355
    static bool msgWait;
  };

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
356 357 358
  // ===========================================================================
  // ===== message macros ======================================================
  // ===========================================================================
359

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
360 361
  /// Should be the first call in every functions. It defines the current 
  /// function name nn for message output via MSG, WARNING, ...
362
#define FUNCNAME(nn) const char *funcName; funcName = nn;
363
  
Praetorius, Simon's avatar
Praetorius, Simon committed
364
#ifdef NDEBUG
365 366 367 368
  #define FUNCNAME_DBG(nn)
#else
  #define FUNCNAME_DBG(nn) const char *funcName; funcName = nn;
#endif
369

370
  /// prints an error message 
371 372 373
#define ERROR Msg::print_error_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_error

374
  /// prints an error message and exits 
375 376 377
#define ERROR_EXIT Msg::print_error_funcname(funcName,__FILE__, __LINE__), \
    Msg::print_error_exit

378
  /// prints a warning
379 380 381
#define WARNING Msg::print_warn_funcname(funcName,__FILE__, __LINE__),	\
    Msg::print_warn

382
  /// if test is false, an error message is printed
383 384
#define TEST(test) if ((test));else ERROR

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

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

397
  /// prints a message
398 399
#define MSG Msg::print_funcname(funcName), Msg::print

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

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

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


Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
415 416 417
  /// 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.
418 419
#define WAIT Msg::wait(Msg::getMsgWait())

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

Thomas Witkowski's avatar
bla  
Thomas Witkowski committed
424 425
  /// internal used indices to represent the different geometrical objects.
  /// Used as parameter for getGeo() and as template parameter for FixVec. 
426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
  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.
		     */
451 452 453
      PROJECTION=-6, /**< index for element and boundary projections */
      
      NO_INDEX =-127
454 455 456 457 458 459 460 461 462 463 464 465 466 467
    } 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
468 469
    /// returns a pointer to \ref referenceElement [dim]. With this pointer you
    /// can get information about the element via Element's getGeo method.
470 471
    static const Element *getReferenceElement(int dim) 
    {
472
      FUNCNAME("Global::getReferenceElement()");
473
      TEST_EXIT(dim > 0 && dim < 4)("invalid dim: %d\n", dim);
474
      return referenceElement[dim];
475
    }
476

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

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

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

501
    /// Inits the Global class with the help of Parameters.
502 503
    static void init();

504 505
    static void clear();

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

  private:
512
    /// Dimension of the simulated world
513 514
    static int dimOfWorld;

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

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

#define COMMA ,

  const int RescheduleErrorCode = 23;

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

540 541 542 543
}

#endif // AMDIS_GLOBAL_H