DOFMatrix.h 19.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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
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
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
// ============================================================================
// ==                                                                        ==
// == 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.
   */
  class MatEntryValueLess : public ::std::unary_function<MatEntry, bool> {
  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.
   */
  class MatEntryColLess : public ::std::unary_function<MatEntry, bool> {
  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).
   */
  struct CmpMatEntryCol : public ::std::binary_function<MatEntry, MatEntry, bool> {
    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.
   */
  struct CmpMatEntryValue : public ::std::binary_function<MatEntry, MatEntry, bool> {
    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.
   */
  struct CmpMatEntryAbsValueLess : public ::std::binary_function<MatEntry, MatEntry, bool> {
    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.
   */
  struct CmpMatEntryAbsValueGreater : public ::std::binary_function<MatEntry, MatEntry, bool> {
    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.
   */
  class DOFMatrix : public DOFIndexed< ::std::vector<MatEntry> >,
                    public Serializable
  {
  public:
    MEMORY_MANAGED(DOFMatrix);

    /** \ingroup DOFAdministration
     * \brief
     * Alias for DOFIterator< ::std::vector<MatEntry<T> > >. So one can access an
     * iterator working on a DOFMatrix via DOFMatrix::Iterator
     */
    class Iterator : public DOFIterator< ::std::vector<MatEntry> > {
    public:
      Iterator(DOFIndexed< ::std::vector<MatEntry> > *c, 
	       DOFIteratorType type)
	: DOFIterator< ::std::vector<MatEntry> >(c, type)
      {};   
    };
    
    /** \brief
     * A MatrixRow represents one row of the sparse matrix. It is realized by
     * a STL vector of MatEntry<T> objects
     */
    typedef ::std::vector<MatEntry> MatrixRow;
    
    /** \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,
	      ::std::string n="");

    /** \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())
     */
    ::std::vector< ::std::vector<MatEntry> >::iterator begin() { 
      return matrix.begin();
    };

    /** \brief
     * Returns an iterator to the end of matrix rows (\ref matrix.end())
     */
    ::std::vector< ::std::vector<MatEntry> >::iterator end() { 
      return matrix.end(); 
    };
 
    /** \brief
     * Used by DOFAdmin to compress the DOFMatrix
     */
    virtual void compressDOFIndexed(int first, int last, 
				    ::std::vector<DegreeOfFreedom> &newDOF);

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

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

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

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

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

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

    inline ::std::vector<Operator*>::iterator getOperatorsEnd() {
      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
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
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
440
441
442
443
444
445
446
447
448
449
450

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

    void setMatrix(::std::vector< ::std::vector<MatEntry> > m) {
      matrix = m;
    };

    /** \brief
     * Returns \ref matrix[n]
     */
    const ::std::vector<MatEntry>& getRow(int n) const { 
      return matrix[n];
    };

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

    /** \brief
     * Returns \ref name
     */
    inline const ::std::string& getName() const { 
      return name; 
    };

    /** \brief
     * Resizes \ref matrix to n rows
     */
    inline void resize(int n) { 
451
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n"); 
452
453
454
455
456
457
458
      matrix.resize(n);
    };

    /** \brief
     * Returns \ref matrix[i] which is the i-th row
     */
    inline const ::std::vector<MatEntry>& operator[](int i) const {
459
460
      TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
      	("Illegal matrix index %d.\n",i); 
461
462
463
464
465
466
467
      return matrix[i];
    };

    /** \brief
     * Returns \ref matrix[i] which is the i-th row
     */
    inline ::std::vector<MatEntry>& operator[](int i) {
468
469
      TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
      	("Illegal vector index %d.\n", i);
470
471
472
473
474
475
      return matrix[i];
    };

    /** \brief
     * Access to \ref matrix[a][b].entry
     */
476
    inline double physAcc(DegreeOfFreedom a, DegreeOfFreedom b) const { 
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
      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
     */
    DegreeOfFreedom logToPhysIndex(DegreeOfFreedom a,DegreeOfFreedom b) const;

    /** \brief
     * Returns value at logical indices a,b
     */
    double logAcc(DegreeOfFreedom a,DegreeOfFreedom b) const;

    /** \brief
     * Changes col at logical indices a,b to c 
     */
    void changeColOfEntry(DegreeOfFreedom a,DegreeOfFreedom b,DegreeOfFreedom c);

    /** \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
     */
518
519
520
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
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

    inline double *hasSparseDOFEntry(int irow, int jcol) {
      ::std::vector<MatEntry>::iterator it;
      ::std::vector<MatEntry>::iterator end = matrix[irow].end();
      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);

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

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

    inline ::std::vector<Operator*> getOperators() { 
      return operators; 
    };
    
    inline ::std::vector<double*> getOperatorFactor() { 
      return operatorFactor; 
    };

    inline ::std::vector<double*> getOperatorEstFactor() { 
      return operatorEstFactor; 
    };

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

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

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

    // ===== Serializable implementation =====
    
    void serialize(::std::ostream &out) 
    {
      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));
      }
    };

    void deserialize(::std::istream &in) {
      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
     */
    ::std::string name;

    /** \brief
     * Sparse matrix stored in an STL vector of MatrixRow objects
     */
    ::std::vector<MatrixRow> matrix;

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

659
660
661
    /** \brief
     *
     */
662
663
    ::std::vector<double*> operatorEstFactor;

664
665
666
    /** \brief
     *
     */
667
668
    BoundaryManager *boundaryManager;

669
670
671
    /** \brief
     *
     */
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    bool coupleMatrix;

    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);
  }

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

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

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

}

#endif  // AMDIS_DOFMATRIX_H