DOFMatrix.cc 13.4 KB
Newer Older
1
2
#include <algorithm>
#include <png.h>
3
#include <boost/numeric/mtl/mtl.hpp>
4
#include "DOFMatrix.h"
5
6
7
8
9
10
11
12
13
14
15
16
#include "QPsiPhi.h"
#include "BasisFunction.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
#include "DOFVector.h"
#include "Operator.h"
#include "BoundaryCondition.h"
#include "BoundaryManager.h"
#include "Assembler.h"
17
18
#include "Utilities.h"

19
20
namespace AMDiS {

21
22
  using namespace mtl;

23
24
25
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
26
27
    : rowFESpace(NULL),
      colFESpace(NULL),
28
29
30
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
31
32
      inserter(NULL)
  {}
33
34
35

  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowFESpace_,
		       const FiniteElemSpace* colFESpace_,
36
		       std::string name_)
37
38
39
    : rowFESpace(rowFESpace_),
      colFESpace(colFESpace_),
      name(name_), 
40
41
      coupleMatrix(false),
      inserter(NULL)
42
43
44
  {
    TEST_EXIT(rowFESpace)("no rowFESpace\n");
  
Thomas Witkowski's avatar
Thomas Witkowski committed
45
    if (!colFESpace)
46
47
48
      colFESpace = rowFESpace;

    if (rowFESpace && rowFESpace->getAdmin())
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
50

Thomas Witkowski's avatar
Thomas Witkowski committed
51
    boundaryManager = NEW BoundaryManager(rowFESpace_);
52
53
54
55
56
57

    nRow = rowFESpace->getBasisFcts()->getNumber();
    nCol = colFESpace->getBasisFcts()->getNumber();
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59

    applyDBCs.clear();
60

61
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
62
63
    applicationOrdering = NULL;
#endif
64
65
66
  }

  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
67
    : name(rhs.name + "copy")
68
  {
69
    *this = rhs;
70
71
    if (rowFESpace && rowFESpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
72
73
74

    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion");
    inserter= 0;
75
76
77
78
  }

  DOFMatrix::~DOFMatrix()
  {
79
    FUNCNAME("DOFMatrix::~DOFMatrix()");
80
    if (rowFESpace && rowFESpace->getAdmin())
81
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
82

83
    if (boundaryManager) delete boundaryManager;
84
    if (inserter) delete inserter;
85
86
87
88
89
90
  }

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");

91
92
93
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type Matrix;
94

95
96
97
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
98

99
100
101
102
103
104
105
    typedef traits::range_generator<major, Matrix>::type      cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
    
    for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) {
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
	    Msg::print(" (%3d,%3d,%20.17lf)", row(*icursor), col(*icursor), value(*icursor));
	Msg::print("\n");
106
107
108
109
110
111
112
    }
  }

  bool DOFMatrix::symmetric()
  {
    FUNCNAME("DOFMatrix::symmetric()");

113
    double tol = 1e-5;
114

115
116
117
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
118

119
120
121
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
122

123
124
125
126
127
128
129
    typedef traits::range_generator<major, Matrix>::type      cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
    
    for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
	// Compare each non-zero entry with its transposed
	if (abs(value(*icursor) - matrix[col(*icursor)][row(*icursor)]) > tol)
130
131
132
133
134
135
136
137
	  return false;
    return true;
  }

  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

138
    int non_symmetric = !symmetric();
139
140
141
142
143
144
145
146
147
148
149
150
151
152

    if (non_symmetric) {
      MSG("matrix `%s' not symmetric.\n", name.data());
    } else {
      MSG("matrix `%s' is symmetric.\n", name.data());
    }
  }

  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
153
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
156
157
158

    /// The matrix values may only be copyed, if there is no active insertion.
    if (rhs.inserter == 0 && inserter == 0)
      matrix = rhs.matrix;

159
    if (rhs.boundaryManager) {
Thomas Witkowski's avatar
Thomas Witkowski committed
160
      boundaryManager = NEW BoundaryManager(*rhs.boundaryManager);
161
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
      boundaryManager = NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
164

165
166
167
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
168
169
170
171
172

    return *this;
  }

  void DOFMatrix::addElementMatrix(double sign, 
173
				   const ElementMatrix& elMat, 
174
				   const BoundaryType *bound,
175
176
				   ElInfo* rowElInfo,
				   ElInfo* colElInfo,
177
178
				   bool add)
  {
179
    FUNCNAME("DOFMatrix::addElementMatrix()");
180

181
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
182
183
    inserter_type &ins= *inserter;

184
    // === Get indices mapping from local to global matrix indices. ===
185

186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
						   rowFESpace->getAdmin(),
						   &rowIndices);
    if (rowFESpace == colFESpace) {
      colIndices = rowIndices;
    } else {
      if (colElInfo) {
	colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &colIndices);
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
	colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &colIndices);
      }
    }
203

204
205
206
207
208
209
210

    // === Add element matrix to the global matrix using the indices mapping. ===

    DegreeOfFreedom row, col;
    double entry;
    for (int i = 0; i < nRow; i++)  {   // for all rows of element matrix
      row = rowIndices[i];
211

212
213
214
215
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
216
	if (!coupleMatrix)
217
	  applyDBCs.insert(static_cast<int>(row));
218
      } else
219
220
	for (int j = 0; j < nCol; j++) {  // for all columns
	  col = colIndices[j];
221
222
223
224
225
226
	  entry = elMat[i][j];
	  if (add)
	    ins[row][col]+= sign * entry;
	  else
	    ins[row][col]= sign * entry;
	}   
227
228
229
    }
  }

230
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
231
  {
232
    return matrix[a][b];
233
234
  }

235
  void DOFMatrix::freeDOFContent(int index)
236
  {
Thomas Witkowski's avatar
nothing    
Thomas Witkowski committed
237
    if (matrix.nnz() == 0) return;
238

239
240
241
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type Matrix;
242

243
244
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
245

246
247
248
249
250
251
    typedef traits::range_generator<major, Matrix>::type      cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
    
    cursor_type cursor = begin<major>(matrix);
    // Jump directly to corresponding row or column
    cursor+= index;
252

253
254
255
256
257
258
    // Requires structural symmetry !!!
    for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) {
      int my_row= row(*icursor), my_col= col(*icursor);
      // Not very efficient (but general)
      matrix.lvalue(my_row, my_col) = 0.0; // Need to call crop somewhere !!!! Peter
      matrix.lvalue(my_col, my_row) = 0.0;
259
260
261
    }
  }

262
  // Should work as before
263
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
264
  {
265
    FUNCNAME("DOFMatrix::assemble()");
266

267
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
268

269
270
271
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
    for (; it != operators.end(); ++it, ++factorIt) {
272
      if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
275
	(*it)->getElementMatrix(elInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);
276
      }
277
    }        
278

279
    addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);   
280
281
  }

282
283
284
285
286
287
288
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

      TEST_EXIT_DBG(op)("No operator!\n");

289
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
290
      op->getElementMatrix(elInfo, elementMatrix, factor);
291
      addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);
292
293
  }

294
295
296
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
297
298
299
300
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

301
    if (!op && operators.size() == 0)
302
303
      return;

304
305
    set_to_zero(elementMatrix);

306
    if (op) {
307
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
308
			   elementMatrix);
309
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311
312
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
      for (; it != operators.end(); ++it, ++factorIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
313
314
315
316
317
318
319
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

320
    addElementMatrix(factor, elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
323
324
325
326
327
328
329
330
331
332
333
  }

  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
			    ElInfo *smallElInfo, ElInfo *largeElInfo,			    
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble2()");

    if (!op && operators.size() == 0) {
      return;
    }

334
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    
336
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339
340
341
342
343
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
344
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
347
348
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
349
350
351
352
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
353
354
      }      
    }
355

356
    addElementMatrix(factor, elementMatrix, bound, mainElInfo, NULL);       
357
358
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
362
363
364
365
366
367
368
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
	 it != operators.end();
	 ++it) {     
      (*it)->finishAssembling();
    }
  }

369
  // Should work as before
370
371
372
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
373
    std::vector<Operator*>::iterator op;
374
375
376
    for(op = operators.begin(); op != operators.end(); ++op) {
      fillFlag |= (*op)->getFillFlag();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
377

378
379
380
    return fillFlag;
  }

381
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
382
  {
383
    matrix+= a * x.matrix + y.matrix;
384
385
  }

386

387
388
  void DOFMatrix::scal(double b) 
  {
389
    matrix*= b;
390
391
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396
397
398
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

399
400
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
401
    matrix= rhs.matrix;
402
403
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
404
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
405
  {      
406
407
408
409
410
411
    inserter_type &ins= *inserter;
   
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
      ins[*it][*it] = 1.0;

    rows->clear();
412
413
  }

414
415
  void DOFMatrix::createPictureFile(const char* filename, int dim)
  {
416
417
418
    TEST_EXIT(0)("Not yet re-implemented.");

#if 0
419
420
    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 
						  NULL, NULL, NULL);
421
422
423
424
425
426
427
428

    if (!png_ptr)
       return;

    png_bytep rowPointers[dim];
    for (int i = 0; i < dim; i++) {
      rowPointers[i] = (png_byte*)png_malloc(png_ptr, dim);

429
      for (int j = 0; j < dim; j++)
430
431
432
433
434
435
436
437
	rowPointers[i][j] = 255;
    }

    double scalFactor = static_cast<double>(dim) / static_cast<double>(matrix.size());

    for (int i = 0; i < static_cast<int>(matrix.size()); i++) {    
      int pi = static_cast<int>(static_cast<double>(i) * scalFactor);

438
      TEST_EXIT_DBG((pi >= 0) && (pi < dim))("PI");
439
440
441
442
443

      for (int j = 0; j < static_cast<int>(matrix[i].size()); j++) {
	
	int pj = static_cast<int>(static_cast<double>(matrix[i][j].col) * scalFactor);

444
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474

	rowPointers[pi][pj] = 0;
      }
    }
   
    FILE *fp = fopen(filename, "wb");
    TEST_EXIT(fp)("Cannot open file for writing matrix picture file!\n");

    png_infop info_ptr = png_create_info_struct(png_ptr);

    if (!info_ptr) {
       png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
       return;
    }

    png_init_io(png_ptr, fp);

    png_set_IHDR(png_ptr, info_ptr, dim, dim, 8,
		 PNG_COLOR_TYPE_GRAY, 
		 PNG_INTERLACE_NONE,
		 PNG_COMPRESSION_TYPE_DEFAULT,
		 PNG_FILTER_TYPE_DEFAULT);

    png_set_rows(png_ptr, info_ptr, rowPointers);

    png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);

    png_destroy_write_struct(&png_ptr, &info_ptr);

    fclose(fp);
475
#endif
476
477
478
479
480
  }


  int DOFMatrix::memsize() 
  {   
481
482
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
483
484
485
  }

}