DOFMatrix.cc 12.1 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),
Thomas Witkowski's avatar
Thomas Witkowski committed
30
      nnzPerRow(0),
31
32
      inserter(NULL)
  {}
33

34
35
36
37
38
39
  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowSpace,
		       const FiniteElemSpace* colSpace,
		       std::string n)
    : rowFESpace(rowSpace),
      colFESpace(colSpace),
      name(n), 
40
      coupleMatrix(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
41
      nnzPerRow(0),
42
      inserter(NULL)
43
  {
44
45
46
    FUNCNAME("DOFMatrix::DOFMatrix()");

    TEST_EXIT(rowFESpace)("No fe space for row!\n");
47
  
Thomas Witkowski's avatar
Thomas Witkowski committed
48
    if (!colFESpace)
49
50
51
      colFESpace = rowFESpace;

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

54
    boundaryManager = new BoundaryManager(rowFESpace);
55
56
57
58
59
60

    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
61
62

    applyDBCs.clear();
63
64
65
  }

  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
66
    : name(rhs.name + "copy")
67
  {
68
69
    FUNCNAME("DOFMatrix::DOFMatrix()");

70
    *this = rhs;
71
72
    if (rowFESpace && rowFESpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
73

74
75
    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
    inserter = 0;
76
77
78
79
  }

  DOFMatrix::~DOFMatrix()
  {
80
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
81

82
    if (rowFESpace && rowFESpace->getAdmin())
83
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
86
87
    if (boundaryManager) 
      delete boundaryManager;
    if (inserter) 
      delete inserter;
88
89
90
91
92
  }

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
93
94
95

    if (inserter) 
      inserter->print();
96
97
98
99
100
101
  }

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

102
    double tol = 1e-5;
103

104
105
106
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
107

108
109
110
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
111

112
113
114
115
116
117
118
    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)
119
120
121
122
123
124
125
126
	  return false;
    return true;
  }

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

127
    int non_symmetric = !symmetric();
128

Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (non_symmetric)
130
      MSG("Matrix `%s' not symmetric.\n", name.data());
Thomas Witkowski's avatar
Thomas Witkowski committed
131
    else
132
      MSG("Matrix `%s' is symmetric.\n", name.data());
133
134
135
136
137
138
139
140
  }

  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
141
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
144
145
146

    /// 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
147
    if (rhs.boundaryManager)
Thomas Witkowski's avatar
Thomas Witkowski committed
148
      boundaryManager = new BoundaryManager(*rhs.boundaryManager);
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
150
      boundaryManager = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
151
    
152
153
154
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
155
156
157
158

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
159
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
160
				   const BoundaryType *bound,
161
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
162
				   ElInfo* colElInfo)
163
  {
164
    FUNCNAME("DOFMatrix::addElementMatrix()");
165

166
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
167
168
    inserter_type &ins= *inserter;
 
169
    // === Get indices mapping from local to global matrix indices. ===
170

Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
    rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFESpace->getAdmin(),
						rowIndices);
174
    if (rowFESpace == colFESpace) {
175
      colIndices = rowIndices;
176
177
    } else {
      if (colElInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
	colFESpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFESpace->getAdmin(),
						    colIndices);
181
182
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
185
	colFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFESpace->getAdmin(),
						    colIndices);
186
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    }
188

189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
    for (int i = 0; i < nRow; i++)  {
191
192
      DegreeOfFreedom row = rowIndices[i];

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

      if (condition && condition->isDirichlet()) {
197
198
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
199
200
 	  if (rankDofs[rowIndices[i]]) 
 	    applyDBCs.insert(static_cast<int>(rowIndices[i]));
201
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
202
	  applyDBCs.insert(static_cast<int>(rowIndices[i]));
203
204
#endif
	}
205
206
207
208
209
210
211
      } else {
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
	  double entry = elMat[i][j];
  
	  ins[row][col] += entry;
	}
212
      }
213
    }
214
  }
215

216
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
217
  {
218
    return matrix[a][b];
219
220
  }

221
  void DOFMatrix::freeDOFContent(int index)
222
  {}
223

224
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
225
  {
226
    FUNCNAME("DOFMatrix::assemble()");
227

228
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
229

230
231
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
232
233
234
    for (; it != operators.end(); ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse() == false)
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
    if (factor != 1.0)
      elementMatrix *= factor;
238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
    addElementMatrix(elementMatrix, bound, elInfo, NULL);   
240
241
  }

242
243
244
245
246
247
248
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

249
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
250
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
253
254
255

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
256
257
  }

258
259
260
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
261
262
263
264
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

265
    if (!op && operators.size() == 0)
266
267
      return;

268
269
    set_to_zero(elementMatrix);

270
    if (op) {
271
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
272
			   elementMatrix);
273
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
      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
277
278
279
280
281
282
283
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
284
285
286
287
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
291
292
293
294
295
296
  }

  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
297
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
      return;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
324
325
    if (factor != 1.0)
      elementMatrix *= factor;

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

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

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

344
345
346
    return fillFlag;
  }

347
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
348
  {
349
    matrix+= a * x.matrix + y.matrix;
350
351
  }

352

353
354
  void DOFMatrix::scal(double b) 
  {
355
    matrix*= b;
356
357
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
360
361
362
363
364
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

365
366
367
368
369
370
371
372
373
374
375
376
377
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
  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;
    }    
  }

423
424
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
425
    matrix = rhs.matrix;
426
427
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
428
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
429
  {      
430
431
432
    inserter_type &ins= *inserter;
   
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
433
434
435
      ins[*it][*it] = 1.0;

    rows->clear();
436
437
  }

438
439
  int DOFMatrix::memsize() 
  {   
440
441
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
442
443
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446
447
448
449
450
451
452
453
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

454
}