DOFMatrix.cc 13 KB
Newer Older
1
2
#include <algorithm>
#include <png.h>
3
4

#include "DOFMatrix.h"
5
6
7
8
9
10
11
12
13
14
15
16
17
#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 "ElementMatrix.h"
#include "Assembler.h"
18
19
20
#include "Utilities.h"

#include <boost/numeric/mtl/mtl.hpp>
21
22
23

namespace AMDiS {

24
25
  using namespace mtl;

26
27
28
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
29
30
31
32
33
    : rowFESpace(NULL),
      colFESpace(NULL),
      elementMatrix(NULL),
      inserter(NULL)
  {}
34
35
36

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
52
    boundaryManager = NEW BoundaryManager(rowFESpace_);
Thomas Witkowski's avatar
Thomas Witkowski committed
53
54
    elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
				      colFESpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
55
56

    applyDBCs.clear();
57
58
59
60

#ifdef HAVE_PARALLEL_AMDIS
    applicationOrdering = NULL;
#endif
61
62
63
  }

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

    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion");
    inserter= 0;
72
73
74
75
  }

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

81
82
83
    if (boundaryManager)
      DELETE boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
    if (elementMatrix)
      DELETE elementMatrix;
86
    if (inserter) delete inserter;
87
88
89
90
91
92
  }

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

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

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

101
102
103
104
105
106
107
    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");
108
109
110
111
112
113
114
    }
  }

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

115
    double tol = 1e-5;
116

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

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

125
126
127
128
129
130
131
    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)
132
133
134
135
136
137
138
139
	  return false;
    return true;
  }

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

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

    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;
155
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
158
159
160

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
170
171
    if (rhs.elementMatrix) {
      elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
					colFESpace->getBasisFcts()->getNumber());
    } else {
      elementMatrix = NULL;
172
173
174
175
176
177
178
179
180
181
182
183
    }

    return *this;
  }

  void DOFMatrix::addElementMatrix(double sign, 
				   const ElementMatrix &elMat, 
				   const BoundaryType *bound,
				   bool add)
  {
    FUNCNAME("DOFMatrix::addElementMatrix");

184
185
186
187
188
    TEST_EXIT(inserter)("DOFMatrix is not in insertion mode");
    inserter_type &ins= *inserter;

    DegreeOfFreedom row, col;
    double entry;
189

190
191
    int n_row = elMat.rowIndices.getSize();
    int n_col = elMat.colIndices.getSize();
192

193
    for (int i = 0; i < n_row; i++)  {   // for all rows of element matrix
194
      row = elMat.rowIndices[i];
195

196
197
198
199
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
200
	if (!coupleMatrix)
201
	  applyDBCs.insert(static_cast<int>(row));
202
203
204
205
206
207
208
209
210
      } else
	for (int j = 0; j < n_col; j++) {  // for all columns
	  col = elMat.colIndices[j];
	  entry = elMat[i][j];
	  if (add)
	    ins[row][col]+= sign * entry;
	  else
	    ins[row][col]= sign * entry;
	}   
211
212
213
    }
  }

214
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
215
  {
216
    return matrix[a][b];
217
218
  }

219
  void DOFMatrix::freeDOFContent(int index)
220
  {
221
222
223
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type Matrix;
224

225
226
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
227

228
229
230
231
232
233
    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;
234

235
236
237
238
239
240
    // 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;
241
242
243
    }
  }

244
  // Should work as before
245
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
246
  {
247
    FUNCNAME("DOFMatrix::assemble()");
248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
    operators[0]->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
251

252
253
254
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
    for (; it != operators.end(); ++it, ++factorIt) {
255
      if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
256
	  (*it)->getElementMatrix(elInfo,
257
258
259
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
      }
260
    }        
261

262
    addElementMatrix(factor, *elementMatrix, bound);   
263
264
  }

265
266
267
268
269
270
271
272
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

      op->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
273
      op->getElementMatrix(elInfo, elementMatrix, factor);
274
275
276
      addElementMatrix(factor, *elementMatrix, bound);
  }

277
278
279
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
280
281
282
283
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

284
    if(!op && operators.size() == 0) {
285
286
287
288
      return;
    }

    Operator *operat = op ? op : operators[0];
289
290
    operat->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, rowElInfo, colElInfo);
291
292
    
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
293
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo,
294
			   elementMatrix);
295
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
      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
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
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

    addElementMatrix(factor, *elementMatrix, bound);       
  }

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

    Operator *operat = op ? op : operators[0];
    operat->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, mainElInfo);
    
324
    if(op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
328
329
330
331
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
332
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
335
336
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
337
338
339
340
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
341
342
      }      
    }
343

344
345
346
    addElementMatrix(factor, *elementMatrix, bound);       
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
350
351
352
353
354
355
356
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
	 it != operators.end();
	 ++it) {     
      (*it)->finishAssembling();
    }
  }

357
  // Should work as before
358
359
360
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
361
    std::vector<Operator*>::iterator op;
362
363
364
    for(op = operators.begin(); op != operators.end(); ++op) {
      fillFlag |= (*op)->getFillFlag();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
365

366
367
368
    return fillFlag;
  }

369
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
370
  {
371
    matrix+= a * x.matrix + y.matrix;
372
373
  }

374

375
376
  void DOFMatrix::scal(double b) 
  {
377
    matrix*= b;
378
379
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
382
383
384
385
386
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

387
388
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
389
    matrix= rhs.matrix;
390
391
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
392
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
393
  {      
394
395
396
397
398
399
    inserter_type &ins= *inserter;
   
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
      ins[*it][*it] = 1.0;

    rows->clear();
400
401
  }

402
403
  void DOFMatrix::createPictureFile(const char* filename, int dim)
  {
404
405
406
    TEST_EXIT(0)("Not yet re-implemented.");

#if 0
407
408
    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 
						  NULL, NULL, NULL);
409
410
411
412
413
414
415
416

    if (!png_ptr)
       return;

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

417
      for (int j = 0; j < dim; j++)
418
419
420
421
422
423
424
425
	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);

426
      TEST_EXIT_DBG((pi >= 0) && (pi < dim))("PI");
427
428
429
430
431

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

432
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462

	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);
463
#endif
464
465
466
467
468
  }


  int DOFMatrix::memsize() 
  {   
469
470
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
471
472
473
  }

}