Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

DOFMatrix.h 14.7 KB
Newer Older
1 2 3 4 5 6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9 10 11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12 13 14 15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17 18 19 20 21 22 23 24 25
// ==                                                                        ==
// ============================================================================

/** \file DOFMatrix.h */

#ifndef AMDIS_DOFMATRIX_H
#define AMDIS_DOFMATRIX_H

#include <vector>
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include <set>
27 28
#include <memory>
#include <list>
29 30
#include <boost/numeric/mtl/mtl.hpp>
#include "AMDiS_fwd.h"
31 32 33 34 35 36 37 38 39 40 41 42 43
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
#include "DOFAdmin.h"
#include "DOFIndexed.h"
#include "Boundary.h"
#include "Serializable.h"

namespace AMDiS {

  /** \ingroup DOFAdministration
   * \brief
   * A DOFMatrix is a sparse matrix representation for matrices that work
44 45
   * on DOFVectors. The underlying matrix type 
   * is a CRS matrix of double.
46
   */
47
  class DOFMatrix : public DOFIndexed<bool>,
48 49 50
                    public Serializable
  {
  public:
51
    /// Type of scalars in the underlying matrix
52
    typedef double value_type;
53 54

    /// Type of underlying matrix
55
    typedef mtl::compressed2D<value_type> base_matrix_type;
56 57

    /// Type of inserter for the base matrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
58
    typedef mtl::matrix::inserter<base_matrix_type, mtl::operations::update_plus<value_type> >  inserter_type;
59 60

  private:
Thomas Witkowski's avatar
Thomas Witkowski committed
61
    /// Symbolic constant for an unused matrix entry
62 63 64 65 66 67 68 69 70 71 72
    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();

73 74
    /// Constructs a DOFMatrix with name n and the given row and olumn FeSpaces.
    DOFMatrix(const FiniteElemSpace* rowFeSpace, const FiniteElemSpace* colFeSpace,
Thomas Witkowski's avatar
Thomas Witkowski committed
75
	      std::string n = "");
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
    /// Copy-Constructor
78 79
    DOFMatrix(const DOFMatrix& rhs);

Thomas Witkowski's avatar
Thomas Witkowski committed
80
    /// Destructor
81 82
    virtual ~DOFMatrix();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    /// Assignment operator.
84 85 86 87
    DOFMatrix& operator=(const DOFMatrix& rhs);

    void copy(const DOFMatrix& rhs);

88 89 90 91
    /// Access underlying matrix directly
    base_matrix_type& getBaseMatrix()
    {
	return matrix;
92
    }
93

94 95 96 97
    /// Access underlying matrix directly (const)
    const base_matrix_type& getBaseMatrix() const
    {
	return matrix;
98
    }
99 100

    // Only to get rid of the abstract functions, I hope they are not used
Thomas Witkowski's avatar
Thomas Witkowski committed
101 102 103 104 105 106 107 108 109 110 111
    std::vector<bool>::iterator begin() 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; 
      return v.begin();
    }
    
    std::vector<bool>::iterator end() 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; 
      return v.end();
    }
112 113
    
    bool dummy; // Must be deleted later
Thomas Witkowski's avatar
Thomas Witkowski committed
114 115 116 117 118 119 120 121 122 123 124
    bool& operator[](int i) 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); 
      return dummy;
    }

    const bool& operator[](int i) const 
    {
      ERROR_EXIT("Shouldn't be used, only fake."); 
      return dummy;
    }
125
 
126
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
Thomas Witkowski's avatar
Thomas Witkowski committed
127 128
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) 
    {}
129

130
    /// Implements DOFIndexedBase::freeDOFContent()
131 132
    virtual void freeDOFContent(int index);

133
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
134 135
    inline bool isCoupleMatrix() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
136
      return coupleMatrix; 
137
    }
138

139
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
140 141
    inline void setCoupleMatrix(bool c) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
142
      coupleMatrix = c; 
143
    }
144

145
    /// a * x + y
Thomas Witkowski's avatar
Thomas Witkowski committed
146
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
147

148
    /// Multiplication with a scalar.
149 150
    void scal(double s);

Thomas Witkowski's avatar
Thomas Witkowski committed
151 152 153 154 155 156
    /** \brief
     * Adds an operator to the DOFMatrix. A factor, that is multipled
     * to the operator, and a multilier factor for the estimator may be
     * also given.
     */
    void addOperator(Operator *op, double* factor = NULL, double* estFactor = NULL);
157

Thomas Witkowski's avatar
Thomas Witkowski committed
158 159
    inline std::vector<double*>::iterator getOperatorFactorBegin() 
    {
160
      return operatorFactor.begin();
161
    }
162

Thomas Witkowski's avatar
Thomas Witkowski committed
163 164
    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
165
      return operatorFactor.end();
166
    }
167

Thomas Witkowski's avatar
Thomas Witkowski committed
168 169
    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
170
      return operatorEstFactor.begin();
171
    }
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173 174
    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
175
      return operatorEstFactor.end();
176
    }
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178 179
    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
180
      return operators.begin();
181
    }
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183 184
    inline std::vector<Operator*>::iterator getOperatorsEnd() 
    {
185
      return operators.end();
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

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

215 216
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
217

218 219 220
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
221 222
		  const BoundaryType *bound,
		  Operator *op = NULL);
223

Thomas Witkowski's avatar
Thomas Witkowski committed
224 225 226 227 228 229 230
    void assemble2(double factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = NULL);

    /// Adds an element matrix to \ref matrix
Thomas Witkowski's avatar
Thomas Witkowski committed
231
    void addElementMatrix(const ElementMatrix& elMat, 
232
			  const BoundaryType *bound,
233
			  ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
234
			  ElInfo* colElInfo);
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236 237 238 239 240 241 242
    /* \brief
     * That function must be called after the matrix assembling has been finished. 
     * This makes it possible to start some cleanup or matrix data compressing 
     * procedures.
     */
    void finishAssembling();

243 244 245 246 247 248
    /** \brief
     * Enable insertion for assembly. You can optionally give an upper limit for
     * entries per row (per column for CCS matrices).  Choosing this parameter
     * too small can induce perceivable overhead for compressed matrices.  Thus,
     * it's better to choose a bit too large than too small.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
249
    void startInsertion(int nnz_per_row = 10);
250

251
    /** \brief
252 253
     * Finishes insertion. For compressed matrix types, this is where the 
     * compression happens.
254 255 256
     */
    void finishInsertion()
    {
257 258 259
      FUNCNAME("DOFMatrix::finishInsertion()");

      TEST_EXIT(inserter)("Inserter wasn't used or is already finished.\n");
260 261 262
      
      delete inserter;
      inserter= 0;
263
    }
264 265 266 267 268

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
269 270
    virtual bool coarseRestrict() 
    {
271
      return false;
272
    }
273

274 275
    /// Returns const \ref rowFeSpace
    const FiniteElemSpace* getRowFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
276
    { 
277
      return rowFeSpace; 
278
    }
279

280 281
    /// Returns const \ref colFeSpace
    const FiniteElemSpace* getColFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    { 
283
      return colFeSpace; 
284
    }
285

286 287
    /// Returns const \ref rowFeSpace
    const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
288
    { 
289
      return rowFeSpace; 
290
    }
291

292
    /// Returns number of rows (\ref matrix.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
293 294
    inline int getSize() const 
    { 
295
      return num_rows(matrix);
296
    }
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298 299 300 301
    /** \brief
     * Returns the number of used rows (equal to number of used DOFs in
     * the row FE space).
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
302 303
    inline int getUsedSize() const 
    {
304
      return rowFeSpace->getAdmin()->getUsedSize();
305
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
306

307
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
308
    inline std::string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
309
    { 
310
      return name; 
311
    }
312

313
    /// Resizes \ref matrix to n rows
Thomas Witkowski's avatar
Thomas Witkowski committed
314 315
    inline void resize(int n) 
    { 
316
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
317
    }
318

319
    /// Returns value at logical indices a,b
320
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
321

322
    /// Changes col at logical indices a,b to c 
323
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
324 325 326 327 328

    /** \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
     */
329 330 331
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
332

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

335
    /// Prints \ref matrix to stdout
336 337
    void print() const;

338
    /// Removes all matrix entries
339 340 341 342
    void clear()
    {
	set_to_zero(matrix);
    }
343

344
    /// Test whether \ref matrix is symmetric. Exits if not.
345 346 347 348
    void test();

    bool symmetric();

349
    inline std::vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    { 
351
      return operators; 
352
    }
353
    
354
    inline std::vector<double*>& getOperatorFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    { 
356
      return operatorFactor; 
357
    }
358

359
    inline std::vector<double*>& getOperatorEstFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    { 
361
      return operatorEstFactor; 
362
    }
363

Thomas Witkowski's avatar
Thomas Witkowski committed
364 365
    inline BoundaryManager* getBoundaryManager() const 
    { 
366
      return boundaryManager; 
367
    }
368

369
    /// Returns a pointer to \ref applyDBCs.
Thomas Witkowski's avatar
Thomas Witkowski committed
370 371
    std::set<int>* getApplyDBCs() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
372 373 374
      return &applyDBCs;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
375 376
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
377
      boundaryManager = bm;
378
    }
379

Thomas Witkowski's avatar
Thomas Witkowski committed
380 381 382 383 384 385 386
    /// Calculates the average of non zero entries per row in matrix.
    void calculateNnz()
    {
      nnzPerRow = 0;

      if (num_rows(matrix) != 0)
	nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2); 
Thomas Witkowski's avatar
Thomas Witkowski committed
387 388
      if (nnzPerRow < 20) 
	nnzPerRow = 20;
Thomas Witkowski's avatar
Thomas Witkowski committed
389 390 391 392 393 394 395 396
    }

    /// Returns \ref nnzPerRow.
    int getNnz()
    {
      return nnzPerRow;
    }

397 398
    /// Writes the matrix to an output stream.
    void serialize(std::ostream &out);
399

400 401
    /// Reads a matrix from an input stream.
    void deserialize(::std::istream &in);
402

403
    ///
404 405
    int memsize();

406
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
407
    void setRankDofs(std::map<DegreeOfFreedom, bool>& dofmap)
Thomas Witkowski's avatar
Thomas Witkowski committed
408
    {
409
      rankDofs = &dofmap;
410 411
    }
#endif
412 413 414 415 416 417

  protected:
    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding row DOFs
     * and basis functions
     */
418
    const FiniteElemSpace *rowFeSpace;
419 420 421 422 423

    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding 
     * column DOFs and basis functions
     */
424
    const FiniteElemSpace *colFeSpace;
425

426
    /// Name of the DOFMatrix
427
    std::string name;
428

429
    /// Sparse matrix, type is a template parameter by default compressed2D<double>
430
    base_matrix_type matrix;
431

432
    /// Used while mesh traversal
433 434 435 436 437 438
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
439
    std::vector<Operator*> operators;
440 441 442 443 444
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
445
    std::vector<double*> operatorFactor;
446

447
    ///
448
    std::vector<double*> operatorEstFactor;
449

450
    ///
451 452
    BoundaryManager *boundaryManager;

453 454 455 456
    /** \brief
     * If false, the matrix is a diagonal matrix within a matrix of DOF matrices.
     * Otherwise the value is true, and the matrix is an off-diagonal matrix.
     */
457 458
    bool coupleMatrix;

459
    /// Temporary variable used in assemble()
460 461 462 463 464 465 466 467 468
    ElementMatrix elementMatrix;

    /// Number of basis functions in the row fe space.
    int nRow;

    /// Number of basis functions in the col fe space.
    int nCol;

    /// Maps local row indices of an element to global matrix indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
469
    std::vector<DegreeOfFreedom> rowIndices;
470 471

    /// Maps local col indices of an element to global matrix indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
472
    std::vector<DegreeOfFreedom> colIndices;
Thomas Witkowski's avatar
Thomas Witkowski committed
473

474 475 476 477 478 479 480
    /* \brief
     * A set of row indices. When assembling the DOF matrix, all rows, that
     * correspond to a dof at a dirichlet boundary, are ignored and the row is
     * left blank. After assembling, the diagonal element of the matrix must be
     * set to 1. The indices of all rows, where this should be done, are stored
     * in this set.
     */    
Thomas Witkowski's avatar
Thomas Witkowski committed
481 482
    std::set<int> applyDBCs;

Thomas Witkowski's avatar
Thomas Witkowski committed
483 484 485 486 487 488 489
    /* \brief
     * Number of non zero entries per row (average). For instationary problems this
     * information may be used in the next timestep to accelerate insertion of
     * elemnts during assembling.
     */
     int nnzPerRow;

490
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
491
    std::map<DegreeOfFreedom, bool> *rankDofs;
492 493
#endif

494
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
495
    inserter_type *inserter;
496
      
497 498 499 500 501 502 503 504 505 506
    friend class DOFAdmin;
    friend class DOFVector<double>;
    friend class DOFVector<unsigned char>;
    friend class DOFVector<int>;
    friend class DOFVector<WorldVector<double> >;
  };

}

#endif  // AMDIS_DOFMATRIX_H