DOFMatrix.h 14.8 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut für Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
// ==                                                                        ==
// ============================================================================

/** \file DOFMatrix.h */

#ifndef AMDIS_DOFMATRIX_H
#define AMDIS_DOFMATRIX_H

#include <vector>
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include <set>
27
28
#include <memory>
#include <list>
29
30
#include <boost/numeric/mtl/mtl.hpp>
#include "AMDiS_fwd.h"
31
32
33
34
35
36
37
38
39
40
41
42
43
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
#include "DOFAdmin.h"
#include "DOFIndexed.h"
#include "Boundary.h"
#include "Serializable.h"

namespace AMDiS {

  /** \ingroup DOFAdministration
   * \brief
   * A DOFMatrix is a sparse matrix representation for matrices that work
44
45
   * on DOFVectors. The underlying matrix type 
   * is a CRS matrix of double.
46
   */
47
  class DOFMatrix : public DOFIndexed<bool>,
48
49
50
                    public Serializable
  {
  public:
51
52
53
54
55
56
57
    /// Type of scalars in the underlying matrix
    typedef double                          value_type;

    /// Type of underlying matrix
    typedef mtl::compressed2D<value_type>       base_matrix_type;

    /// Type of inserter for the base matrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
58
    typedef mtl::matrix::inserter<base_matrix_type, mtl::operations::update_plus<value_type> >  inserter_type;
59
60

  private:
Thomas Witkowski's avatar
Thomas Witkowski committed
61
    /// Symbolic constant for an unused matrix entry
62
63
64
65
66
67
68
69
70
71
72
    static const int UNUSED_ENTRY = -1;
    
    /** \brief
     * Symbolic constant for an unused entry which is not followed by any
     * used entry in this row
     */
    static const int NO_MORE_ENTRIES = -2;

  public:    
    DOFMatrix();

Thomas Witkowski's avatar
Thomas Witkowski committed
73
74
    /// Constructs a DOFMatrix with name n and the given row and olumn FESpaces.
    DOFMatrix(const FiniteElemSpace* rowFESpace, const FiniteElemSpace* colFESpace,
Thomas Witkowski's avatar
Thomas Witkowski committed
75
	      std::string n = "");
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
    /// Copy-Constructor
78
79
    DOFMatrix(const DOFMatrix& rhs);

Thomas Witkowski's avatar
Thomas Witkowski committed
80
    /// Destructor
81
82
    virtual ~DOFMatrix();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    /// Assignment operator.
84
85
86
87
    DOFMatrix& operator=(const DOFMatrix& rhs);

    void copy(const DOFMatrix& rhs);

88
89
90
91
    /// Access underlying matrix directly
    base_matrix_type& getBaseMatrix()
    {
	return matrix;
92
    }
93

94
95
96
97
    /// Access underlying matrix directly (const)
    const base_matrix_type& getBaseMatrix() const
    {
	return matrix;
98
    }
99
100

    // Only to get rid of the abstract functions, I hope they are not used
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
103
104
105
106
107
108
109
110
111
    std::vector<bool>::iterator begin() 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; 
      return v.begin();
    }
    
    std::vector<bool>::iterator end() 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; 
      return v.end();
    }
112
113
    
    bool dummy; // Must be deleted later
Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
116
117
118
119
120
121
122
123
124
    bool& operator[](int i) 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); 
      return dummy;
    }

    const bool& operator[](int i) const 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); 
      return dummy;
    }
125
 
126
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) 
    {}
129

130
    /// Implements DOFIndexedBase::freeDOFContent()
131
132
    virtual void freeDOFContent(int index);

133
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
    inline bool isCoupleMatrix() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
136
      return coupleMatrix; 
137
    }
138

139
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    inline void setCoupleMatrix(bool c) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
142
      coupleMatrix = c; 
143
    }
144

145
    /// a * x + y
Thomas Witkowski's avatar
Thomas Witkowski committed
146
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
147

148
    /// Multiplication with a scalar.
149
150
    void scal(double s);

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
154
155
156
    /** \brief
     * Adds an operator to the DOFMatrix. A factor, that is multipled
     * to the operator, and a multilier factor for the estimator may be
     * also given.
     */
    void addOperator(Operator *op, double* factor = NULL, double* estFactor = NULL);
157

Thomas Witkowski's avatar
Thomas Witkowski committed
158
159
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
160
      return operatorFactor.begin();
161
    }
162

Thomas Witkowski's avatar
Thomas Witkowski committed
163
164
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
165
      return operatorFactor.end();
166
    }
167

Thomas Witkowski's avatar
Thomas Witkowski committed
168
169
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
170
      return operatorEstFactor.begin();
171
    }
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
175
      return operatorEstFactor.end();
176
    }
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
180
      return operators.begin();
181
    }
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
185
      return operators.end();
186
    }
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212

    Flag getAssembleFlag();

    /** \brief
     * Updates the matrix matrix by traversing the underlying mesh and assembling
     * the element contributions into the matrix. Information about the 
     * computation of element matrices and connection of local and global DOFs is
     * stored in minfo; the flags for the mesh traversal are stored at 
     * minfo->fill_flags which specifies the elements to be visited and 
     * information that should be present on the elements for the calculation of 
     * the element matrices and boundary information (if minfo->boundBas is not
     * NULL). On the elements, information about the row DOFs is accessed by 
     * minfo->rowBas using info->row_admin; this vector is also used for the 
     * column DOFs if minfo->nCol is less or equal zero, or minfo->col_admin or 
     * minfo->colBas is a NULL pointer; if row and column DOFs are the same, the 
     * boundary type of the DOFs is accessed by minfo->boundBas if 
     * minfo->boundBas is not a NULL pointer; then the element matrix is 
     * computed by minfo->fillElementMatrix(el info, minfo->myFill); these 
     * contributions, multiplied by minfo->factor, are eventually added to matrix
     * by a call of addElementMatrix() with all information about row and column 
     * DOFs, the element matrix, and boundary types, if available;
     * updateMatrix() only adds element contributions; this makes several calls 
     * for the assemblage of one matrix possible; before the first call, the 
     * matrix should be cleared by calling clear dof matrix().
     */
  
213
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound);
214

215
216
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
217

218
219
220
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
221
222
		  const BoundaryType *bound,
		  Operator *op = NULL);
223

Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
226
227
228
229
230
    void assemble2(double factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = NULL);

    /// Adds an element matrix to \ref matrix
Thomas Witkowski's avatar
Thomas Witkowski committed
231
    void addElementMatrix(const ElementMatrix& elMat, 
232
			  const BoundaryType *bound,
233
			  ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
234
			  ElInfo* colElInfo);
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
238
239
240
241
242
    /* \brief
     * That function must be called after the matrix assembling has been finished. 
     * This makes it possible to start some cleanup or matrix data compressing 
     * procedures.
     */
    void finishAssembling();

243
244
245
246
247
248
    /** \brief
     * Enable insertion for assembly. You can optionally give an upper limit for
     * entries per row (per column for CCS matrices).  Choosing this parameter
     * too small can induce perceivable overhead for compressed matrices.  Thus,
     * it's better to choose a bit too large than too small.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
249
    void startInsertion(int nnz_per_row = 10);
250

251
    /** \brief
252
253
     * Finishes insertion. For compressed matrix types, this is where the 
     * compression happens.
254
255
256
     */
    void finishInsertion()
    {
257
258
259
      FUNCNAME("DOFMatrix::finishInsertion()");

      TEST_EXIT(inserter)("Inserter wasn't used or is already finished.\n");
260
261
262
      
      delete inserter;
      inserter= 0;
263
    }
264
265
266
267
268

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
    virtual bool coarseRestrict() 
    {
271
      return false;
272
    }
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    /// Returns const \ref rowFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
    const FiniteElemSpace* getRowFESpace() const 
    { 
277
      return rowFESpace; 
278
    }
279

280
    /// Returns const \ref colFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
281
282
    const FiniteElemSpace* getColFESpace() const 
    { 
283
      return colFESpace; 
284
    }
285

286
    /// Returns const \ref rowFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
    const FiniteElemSpace* getFESpace() const 
    { 
289
      return rowFESpace; 
290
    }
291

292
    /// Returns number of rows (\ref matrix.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
    inline int getSize() const 
    { 
295
      return num_rows(matrix);
296
    }
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301
    /** \brief
     * Returns the number of used rows (equal to number of used DOFs in
     * the row FE space).
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
    inline int getUsedSize() const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
304
      return rowFESpace->getAdmin()->getUsedSize();
305
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
306

307
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
308
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
309
    { 
310
      return name; 
311
    }
312

313
    /// Resizes \ref matrix to n rows
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
    inline void resize(int n) 
    { 
316
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
317
    }
318

319
    /// Returns value at logical indices a,b
320
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
321

322
    /// Changes col at logical indices a,b to c 
323
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
324
325
326
327
328

    /** \brief
     * Creates an entry with logical indices irow, icol if there is no entry
     * yet. Than sign * entry is added to the value at this logical indices
     */
329
330
331
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
332

Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
    void removeRowsWithDBC(std::set<int> *rows);

335
    /// Prints \ref matrix to stdout
336
337
    void print() const;

338
    /// Removes all matrix entries
339
340
341
342
    void clear()
    {
	set_to_zero(matrix);
    }
343

344
    /// Test whether \ref matrix is symmetric. Exits if not.
345
346
347
348
    void test();

    bool symmetric();

Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
    inline std::vector<Operator*> getOperators() 
    { 
351
      return operators; 
352
    }
353
    
Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
    inline std::vector<double*> getOperatorFactor() 
    { 
356
      return operatorFactor; 
357
    }
358

Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
    inline std::vector<double*> getOperatorEstFactor() 
    { 
361
      return operatorEstFactor; 
362
    }
363

Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
    inline BoundaryManager* getBoundaryManager() const 
    { 
366
      return boundaryManager; 
367
    }
368

369
    /// Returns a pointer to \ref applyDBCs.
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
    std::set<int>* getApplyDBCs() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374
      return &applyDBCs;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
377
      boundaryManager = bm;
378
    }
379

Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
382
383
384
385
386
    /// Calculates the average of non zero entries per row in matrix.
    void calculateNnz()
    {
      nnzPerRow = 0;

      if (num_rows(matrix) != 0)
	nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2); 
Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
      if (nnzPerRow < 20) 
	nnzPerRow = 20;
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
391
392
393
394
395
396
    }

    /// Returns \ref nnzPerRow.
    int getNnz()
    {
      return nnzPerRow;
    }

397
398
    /// Writes the matrix to an output stream.
    void serialize(std::ostream &out);
399

400
401
    /// Reads a matrix from an input stream.
    void deserialize(::std::istream &in);
402

403
    ///
404
405
    int memsize();

406
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
407
    void setRankDofs(std::map<DegreeOfFreedom, bool>& dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
408
    {
409
      rankDofs = dofmap;
410
411
    }
#endif
412
413
414
415
416
417
418
419
420
421
422
423
424
425

  protected:
    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding row DOFs
     * and basis functions
     */
    const FiniteElemSpace *rowFESpace;

    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding 
     * column DOFs and basis functions
     */
    const FiniteElemSpace *colFESpace;

426
    /// Name of the DOFMatrix
427
    std::string name;
428

429
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
430
    base_matrix_type matrix;
431

432
    /// Used while mesh traversal
433
434
435
436
437
438
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
439
    std::vector<Operator*> operators;
440
441
442
443
444
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
445
    std::vector<double*> operatorFactor;
446

447
    ///
448
    std::vector<double*> operatorEstFactor;
449

450
    ///
451
452
    BoundaryManager *boundaryManager;

453
454
455
456
    /** \brief
     * If false, the matrix is a diagonal matrix within a matrix of DOF matrices.
     * Otherwise the value is true, and the matrix is an off-diagonal matrix.
     */
457
458
    bool coupleMatrix;

459
    /// Temporary variable used in assemble()
460
461
462
463
464
465
466
467
468
    ElementMatrix elementMatrix;

    /// Number of basis functions in the row fe space.
    int nRow;

    /// Number of basis functions in the col fe space.
    int nCol;

    /// Maps local row indices of an element to global matrix indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
469
    std::vector<DegreeOfFreedom> rowIndices;
470
471

    /// Maps local col indices of an element to global matrix indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
472
    std::vector<DegreeOfFreedom> colIndices;
Thomas Witkowski's avatar
Thomas Witkowski committed
473

474
475
476
477
478
479
480
    /* \brief
     * A set of row indices. When assembling the DOF matrix, all rows, that
     * correspond to a dof at a dirichlet boundary, are ignored and the row is
     * left blank. After assembling, the diagonal element of the matrix must be
     * set to 1. The indices of all rows, where this should be done, are stored
     * in this set.
     */    
Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
    std::set<int> applyDBCs;

Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
485
486
487
488
489
    /* \brief
     * Number of non zero entries per row (average). For instationary problems this
     * information may be used in the next timestep to accelerate insertion of
     * elemnts during assembling.
     */
     int nnzPerRow;

490
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
491
    std::map<DegreeOfFreedom, bool> rankDofs;
492
493
#endif

494
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
495
    inserter_type *inserter;
496
      
497
498
499
500
501
502
503
504
505
506
    friend class DOFAdmin;
    friend class DOFVector<double>;
    friend class DOFVector<unsigned char>;
    friend class DOFVector<int>;
    friend class DOFVector<WorldVector<double> >;
  };

}

#endif  // AMDIS_DOFMATRIX_H