DOFMatrix.cc 13.1 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
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
59
60
    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
    }

    return *this;
  }

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

184
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
185
186
187
188
    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
  {
Thomas Witkowski's avatar
nothing    
Thomas Witkowski committed
221
    if (matrix.nnz() == 0) return;
222

223
224
225
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type Matrix;
226

227
228
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
229

230
231
232
233
234
235
    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;
236

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

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

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

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

264
    addElementMatrix(factor, *elementMatrix, bound);   
265
266
  }

267
268
269
270
271
272
273
274
  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
275
      op->getElementMatrix(elInfo, elementMatrix, factor);
276
277
278
      addElementMatrix(factor, *elementMatrix, bound);
  }

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

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

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

346
347
348
    addElementMatrix(factor, *elementMatrix, bound);       
  }

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

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

368
369
370
    return fillFlag;
  }

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

376

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

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

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

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

    rows->clear();
402
403
  }

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

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

    if (!png_ptr)
       return;

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

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

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

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

434
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
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
463
464

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


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

}