DOFMatrix.h 20.2 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,
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309

    /** \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 {
      return  (NO_MORE_ENTRIES==((matrix[a])[b]).col);
    };

    /** \brief
     * Returns \ref coupleMatrix.
     */
    inline bool isCoupleMatrix() { return coupleMatrix; };

    /** \brief
     * Returns \ref coupleMatrix.
     */
    inline void setCoupleMatrix(bool c) { coupleMatrix = c; };

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

310
    inline std::vector<double*>::iterator getOperatorFactorBegin() {
311
312
313
      return operatorFactor.begin();
    };

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

318
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
319
320
321
      return operatorEstFactor.begin();
    };

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

326
    inline std::vector<Operator*>::iterator getOperatorsBegin() {
327
328
329
      return operators.begin();
    };

330
    inline std::vector<Operator*>::iterator getOperatorsEnd() {
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
      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().
     */
  
359
360
361
    void assemble(double factor, ElInfo *elInfo, 
		  const BoundaryType *bound, Operator *op = NULL);
    
362
363
364
365
366
367
368
369
370
371
372
373

    /** \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
     */
374
    std::vector< std::vector<MatEntry> >& getMatrix() { 
375
376
377
      return matrix; 
    };

378
    void setMatrix(std::vector< std::vector<MatEntry> > m) {
379
380
381
382
383
384
      matrix = m;
    };

    /** \brief
     * Returns \ref matrix[n]
     */
385
    const std::vector<MatEntry>& getRow(int n) const { 
386
387
388
389
390
391
      return matrix[n];
    };

    /** \brief
     * Returns \ref matrix[n]
     */
392
    std::vector<MatEntry>& getRow(int n) { 
393
394
395
396
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
      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();
    };

440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
    /** \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;
	  }
	}
      }

      return max;
    };

459
460
461
    /** \brief
     * Returns \ref name
     */
462
    inline const std::string& getName() const { 
463
464
465
466
467
468
469
      return name; 
    };

    /** \brief
     * Resizes \ref matrix to n rows
     */
    inline void resize(int n) { 
470
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n"); 
471
472
473
474
475
476
      matrix.resize(n);
    };

    /** \brief
     * Returns \ref matrix[i] which is the i-th row
     */
477
    inline const std::vector<MatEntry>& operator[](int i) const {
478
479
      TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
      	("Illegal matrix index %d.\n",i); 
480
481
482
483
484
485
      return matrix[i];
    };

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

    /** \brief
     * Access to \ref matrix[a][b].entry
     */
495
    inline double physAcc(DegreeOfFreedom a, DegreeOfFreedom b) const { 
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
      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
     */
511
512
513
514
515
516
517
518
519
    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;
    };
520
521
522
523

    /** \brief
     * Returns value at logical indices a,b
     */
524
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
525
526
527
528

    /** \brief
     * Changes col at logical indices a,b to c 
     */
529
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544

    /** \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
     */
545
546
547
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
548
549

    inline double *hasSparseDOFEntry(int irow, int jcol) {
550
551
      std::vector<MatEntry>::iterator it;
      std::vector<MatEntry>::iterator end = matrix[irow].end();
552
553
554
555
556
557
558
559
560
561
562
563
564
      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);

565
    void addRow(std::vector<MatEntry> row);
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588

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

589
    inline std::vector<Operator*> getOperators() { 
590
591
592
      return operators; 
    };
    
593
    inline std::vector<double*> getOperatorFactor() { 
594
595
596
      return operatorFactor; 
    };

597
    inline std::vector<double*> getOperatorEstFactor() { 
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
      return operatorEstFactor; 
    };

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

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

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

    // ===== Serializable implementation =====
    
613
    void serialize(std::ostream &out) 
614
615
616
617
618
619
620
621
622
623
624
625
    {
      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));
      }
    };

626
    void deserialize(std::istream &in) {
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
      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
     */
662
    std::string name;
663
664
665
666

    /** \brief
     * Sparse matrix stored in an STL vector of MatrixRow objects
     */
667
    std::vector<MatrixRow> matrix;
668
669
670
671
672
673
674
675
676
677

    /** \brief
     * Used while mesh traversal
     */
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
678
    std::vector<Operator*> operators;
679
680
681
682
683
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
684
    std::vector<double*> operatorFactor;
685

686
687
688
    /** \brief
     *
     */
689
    std::vector<double*> operatorEstFactor;
690

691
692
693
    /** \brief
     *
     */
694
695
    BoundaryManager *boundaryManager;

696
697
698
    /** \brief
     *
     */
699
700
    bool coupleMatrix;

Thomas Witkowski's avatar
Thomas Witkowski committed
701
702
703
704
705
    /** \brief
     * Temporary variable used in assemble()
     */
    ElementMatrix *elementMatrix;

706
707
708
709
710
711
712
713
714
715
716
717
718
    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);
  }

719
  double norm(std::vector<MatEntry> *row);
720

721
  double min(std::vector<MatEntry> *row);
722

723
  double max(std::vector<MatEntry> *row);
724
725
726
727

}

#endif  // AMDIS_DOFMATRIX_H