FixVec.h 16.7 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
#include "MatrixVector.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "AMDiS_fwd.h"
29
30
31

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
34
35
36
37
38
39
  /// 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 */
  };
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

  /** \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))
    { 
61
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or missing initializer\n");
62
63
64
65
66
67
68
69
70
    }

    /** \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))
    {
71
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
72
      setValues(ini);
73
    }
74
75
76
77
78
79
80
81

    /** \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))
    {
82
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
83
84
85
      this->set(ini);
    }

86
    /// Initialisation for dim.
87
88
    inline void init(int dim) 
    {
89
      this->resize(calcSize(dim));
90
    }
91

92
    /// Initialisation for size
93
94
    inline void initSize(int size) 
    {
95
      this->resize(size);
96
    }
97

98
    /// Returns the \ref size_ of the FixVec.
99
100
    inline int size() const 
    { 
101
      return this->getSize(); 
102
    } 
103
104

  protected:
105
    /// Determines needed vector size.
106
107
108
    static int calcSize(int dim) 
    {
      if (dim < 0)
109
	return Global::getGeo(WORLD);
110
      else
111
	return Global::getGeo(d, dim);
112
    }
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

  public:
    friend class GLWindow;
  };


  /** \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.
     */
137
138
139
    VectorOfFixVecs(int d, int s, InitType initType) 
      : size(s),
	dim(d)
140
    {
141
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
142

143
144
      vec.resize(size);
      for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
145
	vec[i] = new FixVecType(dim, NO_INIT);
146
147
148
149
150
151
152
    }

    /** \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.
     */
153
154
155
    VectorOfFixVecs(int d, int s, InitType initType, const FixVecType* ini)
      : size(s),
	dim(d)
156
    {
157
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
158

159
160
      vec.resize(size);
      for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
161
	vec[i] = new FixVecType(ini[i]);
162
163
164
165
166
167
168
    }

    /** \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.
     */
169
170
171
    VectorOfFixVecs(int d, int s, InitType initType, const FixVecType& ini)
      : size(s),
	dim(d)
172
    {
173
174
175
176
      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
177
	vec[i] = new FixVecType(ini);
178
179
    }

180
    /// Copy constructor
181
182
    VectorOfFixVecs(const VectorOfFixVecs<FixVecType>& rhs)
    {
183
184
      size = rhs.getSize();
      dim = rhs.getDim();
185

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

191
    /// Destructor
192
193
    virtual ~VectorOfFixVecs()
    {
194
      for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
195
	delete vec[i];
196

197
198
199
200
      vec.clear();
    }

    /// Allows the access like in a normal array via []. Used for const objects.
201
202
    inline const FixVecType& operator[](int index) const
    {
203
      TEST_EXIT_DBG(index >= 0 && index < size)("invalid index\n");
204
      return *(vec[index]);
205
    }
206

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

214
    /// Assignment operator
215
216
217
    VectorOfFixVecs<FixVecType>& 
    operator=(const VectorOfFixVecs<FixVecType>& rhs)
    {
218
219
220
221
      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]);
222
223
      }
      return *this;
224
    }
225

226
227
228
229
230
    /// 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
231
	vec[i] = new FixVecType(dim, NO_INIT);
232
233
234
235
      size = newsize;
    }

    /// Returns the \ref size of this VectorOfFixVecs
236
    inline int getSize() const { 
237
238
      return size;
    }
239

240
241
242
243
244
245
    /// Returns \ref dim
    inline int getDim() const {
      return dim;
    }

    /// Returns the size of the contained FixVecs
246
247
248
    inline int getSizeOfFixVec() const { 
      return vec[0]->getSize(); 
    }
249
250

  protected:
251
252
    /// number of contained FixVecs
    int size;
253

254
255
256
257
258
    /// dimension of world
    int dim;

    /// array of pointers to FixVecs
    std::vector<FixVecType*> vec;
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
  };


  /** \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)
    {
279
      TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
      vec = new VectorOfFixVecs<FixVecType>*[rows];
      for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
282
	*i = new VectorOfFixVecs<FixVecType>(dim, columns, NO_INIT);
283
    }
284

285
    /// destructor
286
287
    virtual ~MatrixOfFixVecs()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
288
      for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	delete *i;
Thomas Witkowski's avatar
Thomas Witkowski committed
290
291

      delete [] vec;
292
    }
293

294
    /// assignment operator
295
296
    inline VectorOfFixVecs<FixVecType >& operator[](int index)
    {
297
      TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
298
      return *(vec[index]);
299
    }
300

301
    /// assignment operator
302
303
    inline const VectorOfFixVecs<FixVecType >& operator[](int index) const
    {
304
      TEST_EXIT_DBG(index >= 0 && index < rows)("invalid index\n");
305
      return *(vec[index]);
306
    }
307

308
    /// Returns \ref rows
309
310
    inline int getNumberOfRows() const { 
      return rows; 
311
    }
312

313
    /// Returns \ref columns
314
315
    inline int getNumberOfColumns() const { 
      return columns; 
316
    }
317
318

  private:
319
    /// number of rows of the matrix
320
321
    int rows;

322
    /// number of columns of the matrix
323
324
    int columns;

325
    /// array of pointers to VectorOfFixVecs
326
327
328
329
330
331
332
333
334
335
336
337
338
    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:
339
    DimVec() {}
340

341
342
    /// Calls the corresponding constructor of FixVec
    DimVec(int dim, InitType initType = NO_INIT)
343
      : FixVec<T,PARTS>(dim, initType)
344
    {}
345

346
    /// Calls the corresponding constructor of FixVec
347
348
    DimVec(int dim, InitType initType, T* ini)
      : FixVec<T,PARTS>(dim, initType, ini)
349
    {}
350

351
    /// Calls the corresponding constructor of FixVec
352
353
    DimVec(int dim, InitType initType, const T& ini)
      : FixVec<T,PARTS>(dim, initType, ini)
354
355
    {}

356
    /// Multplication of a matrix with a vector.
357
    void multMatrixVec(DimMat<T> &m, DimVec<T> &v);
358
359
360
361
362
363
364
365
366
367
368
  };


  /** \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
369
    /// Calls the corresponding constructor of VectorOfFixVecs
370
371
    DimMat(int dim, InitType initType = NO_INIT)
      : Matrix<T>(dim + 1, dim + 1)
372
    {}
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
    /// Calls the corresponding constructor of VectorOfFixVecs
375
    DimMat(int dim, InitType initType, const T& ini)
376
      : Matrix<T>(dim + 1, dim + 1)
377
    {
378
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)
379
380
	("wrong initType or wrong initializer\n");    
      this->set(ini);
381
    }
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    /// Calls the corresponding constructor of VectorOfFixVecs
384
    DimMat(int dim, InitType initType, T* ini)
385
      : Matrix<T>(dim + 1, dim + 1)
386
    {
387
      TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
388
      setValues(ini);
389
    }
390
391
392
393
394
395
396
397
398
399
400
401
  };


  /** \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
402
    /// Calls the corresponding constructor of AlgoVec
403
404
    WorldVector() 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), NO_INIT) 
405
    {}
406

Thomas Witkowski's avatar
Thomas Witkowski committed
407
    /// Calls the corresponding constructor of AlgoVec
408
409
    WorldVector(InitType initType, T* ini) 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
410
    {}
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
    /// Calls the corresponding constructor of AlgoVec
413
414
    WorldVector(InitType initType, const T& ini)
      : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
415
    {}
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
    /// Sets all entries to d
418
419
420
421
    inline const WorldVector<T>& operator=(const T& d)
    {
      this->set(d);
      return (*this);
422
    }
423

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

Thomas Witkowski's avatar
Thomas Witkowski committed
429
      for (int i = 0; i < dow; i++)
430
431
432
	this->valArray[i] = 0.5 * (p1[i] + p2[i]);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    /// Multplication of a matrix with a vector.
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
    void multMatrixVec(WorldMatrix<T> &m, WorldVector<T> &v);
  };


  /** \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:
451
    /// Calls the corresponding constructor of FixVec
452
453
    WorldMatrix()
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
454
    {}
455

456
    /// Calls the corresponding constructor of FixVec
457
458
459
    WorldMatrix(InitType initType, T* ini)
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
    {
460
      TEST_EXIT_DBG(initType == VALUE_LIST)("???\n");
461
      setValues(ini);
462
    }
463
464
465
466
467
468
469
470

    /** \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))
    {
471
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
472
      this->set(ini);
473
    }
474
  
475
    /// Returns true if the matrix is a diagonal matrix, returns false otherwise.
476
477
    bool isDiagMatrix() const;

478
    /// Returns true if the matrix is symmetric, returns false otherwise.
479
480
    bool isSymmetric() const;

481
    /// Sets the diagonal entries to the given value.
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
    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);
  };



498
  /// returns the euclidian distance of a and b
499
500
501
  template<typename T, GeoIndex d>
  double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b);

502
  /// FixVec operator for stream output
503
  template<typename T, GeoIndex d>
504
  std::ostream& operator <<(std::ostream& out, const FixVec<T,d>& fixvec)
505
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
506
    for (int i = 0; i < fixvec.getSize() - 1; i++)
507
508
509
510
511
      out << fixvec[i] << " ";
    out << fixvec[fixvec.getSize()-1];
    return out;
  }

512
  /// creates and inits a VectorOfFixVecs<DimVec<double> >
513
514
  VectorOfFixVecs<DimVec<double> > *createAndInit(int dim, int size, ...);

515
  /// creates and inits and double array
516
517
518
519
520
521
  double *createAndInitArray(int size, ...); 

  inline WorldVector<double> operator*(const WorldVector<double>& v, double d) {
    WorldVector<double> result = v;
    result *= d;
    return result;
522
  }
523
524
525
526
527
528
529

  inline WorldVector<double> operator+(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result += v2;
    return result;
530
  }
531
532
533
534
535
536
537

  inline WorldVector<int> operator+(const WorldVector<int>& v1,
				    const WorldVector<int>& v2) 
  {
    WorldVector<int> result = v1;
    result += v2;
    return result;
538
  }
539
540
541
542
543
544
545

  inline WorldVector<double> operator-(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result -= v2;
    return result;
546
  }
547
548
549
550

  inline bool operator<(const WorldVector<double>& v1,
			const WorldVector<double>& v2) 
  {
551
552
553
554
    int dow = Global::getGeo(WORLD);
    for (int i = 0; i < dow; i++) {
      if (abs(v1[i] - v2[i]) < DBL_TOL) 
	continue;
555
556
557
      return v1[i] < v2[i];
    }
    return false;
558
  }
559
560
561
562

  inline bool operator==(const WorldVector<double>& v1,
			 const WorldVector<double>& v2) 
  {
563
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
564
565
566
567
    for (int i = 0; i < dow; i++)
      if (abs(v1[i] - v2[i]) > DBL_TOL) 
	return false;

568
    return true;
569
  }
570
571
572
573
574
575
576
577
578
579
580
581
582
583

  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