DOFMatrix.cc 12.3 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

171
172
173
174
175
176
    std::vector<DegreeOfFreedom> &rowIndices2 = rowElInfo->getLocalIndices(rowFESpace);
    std::vector<DegreeOfFreedom> &colIndices2 = rowIndices2;

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

192
193
    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
194
    for (int i = 0; i < nRow; i++)  {
195
      DegreeOfFreedom row = rowIndices2[i];
196

197
198
199
200
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
201
202
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
203
 	  if ((*rankDofs)[row]) 
204
 	    applyDBCs.insert(static_cast<int>(row));
205
#else
206
	  applyDBCs.insert(static_cast<int>(row));
207
208
#endif
	}
209
      } else {
210
211
	for (int j = 0; j < nCol; j++)
	  ins[row][colIndices2[j]] += elMat[i][j];	
212
      }
213
    }
214
  }
215

216

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

222

223
  void DOFMatrix::freeDOFContent(int index)
224
  {}
225

226

227
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
228
  {
229
    FUNCNAME("DOFMatrix::assemble()");
230

231
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
232

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

Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
    if (factor != 1.0)
      elementMatrix *= factor;
241

242
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
243
244
  }

245

246
247
248
249
250
251
252
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

253
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
254
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
255
256
257
258
259

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
260
261
  }

262

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

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

273
274
    set_to_zero(elementMatrix);

275
    if (op) {
276
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
277
			   elementMatrix);
278
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
281
      for (; it != operators.end(); ++it, ++factorIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
284
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
285
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
291
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
  }

294

Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
297
298
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
302
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
      return;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
329
330
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
331
332
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
333

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

Thomas Witkowski's avatar
Thomas Witkowski committed
342

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

351
352
353
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
354

355
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
356
  {
357
    matrix+= a * x.matrix + y.matrix;
358
359
  }

360

361
362
  void DOFMatrix::scal(double b) 
  {
363
    matrix*= b;
364
365
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
369
370
371
372
373
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
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
423
424
425
426
427
428
429
430
431
432
  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;
    }    
  }

433

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

439

Thomas Witkowski's avatar
Thomas Witkowski committed
440
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
441
  {      
442
443
444
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
445
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
446
      ins[*it][*it] = 1.0;    
447
448

    rows->clear();
449
450
  }

451

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

458

Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
461
462
463
464
465
466
467
468
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

469
}