DOFMatrix.h 16.5 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
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
#include "DOFAdmin.h"
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "Boundary.h"
#include "Serializable.h"

40
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
41
42
43
#include "petscao.h"
#endif

44
45
46
47
48
namespace AMDiS {

  /** \ingroup DOFAdministration
   * \brief
   * A DOFMatrix is a sparse matrix representation for matrices that work
49
50
   * on DOFVectors. The underlying matrix type 
   * is a CRS matrix of double.
51
   */
52
  class DOFMatrix : public DOFIndexed<bool>,
53
54
55
                    public Serializable
  {
  public:
56
57
58
59
60
61
62
63
64
65
    /// 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
66
    /// Symbolic constant for an unused matrix entry
67
68
69
70
71
72
73
74
75
76
77
    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
78
79
    /// 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
80
	      std::string n = "");
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82
    /// Copy-Constructor
83
84
    DOFMatrix(const DOFMatrix& rhs);

Thomas Witkowski's avatar
Thomas Witkowski committed
85
    /// Destructor
86
87
    virtual ~DOFMatrix();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    /// Assignment operator.
89
90
91
92
    DOFMatrix& operator=(const DOFMatrix& rhs);

    void copy(const DOFMatrix& rhs);

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

99
100
101
102
    /// Access underlying matrix directly (const)
    const base_matrix_type& getBaseMatrix() const
    {
	return matrix;
103
    }
104
105
106
107
108
109
110
111

    // Only to get rid of the abstract functions, I hope they are not used
    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();}
    
    bool dummy; // Must be deleted later
    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;}
112
 
113
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
114
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}
115

116
    /// Implements DOFIndexedBase::freeDOFContent()
117
118
    virtual void freeDOFContent(int index);

119
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
    inline bool isCoupleMatrix() { 
      return coupleMatrix; 
122
    }
123

124
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
    inline void setCoupleMatrix(bool c) { 
      coupleMatrix = c; 
127
    }
128

129
    /// a * x + y
Thomas Witkowski's avatar
Thomas Witkowski committed
130
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
131

132
    /// Multiplication with a scalar.
133
134
    void scal(double s);

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
137
138
139
140
    /** \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);
141

142
    inline std::vector<double*>::iterator getOperatorFactorBegin() {
143
      return operatorFactor.begin();
144
    }
145

146
    inline std::vector<double*>::iterator getOperatorFactorEnd() {
147
      return operatorFactor.end();
148
    }
149

150
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
151
      return operatorEstFactor.begin();
152
    }
153

154
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
155
      return operatorEstFactor.end();
156
    }
157

158
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
159
      return operators.begin();
160
    }
161

162
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
163
      return operators.end();
164
    }
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190

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

193
194
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
195

196
197
198
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
		  const BoundaryType *bound,
		  Operator *op = NULL);
201

Thomas Witkowski's avatar
Thomas Witkowski committed
202
203
204
205
206
207
208
    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
209
210
211
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
212
213
			  ElInfo* rowElInfo,
			  ElInfo* colElInfo,
214
215
			  bool add = true);

Thomas Witkowski's avatar
Thomas Witkowski committed
216
217
218
219
220
221
222
    /* \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();

223
224
225
226
227
228
229
230
231
232
    /** \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; 
233

234
      inserter= new inserter_type(matrix, nnz_per_row);
235
    }
236

237
238
239
240
241
242
243
244
245
246
    /** \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;
247
    }
248
249
250
251
252
253
254

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
    virtual bool coarseRestrict() {
      return false;
255
    }
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    /// Returns const \ref rowFESpace
258
259
    const FiniteElemSpace* getRowFESpace() const { 
      return rowFESpace; 
260
    }
261

262
    /// Returns const \ref colFESpace
263
264
    const FiniteElemSpace* getColFESpace() const { 
      return colFESpace; 
265
    }
266

267
    /// Returns const \ref rowFESpace
268
269
    const FiniteElemSpace* getFESpace() const { 
      return rowFESpace; 
270
    }
271

272
    /// Returns number of rows (\ref matrix.size())
273
    inline int getSize() const { 
274
      return num_rows(matrix);
275
    }
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
279
280
281
282
    /** \brief
     * Returns the number of used rows (equal to number of used DOFs in
     * the row FE space).
     */
    inline int getUsedSize() const {
      return rowFESpace->getAdmin()->getUsedSize();
283
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
284

285
    /// Returns \ref name
286
    inline const std::string& getName() const { 
287
      return name; 
288
    }
289

290
    /// Resizes \ref matrix to n rows
291
    inline void resize(int n) { 
292
293
294
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
      //      matrix.change_dim(n, n);
      //WARNING("Obsolete function without effect -- Peter\n");
295
    }
296

297
    /// Returns value at logical indices a,b
298
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
299

300
    /// Changes col at logical indices a,b to c 
301
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
302
303
304
305
306

    /** \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
     */
307
308
309
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
310

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

313
    /// Prints \ref matrix to stdout
314
315
    void print() const;

316
    /// Removes all matrix entries
317
318
319
320
    void clear()
    {
	set_to_zero(matrix);
    }
321

322
    /// Test whether \ref matrix is symmetric. Exits if not.
323
324
325
326
    void test();

    bool symmetric();

327
    inline std::vector<Operator*> getOperators() { 
328
      return operators; 
329
    }
330
    
331
    inline std::vector<double*> getOperatorFactor() { 
332
      return operatorFactor; 
333
    }
334

335
    inline std::vector<double*> getOperatorEstFactor() { 
336
      return operatorEstFactor; 
337
    }
338
339
340

    inline BoundaryManager* getBoundaryManager() const { 
      return boundaryManager; 
341
    }
342

343
    /// Returns a pointer to \ref applyDBCs.
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
347
    std::set<int>* getApplyDBCs() {
      return &applyDBCs;
    }

348
349
    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
350
    }
351
352

    void createPictureFile(const char* filename, int dim);
353
   
354
355
356
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
357
    {
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
      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);
	}
391
    }
392

393
394
395
396
397
398
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
399

400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
  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;
425
      }
426
    }
427

428
    ///
429
430
    int memsize();

431
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
432
433
434
435
436
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif
437
438
439
440
441
442
443
444
445
446
447
448
449
450

  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;

451
    /// Name of the DOFMatrix
452
    std::string name;
453

454
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
455
    base_matrix_type matrix;
456

457
    /// Used while mesh traversal
458
459
460
461
462
463
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
464
    std::vector<Operator*> operators;
465
466
467
468
469
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
470
    std::vector<double*> operatorFactor;
471

472
    ///
473
    std::vector<double*> operatorEstFactor;
474

475
    ///
476
477
    BoundaryManager *boundaryManager;

478
479
480
481
    /** \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.
     */
482
483
    bool coupleMatrix;

484
    /// Temporary variable used in assemble()
485
486
487
488
489
490
491
492
493
494
495
496
497
    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
498

499
500
501
502
503
504
505
    /* \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
506
507
    std::set<int> applyDBCs;

508
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
509
510
511
512
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

513
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
514
    inserter_type *inserter;
515
      
516
517
518
519
520
521
522
523
524
525
    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