SystemVector.h 13.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
// ============================================================================
// ==                                                                        ==
// == 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 {

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

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

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

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

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

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

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

  public:
89
    /// Constructor.
90
91
    SystemVector(const std::string& name_,
		 std::vector<FiniteElemSpace*> feSpace_, 
92
93
94
95
96
97
98
		 int size)
      : name(name_),
	feSpace(feSpace_),
	vectors(size)
      
    {
      vectors.set(NULL);
99
    }
100

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

242
243
    void serialize(std::ostream &out) 
    {
244
      int size = vectors.getSize();
245
      out.write(reinterpret_cast<const char*>(&size), sizeof(int));
246
      for (int i = 0; i < size; i++)
247
	vectors[i]->serialize(out);
248
    }
249

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

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

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

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

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

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

296
297
      return result;
    }
298

299
      
300
301

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

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

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

312

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

}

#endif // AMDIS_SYSTEMVECTOR_H