DOFMatrix.h 17 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \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
58
59
60
    /// 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;
    typedef mtl::matrix::inserter<base_matrix_type>  inserter_type;

  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
230
231
232
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
233
234
			  ElInfo* rowElInfo,
			  ElInfo* colElInfo,
235
236
			  bool add = true);

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
250
251
252
253
    /** \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; 
254

255
      inserter= new inserter_type(matrix, nnz_per_row);
256
    }
257

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

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

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

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

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

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

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

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

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

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

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

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

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

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

342
    /// Prints \ref matrix to stdout
343
344
    void print() const;

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

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

    bool symmetric();

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

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

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

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

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

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

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

404
405
406
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
407
    {
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
      out.write(reinterpret_cast<const char*>(&value), sizeof value);
    }

  public:    
    void serialize(::std::ostream &out) 
    {
      using namespace mtl; 

      typedef traits::range_generator<tag::major, base_matrix_type>::type c_type;
      typedef traits::range_generator<tag::nz, c_type>::type              ic_type;

      typedef Collection<base_matrix_type>::size_type                     size_type;
      typedef Collection<base_matrix_type>::value_type                    value_type;
      
      traits::row<base_matrix_type>::type                                 row(matrix); 
      traits::col<base_matrix_type>::type                                 col(matrix);
      traits::const_value<base_matrix_type>::type                         value(matrix); 

      size_type rows= num_rows(matrix), cols= num_cols(matrix), total= matrix.nnz();
      s_write(out, rows);
      s_write(out, cols);
      s_write(out, total);

      for (c_type cursor(mtl::begin<tag::major>(matrix)), cend(mtl::end<tag::major>(matrix)); 
	   cursor != cend; ++cursor)
	for (ic_type icursor(mtl::begin<tag::nz>(cursor)), icend(mtl::end<tag::nz>(cursor)); 
	     icursor != icend; ++icursor) {
	  size_type   my_row= row(*icursor), my_col= col(*icursor);
	  value_type  my_value= value(*icursor);
	  s_write(out, my_row);
	  s_write(out, my_col);
	  s_write(out, my_value);
	}
441
    }
442

443
444
445
446
447
448
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
449

450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
  public:    
    void deserialize(::std::istream &in) 
    {
      using namespace mtl;

      typedef Collection<base_matrix_type>::size_type                     size_type;
      typedef Collection<base_matrix_type>::value_type                    value_type;
      
      size_type rows, cols, total;
      s_read(in, rows);
      s_read(in, cols);
      s_read(in, total);

      // Prepare matrix insertion
      clear();
      // matrix.change_dim(rows, cols) // do we want this?
      inserter_type ins(matrix);

      for (size_type i= 0; i < total; ++i) {
	size_type   my_row, my_col;
	value_type  my_value;
	s_read(in, my_row);
	s_read(in, my_col);
	s_read(in, my_value);
	ins(my_row, my_col) << my_value;
475
      }
476
    }
477

478
    ///
479
480
    int memsize();

481
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
482
    void setRankDofs(std::map<DegreeOfFreedom, bool>& dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    {
484
      rankDofs = dofmap;
485
486
    }
#endif
487
488
489
490
491
492
493
494
495
496
497
498
499
500

  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;

501
    /// Name of the DOFMatrix
502
    std::string name;
503

504
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
505
    base_matrix_type matrix;
506

507
    /// Used while mesh traversal
508
509
510
511
512
513
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
514
    std::vector<Operator*> operators;
515
516
517
518
519
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
520
    std::vector<double*> operatorFactor;
521

522
    ///
523
    std::vector<double*> operatorEstFactor;
524

525
    ///
526
527
    BoundaryManager *boundaryManager;

528
529
530
531
    /** \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.
     */
532
533
    bool coupleMatrix;

534
    /// Temporary variable used in assemble()
535
536
537
538
539
540
541
542
543
544
545
546
547
    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.
    Vector<DegreeOfFreedom> rowIndices;

    /// Maps local col indices of an element to global matrix indices.
    Vector<DegreeOfFreedom> colIndices;
Thomas Witkowski's avatar
Thomas Witkowski committed
548

549
550
551
552
553
554
555
    /* \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
556
557
    std::set<int> applyDBCs;

Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
560
561
562
563
564
    /* \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;

565
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
566
    std::map<DegreeOfFreedom, bool> rankDofs;
567
568
#endif

569
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
570
    inserter_type *inserter;
571
      
572
573
574
575
576
577
578
579
580
581
    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