SystemVector.h 13.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
26
27
28
29
30
31
32
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
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file SystemVector.h */

#ifndef AMDIS_SYSTEMVECTOR_H
#define AMDIS_SYSTEMVECTOR_H

#include "MemoryManager.h"
#include "MatrixVector.h"
#include "DOFVector.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "OpenMP.h"

namespace AMDiS {

  // ============================================================================
  // ===== class SystemVector ===================================================
  // ============================================================================

  /** \brief
   * A system vector is a vector of dof vectors used for vector valued
   * problems.
   */
  class SystemVector : public Serializable
  {
  public:
    MEMORY_MANAGED(SystemVector);

    /** \brief
     * Creator. 
     */
    class Creator 
      : public CreatorInterface<SystemVector> {
    public:
      MEMORY_MANAGED(Creator);

      /** \brief
       * Creators constructor.
       */
58
59
      Creator(const std::string &name_,
	      std::vector<FiniteElemSpace*> feSpace_, 
60
	      int size_) 
61
62
63
	: name(name_), 
	  feSpace(feSpace_), 
	  size(size_) 
64
65
66
67
68
69
70
71
72
73
74
75
      {};

      /** \brief
       * Destructor
       */
      virtual ~Creator() {};

      /** \brief
       * Implementation of CreatorInterface::create().
       */
      SystemVector *create() { 
	char number[10];
76
	std::string numberedName;
77
	SystemVector *result = NEW SystemVector(name, feSpace, size);
78
	for (int i = 0; i < size; i++) {
79
	  sprintf(number, "[%d]", i);
80
	  numberedName = name + std::string(number);
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
	  result->setDOFVector(i, NEW DOFVector<double>(feSpace[i], numberedName));
	}
	return result; 
      };

      /** \brief
       * Implementation of CreatorInterface::free().
       */ 
      void free(SystemVector *vec) { 
	for (int i = 0; i < size; i++) {
	  DELETE vec->getDOFVector(i);
	  vec->setDOFVector(i, NULL);
	}
	DELETE vec;
      };

    private:
      /** \brief
       * Name of the system vector
       */
101
      std::string name;
102
103
104
105

      /** \brief
       * Finite element space used for creation.
       */
106
      std::vector<FiniteElemSpace*> feSpace;
107
108
109
110
111
112
113
114
115
116
117

      /** \brief
       * Number of component vectors to be created
       */
      int size;
    };

  public:
    /** \brief
     * Constructor.
     */
118
119
    SystemVector(const std::string& name_,
		 std::vector<FiniteElemSpace*> feSpace_, 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
		 int size)
      : name(name_),
	feSpace(feSpace_),
	vectors(size)
      
    {
      vectors.set(NULL);
    };

    /** \brief
     * Copy Constructor.
     */
    SystemVector(const SystemVector& rhs)
      : name(rhs.name),
	feSpace(rhs.feSpace),
	vectors(rhs.vectors.getSize())
      
    {
138
      for (int i = 0; i < vectors.getSize(); i++) {
139
	vectors[i] = new DOFVector<double>(*rhs.vectors[i]);
140
141
142
      }
    };

143
    virtual ~SystemVector() {};
144
145
146
147
148

    /** \brief
     * Sets \ref vectors[index] = vec.
     */ 
    inline void setDOFVector(int index, DOFVector<double> *vec) {
149
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
150
151
152
153
154
155
      vectors[index] = vec;
    };

    /** \brief
     * Returns \ref vectors[index].
     */
156
157
    inline DOFVector<double> *getDOFVector(int index) {      
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
158
159
160
161
162
163
164
      return vectors[index];
    };

    /** \brief
     * Returns \ref vectors[index].
     */
    inline const DOFVector<double> *getDOFVector(int index) const {
165
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
166
167
168
169
170
171
172
      return vectors[index];
    };

    /** \brief
     * Returns sum of used vector sizes.
     */
    inline int getUsedSize() const {
173
174
      int totalSize = 0;
      int size = vectors.getSize();
175
      for (int i = 0; i < size; i++) {
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
	totalSize += vectors[i]->getUsedSize();
      }
      return totalSize;
    };

    /** \brief
     * Returns number of contained vectors.
     */
    inline int getNumVectors() const { 
      return vectors.getSize(); 
    };

    inline int getSize() const {
      return vectors.getSize();
    };

    /** \brief
     * Returns \ref feSpace.
     */
195
196
197
    inline FiniteElemSpace *getFESpace(int i) const { 
      return feSpace[i]; 
    };
198
199
200
201
202
203
204
205
206
207
208

    /** \brief
     * Here the system vector is interpreted as one large vector. The given
     * is used as a global index which indicates a local vector number and
     * a local index on this vector. The operator returns this local vector
     * at the local index.
     */
    inline double& operator[](int index) {
      int localIndex = index;
      int vectorIndex = 0;

209
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
210
211
212
213
214
215
216
217
218
219
220
221
222
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
    };

    /** \brief
     * For const access.
     */
    inline double operator[](int index) const {
      int localIndex = index;
      int vectorIndex = 0;

223
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
    };

    /** \brief
     * Sets all entries in all vectors to value.
     */
    inline void set(double value) {
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
	vectors[i]->set(value);
      }
    };

    /** \brief
     * Sets all entries in all vectors to value.
     */
    inline SystemVector& operator=(double value) {
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
	(*(vectors[i])) = value;
      }
      return *this;
    };

    /** \brief
     * Assignement operator.
     */
    inline SystemVector& operator=(const SystemVector& rhs) {
255
      TEST_EXIT_DBG(rhs.vectors.getSize() == vectors.getSize())("invalied sizes\n");
256
257
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
258
259
260
261
262
	(*(vectors[i])) = (*(rhs.getDOFVector(i)));
      }
      return *this;
    };

263
    void serialize(std::ostream &out) {
264
      int size = vectors.getSize();
265
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
266
      for (int i = 0; i < size; i++) {
267
268
269
270
	vectors[i]->serialize(out);
      }
    };

271
    void deserialize(std::istream &in) {
272
      int size, oldSize = vectors.getSize();
273
274
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      vectors.resize(size);
275
      for (int i = oldSize; i < size; i++) {
276
277
	vectors[i] = NEW DOFVector<double>(feSpace[i], "");
      }
278
      for (int i = 0; i < size; i++) {
279
280
281
282
283
	vectors[i]->deserialize(in);
      }  
    };

    void copy(const SystemVector& rhs) {
284
      int size = vectors.getSize();
285
      TEST_EXIT_DBG(size == rhs.getNumVectors())("invalid sizes\n");
286
      for (int i = 0; i < size; i++) {
287
288
289
290
	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
      }
    };

291
    void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) {
292
293
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
294
295
296
297
298
	vectors[i]->interpol((*f)[i]);
      }
    };

    void print() {
299
300
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
301
302
303
304
305
306
307
308
309
310
	vectors[i]->print();
      }
    };



  protected:
    /** \brief
     * Name of the system vector
     */
311
    std::string name;
312
313
314
315

    /** \brief
     * Finite element space.
     */
316
    std::vector<FiniteElemSpace*> feSpace;
317
318
319
320
321
322
323
324
325
326
327

    /** \brief
     * Local dof vectors.
     */
    Vector<DOFVector<double>*> vectors;
  };

  /** \brief
   * multiplication with scalar
   */
  inline const SystemVector& operator*=(SystemVector& x, double d) {
328
329
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
330
331
332
333
334
335
336
337
338
      *(x.getDOFVector(i)) *= d;
    }
    return x;
  };

  /** \brief
   * scalar product
   */
  inline double operator*(SystemVector& x, SystemVector& y) {
339
340
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
341
    double result = 0.0;
342
343
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
344
345
346
347
348
349
350
351
352
353
354
      result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
    };
    return result;
  };

  /** \brief
   * addition of two system vectors
   */
  inline const SystemVector& operator+=(SystemVector& x, 
					const SystemVector& y) 
  {
355
356
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
357
358
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
359
360
361
362
363
364
365
366
367
368
369
      (*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
    }
    return x;
  };

  /**
   * subtraction of two system vectors.
   */
  inline const SystemVector& operator-=(SystemVector& x, 
					SystemVector& y) 
  {
370
371
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
372
373
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
374
375
376
377
378
379
380
381
382
383
      (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
    }
    return x;
  };

  /** \brief
   * multiplication with a scalar
   */
  inline SystemVector operator*(SystemVector& x, double d) {
    SystemVector result = x;
384
385
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
386
387
388
389
390
391
392
393
394
395
      (*(result.getDOFVector(i))) *= d;
    }
    return result;
  };

  /** \brief
   * multiplication with a scalar
   */
  inline SystemVector operator*(double d, SystemVector& x) {
    SystemVector result = x;
396
397
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
398
399
400
401
402
403
404
405
406
407
408
      (*(result.getDOFVector(i))) *= d;
    }
    return result;
  };

  /** \brief
   * addition of two system vectors
   */
  inline SystemVector operator+(const SystemVector& x,
				const SystemVector& y)
  {
409
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
410
411
      ("invalid size\n");
    SystemVector result = x;
412
413
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
      (*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
    }
    return result;
  };

  /** \brief
   * Calls SystemVector::set(). Used for solving.
   */
  inline void set(SystemVector& x, double value) {
    x.set(value);
  }; 

  /** \brief
   * Calls SystemVector::set(). Used for solving.
   */
  inline void setValue(SystemVector& x, double value) {
    x.set(value);
  }; 

  /** \brief
   * Norm of system vector.
   */
  inline double norm(SystemVector* x) {
    double result = 0.0;
438
439
    int size = x->getNumVectors();
    for (int i = 0; i < size; i++) {
440
441
442
443
444
445
446
447
448
449
      result += x->getDOFVector(i)->squareNrm2();
    }
    return sqrt(result);
  };

  /** \brief
   * L2 norm of system vector.
   */
  inline double L2Norm(SystemVector* x) {
    double result = 0.0;
450
451
    int size = x->getNumVectors();
    for (int i = 0; i < size; i++) {
452
453
454
455
456
457
458
459
460
461
      result += x->getDOFVector(i)->L2NormSquare();
    }
    return sqrt(result);
  };

  /** \brief
   * H1 norm of system vector.
   */
  inline double H1Norm(SystemVector* x) {
    double result = 0.0;
462
463
    int size = x->getNumVectors();
    for (int i = 0; i < size; i++) {
464
465
466
467
468
469
470
471
472
473
474
      result += x->getDOFVector(i)->H1NormSquare();
    }
    return sqrt(result);
  };

  inline void mv(Matrix<DOFMatrix*> &matrix,
		 const SystemVector &x,
		 SystemVector       &result,
		 bool               add = false)
  {
    int size = x.getNumVectors();
475
    int i;
476

477
478
479
    TEST_EXIT_DBG(size == result.getNumVectors())("incompatible sizes\n");
    TEST_EXIT_DBG(size == matrix.getNumRows())("incompatible sizes\n");
    TEST_EXIT_DBG(size == matrix.getNumCols())("incompatible sizes\n");
480

Thomas Witkowski's avatar
Thomas Witkowski committed
481
#ifdef _OPENMP
482
#pragma omp parallel for schedule(static, 1) num_threads(size)
Thomas Witkowski's avatar
Thomas Witkowski committed
483
#endif
484
    for (i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
      if (!add) 
	result.getDOFVector(i)->set(0.0);

      for (int j = 0; j < size; j++) {
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
	if (matrix[i][j]) {
	  mv<double>(NoTranspose, 
		     *(matrix[i][j]), 
		     *(x.getDOFVector(j)), 
		     *(result.getDOFVector(i)),
		     true);
	}
      }
    }
  };

  /** \brief
   * y = a*x + y;
   */
  inline void axpy(double a, SystemVector& x, SystemVector& y)
  {
505
506
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
507

Thomas Witkowski's avatar
Thomas Witkowski committed
508
    int size = x.getNumVectors();
509
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
#ifdef _OPENMP
512
#pragma omp parallel for schedule(static, 1) num_threads(size)
Thomas Witkowski's avatar
Thomas Witkowski committed
513
#endif
514
    for (i = 0; i < size; i++) {
515
516
517
518
519
520
521
522
523
      axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
    }
  };

  /** \brief
   * y = x + a*y
   */
  inline void xpay(double a, SystemVector& x, SystemVector& y)
  {
524
525
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
526
    int size = x.getNumVectors();
527

Thomas Witkowski's avatar
Thomas Witkowski committed
528
    for (int i = 0; i < size; i++) {
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
    }
  }

  /** \brief
   * Returns SystemVector::getUsedSize().
   */
  inline int size(SystemVector* vec) {
    return vec->getUsedSize();
  };

  inline void print(SystemVector* vec) {
    vec->print();
  };

}

#endif // AMDIS_SYSTEMVECTOR_H