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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

265
266
    set_to_zero(elementMatrix);

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

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

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

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

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

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

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

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

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

341
342
343
    return fillFlag;
  }

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

349

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

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

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

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

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

    rows->clear();
433
434
  }

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

}