SystemVector.h 11 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
31
32

/** \file SystemVector.h */

#ifndef AMDIS_SYSTEMVECTOR_H
#define AMDIS_SYSTEMVECTOR_H

#include "MatrixVector.h"
#include "DOFVector.h"
#include "CreatorInterface.h"
#include "Serializable.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
    /// Constructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
38
    SystemVector(std::string name_,
39
		 std::vector<const FiniteElemSpace*> feSpace_, 
40
41
		 int size,
		 bool createVec = false)
42
      : name(name_),
43
	componentSpaces(feSpace_),
Thomas Witkowski's avatar
Thomas Witkowski committed
44
	vectors(size)      
45
46
47
48
49
    {
      if (createVec)
	for (int i = 0; i < size; i++)
	  vectors[i] = new DOFVector<double>(componentSpaces[i], "tmp");
    }
50

51
    /// Copy Constructor.
52
    SystemVector(const SystemVector& rhs)
53
      : name(rhs.getName()),
54
	componentSpaces(rhs.getFeSpaces()),
Thomas Witkowski's avatar
Thomas Witkowski committed
55
	vectors(rhs.getNumVectors())
56
    {
57
      for (unsigned int i = 0; i < vectors.size(); i++)
58
	vectors[i] = new DOFVector<double>(*rhs.getDOFVector(i));
59
60
    }

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

67
68
69
    /// Sets \ref vectors[index] = vec.
    inline void setDOFVector(int index, DOFVector<double> *vec) 
    {
70
71
      TEST_EXIT_DBG(index < static_cast<int>(vectors.size()))
	("Invalid index %d!\n", index);
72
      vectors[index] = vec;
73
    }
74

75
76
77
    /// Returns \ref vectors[index].
    inline DOFVector<double> *getDOFVector(int index) 
    {
78
79
      TEST_EXIT_DBG(index < static_cast<int>(vectors.size()))
	("Invalid index %d!\n", index);
80
      return vectors[index];
81
    }
82

83
84
85
    /// Returns \ref vectors[index].
    inline const DOFVector<double> *getDOFVector(int index) const 
    {
86
87
      TEST_EXIT_DBG(index < static_cast<int>(vectors.size()))
	("Invalid index %d!\n", index);
88
      return vectors[index];
89
    }
90

91
92
93
94
95
    inline std::string getName() const
    {
      return name;
    }

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

105
106
107
    /// Returns number of contained vectors.
    inline int getNumVectors() const 
    { 
108
      return static_cast<int>(vectors.size());
109
    }
110

111
112
    inline int getSize() const 
    {
113
      return static_cast<int>(vectors.size());
114
    }
115

116
    /// Returns the fe space for a given component.
117
    inline const FiniteElemSpace *getFeSpace(int i) const 
118
    { 
119
      return componentSpaces[i]; 
120
    }
121

122
    /// Returns the fe spaces for all components.
123
    inline std::vector<const FiniteElemSpace*> getFeSpaces() const 
124
    {
125
      return componentSpaces;
126
127
    }

128
129
130
131
    /// 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.
132
133
    inline double& operator[](int index) 
    {
134
135
136
      int localIndex = index;
      int vectorIndex = 0;

137
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
138
139
140
141
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
142
    }
143

144
145
146
    /// For const access.
    inline double operator[](int index) const 
    {
147
148
149
      int localIndex = index;
      int vectorIndex = 0;

150
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
151
152
153
154
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
155
    }
156

157
158
159
    /// Sets all entries in all vectors to value.
    inline void set(double value) 
    {
160
      for (unsigned int i = 0; i < vectors.size(); i++)
161
	vectors[i]->set(value);
162
    }
163

164
165
    inline void setCoarsenOperation(CoarsenOperation op) 
    { 
166
      for (int i = 0; i < static_cast<int>(vectors.size()); i++)
167
168
169
	vectors[i]->setCoarsenOperation(op); 
    }

170
171
172
    /// Sets all entries in all vectors to value.
    inline SystemVector& operator=(double value) 
    {
173
      for (int i = 0; i < static_cast<int>(vectors.size()); i++)
174
175
	(*(vectors[i])) = value;
      return *this;
176
    }
177

178
179
180
    /// Assignement operator.
    inline SystemVector& operator=(const SystemVector& rhs) 
    {
181
182
183
      TEST_EXIT_DBG(rhs.vectors.size() == vectors.size())("Invalied sizes!\n");

      for (int i = 0; i < static_cast<int>(vectors.size()); i++)
184
	(*(vectors[i])) = (*(rhs.getDOFVector(i)));
185

186
      return *this;
187
    }
188

189
    void serialize(std::ostream &out);
190

191
    void deserialize(std::istream &in);
192

193
194
    void copy(const SystemVector& rhs) 
    {
195
196
      int size = static_cast<int>(vectors.size());
      TEST_EXIT_DBG(size == rhs.getNumVectors())("Invalid sizes!\n");
197
      for (int i = 0; i < size; i++)
198
	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
199
    }
200

201
202
    void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) 
    {
203
      for (unsigned int i = 0; i < vectors.size(); i++)
204
	vectors[i]->interpol((*f)[i]);
205
206
    }

207
208
209
    void interpol(SystemVector *v, double factor) 
    {
      for (int i = 0; i < v->getSize(); i++)
210
211
	vectors[i]->interpol(v->getDOFVector(i), factor);
    }
212

213
214
    void print() 
    {
215
      for (unsigned int i = 0; i < vectors.size(); i++)
216
	vectors[i]->print();
217
    }
218

219
220
    int calcMemoryUsage() 
    {
221
      int result = 0;
222
      for (unsigned int i = 0; i < vectors.size(); i++)
223
	result += vectors[i]->calcMemoryUsage();
224
225
      result += sizeof(SystemVector);

226
227
      return result;
    }
228

229
      
230
231

  protected:
232
    /// Name of the system vector
233
    std::string name;
234

235
    /// Finite element space.
236
    std::vector<const FiniteElemSpace*> componentSpaces;
237

238
    /// Local dof vectors.
239
    std::vector<DOFVector<double>*> vectors;
240
241
  };

242

243
244
245
  /// multiplication with scalar
  inline const SystemVector& operator*=(SystemVector& x, double d) 
  {
246
    int size = x.getNumVectors();
247
    for (int i = 0; i < size; i++)
248
249
      *(x.getDOFVector(i)) *= d;
    return x;
250
  }
251

252
253
254
255
  /// scalar product
  inline double operator*(SystemVector& x, SystemVector& y) 
  {
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
256
    double result = 0.0;
257
    int size = x.getNumVectors();
258
    for (int i = 0; i < size; i++)
259
260
      result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
    return result;
261
  }
262

263
264
  /// addition of two system vectors
  inline const SystemVector& operator+=(SystemVector& x, const SystemVector& y) 
265
  {
266
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
267
    int size = x.getNumVectors();
268
    for (int i = 0; i < size; i++)
269
270
      (*(x.getDOFVector(i))) += (*(y.getDOFVector(i)));
    return x;
271
  }
272

273
274
  /// subtraction of two system vectors.
  inline const SystemVector& operator-=(SystemVector& x, SystemVector& y) 
275
  {
276
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
277
    int size = x.getNumVectors();
278
    for (int i = 0; i < size; i++)
279
280
      (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
    return x;
281
  }
282

283
284
285
  /// multiplication with a scalar
  inline SystemVector operator*(SystemVector& x, double d) 
  {
286
    SystemVector result = x;
287
    int size = x.getNumVectors();
288
    for (int i = 0; i < size; i++)
289
290
      (*(result.getDOFVector(i))) *= d;
    return result;
291
  }
292

293
294
295
  /// multiplication with a scalar
  inline SystemVector operator*(double d, SystemVector& x) 
  {
296
    SystemVector result = x;
297
    int size = x.getNumVectors();
298
    for (int i = 0; i < size; i++)
299
300
      (*(result.getDOFVector(i))) *= d;
    return result;
301
  }
302

303
304
  /// addition of two system vectors
  inline SystemVector operator+(const SystemVector& x, const SystemVector& y)
305
  {
306
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
307
    SystemVector result = x;
308
    int size = x.getNumVectors();
309
    for (int i = 0; i < size; i++)
310
311
      (*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
    return result;
312
  }
313

314
315
316
  /// Calls SystemVector::set(). Used for solving.
  inline void set(SystemVector& x, double value) 
  {
317
    x.set(value);
318
  } 
319

320
321
322
  /// Calls SystemVector::set(). Used for solving.
  inline void setValue(SystemVector& x, double value) 
  {
323
    x.set(value);
324
  }
325

326
327
328
  /// Norm of system vector.
  inline double norm(SystemVector* x) 
  {
329
    double result = 0.0;
330
    int size = x->getNumVectors();
331
    for (int i = 0; i < size; i++)
332
333
      result += x->getDOFVector(i)->squareNrm2();
    return sqrt(result);
334
  }
335

336
337
338
  /// L2 norm of system vector.
  inline double L2Norm(SystemVector* x) 
  {
339
    double result = 0.0;
340
    int size = x->getNumVectors();
341
    for (int i = 0; i < size; i++)
342
343
      result += x->getDOFVector(i)->L2NormSquare();
    return sqrt(result);
344
  }
345

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

  inline void mv(Matrix<DOFMatrix*> &matrix,
		 const SystemVector &x,
		 SystemVector       &result,
		 bool               add = false)
  {
Naumann, Andreas's avatar
Naumann, Andreas committed
361
    FUNCNAME("mv()");
362

363
    int size = x.getNumVectors();
364
    int i;
365

366
367
368
    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");
369
370

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

374
375
      for (int j = 0; j < size; j++)
	if (matrix[i][j])
376
377
378
379
380
381
	  mv<double>(NoTranspose, 
		     *(matrix[i][j]), 
		     *(x.getDOFVector(j)), 
		     *(result.getDOFVector(i)),
		     true);
    }
382
  }
383

384
  /// y = a*x + y;
385
386
  inline void axpy(double a, SystemVector& x, SystemVector& y)
  {
387
388
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
389

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    int size = x.getNumVectors();
391
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
392

393
    for (i = 0; i < size; i++)
394
      axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
395
  }
396

397
  /// y = x + a*y
398
399
  inline void xpay(double a, SystemVector& x, SystemVector& y)
  {
400
401
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
402
    int size = x.getNumVectors();
403

404
    for (int i = 0; i < size; i++)
405
406
407
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
  }

408
409
410
  /// Returns SystemVector::getUsedSize().
  inline int size(SystemVector* vec) 
  {
411
    return vec->getUsedSize();
412
  }
413

414
415
  inline void print(SystemVector* vec) 
  {
416
    vec->print();
417
  }
418
419
420
421

}

#endif // AMDIS_SYSTEMVECTOR_H