DOFMatrix.h 20 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
 
    /** \brief
260
     * DOFMatrix does not need to be compressed before assembling, when using MTL4.
261
     */
262
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}
263
264
265
266
267
268
269
270
271

    /** \brief
     * Implements DOFIndexedBase::freeDOFContent()
     */ 
    virtual void freeDOFContent(int index);

    /** \brief
     * Returns \ref coupleMatrix.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
    inline bool isCoupleMatrix() { 
      return coupleMatrix; 
274
    }
275
276
277
278

    /** \brief
     * Returns \ref coupleMatrix.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
    inline void setCoupleMatrix(bool c) { 
      coupleMatrix = c; 
281
    }
282

283
      /** \brief
284
285
     * a*x + y
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
286
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
287
288
289
290
291
292

    /** \brief
     * Multiplication with a scalar.
     */
    void scal(double s);

Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
295
296
297
298
    /** \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);
299

300
    inline std::vector<double*>::iterator getOperatorFactorBegin() {
301
      return operatorFactor.begin();
302
    }
303

304
    inline std::vector<double*>::iterator getOperatorFactorEnd() {
305
      return operatorFactor.end();
306
    }
307

308
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
309
      return operatorEstFactor.begin();
310
    }
311

312
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
313
      return operatorEstFactor.end();
314
    }
315

316
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
317
      return operators.begin();
318
    }
319

320
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
321
      return operators.end();
322
    }
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348

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

351
352
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
353

354
355
356
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
		  const BoundaryType *bound,
		  Operator *op = NULL);
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362
363
364
365
366
    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
367
368
369
370
371
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
			  bool add = true);

Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374
375
376
377
378
    /* \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();

379
380
381
382
383
384
385
386
387
388
    /** \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; 
389

390
      inserter= new inserter_type(matrix, nnz_per_row);
391
    }
392

393
394
395
396
397
398
399
400
401
402
    /** \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;
403
    }
404
405
406
407
408
409
410

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

Thomas Witkowski's avatar
Thomas Witkowski committed
413
    /// Returns const \ref rowFESpace
414
415
    const FiniteElemSpace* getRowFESpace() const { 
      return rowFESpace; 
416
    }
417

418
    /// Returns const \ref colFESpace
419
420
    const FiniteElemSpace* getColFESpace() const { 
      return colFESpace; 
421
    }
422

423
    /// Returns const \ref rowFESpace
424
425
    const FiniteElemSpace* getFESpace() const { 
      return rowFESpace; 
426
    }
427

428
    /// Returns number of rows (\ref matrix.size())
429
    inline int getSize() const { 
430
	return num_rows(matrix);
431
    }
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
435
436
437
438
    /** \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();
439
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
440

441
442

    // Only fake, shouldn't be called
443
444
445
446
    /** \brief
     * Returns number of cols. For that, the function iteratos over all
     * rows and searchs for the entry with the highest col number.
     */
447
    int getNumCols() const;
448

449
    /// Returns \ref name
450
    inline const std::string& getName() const { 
451
      return name; 
452
    }
453

454
    /// Resizes \ref matrix to n rows
455
    inline void resize(int n) { 
456
457
458
      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");
459
    }
460

461
    /// Returns value at logical indices a,b
462
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
463

464
    /// Changes col at logical indices a,b to c 
465
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
466
467
468
469
470

    /** \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
     */
471
472
473
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
474
475
476
477
478

    void addMatEntry(int row, MatEntry entry);

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

479
    void addRow(std::vector<MatEntry> row);
480

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

483
    /// Prints \ref matrix to stdout
484
485
    void print() const;

486
    /// Removes all matrix entries
487
488
489
490
    void clear()
    {
	set_to_zero(matrix);
    }
491

492
    /// Test whether \ref matrix is symmetric. Exits if not.
493
494
495
496
    void test();

    bool symmetric();

497
    inline std::vector<Operator*> getOperators() { 
498
      return operators; 
499
    }
500
    
501
    inline std::vector<double*> getOperatorFactor() { 
502
      return operatorFactor; 
503
    }
504

505
    inline std::vector<double*> getOperatorEstFactor() { 
506
      return operatorEstFactor; 
507
    }
508
509
510

    inline BoundaryManager* getBoundaryManager() const { 
      return boundaryManager; 
511
    }
512

Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
515
516
    std::set<int>* getApplyDBCs() {
      return &applyDBCs;
    }

517
518
    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
519
    }
520
521

    void createPictureFile(const char* filename, int dim);
522
   
523
524
525
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
526
    {
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
      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);
	}
560
    }
561

562
563
564
565
566
567
  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }
568

569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
  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;
594
      }
595
    }
596

597
    ///
598
599
    int memsize();

600
601
602
603
604
605
#ifdef HAVE_PARALLEL_AMDIS
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif
606
607
608
609
610
611
612
613
614
615
616
617
618
619

  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;

620
    /// Name of the DOFMatrix
621
    std::string name;
622

623
624
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
    base_matrix_type         matrix;
625

626
    /// Used while mesh traversal
627
628
629
630
631
632
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
633
    std::vector<Operator*> operators;
634
635
636
637
638
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
639
    std::vector<double*> operatorFactor;
640

641
    ///
642
    std::vector<double*> operatorEstFactor;
643

644
    ///
645
646
    BoundaryManager *boundaryManager;

647
    ///
648
649
    bool coupleMatrix;

650
    /// Temporary variable used in assemble()
Thomas Witkowski's avatar
Thomas Witkowski committed
651
652
    ElementMatrix *elementMatrix;

653
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
654
655
    std::set<int> applyDBCs;

656
657
658
659
660
#ifdef HAVE_PARALLEL_AMDIS
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

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