SystemVector.h 14.4 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
// ============================================================================
// ==                                                                        ==
// == 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 "MatrixVector.h"
#include "DOFVector.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "OpenMP.h"

namespace AMDiS {

  /** \brief
   * A system vector is a vector of dof vectors used for vector valued
   * problems.
   */
  class SystemVector : public Serializable
  {
  public:
    /** \brief
     * Creator. 
     */
    class Creator 
      : public CreatorInterface<SystemVector> {
    public:
      /** \brief
       * Creators constructor.
       */
49
50
      Creator(const std::string &name_,
	      std::vector<FiniteElemSpace*> feSpace_, 
51
	      int size_) 
52
53
54
	: name(name_), 
	  feSpace(feSpace_), 
	  size(size_) 
55
      {}
56
57
58
59

      /** \brief
       * Destructor
       */
60
      virtual ~Creator() {}
61
62
63
64
65
66

      /** \brief
       * Implementation of CreatorInterface::create().
       */
      SystemVector *create() { 
	char number[10];
67
	std::string numberedName;
Thomas Witkowski's avatar
Thomas Witkowski committed
68
	SystemVector *result = new SystemVector(name, feSpace, size);
69
	for (int i = 0; i < size; i++) {
70
	  sprintf(number, "[%d]", i);
71
	  numberedName = name + std::string(number);
Thomas Witkowski's avatar
Thomas Witkowski committed
72
	  result->setDOFVector(i, new DOFVector<double>(feSpace[i], numberedName));
73
74
	}
	return result; 
75
      }
76
77
78
79
80
81

      /** \brief
       * Implementation of CreatorInterface::free().
       */ 
      void free(SystemVector *vec) { 
	for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
82
	  delete vec->getDOFVector(i);
83
84
	  vec->setDOFVector(i, NULL);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
85
	delete vec;
86
      }
87
88
89
90
91

    private:
      /** \brief
       * Name of the system vector
       */
92
      std::string name;
93
94
95
96

      /** \brief
       * Finite element space used for creation.
       */
97
      std::vector<FiniteElemSpace*> feSpace;
98
99
100
101
102
103
104
105
106
107
108

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

  public:
    /** \brief
     * Constructor.
     */
109
110
    SystemVector(const std::string& name_,
		 std::vector<FiniteElemSpace*> feSpace_, 
111
112
113
114
115
116
117
		 int size)
      : name(name_),
	feSpace(feSpace_),
	vectors(size)
      
    {
      vectors.set(NULL);
118
    }
119
120
121
122
123
124
125
126
127
128

    /** \brief
     * Copy Constructor.
     */
    SystemVector(const SystemVector& rhs)
      : name(rhs.name),
	feSpace(rhs.feSpace),
	vectors(rhs.vectors.getSize())
      
    {
129
      for (int i = 0; i < vectors.getSize(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
130
	vectors[i] = new DOFVector<double>(*rhs.vectors[i]);
131
      }
132
133
    }

134
135
136
    ~SystemVector() 
    {
      for (int i = 0; i < vectors.getSize(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
137
	delete vectors[i];
138
139
      }
    }
140

141
142
143
    void createNewDOFVectors(std::string name)
    {     
      for (int i = 0; i < vectors.getSize(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
144
	vectors[i] = new DOFVector<double>(feSpace[i], "tmp");
145
146
      }
    }
147
148
149
150
151

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

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

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

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

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

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

    /** \brief
196
     * Returns the fe space for a given component.
197
     */
198
199
    inline FiniteElemSpace *getFESpace(int i) const { 
      return feSpace[i]; 
200
    }
201

202
203
204
    /** \brief
     * Returns the fe spaces for all components.
     */
205
206
207
208
    inline std::vector<FiniteElemSpace*> getFESpaces() const {
      return feSpace;
    }

209
210
211
212
213
214
215
216
217
218
    /** \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;

219
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
220
221
222
223
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
224
    }
225
226
227
228
229
230
231
232

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

233
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
234
235
236
237
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
238
    }
239
240
241
242
243
244
245
246
247

    /** \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);
      }
248
    }
249

250
251
252
253
254
255
256
    inline void setCoarsenOperation(CoarsenOperation op) { 
      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++) {
	vectors[i]->setCoarsenOperation(op); 
      }
    }


257
258
259
260
261
262
263
264
265
    /** \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;
266
    }
267
268
269
270
271

    /** \brief
     * Assignement operator.
     */
    inline SystemVector& operator=(const SystemVector& rhs) {
272
      TEST_EXIT_DBG(rhs.vectors.getSize() == vectors.getSize())("invalied sizes\n");
273
274
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
275
276
277
	(*(vectors[i])) = (*(rhs.getDOFVector(i)));
      }
      return *this;
278
    }
279

280
    void serialize(std::ostream &out) {
281
      int size = vectors.getSize();
282
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
283
      for (int i = 0; i < size; i++) {
284
285
	vectors[i]->serialize(out);
      }
286
    }
287

288
    void deserialize(std::istream &in) {
289
      int size, oldSize = vectors.getSize();
290
291
      in.read(reinterpret_cast<char*>(&size), sizeof(int));
      vectors.resize(size);
292
      for (int i = oldSize; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
293
	vectors[i] = new DOFVector<double>(feSpace[i], "");
294
      }
295
      for (int i = 0; i < size; i++) {
296
297
	vectors[i]->deserialize(in);
      }  
298
    }
299
300

    void copy(const SystemVector& rhs) {
301
      int size = vectors.getSize();
302
      TEST_EXIT_DBG(size == rhs.getNumVectors())("invalid sizes\n");
303
      for (int i = 0; i < size; i++) {
304
305
	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
      }
306
    }
307

308
    void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) {
309
310
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
311
312
	vectors[i]->interpol((*f)[i]);
      }
313
314
315
316
317
318
319
    }

    void interpol(SystemVector *v, double factor) {
      for (int i = 0; i < v->getSize(); i++) {
	vectors[i]->interpol(v->getDOFVector(i), factor);
      }
    }
320
321

    void print() {
322
323
      int size = vectors.getSize();
      for (int i = 0; i < size; i++) {
324
325
	vectors[i]->print();
      }
326
    }
327

328
329
330
331
332
333
334
    int calcMemoryUsage() {
      int result = 0;

      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++) {
	result += vectors[i]->calcMemoryUsage();
      }

335
336
      result += sizeof(SystemVector);

337
338
      return result;
    }
339

340
      
341
342
343
344
345

  protected:
    /** \brief
     * Name of the system vector
     */
346
    std::string name;
347
348
349
350

    /** \brief
     * Finite element space.
     */
351
    std::vector<FiniteElemSpace*> feSpace;
352
353
354
355
356
357
358

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

359

360
361
362
363
  /** \brief
   * multiplication with scalar
   */
  inline const SystemVector& operator*=(SystemVector& x, double d) {
364
365
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
366
367
368
      *(x.getDOFVector(i)) *= d;
    }
    return x;
369
  }
370
371
372
373
374

  /** \brief
   * scalar product
   */
  inline double operator*(SystemVector& x, SystemVector& y) {
375
376
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
377
    double result = 0.0;
378
379
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
380
      result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
381
    }
382
    return result;
383
  }
384
385
386
387
388
389
390

  /** \brief
   * addition of two system vectors
   */
  inline const SystemVector& operator+=(SystemVector& x, 
					const SystemVector& y) 
  {
391
392
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
393
394
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
395
396
397
      (*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
    }
    return x;
398
  }
399
400
401
402
403
404
405

  /**
   * subtraction of two system vectors.
   */
  inline const SystemVector& operator-=(SystemVector& x, 
					SystemVector& y) 
  {
406
407
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
408
409
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
410
411
412
      (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
    }
    return x;
413
  }
414
415
416
417
418
419

  /** \brief
   * multiplication with a scalar
   */
  inline SystemVector operator*(SystemVector& x, double d) {
    SystemVector result = x;
420
421
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
422
423
424
      (*(result.getDOFVector(i))) *= d;
    }
    return result;
425
  }
426
427
428
429
430
431

  /** \brief
   * multiplication with a scalar
   */
  inline SystemVector operator*(double d, SystemVector& x) {
    SystemVector result = x;
432
433
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
434
435
436
      (*(result.getDOFVector(i))) *= d;
    }
    return result;
437
  }
438
439
440
441
442
443
444

  /** \brief
   * addition of two system vectors
   */
  inline SystemVector operator+(const SystemVector& x,
				const SystemVector& y)
  {
445
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
446
447
      ("invalid size\n");
    SystemVector result = x;
448
449
    int size = x.getNumVectors();
    for (int i = 0; i < size; i++) {
450
451
452
      (*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
    }
    return result;
453
  }
454
455
456
457
458
459

  /** \brief
   * Calls SystemVector::set(). Used for solving.
   */
  inline void set(SystemVector& x, double value) {
    x.set(value);
460
  } 
461
462
463
464
465
466

  /** \brief
   * Calls SystemVector::set(). Used for solving.
   */
  inline void setValue(SystemVector& x, double value) {
    x.set(value);
467
  }
468
469
470
471
472
473

  /** \brief
   * Norm of system vector.
   */
  inline double norm(SystemVector* x) {
    double result = 0.0;
474
475
    int size = x->getNumVectors();
    for (int i = 0; i < size; i++) {
476
477
478
      result += x->getDOFVector(i)->squareNrm2();
    }
    return sqrt(result);
479
  }
480
481
482
483
484
485

  /** \brief
   * L2 norm of system vector.
   */
  inline double L2Norm(SystemVector* x) {
    double result = 0.0;
486
487
    int size = x->getNumVectors();
    for (int i = 0; i < size; i++) {
488
489
490
      result += x->getDOFVector(i)->L2NormSquare();
    }
    return sqrt(result);
491
  }
492
493
494
495
496
497

  /** \brief
   * H1 norm of system vector.
   */
  inline double H1Norm(SystemVector* x) {
    double result = 0.0;
498
499
    int size = x->getNumVectors();
    for (int i = 0; i < size; i++) {
500
501
502
      result += x->getDOFVector(i)->H1NormSquare();
    }
    return sqrt(result);
503
  }
504
505
506
507
508
509
510

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

513
514
515
    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");
516

Thomas Witkowski's avatar
Thomas Witkowski committed
517
#ifdef _OPENMP
518
#pragma omp parallel for schedule(static, 1) num_threads(min(size, omp_get_max_threads()))
Thomas Witkowski's avatar
Thomas Witkowski committed
519
#endif
520
    for (i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
523
524
      if (!add) 
	result.getDOFVector(i)->set(0.0);

      for (int j = 0; j < size; j++) {
525
526
527
528
529
530
531
532
533
	if (matrix[i][j]) {
	  mv<double>(NoTranspose, 
		     *(matrix[i][j]), 
		     *(x.getDOFVector(j)), 
		     *(result.getDOFVector(i)),
		     true);
	}
      }
    }
534
  }
535
536
537
538
539
540

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

Thomas Witkowski's avatar
Thomas Witkowski committed
544
    int size = x.getNumVectors();
545
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
546

Thomas Witkowski's avatar
Thomas Witkowski committed
547
#ifdef _OPENMP
548
#pragma omp parallel for schedule(static, 1) num_threads(min(size, omp_get_max_threads()))
Thomas Witkowski's avatar
Thomas Witkowski committed
549
#endif
550
    for (i = 0; i < size; i++) {
551
552
      axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
    }
553
  }
554
555
556
557
558
559

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

Thomas Witkowski's avatar
Thomas Witkowski committed
564
    for (int i = 0; i < size; i++) {
565
566
567
568
569
570
571
572
573
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
    }
  }

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

  inline void print(SystemVector* vec) {
    vec->print();
578
  }
579
580
581
582

}

#endif // AMDIS_SYSTEMVECTOR_H