DOFMatrix.h 16.3 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
44
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
#include "DOFAdmin.h"
#include "DOFIterator.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
53
54
55
56
57
58
59
60
61
    /// 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
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();

Thomas Witkowski's avatar
Thomas Witkowski committed
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
102
103
104
105
106
107

    // 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;}
108
 
109
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
110
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}
111

112
    /// Implements DOFIndexedBase::freeDOFContent()
113
114
    virtual void freeDOFContent(int index);

115
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
    inline bool isCoupleMatrix() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
118
      return coupleMatrix; 
119
    }
120

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

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

130
    /// Multiplication with a scalar.
131
132
    void scal(double s);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
142
      return operatorFactor.begin();
143
    }
144

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

Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
152
      return operatorEstFactor.begin();
153
    }
154

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

Thomas Witkowski's avatar
Thomas Witkowski committed
160
161
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
162
      return operators.begin();
163
    }
164

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

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

197
198
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
199

200
201
202
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
		  const BoundaryType *bound,
		  Operator *op = NULL);
205

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

Thomas Witkowski's avatar
Thomas Witkowski committed
220
221
222
223
224
225
226
    /* \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();

227
228
229
230
231
232
233
234
235
236
    /** \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; 
237

238
      inserter= new inserter_type(matrix, nnz_per_row);
239
    }
240

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

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
257
258
    virtual bool coarseRestrict() 
    {
259
      return false;
260
    }
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    /// Returns const \ref rowFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
    const FiniteElemSpace* getRowFESpace() const 
    { 
265
      return rowFESpace; 
266
    }
267

268
    /// Returns const \ref colFESpace
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
    const FiniteElemSpace* getColFESpace() const 
    { 
271
      return colFESpace; 
272
    }
273

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

280
    /// Returns number of rows (\ref matrix.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
281
282
    inline int getSize() const 
    { 
283
      return num_rows(matrix);
284
    }
285

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

295
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
    inline const std::string& getName() const 
    { 
298
      return name; 
299
    }
300

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

307
    /// Returns value at logical indices a,b
308
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
309

310
    /// Changes col at logical indices a,b to c 
311
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
312
313
314
315
316

    /** \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
     */
317
318
319
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
320

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

323
    /// Prints \ref matrix to stdout
324
325
    void print() const;

326
    /// Removes all matrix entries
327
328
329
330
    void clear()
    {
	set_to_zero(matrix);
    }
331

332
    /// Test whether \ref matrix is symmetric. Exits if not.
333
334
335
336
    void test();

    bool symmetric();

Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
    inline std::vector<Operator*> getOperators() 
    { 
339
      return operators; 
340
    }
341
    
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
    inline std::vector<double*> getOperatorFactor() 
    { 
344
      return operatorFactor; 
345
    }
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
    inline std::vector<double*> getOperatorEstFactor() 
    { 
349
      return operatorEstFactor; 
350
    }
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
    inline BoundaryManager* getBoundaryManager() const 
    { 
354
      return boundaryManager; 
355
    }
356

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

Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
365
      boundaryManager = bm;
366
    }
367

368
369
370
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
371
    {
372
373
374
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
      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);
	}
405
    }
406

407
408
409
410
411
412
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
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
  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;
439
      }
440
    }
441

442
    ///
443
444
    int memsize();

445
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
446
    void setIsRankDOF(std::map<DegreeOfFreedom, bool>& dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
447
    {
448
      isRankDOF = dofmap;
449
450
    }
#endif
451
452
453
454
455
456
457
458
459
460
461
462
463
464

  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;

465
    /// Name of the DOFMatrix
466
    std::string name;
467

468
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
469
    base_matrix_type matrix;
470

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

486
    ///
487
    std::vector<double*> operatorEstFactor;
488

489
    ///
490
491
    BoundaryManager *boundaryManager;

492
493
494
495
    /** \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.
     */
496
497
    bool coupleMatrix;

498
    /// Temporary variable used in assemble()
499
500
501
502
503
504
505
506
507
508
509
510
511
    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
512

513
514
515
516
517
518
519
    /* \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
520
521
    std::set<int> applyDBCs;

522
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
523
    std::map<DegreeOfFreedom, bool> isRankDOF;
524
525
#endif

526
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
527
    inserter_type *inserter;
528
      
529
530
531
532
533
534
535
536
537
538
    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