DOFMatrix.cc 12.8 KB
Newer Older
1
#include <algorithm>
2
#include <boost/numeric/mtl/mtl.hpp>
3
#include "DOFMatrix.h"
4
5
6
7
8
9
10
11
12
13
14
15
#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"
16
#include "Serializer.h"
17

18
19
namespace AMDiS {

20
21
  using namespace mtl;

22
23
24
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
25
26
    : rowFESpace(NULL),
      colFESpace(NULL),
27
28
29
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
30
31
      inserter(NULL)
  {}
32
33
34

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
50
    boundaryManager = new BoundaryManager(rowFESpace_);
51
52
53
54
55
56

    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
57
58

    applyDBCs.clear();
59
60
61
  }

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

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

  DOFMatrix::~DOFMatrix()
  {
74
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
75

76
    if (rowFESpace && rowFESpace->getAdmin())
77
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
78
79
80
81
    if (boundaryManager) 
      delete boundaryManager;
    if (inserter) 
      delete inserter;
82
83
84
85
86
  }

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
87
88
89
90

    if (inserter) 
      inserter->print();

91
  /*  using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
92
93
    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

    std::cout.precision(10);
103
104
    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)
105
	if (value(*icursor) != 0.0)
Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
106
107
	    std::cout << "(" << row(*icursor) << "," << col(*icursor) << "," << value(*icursor) << ") ";
      
108
      std::cout << "\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

Thomas Witkowski's avatar
Thomas Witkowski committed
143
    if (non_symmetric)
144
      MSG("matrix `%s' not symmetric.\n", name.data());
Thomas Witkowski's avatar
Thomas Witkowski committed
145
    else
146
147
148
149
150
151
152
153
154
      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;

Thomas Witkowski's avatar
Thomas Witkowski committed
161
    if (rhs.boundaryManager)
Thomas Witkowski's avatar
Thomas Witkowski committed
162
      boundaryManager = new BoundaryManager(*rhs.boundaryManager);
Thomas Witkowski's avatar
Thomas Witkowski committed
163
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
164
      boundaryManager = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
165
    
166
167
168
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
169
170
171
172

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
173
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
174
				   const BoundaryType *bound,
175
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
176
				   ElInfo* colElInfo)
177
  {
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

Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
187
    rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFESpace->getAdmin(),
						rowIndices);
188
    if (rowFESpace == colFESpace) {
189
      colIndices = rowIndices;
190
191
    } else {
      if (colElInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
192
193
194
	colFESpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFESpace->getAdmin(),
						    colIndices);
195
196
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
197
198
199
	colFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFESpace->getAdmin(),
						    colIndices);
200
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
201
    }
202

Thomas Witkowski's avatar
Thomas Witkowski committed
203
    for (int i = 0; i < nRow; i++)  {
204
205
      DegreeOfFreedom row = rowIndices[i];

206
207
208
209
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
210
211
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
212
213
	  if (rankDofs[rowIndices[i]]) 
	    applyDBCs.insert(static_cast<int>(rowIndices[i]));
214
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
215
	  applyDBCs.insert(static_cast<int>(rowIndices[i]));
216
217
#endif
	}
218
219
220
221
222
223
224
      } else {
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
	  double entry = elMat[i][j];
  
	  ins[row][col] += entry;
	}
225
      }
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

Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
    if (factor != 1.0)
      elementMatrix *= factor;
251

Thomas Witkowski's avatar
Thomas Witkowski committed
252
    addElementMatrix(elementMatrix, bound, elInfo, NULL);   
253
254
  }

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

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

262
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
263
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
266
267
268

      if (factor != 1.0)
	elementMatrix *= factor;

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
269
270
  }

271
272
273
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
274
275
276
277
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

278
    if (!op && operators.size() == 0)
279
280
      return;

281
282
    set_to_zero(elementMatrix);

283
    if (op) {
284
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
285
			   elementMatrix);
286
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
289
      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
290
291
292
293
294
295
296
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
299
300
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
303
304
305
306
307
308
309
  }

  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
			    ElInfo *smallElInfo, ElInfo *largeElInfo,			    
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
310
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
311
312
      return;

313
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
314
    
315
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
318
319
320
321
322
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
323
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
326
327
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
328
329
330
331
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
332
333
      }      
    }
334

Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
338
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
339
340
  }

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

349
  // Should work as before
350
351
352
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
353
354
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
355
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
356

357
358
359
    return fillFlag;
  }

360
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
361
  {
362
    matrix+= a * x.matrix + y.matrix;
363
364
  }

365

366
367
  void DOFMatrix::scal(double b) 
  {
368
    matrix*= b;
369
370
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373
374
375
376
377
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
  void DOFMatrix::serialize(std::ostream &out)
  {
    using namespace mtl; 
    
    typedef traits::range_generator<tag::major, base_matrix_type>::type c_type;
    typedef traits::range_generator<tag::nz, c_type>::type ic_type;
    
    typedef Collection<base_matrix_type>::size_type size_type;
    typedef Collection<base_matrix_type>::value_type value_type;
    
    traits::row<base_matrix_type>::type row(matrix); 
    traits::col<base_matrix_type>::type col(matrix);
    traits::const_value<base_matrix_type>::type value(matrix); 
    
    size_type rows= num_rows(matrix), cols= num_cols(matrix), total= matrix.nnz();
    SerUtil::serialize(out, rows);
    SerUtil::serialize(out, cols);
    SerUtil::serialize(out, total);
    
    for (c_type cursor(mtl::begin<tag::major>(matrix)), 
	   cend(mtl::end<tag::major>(matrix)); cursor != cend; ++cursor)
      for (ic_type icursor(mtl::begin<tag::nz>(cursor)), 
	     icend(mtl::end<tag::nz>(cursor)); icursor != icend; ++icursor) {
	size_type   my_row= row(*icursor), my_col= col(*icursor);
	value_type  my_value= value(*icursor);
	SerUtil::serialize(out, my_row);
	SerUtil::serialize(out, my_col);
	SerUtil::serialize(out, my_value);
      }
  }

  void DOFMatrix::deserialize(std::istream &in)
  {
    using namespace mtl;
    
    typedef Collection<base_matrix_type>::size_type size_type;
    typedef Collection<base_matrix_type>::value_type value_type;
    
    size_type rows, cols, total;
    SerUtil::deserialize(in, rows);
    SerUtil::deserialize(in, cols);
    SerUtil::deserialize(in, total);
    
    // Prepare matrix insertion
    clear();
    // matrix.change_dim(rows, cols) // do we want this?
    inserter_type ins(matrix);
    
    for (size_type i = 0; i < total; ++i) {
      size_type   my_row, my_col;
      value_type  my_value;
      SerUtil::deserialize(in, my_row);
      SerUtil::deserialize(in, my_col);
      SerUtil::deserialize(in, my_value);
      ins(my_row, my_col) << my_value;
    }    
  }

436
437
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
438
    matrix = rhs.matrix;
439
440
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
441
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
442
  {      
443
444
445
    inserter_type &ins= *inserter;
   
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
446
447
448
      ins[*it][*it] = 1.0;

    rows->clear();
449
450
  }

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

}