DOFMatrix.h 20.7 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
    void assemble(double factor, ElInfo *elInfo, 
		  const BoundaryType *bound, Operator *op = NULL);
365
366
367
368


    void assemble(double factor, ElInfo *rowElInfo, ElInfo *colElInfo,
		  const BoundaryType *bound, Operator *op = NULL);
369
    
370
371
372
373
374
375
376
377
378
379
380
381

    /** \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
     */
382
    std::vector< std::vector<MatEntry> >& getMatrix() { 
383
384
385
      return matrix; 
    };

386
    void setMatrix(std::vector< std::vector<MatEntry> > m) {
387
388
389
390
391
392
      matrix = m;
    };

    /** \brief
     * Returns \ref matrix[n]
     */
393
    const std::vector<MatEntry>& getRow(int n) const { 
394
395
396
397
398
399
      return matrix[n];
    };

    /** \brief
     * Returns \ref matrix[n]
     */
400
    std::vector<MatEntry>& getRow(int n) { 
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
444
445
446
447
      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
448
449
450
451
452
453
454
455
    /** \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();
    };

456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
    /** \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
472
473
      // Add one to the maximum value, because the indeces start with 0.
      return max + 1;
474
475
    };

476
477
478
    /** \brief
     * Returns \ref name
     */
479
    inline const std::string& getName() const { 
480
481
482
483
484
485
486
      return name; 
    };

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

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

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

    /** \brief
     * Access to \ref matrix[a][b].entry
     */
512
    inline double physAcc(DegreeOfFreedom a, DegreeOfFreedom b) const { 
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
      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
     */
528
529
530
531
532
533
534
535
536
    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;
    };
537
538
539
540

    /** \brief
     * Returns value at logical indices a,b
     */
541
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
542
543
544
545

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

    /** \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
     */
562
563
564
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
565
566

    inline double *hasSparseDOFEntry(int irow, int jcol) {
567
568
      std::vector<MatEntry>::iterator it;
      std::vector<MatEntry>::iterator end = matrix[irow].end();
569
570
571
572
573
574
575
576
577
578
579
580
581
      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);

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

    /** \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();

606
    inline std::vector<Operator*> getOperators() { 
607
608
609
      return operators; 
    };
    
610
    inline std::vector<double*> getOperatorFactor() { 
611
612
613
      return operatorFactor; 
    };

614
    inline std::vector<double*> getOperatorEstFactor() { 
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
      return operatorEstFactor; 
    };

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

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

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

    // ===== Serializable implementation =====
    
630
    void serialize(std::ostream &out) 
631
632
633
634
635
636
637
638
639
640
641
642
    {
      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));
      }
    };

643
    void deserialize(std::istream &in) {
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
675
676
677
678
      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
     */
679
    std::string name;
680
681
682
683

    /** \brief
     * Sparse matrix stored in an STL vector of MatrixRow objects
     */
684
    std::vector<MatrixRow> matrix;
685
686
687
688
689
690
691
692
693
694

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

703
704
705
    /** \brief
     *
     */
706
    std::vector<double*> operatorEstFactor;
707

708
709
710
    /** \brief
     *
     */
711
712
    BoundaryManager *boundaryManager;

713
714
715
    /** \brief
     *
     */
716
717
    bool coupleMatrix;

Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
720
721
722
    /** \brief
     * Temporary variable used in assemble()
     */
    ElementMatrix *elementMatrix;

723
724
725
726
727
728
729
730
731
732
733
734
735
    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);
  }

736
  double norm(std::vector<MatEntry> *row);
737

738
  double min(std::vector<MatEntry> *row);
739

740
  double max(std::vector<MatEntry> *row);
741

742
743
  void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a, const DOFMatrix *b);

744
745
746
}

#endif  // AMDIS_DOFMATRIX_H