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

    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)
  {
152
153
    TEST_EXIT(rhs.inserter == 0 && inserter == 0)("Cannot copy during insertion");

154
155
156
157
158
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
    matrix = rhs.matrix;
159
    coupleMatrix = rhs.coupleMatrix;
160
    if (rhs.boundaryManager) {
Thomas Witkowski's avatar
Thomas Witkowski committed
161
      boundaryManager = NEW BoundaryManager(*rhs.boundaryManager);
162
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
163
164
165
166
167
168
169
      boundaryManager = NULL;
    }
    if (rhs.elementMatrix) {
      elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
					colFESpace->getBasisFcts()->getNumber());
    } else {
      elementMatrix = NULL;
170
171
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)
  {
    FUNCNAME("DOFMatrix::addElementMatrix");

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

    DegreeOfFreedom row, col;
    double entry;
187

188
189
    int n_row = elMat.rowIndices.getSize();
    int n_col = elMat.colIndices.getSize();
190

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

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

      if (condition && condition->isDirichlet()) {
198
	if (!coupleMatrix)
199
	  applyDBCs.insert(static_cast<int>(row));
200
201
202
203
204
205
206
207
208
      } 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;
	}   
209
210
211
    }
  }

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

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

223
224
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
225

226
227
228
229
230
231
    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;
232

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

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

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

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

260
    addElementMatrix(factor, *elementMatrix, bound);   
261
262
  }

263
264
265
266
267
268
269
270
  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
271
      op->getElementMatrix(elInfo, elementMatrix, factor);
272
273
274
      addElementMatrix(factor, *elementMatrix, bound);
  }

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

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

    Operator *operat = op ? op : operators[0];
287
288
    operat->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, rowElInfo, colElInfo);
289
290
    
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
291
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo,
292
			   elementMatrix);
293
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
296
      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
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
	(*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);
    
322
    if(op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
325
326
327
328
329
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
330
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
333
334
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
335
336
337
338
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
339
340
      }      
    }
341

342
343
344
    addElementMatrix(factor, *elementMatrix, bound);       
  }

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

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

364
365
366
    return fillFlag;
  }

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

372

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

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

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

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

    rows->clear();
398
399
  }

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

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

    if (!png_ptr)
       return;

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

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

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

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

430
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
431
432
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

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


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

}