DOFMatrix.cc 11.9 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
35
36
37
38
  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowSpace,
		       const FiniteElemSpace* colSpace,
		       std::string n)
    : rowFESpace(rowSpace),
      colFESpace(colSpace),
      name(n), 
39
40
      coupleMatrix(false),
      inserter(NULL)
41
  {
42
43
44
    FUNCNAME("DOFMatrix::DOFMatrix()");

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

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

52
    boundaryManager = new BoundaryManager(rowFESpace);
53
54
55
56
57
58

    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
59
60

    applyDBCs.clear();
61
62
63
  }

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

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

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

  DOFMatrix::~DOFMatrix()
  {
78
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
79

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

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
91
92
93

    if (inserter) 
      inserter->print();
94
95
96
97
98
99
  }

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

100
    double tol = 1e-5;
101

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

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

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

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

125
    int non_symmetric = !symmetric();
126

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

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

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

    return *this;
  }

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

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

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

187

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

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

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

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

219
  void DOFMatrix::freeDOFContent(int index)
220
  {}
221

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

226
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
227

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

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

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

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

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

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

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
254
255
  }

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

263
    if (!op && operators.size() == 0)
264
265
      return;

266
267
    set_to_zero(elementMatrix);

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

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

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

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

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

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

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
324
325
  }

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

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

342
343
344
    return fillFlag;
  }

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

350

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

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

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

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

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

    rows->clear();
434
435
  }

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

}