SystemVector.h 11.5 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

/** \file SystemVector.h */

#ifndef AMDIS_SYSTEMVECTOR_H
#define AMDIS_SYSTEMVECTOR_H

#include "MatrixVector.h"
#include "DOFVector.h"
#include "CreatorInterface.h"
#include "Serializable.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
    /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
39
    SystemVector(std::string name_,
40
		 std::vector<FiniteElemSpace*> feSpace_, 
41
42
43
44
45
46
47
		 int size)
      : name(name_),
	feSpace(feSpace_),
	vectors(size)
      
    {
      vectors.set(NULL);
48
    }
49

50
    /// Copy Constructor.
51
    SystemVector(const SystemVector& rhs)
52
53
54
      : name(rhs.getName()),
	feSpace(rhs.getFeSpaces()),
	vectors(rhs.getNumVectors())
55
56
      
    {
57
      for (int i = 0; i < vectors.getSize(); i++)
58
	vectors[i] = new DOFVector<double>(*rhs.getDOFVector(i));
59
60
    }

61
62
    ~SystemVector() 
    {
63
      for (int i = 0; i < vectors.getSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
64
	delete vectors[i];
65
    }
66

67
68
    void createNewDOFVectors(std::string name)
    {     
69
      for (int i = 0; i < vectors.getSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
70
	vectors[i] = new DOFVector<double>(feSpace[i], "tmp");
71
    }
72

73
74
75
    /// Sets \ref vectors[index] = vec.
    inline void setDOFVector(int index, DOFVector<double> *vec) 
    {
76
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
77
      vectors[index] = vec;
78
    }
79

80
81
82
    /// Returns \ref vectors[index].
    inline DOFVector<double> *getDOFVector(int index) 
    {
83
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
84
      return vectors[index];
85
    }
86

87
88
89
    /// Returns \ref vectors[index].
    inline const DOFVector<double> *getDOFVector(int index) const 
    {
90
      TEST_EXIT_DBG(index < vectors.getSize())("invalid index\n");
91
      return vectors[index];
92
    }
93

94
95
96
97
98
    inline std::string getName() const
    {
      return name;
    }

99
100
101
    /// Returns sum of used vector sizes.
    inline int getUsedSize() const 
    {
102
103
      int totalSize = 0;
      int size = vectors.getSize();
104
      for (int i = 0; i < size; i++)
105
106
	totalSize += vectors[i]->getUsedSize();
      return totalSize;
107
    }
108

109
110
111
    /// Returns number of contained vectors.
    inline int getNumVectors() const 
    { 
112
      return vectors.getSize(); 
113
    }
114

115
116
    inline int getSize() const 
    {
117
      return vectors.getSize();
118
    }
119

120
    /// Returns the fe space for a given component.
121
    inline FiniteElemSpace *getFeSpace(int i) const 
122
    { 
123
      return feSpace[i]; 
124
    }
125

126
    /// Returns the fe spaces for all components.
127
    inline std::vector<FiniteElemSpace*> getFeSpaces() const 
128
    {
129
130
131
      return feSpace;
    }

132
133
134
135
136
137
    /** \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.
     */
138
139
    inline double& operator[](int index) 
    {
140
141
142
      int localIndex = index;
      int vectorIndex = 0;

143
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
144
145
146
147
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

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

150
151
152
    /// For const access.
    inline double operator[](int index) const 
    {
153
154
155
      int localIndex = index;
      int vectorIndex = 0;

156
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
157
158
159
160
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

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

163
164
165
    /// Sets all entries in all vectors to value.
    inline void set(double value) 
    {
166
      int size = vectors.getSize();
167
      for (int i = 0; i < size; i++)
168
	vectors[i]->set(value);
169
    }
170

171
172
173
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++)
174
175
176
	vectors[i]->setCoarsenOperation(op); 
    }

177
178
179
    /// Sets all entries in all vectors to value.
    inline SystemVector& operator=(double value) 
    {
180
      int size = vectors.getSize();
181
      for (int i = 0; i < size; i++)
182
183
	(*(vectors[i])) = value;
      return *this;
184
    }
185

186
187
188
    /// Assignement operator.
    inline SystemVector& operator=(const SystemVector& rhs) 
    {
189
      TEST_EXIT_DBG(rhs.vectors.getSize() == vectors.getSize())("invalied sizes\n");
190
      int size = vectors.getSize();
191
      for (int i = 0; i < size; i++)
192
193
	(*(vectors[i])) = (*(rhs.getDOFVector(i)));
      return *this;
194
    }
195

196
197
    void serialize(std::ostream &out) 
    {
198
      int size = vectors.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
199
      SerUtil::serialize(out, size);
200
      for (int i = 0; i < size; i++)
201
	vectors[i]->serialize(out);
202
    }
203

204
205
    void deserialize(std::istream &in) 
    {
206
      int size, oldSize = vectors.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
207
      SerUtil::deserialize(in, size);
208
      vectors.resize(size);
209
      for (int i = oldSize; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
210
	vectors[i] = new DOFVector<double>(feSpace[i], "");
211
      for (int i = 0; i < size; i++)
212
	vectors[i]->deserialize(in);
213
    }
214

215
216
    void copy(const SystemVector& rhs) 
    {
217
      int size = vectors.getSize();
218
      TEST_EXIT_DBG(size == rhs.getNumVectors())("invalid sizes\n");
219
      for (int i = 0; i < size; i++)
220
	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
221
    }
222

223
224
    void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) 
    {
225
      int size = vectors.getSize();
226
      for (int i = 0; i < size; i++)
227
	vectors[i]->interpol((*f)[i]);
228
229
    }

230
231
232
    void interpol(SystemVector *v, double factor) 
    {
      for (int i = 0; i < v->getSize(); i++)
233
234
	vectors[i]->interpol(v->getDOFVector(i), factor);
    }
235

236
237
    void print() 
    {
238
      int size = vectors.getSize();
239
      for (int i = 0; i < size; i++)
240
	vectors[i]->print();
241
    }
242

243
244
    int calcMemoryUsage() 
    {
245
      int result = 0;
246
      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++)
247
	result += vectors[i]->calcMemoryUsage();
248
249
      result += sizeof(SystemVector);

250
251
      return result;
    }
252

253
      
254
255

  protected:
256
    /// Name of the system vector
257
    std::string name;
258

259
    /// Finite element space.
260
    std::vector<FiniteElemSpace*> feSpace;
261

262
    /// Local dof vectors.
263
264
265
    Vector<DOFVector<double>*> vectors;
  };

266

267
268
269
  /// multiplication with scalar
  inline const SystemVector& operator*=(SystemVector& x, double d) 
  {
270
    int size = x.getNumVectors();
271
    for (int i = 0; i < size; i++)
272
273
      *(x.getDOFVector(i)) *= d;
    return x;
274
  }
275

276
277
278
279
  /// scalar product
  inline double operator*(SystemVector& x, SystemVector& y) 
  {
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
280
    double result = 0.0;
281
    int size = x.getNumVectors();
282
    for (int i = 0; i < size; i++)
283
284
      result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
    return result;
285
  }
286

287
288
  /// addition of two system vectors
  inline const SystemVector& operator+=(SystemVector& x, const SystemVector& y) 
289
  {
290
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
291
    int size = x.getNumVectors();
292
    for (int i = 0; i < size; i++)
293
294
      (*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
    return x;
295
  }
296

297
298
  /// subtraction of two system vectors.
  inline const SystemVector& operator-=(SystemVector& x, SystemVector& y) 
299
  {
300
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
301
    int size = x.getNumVectors();
302
    for (int i = 0; i < size; i++)
303
304
      (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
    return x;
305
  }
306

307
308
309
  /// multiplication with a scalar
  inline SystemVector operator*(SystemVector& x, double d) 
  {
310
    SystemVector result = x;
311
    int size = x.getNumVectors();
312
    for (int i = 0; i < size; i++)
313
314
      (*(result.getDOFVector(i))) *= d;
    return result;
315
  }
316

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

327
328
  /// addition of two system vectors
  inline SystemVector operator+(const SystemVector& x, const SystemVector& y)
329
  {
330
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
331
    SystemVector result = x;
332
    int size = x.getNumVectors();
333
    for (int i = 0; i < size; i++)
334
335
      (*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
    return result;
336
  }
337

338
339
340
  /// Calls SystemVector::set(). Used for solving.
  inline void set(SystemVector& x, double value) 
  {
341
    x.set(value);
342
  } 
343

344
345
346
  /// Calls SystemVector::set(). Used for solving.
  inline void setValue(SystemVector& x, double value) 
  {
347
    x.set(value);
348
  }
349

350
351
352
  /// Norm of system vector.
  inline double norm(SystemVector* x) 
  {
353
    double result = 0.0;
354
    int size = x->getNumVectors();
355
    for (int i = 0; i < size; i++)
356
357
      result += x->getDOFVector(i)->squareNrm2();
    return sqrt(result);
358
  }
359

360
361
362
  /// L2 norm of system vector.
  inline double L2Norm(SystemVector* x) 
  {
363
    double result = 0.0;
364
    int size = x->getNumVectors();
365
    for (int i = 0; i < size; i++)
366
367
      result += x->getDOFVector(i)->L2NormSquare();
    return sqrt(result);
368
  }
369

370
371
372
  /// H1 norm of system vector.
  inline double H1Norm(SystemVector* x) 
  {
373
    double result = 0.0;
374
    int size = x->getNumVectors();
375
    for (int i = 0; i < size; i++)
376
377
      result += x->getDOFVector(i)->H1NormSquare();
    return sqrt(result);
378
  }
379
380
381
382
383
384

  inline void mv(Matrix<DOFMatrix*> &matrix,
		 const SystemVector &x,
		 SystemVector       &result,
		 bool               add = false)
  {
Naumann, Andreas's avatar
Naumann, Andreas committed
385
386
387
    FUNCNAME("mv()");
    TEST_EXIT(false)("This function is not supported any more.\n");
#if 0
388
    int size = x.getNumVectors();
389
    int i;
390

391
392
393
    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");
394
395

    for (i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
398
      if (!add) 
	result.getDOFVector(i)->set(0.0);

399
400
      for (int j = 0; j < size; j++)
	if (matrix[i][j])
401
402
403
404
405
406
	  mv<double>(NoTranspose, 
		     *(matrix[i][j]), 
		     *(x.getDOFVector(j)), 
		     *(result.getDOFVector(i)),
		     true);
    }
Naumann, Andreas's avatar
Naumann, Andreas committed
407
#endif
408
  }
409

410
  /// y = a*x + y;
411
412
  inline void axpy(double a, SystemVector& x, SystemVector& y)
  {
413
414
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
415

Thomas Witkowski's avatar
Thomas Witkowski committed
416
    int size = x.getNumVectors();
417
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
418

419
    for (i = 0; i < size; i++)
420
      axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
421
  }
422

423
  /// y = x + a*y
424
425
  inline void xpay(double a, SystemVector& x, SystemVector& y)
  {
426
427
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
428
    int size = x.getNumVectors();
429

430
    for (int i = 0; i < size; i++)
431
432
433
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
  }

434
435
436
  /// Returns SystemVector::getUsedSize().
  inline int size(SystemVector* vec) 
  {
437
    return vec->getUsedSize();
438
  }
439

440
441
  inline void print(SystemVector* vec) 
  {
442
    vec->print();
443
  }
444
445
446
447

}

#endif // AMDIS_SYSTEMVECTOR_H