DOFMatrix.h 14.1 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

/** \file DOFMatrix.h */

#ifndef AMDIS_DOFMATRIX_H
#define AMDIS_DOFMATRIX_H

#include <vector>
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include <set>
28
29
#include <memory>
#include <list>
30
31
#include <boost/numeric/mtl/mtl.hpp>
#include "AMDiS_fwd.h"
32
33
34
35
36
37
38
39
40
41
42
43
44
#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
45
46
   * on DOFVectors. The underlying matrix type 
   * is a CRS matrix of double.
47
   */
48
  class DOFMatrix : public DOFIndexed<bool>,
49
50
51
                    public Serializable
  {
  public:
52
    /// Type of scalars in the underlying matrix
53
    typedef double value_type;
54
55

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

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

  private:
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    /// Symbolic constant for an unused matrix entry
63
64
65
66
67
68
69
70
71
72
73
    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();

74
75
    /// 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
76
	      std::string n = "");
77

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

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

    void copy(const DOFMatrix& rhs);

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

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

    // Only to get rid of the abstract functions, I hope they are not used
Thomas Witkowski's avatar
Thomas Witkowski committed
102
103
104
105
106
107
108
109
110
111
112
    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();
    }
113
114
    
    bool dummy; // Must be deleted later
Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
117
118
119
120
121
122
123
124
125
    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;
    }
126
 
127
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
Thomas Witkowski's avatar
Thomas Witkowski committed
128
129
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) 
    {}
130

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

244
245
246
247
248
249
    /** \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
250
    void startInsertion(int nnz_per_row = 10);
251

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
302
    /** \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
303
304
    inline int getUsedSize() const 
    {
305
      return rowFeSpace->getAdmin()->getUsedSize();
306
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
307

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

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

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

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

    /** \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
     */
330
331
332
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
333

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

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

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

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

    bool symmetric();

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

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

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

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

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

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

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

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

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

404
    ///
405
406
    int memsize();

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

  protected:
    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding row DOFs
     * and basis functions
     */
419
    const FiniteElemSpace *rowFeSpace;
420
421
422
423
424

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

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

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

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

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

451
    ///
452
453
    BoundaryManager *boundaryManager;

454
455
456
457
    /** \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.
     */
458
459
    bool coupleMatrix;

460
    /// Temporary variable used in assemble()
461
462
463
464
465
466
467
468
469
    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
470
    std::vector<DegreeOfFreedom> rowIndices;
471
472

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

475
476
477
478
479
480
481
    /* \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
482
483
    std::set<int> applyDBCs;

Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
486
487
488
489
490
    /* \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;

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

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