SystemVector.h 12.8 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
27
28
29
30

/** \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
31
#include "Serializer.h"
32
33
34

namespace AMDiS {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

298
299
      return result;
    }
300

301
      
302
303

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

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

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

314

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

}

#endif // AMDIS_SYSTEMVECTOR_H