DOFMatrix.h 16.4 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
259
260
261
262
263
264
265
266
267
    /** \brief
     *  Finishes insertion. For compressed matrix types, this is where the 
     *  compression happens.
     */
    void finishInsertion()
    {
      TEST_EXIT(inserter)("Inserter wasn't used or is already finished.");
      
      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
314
    inline const std::string& getName() const 
    { 
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

385
386
387
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
388
    {
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
      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);
	}
422
    }
423

424
425
426
427
428
429
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
430

431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
  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;
456
      }
457
    }
458

459
    ///
460
461
    int memsize();

462
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
463
    void setIsRankDOF(std::map<DegreeOfFreedom, bool>& dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
464
    {
465
      isRankDOF = dofmap;
466
467
    }
#endif
468
469
470
471
472
473
474
475
476
477
478
479
480
481

  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;

482
    /// Name of the DOFMatrix
483
    std::string name;
484

485
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
486
    base_matrix_type matrix;
487

488
    /// Used while mesh traversal
489
490
491
492
493
494
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
495
    std::vector<Operator*> operators;
496
497
498
499
500
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
501
    std::vector<double*> operatorFactor;
502

503
    ///
504
    std::vector<double*> operatorEstFactor;
505

506
    ///
507
508
    BoundaryManager *boundaryManager;

509
510
511
512
    /** \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.
     */
513
514
    bool coupleMatrix;

515
    /// Temporary variable used in assemble()
516
517
518
519
520
521
522
523
524
525
526
527
528
    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
529

530
531
532
533
534
535
536
    /* \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
537
538
    std::set<int> applyDBCs;

539
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
540
    std::map<DegreeOfFreedom, bool> isRankDOF;
541
542
#endif

543
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
544
    inserter_type *inserter;
545
      
546
547
548
549
550
551
552
553
554
555
    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