DOFMatrix.cc 13.4 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");
  
46
    if (!colFESpace) {
47
48
49
50
      colFESpace = rowFESpace;
    }

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

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

    applyDBCs.clear();
58
59
60
61

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

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

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

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

82
83
84
    if (boundaryManager)
      DELETE boundaryManager;

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

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

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

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

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

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

116
    double tol = 1e-5;
117

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

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

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

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

141
    int non_symmetric = !symmetric();
142
143
144
145
146
147
148
149
150

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


151
  // ok
152
153
  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
154
155
    TEST_EXIT(rhs.inserter == 0 && inserter == 0)("Cannot copy during insertion");

156
157
158
159
160
161
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
    matrix = rhs.matrix;
    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
166
167
168
169
170
      boundaryManager = NULL;
    }
    if (rhs.elementMatrix) {
      elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
					colFESpace->getBasisFcts()->getNumber());
    } else {
      elementMatrix = NULL;
171
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
201
202
203
204
205
206
207
208
209
210
	if (!coupleMatrix)
	  ins[row][row]= 1.0;
      } 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
215
    }
  }

  double DOFMatrix::logAcc(DegreeOfFreedom a,DegreeOfFreedom b) const
  {
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
369
370
371
372
    return fillFlag;
  }

  void DOFMatrix::axpy(double a, 
		       const DOFMatrix& x,
		       const DOFMatrix& y)
  {
373
374
375
    FUNCNAME("DOFMatrix::axpy");
    TEST_EXIT(x.getRowFESpace() == y.getRowFESpace() &&
	      rowFESpace == x.getRowFESpace())
376
      ("row fe-spaces not equal\n");
377
378
    TEST_EXIT(x.getColFESpace() == y.getColFESpace() &&
	      colFESpace == x.getColFESpace())
379
      ("col fe-spaces not equal\n");
380
381
382

    matrix+= a * x.matrix + y.matrix;

383
384
  }

385

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
403
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
404
405
  {      
      ERROR_EXIT("TODO: removeRowsWithDBC");
Thomas Witkowski's avatar
Thomas Witkowski committed
406

407
408
409
410
411
412
413
414
415
416
417
//     for (std::set<int>::iterator it = rows->begin();
// 	 it != rows->end(); 
// 	 ++it) {
//       if (coupleMatrix) {
// 	matrix[*it].resize(0);
//       } else {
// 	matrix[*it].resize(1);
// 	matrix[*it][0].col = *it;
// 	matrix[*it][0].entry = 1.0;
//       }
//     }
Thomas Witkowski's avatar
Thomas Witkowski committed
418

419
//     rows->clear();
420
421
  }

422
423
  void DOFMatrix::createPictureFile(const char* filename, int dim)
  {
424
425
426
    TEST_EXIT(0)("Not yet re-implemented.");

#if 0
427
428
    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 
						  NULL, NULL, NULL);
429
430
431
432
433
434
435
436

    if (!png_ptr)
       return;

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

437
      for (int j = 0; j < dim; j++)
438
439
440
441
442
443
444
445
	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);

446
      TEST_EXIT_DBG((pi >= 0) && (pi < dim))("PI");
447
448
449
450
451

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

452
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482

	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);
483
#endif
484
485
486
487
488
  }


  int DOFMatrix::memsize() 
  {   
489
490
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
491
492
493
  }

}