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
Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
    inline int getSize() const 
    { 
238
239
      return size;
    }
240

241
    /// Returns \ref dim
Thomas Witkowski's avatar
Thomas Witkowski committed
242
243
    inline int getDim() const 
    {
244
245
246
247
      return dim;
    }

    /// Returns the size of the contained FixVecs
Thomas Witkowski's avatar
Thomas Witkowski committed
248
249
    inline int getSizeOfFixVec() const 
    { 
250
251
      return vec[0]->getSize(); 
    }
252
253

  protected:
254
255
    /// number of contained FixVecs
    int size;
256

257
258
259
260
261
    /// dimension of world
    int dim;

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


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

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

      delete [] vec;
295
    }
296

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

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

311
    /// Returns \ref rows
312
313
    inline int getNumberOfRows() const { 
      return rows; 
314
    }
315

316
    /// Returns \ref columns
317
318
    inline int getNumberOfColumns() const { 
      return columns; 
319
    }
320
321

  private:
322
    /// number of rows of the matrix
323
324
    int rows;

325
    /// number of columns of the matrix
326
327
    int columns;

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

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

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

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

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


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

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

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


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
436
    /// Multplication of a matrix with a vector.
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
    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:
454
    /// Calls the corresponding constructor of FixVec
455
456
    WorldMatrix()
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
457
    {}
458

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

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

481
    /// Returns true if the matrix is symmetric, returns false otherwise.
482
483
    bool isSymmetric() const;

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



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

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

515
  /// creates and inits a VectorOfFixVecs<DimVec<double> >
516
517
  VectorOfFixVecs<DimVec<double> > *createAndInit(int dim, int size, ...);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
  inline WorldVector<double> operator*(const WorldVector<double>& v, double d) 
  {
523
524
525
    WorldVector<double> result = v;
    result *= d;
    return result;
526
  }
527
528
529
530
531
532
533

  inline WorldVector<double> operator+(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result += v2;
    return result;
534
  }
535
536
537
538
539
540
541

  inline WorldVector<int> operator+(const WorldVector<int>& v1,
				    const WorldVector<int>& v2) 
  {
    WorldVector<int> result = v1;
    result += v2;
    return result;
542
  }
543
544
545
546
547
548
549

  inline WorldVector<double> operator-(const WorldVector<double>& v1,
				       const WorldVector<double>& v2) 
  {
    WorldVector<double> result = v1;
    result -= v2;
    return result;
550
  }
551
552
553
554

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

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

572
    return true;
573
  }
574
575
576
577
578
579
580
581
582
583
584
585
586
587

  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