SystemVector.h 10.9 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
42
		 int size)
      : name(name_),
	feSpace(feSpace_),
Thomas Witkowski's avatar
Thomas Witkowski committed
43
	vectors(size)      
44
    {}
45

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

56
57
    ~SystemVector() 
    {
58
      for (unsigned int i = 0; i < vectors.size(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
59
	delete vectors[i];
60
    }
61

62
63
64
    /// Sets \ref vectors[index] = vec.
    inline void setDOFVector(int index, DOFVector<double> *vec) 
    {
65
66
      TEST_EXIT_DBG(index < static_cast<int>(vectors.size()))
	("Invalid index %d!\n", index);
67
      vectors[index] = vec;
68
    }
69

70
71
72
    /// Returns \ref vectors[index].
    inline DOFVector<double> *getDOFVector(int index) 
    {
73
74
      TEST_EXIT_DBG(index < static_cast<int>(vectors.size()))
	("Invalid index %d!\n", index);
75
      return vectors[index];
76
    }
77

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

86
87
88
89
90
    inline std::string getName() const
    {
      return name;
    }

91
92
93
    /// Returns sum of used vector sizes.
    inline int getUsedSize() const 
    {
94
      int totalSize = 0;
95
      for (unsigned int i = 0; i < vectors.size(); i++)
96
97
	totalSize += vectors[i]->getUsedSize();
      return totalSize;
98
    }
99

100
101
102
    /// Returns number of contained vectors.
    inline int getNumVectors() const 
    { 
103
      return static_cast<int>(vectors.size());
104
    }
105

106
107
    inline int getSize() const 
    {
108
      return static_cast<int>(vectors.size());
109
    }
110

111
    /// Returns the fe space for a given component.
112
    inline const FiniteElemSpace *getFeSpace(int i) const 
113
    { 
114
      return feSpace[i]; 
115
    }
116

117
    /// Returns the fe spaces for all components.
118
    inline std::vector<const FiniteElemSpace*> getFeSpaces() const 
119
    {
120
121
122
      return feSpace;
    }

123
124
125
126
127
128
    vector<const FiniteElemSpace*> getComponentFeSpaces() const;

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

134
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
135
136
137
138
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
139
    }
140

141
142
143
    /// For const access.
    inline double operator[](int index) const 
    {
144
145
146
      int localIndex = index;
      int vectorIndex = 0;

147
      while (localIndex >= vectors[vectorIndex]->getUsedSize()) {
148
149
150
151
	localIndex -= vectors[vectorIndex++]->getUsedSize();
      }

      return (*(vectors[vectorIndex]))[localIndex];
152
    }
153

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

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

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

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

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

183
      return *this;
184
    }
185

186
    void serialize(std::ostream &out);
187

188
    void deserialize(std::istream &in);
189

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

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

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

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

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

223
224
      return result;
    }
225

226
      
227
228

  protected:
229
    /// Name of the system vector
230
    std::string name;
231

232
    /// Finite element space.
233
    std::vector<const FiniteElemSpace*> feSpace;
234

235
    /// Local dof vectors.
236
    std::vector<DOFVector<double>*> vectors;
237
238
  };

239

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

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

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

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

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

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

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

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

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

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

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

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

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

360
    int size = x.getNumVectors();
361
    int i;
362

363
364
365
    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");
366
367

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    int size = x.getNumVectors();
388
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
389

390
    for (i = 0; i < size; i++)
391
      axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
392
  }
393

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

401
    for (int i = 0; i < size; i++)
402
403
404
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
  }

405
406
407
  /// Returns SystemVector::getUsedSize().
  inline int size(SystemVector* vec) 
  {
408
    return vec->getUsedSize();
409
  }
410

411
412
  inline void print(SystemVector* vec) 
  {
413
    vec->print();
414
  }
415
416
417
418

}

#endif // AMDIS_SYSTEMVECTOR_H