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() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
122
      return coupleMatrix; 
123
    }
124

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

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

134
    /// Multiplication with a scalar.
135
136
    void scal(double s);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
146
      return operatorFactor.begin();
147
    }
148

Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
151
      return operatorFactor.end();
152
    }
153

Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
156
      return operatorEstFactor.begin();
157
    }
158

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

Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
166
      return operators.begin();
167
    }
168

Thomas Witkowski's avatar
Thomas Witkowski committed
169
170
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
171
      return operators.end();
172
    }
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198

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

201
202
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
203

204
205
206
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
		  const BoundaryType *bound,
		  Operator *op = NULL);
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
211
212
213
214
215
216
    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
217
218
219
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
220
221
			  ElInfo* rowElInfo,
			  ElInfo* colElInfo,
222
223
			  bool add = true);

Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
226
227
228
229
230
    /* \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();

231
232
233
234
235
236
237
238
239
240
    /** \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; 
241

242
      inserter= new inserter_type(matrix, nnz_per_row);
243
    }
244

245
246
247
248
249
250
251
252
253
254
    /** \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;
255
    }
256
257
258
259
260

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
    virtual bool coarseRestrict() 
    {
263
      return false;
264
    }
265

Thomas Witkowski's avatar
Thomas Witkowski committed
266
    /// Returns const \ref rowFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
267
268
    const FiniteElemSpace* getRowFESpace() const 
    { 
269
      return rowFESpace; 
270
    }
271

272
    /// Returns const \ref colFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
    const FiniteElemSpace* getColFESpace() const 
    { 
275
      return colFESpace; 
276
    }
277

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

284
    /// Returns number of rows (\ref matrix.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
    inline int getSize() const 
    { 
287
      return num_rows(matrix);
288
    }
289

Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
292
293
    /** \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
294
295
    inline int getUsedSize() const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
296
      return rowFESpace->getAdmin()->getUsedSize();
297
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
298

299
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
    inline const std::string& getName() const 
    { 
302
      return name; 
303
    }
304

305
    /// Resizes \ref matrix to n rows
Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
    inline void resize(int n) 
    { 
308
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
309
    }
310

311
    /// Returns value at logical indices a,b
312
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
313

314
    /// Changes col at logical indices a,b to c 
315
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
316
317
318
319
320

    /** \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
     */
321
322
323
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
324

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

327
    /// Prints \ref matrix to stdout
328
329
    void print() const;

330
    /// Removes all matrix entries
331
332
333
334
    void clear()
    {
	set_to_zero(matrix);
    }
335

336
    /// Test whether \ref matrix is symmetric. Exits if not.
337
338
339
340
    void test();

    bool symmetric();

Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
    inline std::vector<Operator*> getOperators() 
    { 
343
      return operators; 
344
    }
345
    
Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
    inline std::vector<double*> getOperatorFactor() 
    { 
348
      return operatorFactor; 
349
    }
350

Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
    inline std::vector<double*> getOperatorEstFactor() 
    { 
353
      return operatorEstFactor; 
354
    }
355

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
    inline BoundaryManager* getBoundaryManager() const 
    { 
358
      return boundaryManager; 
359
    }
360

361
    /// Returns a pointer to \ref applyDBCs.
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
    std::set<int>* getApplyDBCs() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
      return &applyDBCs;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
369
      boundaryManager = bm;
370
    }
371

372
373
374
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
375
    {
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
      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);
	}
409
    }
410

411
412
413
414
415
416
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
417

418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
  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;
443
      }
444
    }
445

446
    ///
447
448
    int memsize();

449
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
450
    /// Sets the petsc application ordering object to map dof indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
    void useApplicationOrdering(AO *ao) 
    {
453
454
455
      applicationOrdering = ao;
    }
#endif
456
457
458
459
460
461
462
463
464
465
466
467
468
469

  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;

470
    /// Name of the DOFMatrix
471
    std::string name;
472

473
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
474
    base_matrix_type matrix;
475

476
    /// Used while mesh traversal
477
478
479
480
481
482
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
483
    std::vector<Operator*> operators;
484
485
486
487
488
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
489
    std::vector<double*> operatorFactor;
490

491
    ///
492
    std::vector<double*> operatorEstFactor;
493

494
    ///
495
496
    BoundaryManager *boundaryManager;

497
498
499
500
    /** \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.
     */
501
502
    bool coupleMatrix;

503
    /// Temporary variable used in assemble()
504
505
506
507
508
509
510
511
512
513
514
515
516
    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
517

518
519
520
521
522
523
524
    /* \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
525
526
    std::set<int> applyDBCs;

527
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
528
529
530
531
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

532
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
533
    inserter_type *inserter;
534
      
535
536
537
538
539
540
541
542
543
544
    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