DOFMatrix.h 14.9 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.
127
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}
128

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
152
153
154
155
    /** \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);
156

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
184
      return operators.end();
185
    }
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

    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().
     */
  
212
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound);
213

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
225
226
227
228
229
    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
230
    void addElementMatrix(const ElementMatrix& elMat, 
231
			  const BoundaryType *bound,
232
			  ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
233
			  ElInfo* colElInfo);
234

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
239
240
241
    /* \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();

242
243
244
245
246
247
248
249
250
251
    /** \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.
     */
    void startInsertion(int nnz_per_row= 10)
    {
      if (inserter) 
	delete inserter, inserter= 0; 
252

253
      inserter= new inserter_type(matrix, nnz_per_row);
254
    }
255

256
    /** \brief
257
258
     * Finishes insertion. For compressed matrix types, this is where the 
     * compression happens.
259
260
261
     */
    void finishInsertion()
    {
262
263
264
      FUNCNAME("DOFMatrix::finishInsertion()");

      TEST_EXIT(inserter)("Inserter wasn't used or is already finished.\n");
265
266
267
      
      delete inserter;
      inserter= 0;
268
    }
269
270
271
272
273

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
    virtual bool coarseRestrict() 
    {
276
      return false;
277
    }
278

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

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

291
    /// Returns const \ref rowFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
    const FiniteElemSpace* getFESpace() const 
    { 
294
      return rowFESpace; 
295
    }
296

297
    /// Returns number of rows (\ref matrix.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
    inline int getSize() const 
    { 
300
      return num_rows(matrix);
301
    }
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
306
    /** \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
307
308
    inline int getUsedSize() const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
309
      return rowFESpace->getAdmin()->getUsedSize();
310
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
311

312
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
313
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
314
    { 
315
      return name; 
316
    }
317

318
    /// Resizes \ref matrix to n rows
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
    inline void resize(int n) 
    { 
321
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
322
    }
323

324
    /// Returns value at logical indices a,b
325
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
326

327
    /// Changes col at logical indices a,b to c 
328
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
329
330
331
332
333

    /** \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
     */
334
335
336
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
337

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

340
    /// Prints \ref matrix to stdout
341
342
    void print() const;

343
    /// Removes all matrix entries
344
345
346
347
    void clear()
    {
	set_to_zero(matrix);
    }
348

349
    /// Test whether \ref matrix is symmetric. Exits if not.
350
351
352
353
    void test();

    bool symmetric();

Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
    inline std::vector<Operator*> getOperators() 
    { 
356
      return operators; 
357
    }
358
    
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
    inline std::vector<double*> getOperatorFactor() 
    { 
361
      return operatorFactor; 
362
    }
363

Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
    inline std::vector<double*> getOperatorEstFactor() 
    { 
366
      return operatorEstFactor; 
367
    }
368

Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
    inline BoundaryManager* getBoundaryManager() const 
    { 
371
      return boundaryManager; 
372
    }
373

374
    /// Returns a pointer to \ref applyDBCs.
Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
    std::set<int>* getApplyDBCs() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
      return &applyDBCs;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
382
      boundaryManager = bm;
383
    }
384

Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
387
388
389
390
391
    /// 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
392
393
      if (nnzPerRow < 20) 
	nnzPerRow = 20;
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
396
397
398
399
400
401
    }

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

402
403
    /// Writes the matrix to an output stream.
    void serialize(std::ostream &out);
404

405
406
    /// Reads a matrix from an input stream.
    void deserialize(::std::istream &in);
407

408
    ///
409
410
    int memsize();

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

  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;

431
    /// Name of the DOFMatrix
432
    std::string name;
433

434
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
435
    base_matrix_type matrix;
436

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

452
    ///
453
    std::vector<double*> operatorEstFactor;
454

455
    ///
456
457
    BoundaryManager *boundaryManager;

458
459
460
461
    /** \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.
     */
462
463
    bool coupleMatrix;

464
    /// Temporary variable used in assemble()
465
466
467
468
469
470
471
472
473
    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
474
    std::vector<DegreeOfFreedom> rowIndices;
475
476

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

479
480
481
482
483
484
485
    /* \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
486
487
    std::set<int> applyDBCs;

Thomas Witkowski's avatar
Thomas Witkowski committed
488
489
490
491
492
493
494
    /* \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;

495
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
496
    std::map<DegreeOfFreedom, bool> rankDofs;
497
498
#endif

499
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
500
    inserter_type *inserter;
501
      
502
503
504
505
506
507
508
509
510
511
    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