DOFMatrix.cc 12.3 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
    typedef traits::range_generator<major, Matrix>::type      cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
101
102
103

    std::cout.precision(10);
    for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
104
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
105
	  std::cout << row(*icursor) << " " << col(*icursor) << " " << value(*icursor) << "\n";
106
107
108
109
110
111
  }

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

112
    double tol = 1e-5;
113

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

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

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

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

137
    int non_symmetric = !symmetric();
138
139
140
141
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)
  {
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
152
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
153
154
155
156
157

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

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

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

    return *this;
  }

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

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

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

185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    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);
      }
    }
202

203
204
205
206
207
208
    // === 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];
209

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

      if (condition && condition->isDirichlet()) {
214
	if (!coupleMatrix)
215
	  applyDBCs.insert(static_cast<int>(row));
216
      } else
217
218
	for (int j = 0; j < nCol; j++) {  // for all columns
	  col = colIndices[j];
219
	  entry = elMat[i][j];
220

221
222
223
224
225
	  if (add)
	    ins[row][col]+= sign * entry;
	  else
	    ins[row][col]= sign * entry;
	}   
226
227
228
    }
  }

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

234
  void DOFMatrix::freeDOFContent(int index)
235
  {}
236

237
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
238
  {
239
    FUNCNAME("DOFMatrix::assemble()");
240

241
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
242

243
244
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
245
246
247
    for (; it != operators.end(); ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse() == false)
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
248

249
    addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);   
250
251
  }

252
253
254
255
256
257
258
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

259
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
260
      op->getElementMatrix(elInfo, elementMatrix, factor);
261
      addElementMatrix(factor, elementMatrix, bound, elInfo, NULL);
262
263
  }

264
265
266
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
267
268
269
270
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

271
    if (!op && operators.size() == 0)
272
273
      return;

274
275
    set_to_zero(elementMatrix);

276
    if (op) {
277
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
278
			   elementMatrix);
279
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
      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
283
284
285
286
287
288
289
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

290
    addElementMatrix(factor, elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
291
292
293
294
295
296
297
298
299
300
301
302
303
  }

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

304
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
305
    
306
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
309
310
311
312
313
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
314
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
317
318
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
319
320
321
322
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
323
324
      }      
    }
325

326
    addElementMatrix(factor, elementMatrix, bound, mainElInfo, NULL);       
327
328
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
332
333
334
335
336
337
338
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
	 it != operators.end();
	 ++it) {     
      (*it)->finishAssembling();
    }
  }

339
  // Should work as before
340
341
342
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
343
    std::vector<Operator*>::iterator op;
344
345
346
    for(op = operators.begin(); op != operators.end(); ++op) {
      fillFlag |= (*op)->getFillFlag();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
347

348
349
350
    return fillFlag;
  }

351
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
352
  {
353
    matrix+= a * x.matrix + y.matrix;
354
355
  }

356

357
358
  void DOFMatrix::scal(double b) 
  {
359
    matrix*= b;
360
361
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
364
365
366
367
368
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

369
370
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
371
    matrix= rhs.matrix;
372
373
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
374
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
375
  {      
376
377
378
379
380
381
    inserter_type &ins= *inserter;
   
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
      ins[*it][*it] = 1.0;

    rows->clear();
382
383
  }

384
385
  void DOFMatrix::createPictureFile(const char* filename, int dim)
  {
386
387
388
    TEST_EXIT(0)("Not yet re-implemented.");

#if 0
389
390
    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 
						  NULL, NULL, NULL);
391
392
393
394
395
396
397
398

    if (!png_ptr)
       return;

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

399
      for (int j = 0; j < dim; j++)
400
401
402
403
404
405
406
407
	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);

408
      TEST_EXIT_DBG((pi >= 0) && (pi < dim))("PI");
409
410
411
412
413

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

414
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444

	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);
445
#endif
446
447
448
449
450
  }


  int DOFMatrix::memsize() 
  {   
451
452
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
453
454
455
  }

}