DOFMatrix.h 20.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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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
78
79
80
81
82
83
// ============================================================================
// ==                                                                        ==
// == 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>
#include <memory>
#include <list>
#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"

namespace AMDiS {

  // ============================================================================
  // ===== forward declarations =================================================
  // ============================================================================

  class Mesh;
  class FiniteElemSpace;
  class ElInfo;
  class BasisFunction;
  class FillInfo;
  class Operator;
  class ElementMatrix;
  class BoundaryManager;

  template<typename T> class DOFVector;

  // ===========================================================================
  // ===== 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.
   */
84
  class MatEntryValueLess : public std::unary_function<MatEntry, bool> {
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
  private:
    const double value_;
  public:
    MatEntryValueLess(const double& value) 
      : value_(value) 
    {};

    bool operator()(const MatEntry& m) {
      return (fabs(m.entry) < value_);
    };      
  };


  /** \brief
   * Is used to search for all entries of a matrix which value is smaller 
   * than a given value.
   */
102
  class MatEntryColLess : public std::unary_function<MatEntry, bool> {
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
  private:
    const int col_;
  public:
    MatEntryColLess(const int& col) 
      : col_(col) 
    {};

    bool operator()(const MatEntry& m) {
      return (fabs(m.col) < col_);
    };      
  };


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

  /** \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.
   */
131
  struct CmpMatEntryValue : public std::binary_function<MatEntry, MatEntry, bool> {
132
133
134
135
136
137
138
139
140
141
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return m1.entry < m2.entry;
    };
  };


  /** \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.
   */
142
  struct CmpMatEntryAbsValueLess : public std::binary_function<MatEntry, MatEntry, bool> {
143
144
145
146
147
148
149
150
151
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return fabs(m1.entry) < fabs(m2.entry);
    };
  };
  
  /** \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.
   */
152
  struct CmpMatEntryAbsValueGreater : public std::binary_function<MatEntry, MatEntry, bool> {
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return fabs(m1.entry) > fabs(m2.entry);
    };
  };

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

  /** \ingroup DOFAdministration
   * \brief
   * A DOFMatrix is a sparse matrix representation for matrices that work
   * on DOFVectors. Every row of a matrix is realized as a vector of MatEntry 
   * objects. Each entry consists of a column DOF index and the corresponding 
   * double matrix entry. Unused entries are marked with a negative column index.
   */
169
  class DOFMatrix : public DOFIndexed< std::vector<MatEntry> >,
170
171
172
173
174
175
176
                    public Serializable
  {
  public:
    MEMORY_MANAGED(DOFMatrix);

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

    /** \brief
     * Returns an iterator to the begin of matrix rows (\ref matrix.begin())
     */
237
    std::vector< std::vector<MatEntry> >::iterator begin() { 
238
239
240
241
242
243
      return matrix.begin();
    };

    /** \brief
     * Returns an iterator to the end of matrix rows (\ref matrix.end())
     */
244
    std::vector< std::vector<MatEntry> >::iterator end() { 
245
246
247
248
249
250
251
      return matrix.end(); 
    };
 
    /** \brief
     * Used by DOFAdmin to compress the DOFMatrix
     */
    virtual void compressDOFIndexed(int first, int last, 
252
				    std::vector<DegreeOfFreedom> &newDOF);
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

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

    /** \brief
     * Returns whether entry matrix[a][b] is used
     */    
    inline bool entryUsed(DegreeOfFreedom a,int b) const {
      return (((matrix[a])[b]).col >= 0);
    };

    /** \brief
     * Returns true if matrix[a][b] has entry \ref NO_MORE_ENTRIES
     */
    inline bool noMoreEntries(DegreeOfFreedom a,int b) const {
Thomas Witkowski's avatar
Thomas Witkowski committed
270
      return (NO_MORE_ENTRIES == ((matrix[a])[b]).col);
271
272
273
274
275
    };

    /** \brief
     * Returns \ref coupleMatrix.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
    inline bool isCoupleMatrix() { 
      return coupleMatrix; 
    };
279
280
281
282

    /** \brief
     * Returns \ref coupleMatrix.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
285
    inline void setCoupleMatrix(bool c) { 
      coupleMatrix = c; 
    };
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313

    /** \brief
     * Matrix-matrix multiplication.
     */
    void mm(MatrixTranspose aTranspose, DOFMatrix& a, 
	    MatrixTranspose bTranspose, DOFMatrix& b);
  
    /** \brief
     * a*x + y
     */
    void axpy(double a, 
	      const DOFMatrix& x,
	      const DOFMatrix& y);

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

    inline void addOperator(Operator *op, 
			    double* factor = NULL,
			    double* estFactor = NULL) 
    { 
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
    };

314
    inline std::vector<double*>::iterator getOperatorFactorBegin() {
315
316
317
      return operatorFactor.begin();
    };

318
    inline std::vector<double*>::iterator getOperatorFactorEnd() {
319
320
321
      return operatorFactor.end();
    };

322
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
323
324
325
      return operatorEstFactor.begin();
    };

326
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
327
328
329
      return operatorEstFactor.end();
    };

330
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
331
332
333
      return operators.begin();
    };

334
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
      return operators.end();
    };

    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().
     */
  
363
364
365
    void assemble(double factor, ElInfo *elInfo, 
		  const BoundaryType *bound, Operator *op = NULL);
    
366
367
368
369
370
371
372
373
374
375
376
377

    /** \brief
     * Adds an element matrix to \ref matrix
     */
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
			  bool add = true);

    /** \brief
     * Returns \ref matrix
     */
378
    std::vector< std::vector<MatEntry> >& getMatrix() { 
379
380
381
      return matrix; 
    };

382
    void setMatrix(std::vector< std::vector<MatEntry> > m) {
383
384
385
386
387
388
      matrix = m;
    };

    /** \brief
     * Returns \ref matrix[n]
     */
389
    const std::vector<MatEntry>& getRow(int n) const { 
390
391
392
393
394
395
      return matrix[n];
    };

    /** \brief
     * Returns \ref matrix[n]
     */
396
    std::vector<MatEntry>& getRow(int n) { 
397
398
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
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
      return matrix[n];
    };

    /** \brief
     * Returns whether interpolation should be performed after refinement
     * (false by default)
     */
    virtual bool refineInterpol() {
      return false;
    };

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

    /** \brief
     * Returns const \ref rowFESpace
     */
    const FiniteElemSpace* getRowFESpace() const { 
      return rowFESpace; 
    };

    /** \brief
     * Returns const \ref colFESpace
     */
    const FiniteElemSpace* getColFESpace() const { 
      return colFESpace; 
    };

    /** \brief
     * Returns const \ref rowFESpace
     */
    const FiniteElemSpace* getFESpace() const { 
      return rowFESpace; 
    };

    /** \brief
     * Returns number of rows (\ref matrix.size())
     */
    inline int getSize() const { 
      return matrix.size();
    };

Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446
447
448
449
450
451
    /** \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();
    };

452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
    /** \brief
     * Returns number of cols. For that, the function iteratos over all
     * rows and searchs for the entry with the highest col number.
     */
    int getNumCols() const {
      int max = 0;
      int rowSize = static_cast<int>(matrix.size());
      for (int i = 0; i < rowSize; i++) {
	int colSize = static_cast<int>(matrix[i].size());
	for (int j = 0; j < colSize; j++) {
	  if (matrix[i][j].col > max) {
	    max = matrix[i][j].col;
	  }
	}
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
      // Add one to the maximum value, because the indeces start with 0.
      return max + 1;
470
471
    };

472
473
474
    /** \brief
     * Returns \ref name
     */
475
    inline const std::string& getName() const { 
476
477
478
479
480
481
482
      return name; 
    };

    /** \brief
     * Resizes \ref matrix to n rows
     */
    inline void resize(int n) { 
483
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n"); 
484
485
486
487
488
489
      matrix.resize(n);
    };

    /** \brief
     * Returns \ref matrix[i] which is the i-th row
     */
490
    inline const std::vector<MatEntry>& operator[](int i) const {
491
492
      TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
      	("Illegal matrix index %d.\n",i); 
493
494
495
496
497
498
      return matrix[i];
    };

    /** \brief
     * Returns \ref matrix[i] which is the i-th row
     */
499
    inline std::vector<MatEntry>& operator[](int i) {
500
501
      TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
      	("Illegal vector index %d.\n", i);
502
503
504
505
506
507
      return matrix[i];
    };

    /** \brief
     * Access to \ref matrix[a][b].entry
     */
508
    inline double physAcc(DegreeOfFreedom a, DegreeOfFreedom b) const { 
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
      return matrix[a][b].entry; 
    };

    /** \brief
     * Access to \ref matrix[a][b].col
     */
    inline DegreeOfFreedom physToLogIndex(DegreeOfFreedom a,
					  DegreeOfFreedom b) const 
    {
      return matrix[a][b].col;
    };

    /** \brief
     * Returns physical column index of logical index b in row a
     */
524
525
526
527
528
529
530
531
532
    inline DegreeOfFreedom logToPhysIndex(DegreeOfFreedom a, DegreeOfFreedom b) const
    {      
      int size = static_cast<int>(matrix[a].size());
      for (int j = 0; j < size; j++) 
	if (b == matrix[a][j].col) 
	  return j;

      return -1;
    };
533
534
535
536

    /** \brief
     * Returns value at logical indices a,b
     */
537
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
538
539
540
541

    /** \brief
     * Changes col at logical indices a,b to c 
     */
542
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557

    /** \brief
     * Changes col of \ref matrix[a][b] to c
     */
    inline void changePhysColOfEntry(DegreeOfFreedom a,
				     DegreeOfFreedom b,
				     DegreeOfFreedom c) 
    { 
      matrix[a][b].col = c;
    };

    /** \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
     */
558
559
560
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
561
562

    inline double *hasSparseDOFEntry(int irow, int jcol) {
563
564
      std::vector<MatEntry>::iterator it;
      std::vector<MatEntry>::iterator end = matrix[irow].end();
565
566
567
568
569
570
571
572
573
574
575
576
577
      for (it = matrix[irow].begin(); it != end; ++it) {
	if (it->col == NO_MORE_ENTRIES) 
	  return NULL;
	if (it->col == jcol) 
	  return &(it->entry);
      };
      return NULL;
    };

    void addMatEntry(int row, MatEntry entry);

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

578
    void addRow(std::vector<MatEntry> row);
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601

    /** \brief
     * Prints \ref matrix to stdout
     */
    void print() const;

    /** \brief
     * Prints a row of \ref matrix to stdout
     */
    void printRow(int i) const;

    /** \brief
     * Removes all matrix entries
     */
    void clear();

    /** \brief
     * Test whether \ref matrix is symmetric. Exits if not.
     */
    void test();

    bool symmetric();

602
    inline std::vector<Operator*> getOperators() { 
603
604
605
      return operators; 
    };
    
606
    inline std::vector<double*> getOperatorFactor() { 
607
608
609
      return operatorFactor; 
    };

610
    inline std::vector<double*> getOperatorEstFactor() { 
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
      return operatorEstFactor; 
    };

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

    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
    };

    void createPictureFile(const char* filename, int dim);

    // ===== Serializable implementation =====
    
626
    void serialize(std::ostream &out) 
627
628
629
630
631
632
633
634
635
636
637
638
    {
      unsigned int matrixSize = matrix.size();
      unsigned int vecSize = 0;

      out.write(reinterpret_cast<const char*>(&matrixSize), sizeof(unsigned int));
      for(unsigned int i = 0; i < matrixSize; i++) {
	vecSize = matrix[i].size();
	out.write(reinterpret_cast<const char*>(&vecSize), sizeof(unsigned int));
	out.write(reinterpret_cast<const char*>(&(matrix[i][0])), vecSize * sizeof(MatEntry));
      }
    };

639
    void deserialize(std::istream &in) {
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
      unsigned int matrixSize, vecSize;

      in.read(reinterpret_cast<char*>(&matrixSize), sizeof(unsigned int));
      matrix.resize(matrixSize);
      for(unsigned int i = 0; i < matrixSize; i++) {
	in.read(reinterpret_cast<char*>(&vecSize), sizeof(unsigned int));
	matrix[i].resize(vecSize);
	in.read(reinterpret_cast<char*>(&(matrix[i][0])), vecSize * sizeof(MatEntry));
      }
    };

    int memsize();

  protected:
    /** \brief
     * Used by updateMatrix
     */
    //static int assembleFct(ElInfo *el_info);

  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;

    /** \brief
     * Name of the DOFMatrix
     */
675
    std::string name;
676
677
678
679

    /** \brief
     * Sparse matrix stored in an STL vector of MatrixRow objects
     */
680
    std::vector<MatrixRow> matrix;
681
682
683
684
685
686
687
688
689
690

    /** \brief
     * Used while mesh traversal
     */
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
691
    std::vector<Operator*> operators;
692
693
694
695
696
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
697
    std::vector<double*> operatorFactor;
698

699
700
701
    /** \brief
     *
     */
702
    std::vector<double*> operatorEstFactor;
703

704
705
706
    /** \brief
     *
     */
707
708
    BoundaryManager *boundaryManager;

709
710
711
    /** \brief
     *
     */
712
713
    bool coupleMatrix;

Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
716
717
718
    /** \brief
     * Temporary variable used in assemble()
     */
    ElementMatrix *elementMatrix;

719
720
721
722
723
724
725
726
727
728
729
730
731
    friend class DOFAdmin;
    friend class DOFVector<double>;
    friend class DOFVector<unsigned char>;
    friend class DOFVector<int>;
    friend class DOFVector<WorldVector<double> >;
  };


  inline DegreeOfFreedom logToPhysIndex(DOFMatrix *m, DegreeOfFreedom a, DegreeOfFreedom b)
  {
    return m->logToPhysIndex(a, b);
  }

732
  double norm(std::vector<MatEntry> *row);
733

734
  double min(std::vector<MatEntry> *row);
735

736
  double max(std::vector<MatEntry> *row);
737
738
739
740

}

#endif  // AMDIS_DOFMATRIX_H