DOFMatrix.cc 13.3 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
29
30
31
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/ParallelDofMapping.h"
#endif
32

33
34
namespace AMDiS {

35
36
  using namespace mtl;

37
38
39
  DOFMatrix *DOFMatrix::traversePtr = NULL;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
49

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

62
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
63

64
65
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
66

67
68
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->addDOFIndexed(this);
69

70
    boundaryManager = new BoundaryManager(rowFeSpace);
71

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

    applyDBCs.clear();
79
80
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
81

82
  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    : name(rhs.name + "copy")
84
  {
85
86
    FUNCNAME("DOFMatrix::DOFMatrix()");

87
    *this = rhs;
88
89
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFeSpace->getAdmin()))->addDOFIndexed(this);
90

91
92
    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
    inserter = 0;
93
94
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
95

96
97
  DOFMatrix::~DOFMatrix()
  {
98
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
99

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

Thomas Witkowski's avatar
Thomas Witkowski committed
108

109
110
111
  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
112
113
114

    if (inserter) 
      inserter->print();
115
116
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
117

118
119
120
121
  bool DOFMatrix::symmetric()
  {
    FUNCNAME("DOFMatrix::symmetric()");

122
    double tol = 1e-5;
123

124
125
126
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
127

128
129
130
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
131

132
133
134
135
136
137
138
    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)
139
140
141
142
	  return false;
    return true;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
143

144
145
146
147
  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

148
    int non_symmetric = !symmetric();
149

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

Thomas Witkowski's avatar
Thomas Witkowski committed
156

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

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
181

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

189
190
191
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode\n");
    TEST_EXIT_DBG(rowFeSpace)("Have now rowFeSpace!\n");

192
193
    inserter_type &ins= *inserter;
 
194
    // === Get indices mapping from local to global matrix indices. ===
195

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

214
215
    using namespace mtl;

216
#if 0
217
218
219
220
221
222
223
224
225
226
227
228
    if (MPI::COMM_WORLD.Get_rank() == 0) {
      std::cout << "----- PRINT MAT " << rowElInfo->getElement()->getIndex() << "--------" << 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;
    }
229
230
#endif
         
Thomas Witkowski's avatar
Thomas Witkowski committed
231
    for (int i = 0; i < nRow; i++)  {
232
      DegreeOfFreedom row = rowIndices[i];
233

234
235
236
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

237
      if (condition && condition->isDirichlet()) {	
238
	if (condition->applyBoundaryCondition()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
239

240
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
 	  if (dofMap->isRankDof(rowIndices[i])) 
 	    applyDBCs.insert(static_cast<int>(row));
243
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
244
 	  applyDBCs.insert(static_cast<int>(row));
245
246
#endif
	}
247
      } else {
248
249
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
250
	  if (fabs(elMat[i][j]) > 1e-10) {
251
	    ins[row][col] += elMat[i][j];
252
	  }
253
	}
254
      }
255
    }
256
  }
257

258

259
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
260
  {
261
    return matrix[a][b];
262
263
  }

264

265
  void DOFMatrix::freeDOFContent(int index)
266
  {}
267

268

269
270
271
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
272
  {
273
    FUNCNAME("DOFMatrix::assemble()");
274

275
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
276

277
278
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
279
    for (; it != operators.end(); ++it, ++factorIt)
280
281
      if ((*it)->getNeedDualTraverse() == false) 
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
282

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

286
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
287
288
  }

289

290
291
292
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
293
294
295
296
297
298
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

299
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
300
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
303
304
305

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
306
307
  }

308

309
310
311
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
312
313
314
315
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

316
    if (!op && operators.size() == 0)
317
318
      return;

319
320
    set_to_zero(elementMatrix);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339
340
    if (factor != 1.0)
      elementMatrix *= factor;

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

343

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
346
347
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
350
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
351
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
      return;

354
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    
356
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
359
360
361
362
363
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
364
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
368
369
370
371
372
	  (*it)->getElementMatrix(mainElInfo, 
				  auxElInfo,
				  smallElInfo, 
				  largeElInfo,
				  rowFeSpace == colFeSpace,
373
374
375
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
376
377
      }      
    }
378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
382
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
383
384
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
385

Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
393
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
394

395
  // Should work as before
396
397
398
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
399
400
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
401
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
402

403
404
405
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
406

407
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
408
  {
409
    matrix+= a * x.matrix + y.matrix;
410
411
  }

412

413
414
  void DOFMatrix::scal(double b) 
  {
415
    matrix*= b;
416
417
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
421
422
423
424
425
  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
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
471
472
473
474
475
476
477
478
479
480
481
482
483
484
  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;
    }    
  }

485

486
487
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
488
    matrix = rhs.matrix;
489
490
  }

491

Thomas Witkowski's avatar
Thomas Witkowski committed
492
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
493
  {      
494
495
496
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
497
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
498
      ins[*it][*it] = 1.0;    
499
500

    rows->clear();
501
502
  }

503

504
505
  int DOFMatrix::memsize() 
  {   
506
507
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
508
509
  }

510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
512
513
514
515
516
517
518
519
520
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
523
524
525
526

  void DOFMatrix::setDofMap(FeSpaceDofMap& m)
  {
    dofMap = &m; 
  }

527
}