DOFMatrix.cc 12.6 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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
    using namespace mtl;

#if 0
    std::cout << "----- PRINT MAT --------" << std::endl;
    std::cout << elMat << std::endl;
    std::cout << "rows: ";
    for (int i = 0; i < rowIndices.size(); i++)
      std::cout << rowIndices[i] << " ";
    std::cout << std::endl;
    std::cout << "cols: ";
    for (int i = 0; i < colIndices.size(); i++)
      std::cout << colIndices[i] << " ";
    std::cout << std::endl;
#endif
         
Thomas Witkowski's avatar
Thomas Witkowski committed
204
    for (int i = 0; i < nRow; i++)  {
205
206
      DegreeOfFreedom row = rowIndices[i];

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

      if (condition && condition->isDirichlet()) {
211
212
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
213
 	  if ((*rankDofs)[rowIndices[i]]) 
214
 	    applyDBCs.insert(static_cast<int>(row));
215
#else
216
	  applyDBCs.insert(static_cast<int>(row));
217
218
#endif
	}
219
220
221
222
223
224
225
      } else {
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
	  double entry = elMat[i][j];
  
	  ins[row][col] += entry;
	}
226
      }
227
    }
228
  }
229

230

231
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
232
  {
233
    return matrix[a][b];
234
235
  }

236

237
  void DOFMatrix::freeDOFContent(int index)
238
  {}
239

240

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

245
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
246

247
248
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
249
250
251
    for (; it != operators.end(); ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse() == false)
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
252

Thomas Witkowski's avatar
Thomas Witkowski committed
253
254
    if (factor != 1.0)
      elementMatrix *= factor;
255

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

259

260
261
262
263
264
265
266
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

267
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
271
272
273

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
274
275
  }

276

277
278
279
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
280
281
282
283
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

284
    if (!op && operators.size() == 0)
285
286
      return;

287
288
    set_to_zero(elementMatrix);

289
    if (op) {
290
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
291
			   elementMatrix);
292
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
295
      for (; it != operators.end(); ++it, ++factorIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
299
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
305
    if (factor != 1.0)
      elementMatrix *= factor;

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

308

Thomas Witkowski's avatar
Thomas Witkowski committed
309
310
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
311
312
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
313
314
315
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
316
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
      return;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
344
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
345
346
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
350
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
354
      (*it)->finishAssembling();
  }

355
  // Should work as before
356
357
358
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
359
360
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
361
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
362

363
364
365
    return fillFlag;
  }

366
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
367
  {
368
    matrix+= a * x.matrix + y.matrix;
369
370
  }

371

372
373
  void DOFMatrix::scal(double b) 
  {
374
    matrix*= b;
375
376
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
380
381
382
383
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

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
436
437
438
439
440
441
  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;
    }    
  }

442

443
444
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
445
    matrix = rhs.matrix;
446
447
  }

448

Thomas Witkowski's avatar
Thomas Witkowski committed
449
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
450
  {      
451
452
453
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
454
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
455
      ins[*it][*it] = 1.0;    
456
457

    rows->clear();
458
459
  }

460

461
462
  int DOFMatrix::memsize() 
  {   
463
464
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
465
466
  }

467

Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
470
471
472
473
474
475
476
477
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

478
}