DOFMatrix.cc 13.4 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);
77
78
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
79

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
93

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
106

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
115

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

120
    double tol = 1e-5;
121

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
141

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

146
    int non_symmetric = !symmetric();
147

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

Thomas Witkowski's avatar
Thomas Witkowski committed
154

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

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
179

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

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

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

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

212
213
    using namespace mtl;

214
#if 0
215
216
217
218
219
220
221
222
223
224
225
226
    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;
    }
227
228
#endif
         
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    for (int i = 0; i < nRow; i++)  {
230
      DegreeOfFreedom row = rowIndices[i];
231

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

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

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

255

256
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
257
  {
258
    return matrix[a][b];
259
260
  }

261

262
  void DOFMatrix::freeDOFContent(int index)
263
  {}
264

265

266
267
268
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
269
  {
270
    FUNCNAME("DOFMatrix::assemble()");
271

272
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
273

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

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

283
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
284
285
  }

286

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

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

296
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
297
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301
302

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
303
304
  }

305

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

313
    if (!op && operators.size() == 0)
314
315
      return;

316
317
    set_to_zero(elementMatrix);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
337
    if (factor != 1.0)
      elementMatrix *= factor;

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

340

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

Thomas Witkowski's avatar
Thomas Witkowski committed
348
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
      return;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
378
379
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
380
381
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
382

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

Thomas Witkowski's avatar
Thomas Witkowski committed
391

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

400
401
402
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
403

404
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
405
  {
406
    matrix+= a * x.matrix + y.matrix;
407
408
  }

409

410
411
  void DOFMatrix::scal(double b) 
  {
412
    matrix*= b;
413
414
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
415

Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
418
419
420
421
422
  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
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
471
472
473
474
475
476
477
478
479
480
481
  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;
    }    
  }

482

483
484
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
485
    matrix = rhs.matrix;
486
487
  }

488

Thomas Witkowski's avatar
Thomas Witkowski committed
489
  void DOFMatrix::clearDirichletRows()
490
  {      
Thomas Witkowski's avatar
Thomas Witkowski committed
491
    FUNCNAME("DOFMatrix::clearDirichletRows()");
492

493
494
495
    // Do the following only in sequential code. In parallel mode, the specific
    // solver method must care about dirichlet boundary conditions.

496
    inserter_type &ins = *inserter;  
Thomas Witkowski's avatar
Thomas Witkowski committed
497
498
    for (std::set<int>::iterator it = dirichletDofs.begin(); 
	 it != dirichletDofs.end(); ++it)
499
      ins[*it][*it] = 1.0;    
500
501
  }

502

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

509

Thomas Witkowski's avatar
Thomas Witkowski committed
510
511
512
513
514
515
516
517
  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
518
519

    dirichletDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
522

523
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
524
525
526
527
  void DOFMatrix::setDofMap(FeSpaceDofMap& m)
  {
    dofMap = &m; 
  }
528
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
529

530
}