SystemVector.h 13.4 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
26
27
28
29
// ==                                                                        ==
// ============================================================================

/** \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"
Thomas Witkowski's avatar
Thomas Witkowski committed
30
#include "Serializer.h"
31
32
33

namespace AMDiS {

34
  /// A system vector is a vector of dof vectors used for vector valued problems.
35
36
37
  class SystemVector : public Serializable
  {
  public:
38
    /// Creator. 
39
40
41
    class Creator 
      : public CreatorInterface<SystemVector> {
    public:
42
      /// Creators constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
43
      Creator(std::string name_,
44
	      std::vector<FiniteElemSpace*> feSpace_, 
45
	      int size_) 
46
47
48
	: name(name_), 
	  feSpace(feSpace_), 
	  size(size_) 
49
      {}
50

51
      /// Destructor
52
      virtual ~Creator() {}
53

54
55
56
      /// Implementation of CreatorInterface::create().
      SystemVector *create() 
      { 
57
	char number[10];
58
	std::string numberedName;
Thomas Witkowski's avatar
Thomas Witkowski committed
59
	SystemVector *result = new SystemVector(name, feSpace, size);
60
	for (int i = 0; i < size; i++) {
61
	  sprintf(number, "[%d]", i);
62
	  numberedName = name + std::string(number);
Thomas Witkowski's avatar
Thomas Witkowski committed
63
	  result->setDOFVector(i, new DOFVector<double>(feSpace[i], numberedName));
64
65
	}
	return result; 
66
      }
67

68
69
70
      /// Implementation of CreatorInterface::free().
      void free(SystemVector *vec) 
      { 
71
	for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
72
	  delete vec->getDOFVector(i);
73
74
	  vec->setDOFVector(i, NULL);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
75
	delete vec;
76
      }
77
78

    private:
79
      /// Name of the system vector
80
      std::string name;
81

82
      /// Finite element space used for creation.
83
      std::vector<FiniteElemSpace*> feSpace;
84

85
      /// Number of component vectors to be created
86
87
88
89
      int size;
    };

  public:
90
    /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
91
    SystemVector(std::string name_,
92
		 std::vector<FiniteElemSpace*> feSpace_, 
93
94
95
96
97
98
99
		 int size)
      : name(name_),
	feSpace(feSpace_),
	vectors(size)
      
    {
      vectors.set(NULL);
100
    }
101

102
    /// Copy Constructor.
103
104
105
106
107
108
    SystemVector(const SystemVector& rhs)
      : name(rhs.name),
	feSpace(rhs.feSpace),
	vectors(rhs.vectors.getSize())
      
    {
109
      for (int i = 0; i < vectors.getSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
110
	vectors[i] = new DOFVector<double>(*rhs.vectors[i]);
111
112
    }

113
114
    ~SystemVector() 
    {
115
      for (int i = 0; i < vectors.getSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
116
	delete vectors[i];
117
    }
118

119
120
    void createNewDOFVectors(std::string name)
    {     
121
      for (int i = 0; i < vectors.getSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
122
	vectors[i] = new DOFVector<double>(feSpace[i], "tmp");
123
    }
124

125
126
127
    /// Sets \ref vectors[index] = vec.
    inline void setDOFVector(int index, DOFVector<double> *vec) 
    {
128
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
129
      vectors[index] = vec;
130
    }
131

132
133
134
    /// Returns \ref vectors[index].
    inline DOFVector<double> *getDOFVector(int index) 
    {
135
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
136
      return vectors[index];
137
    }
138

139
140
141
    /// Returns \ref vectors[index].
    inline const DOFVector<double> *getDOFVector(int index) const 
    {
142
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
143
      return vectors[index];
144
    }
145

146
147
148
    /// Returns sum of used vector sizes.
    inline int getUsedSize() const 
    {
149
150
      int totalSize = 0;
      int size = vectors.getSize();
151
      for (int i = 0; i < size; i++)
152
153
	totalSize += vectors[i]->getUsedSize();
      return totalSize;
154
    }
155

156
157
158
    /// Returns number of contained vectors.
    inline int getNumVectors() const 
    { 
159
      return vectors.getSize(); 
160
    }
161

162
163
    inline int getSize() const 
    {
164
      return vectors.getSize();
165
    }
166

167
    /// Returns the fe space for a given component.
168
    inline FiniteElemSpace *getFeSpace(int i) const 
169
    { 
170
      return feSpace[i]; 
171
    }
172

173
    /// Returns the fe spaces for all components.
174
    inline std::vector<FiniteElemSpace*> getFeSpaces() const 
175
    {
176
177
178
      return feSpace;
    }

179
180
181
182
183
184
    /** \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.
     */
185
186
    inline double& operator[](int index) 
    {
187
188
189
      int localIndex = index;
      int vectorIndex = 0;

190
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
191
192
193
194
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

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

197
198
199
    /// For const access.
    inline double operator[](int index) const 
    {
200
201
202
      int localIndex = index;
      int vectorIndex = 0;

203
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
204
205
206
207
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

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

210
211
212
    /// Sets all entries in all vectors to value.
    inline void set(double value) 
    {
213
      int size = vectors.getSize();
214
      for (int i = 0; i < size; i++)
215
	vectors[i]->set(value);
216
    }
217

218
219
220
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++)
221
222
223
	vectors[i]->setCoarsenOperation(op); 
    }

224
225
226
    /// Sets all entries in all vectors to value.
    inline SystemVector& operator=(double value) 
    {
227
      int size = vectors.getSize();
228
      for (int i = 0; i < size; i++)
229
230
	(*(vectors[i])) = value;
      return *this;
231
    }
232

233
234
235
    /// Assignement operator.
    inline SystemVector& operator=(const SystemVector& rhs) 
    {
236
      TEST_EXIT_DBG(rhs.vectors.getSize() == vectors.getSize())("invalied sizes\n");
237
      int size = vectors.getSize();
238
      for (int i = 0; i < size; i++)
239
240
	(*(vectors[i])) = (*(rhs.getDOFVector(i)));
      return *this;
241
    }
242

243
244
    void serialize(std::ostream &out) 
    {
245
      int size = vectors.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
246
      SerUtil::serialize(out, size);
247
      for (int i = 0; i < size; i++)
248
	vectors[i]->serialize(out);
249
    }
250

251
252
    void deserialize(std::istream &in) 
    {
253
      int size, oldSize = vectors.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
254
      SerUtil::deserialize(in, size);
255
      vectors.resize(size);
256
      for (int i = oldSize; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
257
	vectors[i] = new DOFVector<double>(feSpace[i], "");
258
      for (int i = 0; i < size; i++)
259
	vectors[i]->deserialize(in);
260
    }
261

262
263
    void copy(const SystemVector& rhs) 
    {
264
      int size = vectors.getSize();
265
      TEST_EXIT_DBG(size == rhs.getNumVectors())("invalid sizes\n");
266
      for (int i = 0; i < size; i++)
267
	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
268
    }
269

270
271
    void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) 
    {
272
      int size = vectors.getSize();
273
      for (int i = 0; i < size; i++)
274
	vectors[i]->interpol((*f)[i]);
275
276
    }

277
278
279
    void interpol(SystemVector *v, double factor) 
    {
      for (int i = 0; i < v->getSize(); i++)
280
281
	vectors[i]->interpol(v->getDOFVector(i), factor);
    }
282

283
284
    void print() 
    {
285
      int size = vectors.getSize();
286
      for (int i = 0; i < size; i++)
287
	vectors[i]->print();
288
    }
289

290
291
    int calcMemoryUsage() 
    {
292
      int result = 0;
293
      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++)
294
	result += vectors[i]->calcMemoryUsage();
295
296
      result += sizeof(SystemVector);

297
298
      return result;
    }
299

300
      
301
302

  protected:
303
    /// Name of the system vector
304
    std::string name;
305

306
    /// Finite element space.
307
    std::vector<FiniteElemSpace*> feSpace;
308

309
    /// Local dof vectors.
310
311
312
    Vector<DOFVector<double>*> vectors;
  };

313

314
315
316
  /// multiplication with scalar
  inline const SystemVector& operator*=(SystemVector& x, double d) 
  {
317
    int size = x.getNumVectors();
318
    for (int i = 0; i < size; i++)
319
320
      *(x.getDOFVector(i)) *= d;
    return x;
321
  }
322

323
324
325
326
  /// scalar product
  inline double operator*(SystemVector& x, SystemVector& y) 
  {
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
327
    double result = 0.0;
328
    int size = x.getNumVectors();
329
    for (int i = 0; i < size; i++)
330
331
      result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
    return result;
332
  }
333

334
335
  /// addition of two system vectors
  inline const SystemVector& operator+=(SystemVector& x, const SystemVector& y) 
336
  {
337
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
338
    int size = x.getNumVectors();
339
    for (int i = 0; i < size; i++)
340
341
      (*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
    return x;
342
  }
343

344
345
  /// subtraction of two system vectors.
  inline const SystemVector& operator-=(SystemVector& x, SystemVector& y) 
346
  {
347
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
348
    int size = x.getNumVectors();
349
    for (int i = 0; i < size; i++)
350
351
      (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
    return x;
352
  }
353

354
355
356
  /// multiplication with a scalar
  inline SystemVector operator*(SystemVector& x, double d) 
  {
357
    SystemVector result = x;
358
    int size = x.getNumVectors();
359
    for (int i = 0; i < size; i++)
360
361
      (*(result.getDOFVector(i))) *= d;
    return result;
362
  }
363

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

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

385
386
387
  /// Calls SystemVector::set(). Used for solving.
  inline void set(SystemVector& x, double value) 
  {
388
    x.set(value);
389
  } 
390

391
392
393
  /// Calls SystemVector::set(). Used for solving.
  inline void setValue(SystemVector& x, double value) 
  {
394
    x.set(value);
395
  }
396

397
398
399
  /// Norm of system vector.
  inline double norm(SystemVector* x) 
  {
400
    double result = 0.0;
401
    int size = x->getNumVectors();
402
    for (int i = 0; i < size; i++)
403
404
      result += x->getDOFVector(i)->squareNrm2();
    return sqrt(result);
405
  }
406

407
408
409
  /// L2 norm of system vector.
  inline double L2Norm(SystemVector* x) 
  {
410
    double result = 0.0;
411
    int size = x->getNumVectors();
412
    for (int i = 0; i < size; i++)
413
414
      result += x->getDOFVector(i)->L2NormSquare();
    return sqrt(result);
415
  }
416

417
418
419
  /// H1 norm of system vector.
  inline double H1Norm(SystemVector* x) 
  {
420
    double result = 0.0;
421
    int size = x->getNumVectors();
422
    for (int i = 0; i < size; i++)
423
424
      result += x->getDOFVector(i)->H1NormSquare();
    return sqrt(result);
425
  }
426
427
428
429
430
431
432

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

435
436
437
    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");
438

Thomas Witkowski's avatar
Thomas Witkowski committed
439
#ifdef _OPENMP
440
#pragma omp parallel for schedule(static, 1) num_threads(min(size, omp_get_max_threads()))
Thomas Witkowski's avatar
Thomas Witkowski committed
441
#endif
442
    for (i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
443
444
445
      if (!add) 
	result.getDOFVector(i)->set(0.0);

446
447
      for (int j = 0; j < size; j++)
	if (matrix[i][j])
448
449
450
451
452
453
	  mv<double>(NoTranspose, 
		     *(matrix[i][j]), 
		     *(x.getDOFVector(j)), 
		     *(result.getDOFVector(i)),
		     true);
    }
454
  }
455

456
  /// y = a*x + y;
457
458
  inline void axpy(double a, SystemVector& x, SystemVector& y)
  {
459
460
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
461

Thomas Witkowski's avatar
Thomas Witkowski committed
462
    int size = x.getNumVectors();
463
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
464

Thomas Witkowski's avatar
Thomas Witkowski committed
465
#ifdef _OPENMP
466
#pragma omp parallel for schedule(static, 1) num_threads(min(size, omp_get_max_threads()))
Thomas Witkowski's avatar
Thomas Witkowski committed
467
#endif
468
    for (i = 0; i < size; i++)
469
      axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
470
  }
471

472
  /// y = x + a*y
473
474
  inline void xpay(double a, SystemVector& x, SystemVector& y)
  {
475
476
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
477
    int size = x.getNumVectors();
478

479
    for (int i = 0; i < size; i++)
480
481
482
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
  }

483
484
485
  /// Returns SystemVector::getUsedSize().
  inline int size(SystemVector* vec) 
  {
486
    return vec->getUsedSize();
487
  }
488

489
490
  inline void print(SystemVector* vec) 
  {
491
    vec->print();
492
  }
493
494
495
496

}

#endif // AMDIS_SYSTEMVECTOR_H