DOFMatrix.h 14.3 KB
Newer Older
1 2 3 4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6 7
// ==                                                                        ==
// ============================================================================
8 9 10 11 12 13 14 15 16 17 18 19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20 21 22 23 24 25 26

/** \file DOFMatrix.h */

#ifndef AMDIS_DOFMATRIX_H
#define AMDIS_DOFMATRIX_H

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

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
42 43
  using namespace std;

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

57 58 59
    typedef unsigned size_type;
    typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;

60
    /// Type of underlying matrix
61
    typedef mtl::compressed2D<value_type, para> base_matrix_type;
62 63

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

  private:
Thomas Witkowski's avatar
Thomas Witkowski committed
67
    /// Symbolic constant for an unused matrix entry
68 69
    static const int UNUSED_ENTRY = -1;
    
Thomas Witkowski's avatar
Thomas Witkowski committed
70 71
    /// Symbolic constant for an unused entry which is not followed by any
    /// used entry in this row
72 73 74 75 76
    static const int NO_MORE_ENTRIES = -2;

  public:    
    DOFMatrix();

77
    /// Constructs a DOFMatrix with name n and the given row and olumn FeSpaces.
Thomas Witkowski's avatar
Thomas Witkowski committed
78 79 80
    DOFMatrix(const FiniteElemSpace* rowFeSpace, 
	      const FiniteElemSpace* colFeSpace,
	      string n = "");
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82
    /// Copy-Constructor
83 84
    DOFMatrix(const DOFMatrix& rhs);

Thomas Witkowski's avatar
Thomas Witkowski committed
85
    /// Destructor
86 87
    virtual ~DOFMatrix();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    /// Assignment operator.
89 90 91 92
    DOFMatrix& operator=(const DOFMatrix& rhs);

    void copy(const DOFMatrix& rhs);

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

99 100 101 102
    /// Access underlying matrix directly (const)
    const base_matrix_type& getBaseMatrix() const
    {
	return matrix;
103
    }
104 105

    // Only to get rid of the abstract functions, I hope they are not used
Thomas Witkowski's avatar
Thomas Witkowski committed
106
    vector<bool>::iterator begin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
107
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
108
      ERROR_EXIT("Shouldn't be used, only fake."); vector<bool> v; 
Thomas Witkowski's avatar
Thomas Witkowski committed
109 110 111
      return v.begin();
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
112
    vector<bool>::iterator end() 
Thomas Witkowski's avatar
Thomas Witkowski committed
113
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
114
      ERROR_EXIT("Shouldn't be used, only fake."); vector<bool> v; 
Thomas Witkowski's avatar
Thomas Witkowski committed
115 116
      return v.end();
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
117 118 119 120 121 122 123 124 125 126 127 128 129

    // Only to get rid of the abstract functions, I hope they are not used
    vector<bool>::const_iterator begin() const
    {
      ERROR_EXIT("Shouldn't be used, only fake."); vector<bool> v;
      return v.begin();
    }

    vector<bool>::const_iterator end() const
    {
      ERROR_EXIT("Shouldn't be used, only fake."); vector<bool> v;
      return v.end();
    }
130 131
    
    bool dummy; // Must be deleted later
Thomas Witkowski's avatar
Thomas Witkowski committed
132 133 134 135 136 137 138 139 140 141 142
    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;
    }
143
 
144
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
Thomas Witkowski's avatar
Thomas Witkowski committed
145
    void compressDOFIndexed(int first, int last, vector<DegreeOfFreedom> &newDOF) 
Thomas Witkowski's avatar
Thomas Witkowski committed
146
    {}
147

148
    /// Implements DOFIndexedBase::freeDOFContent()
149 150
    virtual void freeDOFContent(int index);

151
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
152 153
    inline bool isCoupleMatrix() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
154
      return coupleMatrix; 
155
    }
156

157
    /// Returns \ref coupleMatrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
158 159
    inline void setCoupleMatrix(bool c) 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
160
      coupleMatrix = c; 
161
    }
162

163
    /// a * x + y
Thomas Witkowski's avatar
Thomas Witkowski committed
164
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
165

166
    /// Multiplication with a scalar.
167 168
    void scal(double s);

Thomas Witkowski's avatar
Thomas Witkowski committed
169 170 171 172 173
    /// 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);
174

Thomas Witkowski's avatar
Thomas Witkowski committed
175
    inline vector<double*>::iterator getOperatorFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
176
    {
177
      return operatorFactor.begin();
178
    }
179

Thomas Witkowski's avatar
Thomas Witkowski committed
180
    inline vector<double*>::iterator getOperatorFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    {
182
      return operatorFactor.end();
183
    }
184

Thomas Witkowski's avatar
Thomas Witkowski committed
185
    inline vector<double*>::iterator getOperatorEstFactorBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
186
    {
187
      return operatorEstFactor.begin();
188
    }
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
    inline vector<double*>::iterator getOperatorEstFactorEnd() 
Thomas Witkowski's avatar
Thomas Witkowski committed
191
    {
192
      return operatorEstFactor.end();
193
    }
194

Thomas Witkowski's avatar
Thomas Witkowski committed
195
    inline vector<Operator*>::iterator getOperatorsBegin() 
Thomas Witkowski's avatar
Thomas Witkowski committed
196
    {
197
      return operators.begin();
198
    }
199

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

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

232 233
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);
234

235 236 237
    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
238 239
		  const BoundaryType *bound,
		  Operator *op = NULL);
240

Thomas Witkowski's avatar
Thomas Witkowski committed
241 242 243 244 245 246 247
    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
248
    void addElementMatrix(const ElementMatrix& elMat, 
249
			  const BoundaryType *bound,
250
			  ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
251
			  ElInfo* colElInfo);
252

Thomas Witkowski's avatar
Thomas Witkowski committed
253 254 255
    /// 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.
Thomas Witkowski's avatar
Thomas Witkowski committed
256 257
    void finishAssembling();

Thomas Witkowski's avatar
Thomas Witkowski committed
258 259 260 261
    /// 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
262
    void startInsertion(int nnz_per_row = 10);
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264 265
    /// Finishes insertion. For compressed matrix types, this is where the
    /// compression happens.
266 267
    void finishInsertion()
    {
268 269 270
      FUNCNAME("DOFMatrix::finishInsertion()");

      TEST_EXIT(inserter)("Inserter wasn't used or is already finished.\n");
271 272 273
      
      delete inserter;
      inserter= 0;
274
    }
275

276
#if 0
Thomas Witkowski's avatar
Thomas Witkowski committed
277 278
    /// Returns whether restriction should be performed after coarsening
    /// (false by default)
Thomas Witkowski's avatar
Thomas Witkowski committed
279 280
    virtual bool coarseRestrict() 
    {
281
      return false;
282
    }
283
#endif
284

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

291 292
    /// Returns const \ref colFeSpace
    const FiniteElemSpace* getColFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    { 
294
      return colFeSpace; 
295
    }
296

297 298
    /// Returns const \ref rowFeSpace
    const FiniteElemSpace* getFeSpace() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
299
    { 
300
      return rowFeSpace; 
301
    }
302

303
    /// Returns number of rows (\ref matrix.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
304 305
    inline int getSize() const 
    { 
306
      return num_rows(matrix);
307
    }
308

Thomas Witkowski's avatar
Thomas Witkowski committed
309 310
    /// Returns the number of used rows (equal to number of used DOFs in
    /// the row FE space).
Thomas Witkowski's avatar
Thomas Witkowski committed
311 312
    inline int getUsedSize() const 
    {
313
      return rowFeSpace->getAdmin()->getUsedSize();
314
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
315

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
316 317 318 319 320
    std::set<DegreeOfFreedom>& getDirichletRows()
    {
      return dirichletDofs;
    }

321
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
322
    inline string getName() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
323
    { 
324
      return name; 
325
    }
326

327
    /// Resizes \ref matrix to n rows
Thomas Witkowski's avatar
Thomas Witkowski committed
328 329
    inline void resize(int n) 
    { 
330
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
331
    }
332

333
    /// Returns value at logical indices a,b
334
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
335

336
    /// Changes col at logical indices a,b to c 
337
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);
338

Thomas Witkowski's avatar
Thomas Witkowski committed
339 340
    /// 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
341 342 343
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);
344

Thomas Witkowski's avatar
Thomas Witkowski committed
345
    void clearDirichletRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
346

347
    /// Prints \ref matrix to stdout
348 349
    void print() const;

350
    /// Removes all matrix entries
351 352 353 354
    void clear()
    {
	set_to_zero(matrix);
    }
355

356
    /// Test whether \ref matrix is symmetric. Exits if not.
357 358 359 360
    void test();

    bool symmetric();

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    inline vector<Operator*>& getOperators() 
Thomas Witkowski's avatar
Thomas Witkowski committed
362
    { 
363
      return operators; 
364
    }
365
    
Thomas Witkowski's avatar
Thomas Witkowski committed
366
    inline vector<double*>& getOperatorFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
367
    { 
368
      return operatorFactor; 
369
    }
370

Thomas Witkowski's avatar
Thomas Witkowski committed
371
    inline vector<double*>& getOperatorEstFactor() 
Thomas Witkowski's avatar
Thomas Witkowski committed
372
    { 
373
      return operatorEstFactor; 
374
    }
375

Thomas Witkowski's avatar
Thomas Witkowski committed
376 377
    inline BoundaryManager* getBoundaryManager() const 
    { 
378
      return boundaryManager; 
379
    }
380

Thomas Witkowski's avatar
Thomas Witkowski committed
381 382
    inline void setBoundaryManager(BoundaryManager *bm) 
    {
383
      boundaryManager = bm;
384
    }
385

Thomas Witkowski's avatar
Thomas Witkowski committed
386 387 388 389 390 391 392
    /// 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
393 394
      if (nnzPerRow < 20) 
	nnzPerRow = 20;
Thomas Witkowski's avatar
Thomas Witkowski committed
395 396 397 398 399 400 401 402
    }

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

403
    /// Writes the matrix to an output stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
404
    void serialize(ostream &out);
405

406
    /// Reads a matrix from an input stream.
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    void deserialize(istream &in);
408

409
    ///
410 411
    int memsize();

412
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
413
    void setDofMap(FeSpaceDofMap& m);
414
#endif
415 416

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
417 418
    /// Pointer to a FiniteElemSpace with information about corresponding row DOFs
    /// and basis functions
419
    const FiniteElemSpace *rowFeSpace;
420

Thomas Witkowski's avatar
Thomas Witkowski committed
421 422
    /// Pointer to a FiniteElemSpace with information about corresponding 
    /// column DOFs and basis functions
423
    const FiniteElemSpace *colFeSpace;
424

425
    /// Name of the DOFMatrix
Thomas Witkowski's avatar
Thomas Witkowski committed
426
    string name;
427

Thomas Witkowski's avatar
Thomas Witkowski committed
428 429
    /// Sparse matrix, type is a template parameter by 
    /// default compressed2D<double>
430
    base_matrix_type matrix;
431

Thomas Witkowski's avatar
Thomas Witkowski committed
432 433 434
    /// Pointers to all operators of the equation systems. Are used in the
    /// assembling process.
    vector<Operator*> operators;
435
    
Thomas Witkowski's avatar
Thomas Witkowski committed
436 437 438
    /// Defines for each operator a factor which is used to scal the element
    /// matrix after the assembling process of the operator.
    vector<double*> operatorFactor;
439

440
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
441
    vector<double*> operatorEstFactor;
442

443
    ///
444 445
    BoundaryManager *boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
446 447
    /// 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.
448 449
    bool coupleMatrix;

450
    /// Temporary variable used in assemble()
451 452 453 454 455 456 457 458 459
    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
460
    vector<DegreeOfFreedom> rowIndices;
461 462

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

465 466 467 468
    /// 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
Thomas Witkowski's avatar
Thomas Witkowski committed
469 470
    /// in this set.
    std::set<DegreeOfFreedom> dirichletDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
471

472 473 474 475
    /// 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
476

477
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
478 479 480 481
     /// Is used in parallel computations to check whether specific DOFs in the
     /// row FE spaces are owned by the rank or not. This is used to ensure that 
     /// Dirichlet BC is handled correctly in parallel computations.
     FeSpaceDofMap *dofMap;
482 483
#endif

484
    /// Inserter object: implemented as pointer, allocated and deallocated as needed
485
    inserter_type *inserter;
486
      
487 488 489 490 491 492 493 494 495 496
    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