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<FiniteElemSpace*> feSpace_, 
40
41
42
43
44
45
46
		 int size)
      : name(name_),
	feSpace(feSpace_),
	vectors(size)
      
    {
      vectors.set(NULL);
47
    }
48

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

195
    void serialize(std::ostream &out);
196

197
    void deserialize(std::istream &in);
198

199
200
    void copy(const SystemVector& rhs) 
    {
201
      int size = vectors.getSize();
202
      TEST_EXIT_DBG(size == rhs.getNumVectors())("invalid sizes\n");
203
      for (int i = 0; i < size; i++)
204
	vectors[i]->copy(*(const_cast<SystemVector&>(rhs).getDOFVector(i)));
205
    }
206

207
208
    void interpol(std::vector<AbstractFunction<double, WorldVector<double> >*> *f) 
    {
209
      int size = vectors.getSize();
210
      for (int i = 0; i < size; i++)
211
	vectors[i]->interpol((*f)[i]);
212
213
    }

214
215
216
    void interpol(SystemVector *v, double factor) 
    {
      for (int i = 0; i < v->getSize(); i++)
217
218
	vectors[i]->interpol(v->getDOFVector(i), factor);
    }
219

220
221
    void print() 
    {
222
      int size = vectors.getSize();
223
      for (int i = 0; i < size; i++)
224
	vectors[i]->print();
225
    }
226

227
228
    int calcMemoryUsage() 
    {
229
      int result = 0;
230
      for (int i = 0; i < static_cast<int>(vectors.getSize()); i++)
231
	result += vectors[i]->calcMemoryUsage();
232
233
      result += sizeof(SystemVector);

234
235
      return result;
    }
236

237
      
238
239

  protected:
240
    /// Name of the system vector
241
    std::string name;
242

243
    /// Finite element space.
244
    std::vector<FiniteElemSpace*> feSpace;
245

246
    /// Local dof vectors.
247
248
249
    Vector<DOFVector<double>*> vectors;
  };

250

251
252
253
  /// multiplication with scalar
  inline const SystemVector& operator*=(SystemVector& x, double d) 
  {
254
    int size = x.getNumVectors();
255
    for (int i = 0; i < size; i++)
256
257
      *(x.getDOFVector(i)) *= d;
    return x;
258
  }
259

260
261
262
263
  /// scalar product
  inline double operator*(SystemVector& x, SystemVector& y) 
  {
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
264
    double result = 0.0;
265
    int size = x.getNumVectors();
266
    for (int i = 0; i < size; i++)
267
268
      result += (*x.getDOFVector(i)) * (*y.getDOFVector(i));
    return result;
269
  }
270

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

281
282
  /// subtraction of two system vectors.
  inline const SystemVector& operator-=(SystemVector& x, SystemVector& y) 
283
  {
284
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
285
    int size = x.getNumVectors();
286
    for (int i = 0; i < size; i++)
287
288
      (*(x.getDOFVector(i))) -= (*(y.getDOFVector(i)));
    return x;
289
  }
290

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

301
302
303
  /// multiplication with a scalar
  inline SystemVector operator*(double d, SystemVector& x) 
  {
304
    SystemVector result = x;
305
    int size = x.getNumVectors();
306
    for (int i = 0; i < size; i++)
307
308
      (*(result.getDOFVector(i))) *= d;
    return result;
309
  }
310

311
312
  /// addition of two system vectors
  inline SystemVector operator+(const SystemVector& x, const SystemVector& y)
313
  {
314
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())("invalid size\n");
315
    SystemVector result = x;
316
    int size = x.getNumVectors();
317
    for (int i = 0; i < size; i++)
318
319
      (*(result.getDOFVector(i))) += (*(y.getDOFVector(i)));
    return result;
320
  }
321

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

328
329
330
  /// Calls SystemVector::set(). Used for solving.
  inline void setValue(SystemVector& x, double value) 
  {
331
    x.set(value);
332
  }
333

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

344
345
346
  /// L2 norm of system vector.
  inline double L2Norm(SystemVector* x) 
  {
347
    double result = 0.0;
348
    int size = x->getNumVectors();
349
    for (int i = 0; i < size; i++)
350
351
      result += x->getDOFVector(i)->L2NormSquare();
    return sqrt(result);
352
  }
353

354
355
356
  /// H1 norm of system vector.
  inline double H1Norm(SystemVector* x) 
  {
357
    double result = 0.0;
358
    int size = x->getNumVectors();
359
    for (int i = 0; i < size; i++)
360
361
      result += x->getDOFVector(i)->H1NormSquare();
    return sqrt(result);
362
  }
363
364
365
366
367
368

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

371
    int size = x.getNumVectors();
372
    int i;
373

374
375
376
    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");
377
378

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

382
383
      for (int j = 0; j < size; j++)
	if (matrix[i][j])
384
385
386
387
388
389
	  mv<double>(NoTranspose, 
		     *(matrix[i][j]), 
		     *(x.getDOFVector(j)), 
		     *(result.getDOFVector(i)),
		     true);
    }
390
  }
391

392
  /// y = a*x + y;
393
394
  inline void axpy(double a, SystemVector& x, SystemVector& y)
  {
395
396
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
397

Thomas Witkowski's avatar
Thomas Witkowski committed
398
    int size = x.getNumVectors();
399
    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
400

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

405
  /// y = x + a*y
406
407
  inline void xpay(double a, SystemVector& x, SystemVector& y)
  {
408
409
    TEST_EXIT_DBG(x.getNumVectors() == y.getNumVectors())
          ("invalid size\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
410
    int size = x.getNumVectors();
411

412
    for (int i = 0; i < size; i++)
413
414
415
      xpay(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
  }

416
417
418
  /// Returns SystemVector::getUsedSize().
  inline int size(SystemVector* vec) 
  {
419
    return vec->getUsedSize();
420
  }
421

422
423
  inline void print(SystemVector* vec) 
  {
424
    vec->print();
425
  }
426
427
428
429

}

#endif // AMDIS_SYSTEMVECTOR_H