DOFMatrix.h 20.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
26
27
28
29
// ============================================================================
// ==                                                                        ==
// == 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

// ============================================================================
// ===== includes =============================================================
// ============================================================================

#include <vector>
Thomas Witkowski's avatar
Thomas Witkowski committed
30
#include <set>
31
32
#include <memory>
#include <list>
33
#include <boost/numeric/mtl/mtl.hpp>
34

35
#include "AMDiS_fwd.h"
36
37
38
39
40
41
42
43
44
45
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
#include "DOFAdmin.h"
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "MemoryManager.h"
#include "Boundary.h"
#include "Serializable.h"

46
47
48
49
#ifdef HAVE_PARALLEL_AMDIS
#include "petscao.h"
#endif

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
namespace AMDiS {


  // ===========================================================================
  // ===== struct MatEntry =====================================================
  // ===========================================================================

  /** \brief
   * Represents one entry of a DOFMatrix
   */
  struct MatEntry
  {
    /** \brief
     * column index of entry (if >= 0; else unused)
     */
    DegreeOfFreedom col;

    /** \brief
     * matrix entry
     */
    double entry;
  };


  /** \brief
   * Is used to search for all entries of a matrix which column index is 
   * smaller than a given value.
   */
78
  class MatEntryValueLess : public std::unary_function<MatEntry, bool> {
79
80
81
82
83
  private:
    const double value_;
  public:
    MatEntryValueLess(const double& value) 
      : value_(value) 
84
    {}
85
86
87

    bool operator()(const MatEntry& m) {
      return (fabs(m.entry) < value_);
88
    }      
89
90
91
92
93
94
95
  };


  /** \brief
   * Is used to search for all entries of a matrix which value is smaller 
   * than a given value.
   */
96
  class MatEntryColLess : public std::unary_function<MatEntry, bool> {
97
98
99
100
101
  private:
    const int col_;
  public:
    MatEntryColLess(const int& col) 
      : col_(col) 
102
    {}
103
104
105

    bool operator()(const MatEntry& m) {
      return (fabs(m.col) < col_);
106
    }      
107
108
109
110
111
112
113
  };


  /** \brief
   * This function is required if two MatEntries are compared by their col
   * entry (for example when sorting a vector of MatEntries).
   */
114
  struct CmpMatEntryCol : public std::binary_function<MatEntry, MatEntry, bool> {
115
116
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return m1.col < m2.col;
117
    }
118
119
120
121
122
123
124
  };
  

  /** \brief
   * This function compares two matrix entries by their values. It returns true,
   * if the value of m2 is greater than the value of m1.
   */
125
  struct CmpMatEntryValue : public std::binary_function<MatEntry, MatEntry, bool> {
126
127
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return m1.entry < m2.entry;
128
    }
129
130
131
132
133
134
135
  };


  /** \brief
   * This function compares two matrix entries by their values. It returns true,
   * if the value of m2 is greater than the value of m1.
   */
136
  struct CmpMatEntryAbsValueLess : public std::binary_function<MatEntry, MatEntry, bool> {
137
138
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return fabs(m1.entry) < fabs(m2.entry);
139
    }
140
141
142
143
144
145
  };
  
  /** \brief
   * This function compares two matrix entries by their values. It returns true,
   * if the value of m1 is greater than the value of m2.
   */
146
  struct CmpMatEntryAbsValueGreater : public std::binary_function<MatEntry, MatEntry, bool> {
147
148
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return fabs(m1.entry) > fabs(m2.entry);
149
    }
150
151
152
153
154
155
156
157
158
  };

  // ============================================================================
  // ===== class DOFMatrix ======================================================
  // ============================================================================

  /** \ingroup DOFAdministration
   * \brief
   * A DOFMatrix is a sparse matrix representation for matrices that work
159
160
   * on DOFVectors. The underlying matrix type 
   * is a CRS matrix of double.
161
   */
162
  class DOFMatrix : public DOFIndexed<bool>,
163
164
165
166
167
                    public Serializable
  {
  public:
    MEMORY_MANAGED(DOFMatrix);

168
169
170
171
172
173
174
175
176
177
    /// 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:
178
179
    /** \ingroup DOFAdministration
     * \brief
180
     * Alias for DOFIterator< std::vector<MatEntry<T> > >. So one can access an
181
182
     * iterator working on a DOFMatrix via DOFMatrix::Iterator
     */
183
    class Iterator : public DOFIterator< std::vector<MatEntry> > {
184
    public:
185
      Iterator(DOFIndexed< std::vector<MatEntry> > *c, 
186
	       DOFIteratorType type)
187
	: DOFIterator< std::vector<MatEntry> >(c, type)
188
      {}
189
190
191
192
193
194
    };
    
    /** \brief
     * A MatrixRow represents one row of the sparse matrix. It is realized by
     * a STL vector of MatEntry<T> objects
     */
195
    typedef std::vector<MatEntry> MatrixRow;
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
    
    /** \brief
     * Symbolic constant for an unused matrix entry
     */
    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();

    /** \brief
     * Constructs a DOFMatrix with name n and the given row and 
     * column FESpaces.
     */
    DOFMatrix(const FiniteElemSpace* rowFESpace,
	      const FiniteElemSpace* colFESpace,
Thomas Witkowski's avatar
Thomas Witkowski committed
217
	      std::string n = "");
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236

    /** \brief
     * Copy-Constructor
     */
    DOFMatrix(const DOFMatrix& rhs);


    /** \brief
     * Destructor
     */
    virtual ~DOFMatrix();
  
    /** \brief
     * operator=
     */
    DOFMatrix& operator=(const DOFMatrix& rhs);

    void copy(const DOFMatrix& rhs);

237
238
239
240
    /// Access underlying matrix directly
    base_matrix_type& getBaseMatrix()
    {
	return matrix;
241
    }
242

243
244
245
246
    /// Access underlying matrix directly (const)
    const base_matrix_type& getBaseMatrix() const
    {
	return matrix;
247
    }
248
249
250
251
252
253
254
255
256
257

    // 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();}
#if 1
    
    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;}
#endif
258
 
259
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
260
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}
261

262
    /// Implements DOFIndexedBase::freeDOFContent()
263
264
    virtual void freeDOFContent(int index);

265
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
    inline bool isCoupleMatrix() { 
      return coupleMatrix; 
268
    }
269

270
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
    inline void setCoupleMatrix(bool c) { 
      coupleMatrix = c; 
273
    }
274

275
    /// a * x + y
Thomas Witkowski's avatar
Thomas Witkowski committed
276
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
277

278
    /// Multiplication with a scalar.
279
280
    void scal(double s);

Thomas Witkowski's avatar
Thomas Witkowski committed
281
282
283
284
285
286
    /** \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);
287

288
    inline std::vector<double*>::iterator getOperatorFactorBegin() {
289
      return operatorFactor.begin();
290
    }
291

292
    inline std::vector<double*>::iterator getOperatorFactorEnd() {
293
      return operatorFactor.end();
294
    }
295

296
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
297
      return operatorEstFactor.begin();
298
    }
299

300
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
301
      return operatorEstFactor.end();
302
    }
303

304
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
305
      return operators.begin();
306
    }
307

308
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
309
      return operators.end();
310
    }
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336

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

339
340
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
341

342
343
344
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
		  const BoundaryType *bound,
		  Operator *op = NULL);
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
350
351
352
353
354
    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
355
356
357
358
359
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
			  bool add = true);

Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362
363
364
365
366
    /* \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();

367
368
369
370
371
372
373
374
375
376
    /** \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; 
377

378
      inserter= new inserter_type(matrix, nnz_per_row);
379
    }
380

381
382
383
384
385
386
387
388
389
390
    /** \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;
391
    }
392
393
394
395
396
397
398

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

Thomas Witkowski's avatar
Thomas Witkowski committed
401
    /// Returns const \ref rowFESpace
402
403
    const FiniteElemSpace* getRowFESpace() const { 
      return rowFESpace; 
404
    }
405

406
    /// Returns const \ref colFESpace
407
408
    const FiniteElemSpace* getColFESpace() const { 
      return colFESpace; 
409
    }
410

411
    /// Returns const \ref rowFESpace
412
413
    const FiniteElemSpace* getFESpace() const { 
      return rowFESpace; 
414
    }
415

416
    /// Returns number of rows (\ref matrix.size())
417
    inline int getSize() const { 
418
      return num_rows(matrix);
419
    }
420

Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
423
424
425
426
    /** \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();
427
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
428

429
    // Only fake, shouldn't be called
430
431
432
433
    /** \brief
     * Returns number of cols. For that, the function iteratos over all
     * rows and searchs for the entry with the highest col number.
     */
434
    int getNumCols() const;
435

436
    /// Returns \ref name
437
    inline const std::string& getName() const { 
438
      return name; 
439
    }
440

441
    /// Resizes \ref matrix to n rows
442
    inline void resize(int n) { 
443
444
445
      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");
446
    }
447

448
    /// Returns value at logical indices a,b
449
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
450

451
    /// Changes col at logical indices a,b to c 
452
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
453
454
455
456
457

    /** \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
     */
458
459
460
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
461
462
463
464
465

    void addMatEntry(int row, MatEntry entry);

    void addMatEntry(int row, int DegreeOfFreedom, double value);

466
    void addRow(std::vector<MatEntry> row);
467

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

470
    /// Prints \ref matrix to stdout
471
472
    void print() const;

473
    /// Removes all matrix entries
474
475
476
477
    void clear()
    {
	set_to_zero(matrix);
    }
478

479
    /// Test whether \ref matrix is symmetric. Exits if not.
480
481
482
483
    void test();

    bool symmetric();

484
    inline std::vector<Operator*> getOperators() { 
485
      return operators; 
486
    }
487
    
488
    inline std::vector<double*> getOperatorFactor() { 
489
      return operatorFactor; 
490
    }
491

492
    inline std::vector<double*> getOperatorEstFactor() { 
493
      return operatorEstFactor; 
494
    }
495
496
497

    inline BoundaryManager* getBoundaryManager() const { 
      return boundaryManager; 
498
    }
499

500
    /// Returns a pointer to \ref applyDBCs.
Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
503
504
    std::set<int>* getApplyDBCs() {
      return &applyDBCs;
    }

505
506
    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
507
    }
508
509

    void createPictureFile(const char* filename, int dim);
510
   
511
512
513
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
514
    {
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
      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);
	}
548
    }
549

550
551
552
553
554
555
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
556

557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
  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;
582
      }
583
    }
584

585
    ///
586
587
    int memsize();

588
589
590
591
592
593
#ifdef HAVE_PARALLEL_AMDIS
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif
594
595
596
597
598
599
600
601
602
603
604
605
606
607

  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;

608
    /// Name of the DOFMatrix
609
    std::string name;
610

611
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
612
    base_matrix_type matrix;
613

614
    /// Used while mesh traversal
615
616
617
618
619
620
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
621
    std::vector<Operator*> operators;
622
623
624
625
626
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
627
    std::vector<double*> operatorFactor;
628

629
    ///
630
    std::vector<double*> operatorEstFactor;
631

632
    ///
633
634
    BoundaryManager *boundaryManager;

635
636
637
638
    /** \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.
     */
639
640
    bool coupleMatrix;

641
    /// Temporary variable used in assemble()
Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
    ElementMatrix *elementMatrix;

644
645
646
647
648
649
650
    /* \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
651
652
    std::set<int> applyDBCs;

653
654
655
656
657
#ifdef HAVE_PARALLEL_AMDIS
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

658
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
659
    inserter_type *inserter;
660
      
661
662
663
664
665
666
667
668
669
670
    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