FixVec.h 20 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file FixVec.h */

#ifndef AMDIS_FIXVEC_H
#define AMDIS_FIXVEC_H

// ===========================================================================
// ===== includes ============================================================
// ===========================================================================

#include "Global.h"
#include "MemoryManager.h"
#include <iostream>
#include "MatrixVector.h"

namespace AMDiS {

  // ===========================================================================
  // ===== forward declarations ================================================
  // ===========================================================================

  class Mesh;
41
  template<typename T> class DimMat;
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
  template<typename T> class WorldVector;
  template<typename T> class WorldMatrix;


  // ===========================================================================
  // ===== definitions =========================================================
  // ===========================================================================

  /** determines how to initialize a FixVec when constructed */
  enum InitType 
    {
      NO_INIT = 0,       /**< no initialisation */
      VALUE_LIST = 1,    /**< a complete value list is given */
      DEFAULT_VALUE = 2, /**< all values ar set to a given default value */
      UNIT_VECTOR = 3,   /**< the i-th value is 1, all other values are 0 */
      UNIT_MATRIX = 4    /**< values at the diagonal of a matrix are set to one */
    };

  // ===========================================================================
  // ===== class FixVec ========================================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * A FixVec is a template vector of a fixed size. 
   *
   * The size is determined at construction time and depends on the dimension
   * and the template parameter d. So a FixVec<int, VERTEX> is a integer vector
   * with 3 entries in 2d and with 4 entries in 3d. The dimension and the way
   * the vector should be initialized are specified by the constructor call.
   */
  template<typename T,GeoIndex d>
  class FixVec : public Vector<T>
  {
  public:
    MEMORY_MANAGED(FixVec<T COMMA d>);

    /** \brief
     * constructor without initialisation. initType must be NO_INIT. If dim is
     * not spezified, a FixVec for DIM_OF_WORLD is created.
     */
    FixVec(int dim = -1, InitType initType = NO_INIT)
      : Vector<T>(calcSize(dim))
    { 
86
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or missing initializer\n");
87
88
89
90
91
92
93
94
95
    }

    /** \brief
     * constructor with value list initialisation. initType must be VALUE_LIST.
     * ini is an array which contains the initialisation values.
     */
    FixVec(int dim, InitType initType, const T* ini)
      : Vector<T>(calcSize(dim))
    {
96
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
97
      setValues(ini);
98
    }
99
100
101
102
103
104
105
106

    /** \brief
     * constructor with default value initialisation. initType must be
     * DEFAULT_VALUE. All vector entries are set to ini.
     */
    FixVec(int dim, InitType initType, const T& ini)
      : Vector<T>(calcSize(dim))
    {
107
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
108
109
110
      this->set(ini);
    }

111
112
    /// Initialisation for dim.
    inline void init(int dim) {
113
      this->resize(calcSize(dim));
114
    }
115

116
117
    /// Initialisation for size
    inline void initSize(int size) {
118
      this->resize(size);
119
    }
120

121
    /// Returns the \ref size_ of the FixVec.
122
123
    inline int size() const { 
      return this->getSize(); 
124
    } 
125
126

  protected:
127
    /// Determines needed vector size.
128
129
130
131
132
133
    static int calcSize(int dim) {
      if (dim < 0) {
	return Global::getGeo(WORLD);
      } else {
	return Global::getGeo(d, dim);
      }
134
    }
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

  public:
    friend class GLWindow;
  };


  // ===========================================================================
  // ===== class VectorOfFixVecs ===============================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * Contains an vector of FixVecs of same type. To simply allocate an array of 
   * FixVecs
   * isn't possible, because the FixVec constructor normally needs at least
   * the corresponding dimension. So you must create an array of FixVec pointers
   * and call the constructor of each FixVec manually. When you use 
   * VectorOfFixVecs, this work is done by VectorOfFixVecs's constructor.
   */
  template<typename FixVecType>
  class VectorOfFixVecs
  {
  public:
    MEMORY_MANAGED(VectorOfFixVecs<FixVecType>);

    /** \brief
     * constructs a VectorOfFixVecs without initialisation. dim is passed to 
     * FixVec's constructors. size_ is the number of contained FixVecs. initType
     * must be NO_INIT.
     */
165
166
167
    VectorOfFixVecs(int d, int s, InitType initType) 
      : size(s),
	dim(d)
168
    {
169
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
170

171
172
173
      vec.resize(size);
      for (int i = 0; i < size; i++)
	vec[i] = NEW FixVecType(dim, NO_INIT);
174
175
176
177
178
179
180
    }

    /** \brief
     * constructs a VectorOfFixVecs via an value list.  dim is passed to 
     * FixVec's constructors. size_ is the number of contained FixVecs. initType
     * must be VALUE_LIST. ini contains the initialisation values.
     */
181
182
183
    VectorOfFixVecs(int d, int s, InitType initType, const FixVecType* ini)
      : size(s),
	dim(d)
184
    {
185
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
186

187
188
189
      vec.resize(size);
      for (int i = 0; i < size; i++)
	vec[i] = NEW FixVecType(ini[i]);
190
191
192
193
194
195
196
    }

    /** \brief
     * constructs a VectorOfFixVecs with an default value.  dim is passed to 
     * FixVec's constructors. size_ is the number of contained FixVecs. initType
     * must be DEFAULT_VALUE. All entries are set to ini.
     */
197
198
199
    VectorOfFixVecs(int d, int s, InitType initType, const FixVecType& ini)
      : size(s),
	dim(d)
200
    {
201
202
203
204
205
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");

      vec.resize(size);
      for (int i = 0; i < size; i++) 
	vec[i] = NEW FixVecType(ini);
206
207
    }

208
    /// Copy constructor
209
210
    VectorOfFixVecs(const VectorOfFixVecs<FixVecType>& rhs)
    {
211
212
      size = rhs.getSize();
      dim = rhs.getDim();
213

214
215
      vec.resize(size);
      for (int i = 0; i < size; i++) 
216
	vec[i] = NEW FixVecType(*(rhs.vec[i]));
217
    }
218

219
    /// Destructor
220
221
    virtual ~VectorOfFixVecs()
    {
222
223
      for (int i = 0; i < size; i++)
	DELETE vec[i];
224

225
226
227
228
      vec.clear();
    }

    /// Allows the access like in a normal array via []. Used for const objects.
229
230
    inline const FixVecType& operator[](int index) const
    {
231
      TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
232
      return *(vec[index]);
233
    }
234

235
    /// Allows the access like in a normal array via []. 
236
237
    inline FixVecType& operator[](int index)
    {
238
      TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
239
      return *(vec[index]);
240
    }
241

242
    /// Assignment operator
243
244
245
    VectorOfFixVecs<FixVecType>& 
    operator=(const VectorOfFixVecs<FixVecType>& rhs)
    {
246
247
248
249
      TEST_EXIT_DBG(size == rhs.size)("vectors of different size\n");
      if (this != &rhs) {
	for (int i = 0; i < size; i++)
	  *(vec[i]) = *(rhs.vec[i]);
250
251
      }
      return *this;
252
    }
253

254
255
256
257
258
259
260
261
262
263
    /// Resize the vector
    inline void resize(int newsize)
    {
      vec.resize(newsize);
      for (int i = size; i < newsize; i++)
	vec[i] = NEW FixVecType(dim, NO_INIT);
      size = newsize;
    }

    /// Returns the \ref size of this VectorOfFixVecs
264
    inline int getSize() const { 
265
266
      return size;
    }
267

268
269
270
271
272
273
    /// Returns \ref dim
    inline int getDim() const {
      return dim;
    }

    /// Returns the size of the contained FixVecs
274
275
276
    inline int getSizeOfFixVec() const { 
      return vec[0]->getSize(); 
    }
277
278

  protected:
279
280
    /// number of contained FixVecs
    int size;
281

282
283
284
285
286
    /// dimension of world
    int dim;

    /// array of pointers to FixVecs
    std::vector<FixVecType*> vec;
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
  };

  // ===========================================================================
  // ===== class MatrixOfFixVecs ===============================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * Like the class VectorOfFixVecs contains a vector of FixVecs, this class
   * contains a matrix of FixVecs of same type.
   */
  template<typename FixVecType>
  class MatrixOfFixVecs
  {
  public:
    MEMORY_MANAGED(MatrixOfFixVecs<FixVecType>);

    /** \brief
     * Constructs the matrix without initialisation. r is the number of rows,
     * c is the number of columns. The other parameters are analog to the
     * VectorOfFixVecs constructors.
     */
    MatrixOfFixVecs(int dim, int r, int c, InitType initType)
      : rows(r), columns(c)
    {
312
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
313
      vec = GET_MEMORY(VectorOfFixVecs<FixVecType>*, rows);
314
      for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) {
315
316
	*i = NEW VectorOfFixVecs<FixVecType>(dim, columns, NO_INIT);
      }
317
    }
318

319
    /// destructor
320
321
    virtual ~MatrixOfFixVecs()
    {
322
      for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) {
323
324
325
	DELETE *i;
      }
      FREE_MEMORY(vec, VectorOfFixVecs<FixVecType>*, rows);
326
    }
327

328
    /// assignment operator
329
330
    inline VectorOfFixVecs<FixVecType >& operator[](int index)
    {
331
      TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
332
      return *(vec[index]);
333
    }
334

335
    /// assignment operator
336
337
    inline const VectorOfFixVecs<FixVecType >& operator[](int index) const
    {
338
      TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
339
      return *(vec[index]);
340
    }
341

342
    /// Returns \ref rows
343
344
    inline int getNumberOfRows() const { 
      return rows; 
345
    }
346

347
    /// Returns \ref columns
348
349
    inline int getNumberOfColumns() const { 
      return columns; 
350
    }
351
352

  private:
353
    /// number of rows of the matrix
354
355
    int rows;

356
    /// number of columns of the matrix
357
358
    int columns;

359
    /// array of pointers to VectorOfFixVecs
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
    VectorOfFixVecs<FixVecType> **vec;
  };


  // ===========================================================================
  // ===== class DimVec ========================================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * A DimVec is a FixVec with dim + 1 entries. It can be used for storing
   * barycentric coordinates or information associated to vertices or
   * parts of an element.
   */
  template<typename T>
  class DimVec : public FixVec<T,PARTS> {
  public:
    MEMORY_MANAGED(DimVec<T>);

379
    DimVec() {}
380

381
382
    /// Calls the corresponding constructor of FixVec
    DimVec(int dim, InitType initType = NO_INIT)
383
      : FixVec<T,PARTS>(dim, initType)
384
    {}
385

386
    /// Calls the corresponding constructor of FixVec
387
388
    DimVec(int dim, InitType initType, T* ini)
      : FixVec<T,PARTS>(dim, initType, ini)
389
    {}
390

391
    /// Calls the corresponding constructor of FixVec
392
393
    DimVec(int dim, InitType initType, const T& ini)
      : FixVec<T,PARTS>(dim, initType, ini)
394
395
    {}

396
    /// Multplication of a matrix with a vector.
397
    void multMatrixVec(DimMat<T> &m, DimVec<T> &v);
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
  };

  // ===========================================================================
  // ===== class DimMat ========================================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * A DimMat is a VectorOfFixVecs of dim+1 DimVecs. 
   */
  template<typename T>
  class DimMat : public Matrix<T>
  {
  public:
    MEMORY_MANAGED(DimMat<T>);

    /** \brief
     * Calls the corresponding constructor of VectorOfFixVecs
     */
417
418
    DimMat(int dim, InitType initType = NO_INIT)
      : Matrix<T>(dim + 1, dim + 1)
419
    {}
420
421
422
423
424

    /** \brief
     * Calls the corresponding constructor of VectorOfFixVecs
     */
    DimMat(int dim, InitType initType, const T& ini)
425
      : Matrix<T>(dim + 1, dim + 1)
426
    {
427
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)
428
429
	("wrong initType or wrong initializer\n");    
      this->set(ini);
430
    }
431
432
433
434
435

    /** \brief
     * Calls the corresponding constructor of VectorOfFixVecs
     */
    DimMat(int dim, InitType initType, T* ini)
436
      : Matrix<T>(dim + 1, dim + 1)
437
    {
438
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
439
      setValues(ini);
440
    }
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
  };

  // ===========================================================================
  // ===== class WorldVector ===================================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * A WorldVector is an AlgoVec with DIM_OF_WORLD entries of type double.
   * Can be used for storing world coordinates.
   */
  template<typename T>
  class WorldVector : public FixVec<T, WORLD>
  {
  public:
    MEMORY_MANAGED(WorldVector<T>);

    /** \brief
     * Calls the corresponding constructor of AlgoVec
     */
    WorldVector() 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), NO_INIT) 
463
    {}
464
465
466
467
468
469

    /** \brief
     * Calls the corresponding constructor of AlgoVec
     */
    WorldVector(InitType initType, T* ini) 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
470
    {}
471
472
473
474
475
476

    /** \brief
     * Calls the corresponding constructor of AlgoVec
     */
    WorldVector(InitType initType, const T& ini)
      : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
477
    {}
478
479
480
481
482
483
484
485

    /** \brief
     * Sets all entries to d
     */
    inline const WorldVector<T>& operator=(const T& d)
    {
      this->set(d);
      return (*this);
486
    }
487

488
489
490
491
492
493
494
495
496
497
498
499
500
    /** \brief
     * Sets the arrays value to the geometric midpoint of the points
     * p1 and p2;
     */
    inline void setMidpoint(const WorldVector<T> &p1, const WorldVector<T> &p2)
    {
      int dow = Global::getGeo(WORLD);

      for (int i = 0; i < dow; i++) {
	this->valArray[i] = 0.5 * (p1[i] + p2[i]);
      }
    }

501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
    /** \brief
     * Multplication of a matrix with a vector.
     */
    void multMatrixVec(WorldMatrix<T> &m, WorldVector<T> &v);
  };


  // ===========================================================================
  // ===== class WorldMatrix ===================================================
  // ===========================================================================

  /** \ingroup Common
   * \brief
   * A WorldMatrix is a FixVec with DIM_OF_WORLD WorldVector entries.
   * So it can be seen as matrix with DIM_OF_WORLD x DIM_OF_WORLD double 
   * entries.
   * Here it's not necessary to use the VectorOfFixVecs class, because the 
   * FixVec constructor assumes dim = DIM_OF_WORLD, if no dim is spezified, 
   * what is the right assumption in this case.
   */
  template<typename T>
  class WorldMatrix : public Matrix<T>
  {
  public:
    MEMORY_MANAGED(WorldMatrix<T>);

527
    /// Calls the corresponding constructor of FixVec
528
529
    WorldMatrix()
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
530
    {}
531

532
    /// Calls the corresponding constructor of FixVec
533
534
535
    WorldMatrix(InitType initType, T* ini)
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
    {
536
      TEST_EXIT_DBG(initType == VALUE_LIST)("???\n");
537
      setValues(ini);
538
    }
539
540
541
542
543
544
545
546

    /** \brief
     * Calls the corresponding constructor of FixVec and sets all matrix entries
     * to ini
     */
    WorldMatrix(InitType initType, const T& ini)
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
    {
547
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
548
      this->set(ini);
549
    }
550
  
551
    /// Returns true if the matrix is a diagonal matrix, returns false otherwise.
552
553
    bool isDiagMatrix() const;

554
    /// Returns true if the matrix is symmetric, returns false otherwise.
555
556
    bool isSymmetric() const;

557
    /// Sets the diagonal entries to the given value.
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
    void setDiag(T value);

    /** \brief
     * Creates a matrix from the "multiplication" of two vectors.
     *
     * 2d-Example:
     *
     * /a\   /c\   /ac ad\
     * | | * | | = |     |
     * \b/   \d/   \bc bd/
     */
    void vecProduct(const WorldVector<T>& v1, const WorldVector<T>& v2);
  };


  // ===========================================================================
  // ===== global functions ====================================================
  // ===========================================================================


578
  /// returns the euclidian distance of a and b
579
580
581
  template<typename T, GeoIndex d>
  double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b);

582
  /// FixVec operator for stream output
583
  template<typename T, GeoIndex d>
584
  std::ostream& operator <<(std::ostream& out, const FixVec<T,d>& fixvec)
585
  {
586
    for (int i = 0; i < fixvec.getSize()-1; i++) {
587
588
589
590
591
592
      out << fixvec[i] << " ";
    }
    out << fixvec[fixvec.getSize()-1];
    return out;
  }

593
  /// creates and inits a VectorOfFixVecs<DimVec<double> >
594
595
  VectorOfFixVecs<DimVec<double> > *createAndInit(int dim, int size, ...);

596
  /// creates and inits and double array
597
598
599
600
601
602
  double *createAndInitArray(int size, ...); 

  inline WorldVector<double> operator*(const WorldVector<double>& v, double d) {
    WorldVector<double> result = v;
    result *= d;
    return result;
603
  }
604
605
606
607
608
609
610

  inline WorldVector<double> operator+(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result += v2;
    return result;
611
  }
612
613
614
615
616
617
618

  inline WorldVector<int> operator+(const WorldVector<int>& v1,
				    const WorldVector<int>& v2) 
  {
    WorldVector<int> result = v1;
    result += v2;
    return result;
619
  }
620
621
622
623
624
625
626

  inline WorldVector<double> operator-(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result -= v2;
    return result;
627
  }
628
629
630
631

  inline bool operator<(const WorldVector<double>& v1,
			const WorldVector<double>& v2) 
  {
632
633
634
635
    int dow = Global::getGeo(WORLD);
    for (int i = 0; i < dow; i++) {
      if (abs(v1[i] - v2[i]) < DBL_TOL) 
	continue;
636
637
638
      return v1[i] < v2[i];
    }
    return false;
639
  }
640
641
642
643

  inline bool operator==(const WorldVector<double>& v1,
			 const WorldVector<double>& v2) 
  {
644
645
646
    int dow = Global::getGeo(WORLD);
    for (int i = 0; i < dow; i++) {
      if (abs(v1[i] - v2[i]) > DBL_TOL) return false;
647
648
    }
    return true;
649
  }
650
651
652
653
654
655
656
657
658
659
660
661
662
663

  template<typename T>
  const WorldMatrix<T>& operator*=(WorldMatrix<T>& m, T scal);

  template<typename T>
    const WorldMatrix<T>& operator-=(WorldMatrix<T>& m1, const WorldMatrix<T>& m2);

  template<typename T>
    const WorldMatrix<T>& operator+=(WorldMatrix<T>& m1, const WorldMatrix<T>& m2);
}

#include "FixVec.hh"

#endif // AMDIS_FIXVEC_H