FixVec.h 18 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
// ============================================================================
// ==                                                                        ==
// == 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

#include <iostream>
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include "Global.h"
27
28
29
30
31
#include "MatrixVector.h"

namespace AMDiS {

  class Mesh;
32
  template<typename T> class DimMat;
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
  template<typename T> class WorldVector;
  template<typename T> class WorldMatrix;

  /** 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 */
    };

  /** \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:
    /** \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))
    { 
66
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or missing initializer\n");
67
68
69
70
71
72
73
74
75
    }

    /** \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))
    {
76
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
77
      setValues(ini);
78
    }
79
80
81
82
83
84
85
86

    /** \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))
    {
87
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
88
89
90
      this->set(ini);
    }

91
    /// Initialisation for dim.
92
93
    inline void init(int dim) 
    {
94
      this->resize(calcSize(dim));
95
    }
96

97
    /// Initialisation for size
98
99
    inline void initSize(int size) 
    {
100
      this->resize(size);
101
    }
102

103
    /// Returns the \ref size_ of the FixVec.
104
105
    inline int size() const 
    { 
106
      return this->getSize(); 
107
    } 
108
109

  protected:
110
    /// Determines needed vector size.
111
112
113
    static int calcSize(int dim) 
    {
      if (dim < 0)
114
	return Global::getGeo(WORLD);
115
      else
116
	return Global::getGeo(d, dim);
117
    }
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

  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:
    /** \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.
     */
146
147
148
    VectorOfFixVecs(int d, int s, InitType initType) 
      : size(s),
	dim(d)
149
    {
150
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
151

152
153
      vec.resize(size);
      for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
154
	vec[i] = new FixVecType(dim, NO_INIT);
155
156
157
158
159
160
161
    }

    /** \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.
     */
162
163
164
    VectorOfFixVecs(int d, int s, InitType initType, const FixVecType* ini)
      : size(s),
	dim(d)
165
    {
166
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
167

168
169
      vec.resize(size);
      for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
170
	vec[i] = new FixVecType(ini[i]);
171
172
173
174
175
176
177
    }

    /** \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.
     */
178
179
180
    VectorOfFixVecs(int d, int s, InitType initType, const FixVecType& ini)
      : size(s),
	dim(d)
181
    {
182
183
184
185
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");

      vec.resize(size);
      for (int i = 0; i < size; i++) 
Thomas Witkowski's avatar
Thomas Witkowski committed
186
	vec[i] = new FixVecType(ini);
187
188
    }

189
    /// Copy constructor
190
191
    VectorOfFixVecs(const VectorOfFixVecs<FixVecType>& rhs)
    {
192
193
      size = rhs.getSize();
      dim = rhs.getDim();
194

195
196
      vec.resize(size);
      for (int i = 0; i < size; i++) 
Thomas Witkowski's avatar
Thomas Witkowski committed
197
	vec[i] = new FixVecType(*(rhs.vec[i]));
198
    }
199

200
    /// Destructor
201
202
    virtual ~VectorOfFixVecs()
    {
203
      for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
204
	delete vec[i];
205

206
207
208
209
      vec.clear();
    }

    /// Allows the access like in a normal array via []. Used for const objects.
210
211
    inline const FixVecType& operator[](int index) const
    {
212
      TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
213
      return *(vec[index]);
214
    }
215

216
    /// Allows the access like in a normal array via []. 
217
218
    inline FixVecType& operator[](int index)
    {
219
      TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
220
      return *(vec[index]);
221
    }
222

223
    /// Assignment operator
224
225
226
    VectorOfFixVecs<FixVecType>& 
    operator=(const VectorOfFixVecs<FixVecType>& rhs)
    {
227
228
229
230
      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]);
231
232
      }
      return *this;
233
    }
234

235
236
237
238
239
    /// Resize the vector
    inline void resize(int newsize)
    {
      vec.resize(newsize);
      for (int i = size; i < newsize; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
240
	vec[i] = new FixVecType(dim, NO_INIT);
241
242
243
244
      size = newsize;
    }

    /// Returns the \ref size of this VectorOfFixVecs
245
    inline int getSize() const { 
246
247
      return size;
    }
248

249
250
251
252
253
254
    /// Returns \ref dim
    inline int getDim() const {
      return dim;
    }

    /// Returns the size of the contained FixVecs
255
256
257
    inline int getSizeOfFixVec() const { 
      return vec[0]->getSize(); 
    }
258
259

  protected:
260
261
    /// number of contained FixVecs
    int size;
262

263
264
265
266
267
    /// dimension of world
    int dim;

    /// array of pointers to FixVecs
    std::vector<FixVecType*> vec;
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
  };

  // ===========================================================================
  // ===== 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:
    /** \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)
    {
291
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
      vec = new VectorOfFixVecs<FixVecType>*[rows];
      for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
294
	*i = new VectorOfFixVecs<FixVecType>(dim, columns, NO_INIT);
295
    }
296

297
    /// destructor
298
299
    virtual ~MatrixOfFixVecs()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
300
      for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
301
	delete *i;
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303

      delete [] vec;
304
    }
305

306
    /// assignment operator
307
308
    inline VectorOfFixVecs<FixVecType >& operator[](int index)
    {
309
      TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
310
      return *(vec[index]);
311
    }
312

313
    /// assignment operator
314
315
    inline const VectorOfFixVecs<FixVecType >& operator[](int index) const
    {
316
      TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
317
      return *(vec[index]);
318
    }
319

320
    /// Returns \ref rows
321
322
    inline int getNumberOfRows() const { 
      return rows; 
323
    }
324

325
    /// Returns \ref columns
326
327
    inline int getNumberOfColumns() const { 
      return columns; 
328
    }
329
330

  private:
331
    /// number of rows of the matrix
332
333
    int rows;

334
    /// number of columns of the matrix
335
336
    int columns;

337
    /// array of pointers to VectorOfFixVecs
338
339
340
341
342
343
344
345
346
347
348
349
350
    VectorOfFixVecs<FixVecType> **vec;
  };


  /** \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:
351
    DimVec() {}
352

353
354
    /// Calls the corresponding constructor of FixVec
    DimVec(int dim, InitType initType = NO_INIT)
355
      : FixVec<T,PARTS>(dim, initType)
356
    {}
357

358
    /// Calls the corresponding constructor of FixVec
359
360
    DimVec(int dim, InitType initType, T* ini)
      : FixVec<T,PARTS>(dim, initType, ini)
361
    {}
362

363
    /// Calls the corresponding constructor of FixVec
364
365
    DimVec(int dim, InitType initType, const T& ini)
      : FixVec<T,PARTS>(dim, initType, ini)
366
367
    {}

368
    /// Multplication of a matrix with a vector.
369
    void multMatrixVec(DimMat<T> &m, DimVec<T> &v);
370
371
372
373
374
375
376
377
378
379
380
  };


  /** \ingroup Common
   * \brief
   * A DimMat is a VectorOfFixVecs of dim+1 DimVecs. 
   */
  template<typename T>
  class DimMat : public Matrix<T>
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    /// Calls the corresponding constructor of VectorOfFixVecs
382
383
    DimMat(int dim, InitType initType = NO_INIT)
      : Matrix<T>(dim + 1, dim + 1)
384
    {}
385

Thomas Witkowski's avatar
Thomas Witkowski committed
386
    /// Calls the corresponding constructor of VectorOfFixVecs
387
    DimMat(int dim, InitType initType, const T& ini)
388
      : Matrix<T>(dim + 1, dim + 1)
389
    {
390
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)
391
392
	("wrong initType or wrong initializer\n");    
      this->set(ini);
393
    }
394

Thomas Witkowski's avatar
Thomas Witkowski committed
395
    /// Calls the corresponding constructor of VectorOfFixVecs
396
    DimMat(int dim, InitType initType, T* ini)
397
      : Matrix<T>(dim + 1, dim + 1)
398
    {
399
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
400
      setValues(ini);
401
    }
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
  };

  // ===========================================================================
  // ===== 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:
Thomas Witkowski's avatar
Thomas Witkowski committed
417
    /// Calls the corresponding constructor of AlgoVec
418
419
    WorldVector() 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), NO_INIT) 
420
    {}
421

Thomas Witkowski's avatar
Thomas Witkowski committed
422
    /// Calls the corresponding constructor of AlgoVec
423
424
    WorldVector(InitType initType, T* ini) 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
425
    {}
426

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    /// Calls the corresponding constructor of AlgoVec
428
429
    WorldVector(InitType initType, const T& ini)
      : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
430
    {}
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    /// Sets all entries to d
433
434
435
436
    inline const WorldVector<T>& operator=(const T& d)
    {
      this->set(d);
      return (*this);
437
    }
438

Thomas Witkowski's avatar
Thomas Witkowski committed
439
    /// Sets the arrays value to the geometric midpoint of the points  p1 and p2.
440
441
442
443
    inline void setMidpoint(const WorldVector<T> &p1, const WorldVector<T> &p2)
    {
      int dow = Global::getGeo(WORLD);

Thomas Witkowski's avatar
Thomas Witkowski committed
444
      for (int i = 0; i < dow; i++)
445
446
447
	this->valArray[i] = 0.5 * (p1[i] + p2[i]);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
448
    /// Multplication of a matrix with a vector.
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
    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:
470
    /// Calls the corresponding constructor of FixVec
471
472
    WorldMatrix()
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
473
    {}
474

475
    /// Calls the corresponding constructor of FixVec
476
477
478
    WorldMatrix(InitType initType, T* ini)
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
    {
479
      TEST_EXIT_DBG(initType == VALUE_LIST)("???\n");
480
      setValues(ini);
481
    }
482
483
484
485
486
487
488
489

    /** \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))
    {
490
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
491
      this->set(ini);
492
    }
493
  
494
    /// Returns true if the matrix is a diagonal matrix, returns false otherwise.
495
496
    bool isDiagMatrix() const;

497
    /// Returns true if the matrix is symmetric, returns false otherwise.
498
499
    bool isSymmetric() const;

500
    /// Sets the diagonal entries to the given value.
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
    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 ====================================================
  // ===========================================================================


521
  /// returns the euclidian distance of a and b
522
523
524
  template<typename T, GeoIndex d>
  double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b);

525
  /// FixVec operator for stream output
526
  template<typename T, GeoIndex d>
527
  std::ostream& operator <<(std::ostream& out, const FixVec<T,d>& fixvec)
528
  {
529
    for (int i = 0; i < fixvec.getSize()-1; i++) {
530
531
532
533
534
535
      out << fixvec[i] << " ";
    }
    out << fixvec[fixvec.getSize()-1];
    return out;
  }

536
  /// creates and inits a VectorOfFixVecs<DimVec<double> >
537
538
  VectorOfFixVecs<DimVec<double> > *createAndInit(int dim, int size, ...);

539
  /// creates and inits and double array
540
541
542
543
544
545
  double *createAndInitArray(int size, ...); 

  inline WorldVector<double> operator*(const WorldVector<double>& v, double d) {
    WorldVector<double> result = v;
    result *= d;
    return result;
546
  }
547
548
549
550
551
552
553

  inline WorldVector<double> operator+(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result += v2;
    return result;
554
  }
555
556
557
558
559
560
561

  inline WorldVector<int> operator+(const WorldVector<int>& v1,
				    const WorldVector<int>& v2) 
  {
    WorldVector<int> result = v1;
    result += v2;
    return result;
562
  }
563
564
565
566
567
568
569

  inline WorldVector<double> operator-(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result -= v2;
    return result;
570
  }
571
572
573
574

  inline bool operator<(const WorldVector<double>& v1,
			const WorldVector<double>& v2) 
  {
575
576
577
578
    int dow = Global::getGeo(WORLD);
    for (int i = 0; i < dow; i++) {
      if (abs(v1[i] - v2[i]) < DBL_TOL) 
	continue;
579
580
581
      return v1[i] < v2[i];
    }
    return false;
582
  }
583
584
585
586

  inline bool operator==(const WorldVector<double>& v1,
			 const WorldVector<double>& v2) 
  {
587
588
589
    int dow = Global::getGeo(WORLD);
    for (int i = 0; i < dow; i++) {
      if (abs(v1[i] - v2[i]) > DBL_TOL) return false;
590
591
    }
    return true;
592
  }
593
594
595
596
597
598
599
600
601
602
603
604
605
606

  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