FixVec.h 18.1 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25
26

/** \file FixVec.h */

#ifndef AMDIS_FIXVEC_H
#define AMDIS_FIXVEC_H

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

namespace AMDiS {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

199
200
201
202
      vec.clear();
    }

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

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

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

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

    /// Returns the \ref size of this VectorOfFixVecs
Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
    inline int getSize() const 
    { 
240
241
      return size;
    }
242

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

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

  protected:
256
257
    /// number of contained FixVecs
    int size;
258

259
260
261
262
263
    /// dimension of world
    int dim;

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

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

      delete [] vec;
298
    }
299

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

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

314
    /// Returns \ref rows
315
316
    inline int getNumberOfRows() const 
    { 
317
      return rows; 
318
    }
319

320
    /// Returns \ref columns
321
322
    inline int getNumberOfColumns() const 
    { 
323
      return columns; 
324
    }
325
326

  private:
327
    /// number of rows of the matrix
328
329
    int rows;

330
    /// number of columns of the matrix
331
332
    int columns;

333
    /// array of pointers to VectorOfFixVecs
334
335
336
337
338
339
340
341
342
343
344
    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>
345
346
  class DimVec : public FixVec<T,PARTS> 
  {
347
  public:
348
    DimVec() {}
349

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

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

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

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


  /** \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
378
    /// Calls the corresponding constructor of VectorOfFixVecs
379
380
    DimMat(int dim, InitType initType = NO_INIT)
      : Matrix<T>(dim + 1, dim + 1)
381
    {}
382

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

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


  /** \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
411
    /// Calls the corresponding constructor of AlgoVec
412
413
    WorldVector() 
      : FixVec<T, WORLD>(Global::getGeo(WORLD), NO_INIT) 
414
    {}
415

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
426
    /// Sets all entries to d
427
428
429
430
    inline const WorldVector<T>& operator=(const T& d)
    {
      this->set(d);
      return (*this);
431
    }
432

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

Thomas Witkowski's avatar
Thomas Witkowski committed
438
      for (int i = 0; i < dow; i++)
439
440
441
	this->valArray[i] = 0.5 * (p1[i] + p2[i]);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
442
    /// Multplication of a matrix with a vector.
443
    void multMatrixVec(const WorldMatrix<T> &m, const WorldVector<T> &v);
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
  };


  /** \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:
460
    /// Calls the corresponding constructor of FixVec
461
462
    WorldMatrix()
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
463
    {}
464

465
    /// Calls the corresponding constructor of FixVec
466
467
468
    WorldMatrix(InitType initType, T* ini)
      : Matrix<T>(Global::getGeo(WORLD), Global::getGeo(WORLD))
    {
469
      TEST_EXIT_DBG(initType == VALUE_LIST)("???\n");
470
      setValues(ini);
471
    }
472
473
474
475
476
477
478
479

    /** \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))
    {
480
      TEST_EXIT_DBG(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
481
      this->set(ini);
482
    }
483
  
484
    /// Returns true if the matrix is a diagonal matrix, returns false otherwise.
485
486
    bool isDiagMatrix() const;

487
    /// Returns true if the matrix is symmetric, returns false otherwise.
488
489
    bool isSymmetric() const;

490
    /// Sets the diagonal entries to the given value.
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
    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);
  };



507
  /// returns the euclidian distance of a and b
508
  template<typename T, GeoIndex d>
509
  double absteukl(const FixVec<T,d>& a, const FixVec<T,d>& b);
510

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

521
  /// creates and inits a VectorOfFixVecs<DimVec<double> >
522
523
  VectorOfFixVecs<DimVec<double> > *createAndInit(int dim, int size, ...);

524
  /// creates and inits and double array
525
526
  double *createAndInitArray(int size, ...); 

527
528
  template<typename T>
  WorldVector<T> operator*(const WorldVector<T>& v, double d)
Thomas Witkowski's avatar
Thomas Witkowski committed
529
  {
530
    WorldVector<T> result = v;
531
532
    result *= d;
    return result;
533
  }
534

535
536
  template<typename T>
  WorldVector<T> operator*(double d, const WorldVector<T>& v)
537
  {
538
    return v * d;
539
  }
540

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

550
551
552
553
554
555
556
//   inline WorldVector<int> operator+(const WorldVector<int>& v1,
// 				    const WorldVector<int>& v2) 
//   {
//     WorldVector<int> result = v1;
//     result += v2;
//     return result;
//   }
557

558
559
560
  template<typename T>
  WorldVector<T> operator-(const WorldVector<T>& v1,
			   const WorldVector<T>& v2)
561
  {
562
    WorldVector<T> result = v1;
563
564
    result -= v2;
    return result;
565
  }
566

567
  inline bool operator<(const WorldVector<double>& v1, const WorldVector<double>& v2) 
568
  {
569
570
571
572
    int dow = Global::getGeo(WORLD);
    for (int i = 0; i < dow; i++) {
      if (abs(v1[i] - v2[i]) < DBL_TOL) 
	continue;
573
574
575
      return v1[i] < v2[i];
    }
    return false;
576
  }
577

578
  inline bool operator==(const WorldVector<double>& v1, const WorldVector<double>& v2) 
579
  {
580
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
581
582
583
584
    for (int i = 0; i < dow; i++)
      if (abs(v1[i] - v2[i]) > DBL_TOL) 
	return false;

585
    return true;
586
  }
587
588
589
590
591
592
593
594
595

  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);
596
597
598
599
600
601
602
603
604


  inline double norm(const WorldVector<double>& v)
  {
    double val = 0.0;
    for (int i = 0; i < Global::getGeo(WORLD); i++)
      val += v[i] * v[i];
    return sqrt(val);
  }
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647


  template<typename T>
  struct GradientType
  {
    typedef WorldVector<T> type;
    
    static mtl::dense_vector<double> getValues(type &t) {
      mtl::dense_vector<double> result(t.getSize());
      for (size_t i = 0; i < t.getSize(); i++)
	result[i] = static_cast<T>(t[i]);
      return result;
    }
  };

  template<typename T>
  struct GradientType<WorldVector<T> >
  {
    typedef WorldVector<WorldVector<T> > type;

    static mtl::dense_vector<double> getValues(type &t) {
      mtl::dense_vector<double> result(sqr(t.getSize()));
      for (size_t i = 0; i < t.getSize(); i++)
	for (size_t j = 0; j < t.getSize(); j++)
	  result[i*(t.getSize()) + j] = static_cast<T>(t[i][j]);
      return result;
    }
  };

  template<typename T>
  struct D2Type
  {
    typedef WorldMatrix<T> type;
    
    static mtl::dense_vector<double> getValues(type &t) {
      mtl::dense_vector<double> result(t.numRows() * t.numCols());
      for (size_t i = 0; i < t.numRows(); i++)
	for (size_t j = 0; j < t.numCols(); j++)
	  result[i*(t.numCols()) + j] = static_cast<T>(t[i][j]);
      return result;
    }
  };

648
649
650
651
652
653
654
655
656
657
658
659
660
  /// defines of wich type the product T1*T2 is
  template<typename T1, typename T2>
  struct ProductType
  {
    typedef T1 type;
  };
  
  // scalar times vector
  template<typename T1>
  struct ProductType<T1, WorldVector<T1> >
  {
    typedef WorldVector<typename ProductType<T1,T1>::type > type;
  };
661
  
662
663
664
665
666
667
668
669
670
671
672
673
674
  // vector times scalar
  template<typename T1>
  struct ProductType<WorldVector<T1>, T1>
  {
    typedef WorldVector<typename ProductType<T1,T1>::type > type;
  };
  
  // vector times vector
  template<typename T1, typename T2>
  struct ProductType<WorldVector<T1>, WorldVector<T2> >
  {
    typedef typename ProductType<T1,T2>::type type;
  };
675
676
}

677

678
679
680
#include "FixVec.hh"

#endif // AMDIS_FIXVEC_H