DOFMatrix.cc 12.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include <algorithm>
14
#include <boost/numeric/mtl/mtl.hpp>
15
#include "DOFMatrix.h"
16
17
18
19
20
21
22
23
24
25
26
27
#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"
28
#include "Serializer.h"
29

30
31
namespace AMDiS {

32
33
  using namespace mtl;

34
35
36
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
37
38
    : rowFeSpace(NULL),
      colFeSpace(NULL),
39
40
41
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      nnzPerRow(0),
43
44
      inserter(NULL)
  {}
45

Thomas Witkowski's avatar
Thomas Witkowski committed
46

47
48
49
  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowSpace,
		       const FiniteElemSpace* colSpace,
		       std::string n)
50
51
    : rowFeSpace(rowSpace),
      colFeSpace(colSpace),
52
      name(n), 
53
      coupleMatrix(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
54
      nnzPerRow(0),
55
      inserter(NULL)
56
  {
57
58
    FUNCNAME("DOFMatrix::DOFMatrix()");

59
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
60
  
61
62
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
63

64
65
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->addDOFIndexed(this);
66

67
    boundaryManager = new BoundaryManager(rowFeSpace);
68

69
70
    nRow = rowFeSpace->getBasisFcts()->getNumber();
    nCol = colFeSpace->getBasisFcts()->getNumber();
71
72
73
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
74
75

    applyDBCs.clear();
76
77
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
78

79
  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
80
    : name(rhs.name + "copy")
81
  {
82
83
    FUNCNAME("DOFMatrix::DOFMatrix()");

84
    *this = rhs;
85
86
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFeSpace->getAdmin()))->addDOFIndexed(this);
87

88
89
    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
    inserter = 0;
90
91
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
92

93
94
  DOFMatrix::~DOFMatrix()
  {
95
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
96

97
98
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
99
100
101
102
    if (boundaryManager) 
      delete boundaryManager;
    if (inserter) 
      delete inserter;
103
104
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
105

106
107
108
  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
109
110
111

    if (inserter) 
      inserter->print();
112
113
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
114

115
116
117
118
  bool DOFMatrix::symmetric()
  {
    FUNCNAME("DOFMatrix::symmetric()");

119
    double tol = 1e-5;
120

121
122
123
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
124

125
126
127
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
128

129
130
131
132
133
134
135
    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)
136
137
138
139
	  return false;
    return true;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
140

141
142
143
144
  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

145
    int non_symmetric = !symmetric();
146

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (non_symmetric)
148
      MSG("Matrix `%s' not symmetric.\n", name.data());
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    else
150
      MSG("Matrix `%s' is symmetric.\n", name.data());
151
152
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
153

154
155
  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
156
157
    rowFeSpace = rhs.rowFeSpace;
    colFeSpace = rhs.colFeSpace;
158
159
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
160
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
161
162
163
164
165

    /// 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
166
    if (rhs.boundaryManager)
Thomas Witkowski's avatar
Thomas Witkowski committed
167
      boundaryManager = new BoundaryManager(*rhs.boundaryManager);
Thomas Witkowski's avatar
Thomas Witkowski committed
168
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
169
      boundaryManager = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
170
    
171
172
173
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
174
175
176
177

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
178

Thomas Witkowski's avatar
Thomas Witkowski committed
179
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
180
				   const BoundaryType *bound,
181
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
182
				   ElInfo* colElInfo)
183
  {
184
    FUNCNAME("DOFMatrix::addElementMatrix()");
185

186
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
187
188
    inserter_type &ins= *inserter;
 
189
    // === Get indices mapping from local to global matrix indices. ===
190

191
192
    rowFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFeSpace->getAdmin(),
193
						rowIndices);
194
    if (rowFeSpace == colFeSpace) {
195
      colIndices = rowIndices;
196
197
    } else {
      if (colElInfo) {
198
199
	colFeSpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
200
						    colIndices);
201
202
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
203
204
	colFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
205
						    colIndices);
206
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
207
    }
208

209
210
    using namespace mtl;

211
212
213
214
215
216
217
218
219
220
221
222
223
#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
224
    for (int i = 0; i < nRow; i++)  {
225
      DegreeOfFreedom row = rowIndices[i];
226

227
228
229
230
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
231
232
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
233
234
	  if ((*rankDofs)[rowIndices[i]]) 
	    applyDBCs.insert(static_cast<int>(row));
235
#else
236
	  applyDBCs.insert(static_cast<int>(row));
237
238
#endif
	}
239
      } else {
240
241
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
242
	  ins[row][col] += elMat[i][j];
243
	}
244
      }
245
    }
246
  }
247

248

249
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
250
  {
251
    return matrix[a][b];
252
253
  }

254

255
  void DOFMatrix::freeDOFContent(int index)
256
  {}
257

258

259
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
260
  {
261
    FUNCNAME("DOFMatrix::assemble()");
262

263
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
264

265
266
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
267
    for (; it != operators.end(); ++it, ++factorIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      if ((*it)->getNeedDualTraverse() == false)
269
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);      
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
    if (factor != 1.0)
      elementMatrix *= factor;
273

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

277

278
279
280
281
282
283
284
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

285
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
286
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
289
290
291

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
292
293
  }

294

295
296
297
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
298
299
300
301
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

302
    if (!op && operators.size() == 0)
303
304
      return;

305
306
    set_to_zero(elementMatrix);

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

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

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

329

Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
332
333
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
337
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
      return;

340
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
341
    
342
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
346
347
348
349
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
350
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
354
355
356
357
358
	  (*it)->getElementMatrix(mainElInfo, 
				  auxElInfo,
				  smallElInfo, 
				  largeElInfo,
				  rowFeSpace == colFeSpace,
359
360
361
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
362
363
      }      
    }
364

Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
368
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
369
370
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
371

Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374
375
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
376
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
380

381
  // Should work as before
382
383
384
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
385
386
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
387
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
388

389
390
391
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
392

393
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
394
  {
395
    matrix+= a * x.matrix + y.matrix;
396
397
  }

398

399
400
  void DOFMatrix::scal(double b) 
  {
401
    matrix*= b;
402
403
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
409
410
411
  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
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
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
  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;
    }    
  }

471

472
473
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
474
    matrix = rhs.matrix;
475
476
  }

477

Thomas Witkowski's avatar
Thomas Witkowski committed
478
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
479
  {      
480
481
482
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
483
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
484
      ins[*it][*it] = 1.0;    
485
486

    rows->clear();
487
488
  }

489

490
491
  int DOFMatrix::memsize() 
  {   
492
493
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
494
495
  }

496

Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
499
500
501
502
503
504
505
506
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

507
}