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
  DOFMatrix::DOFMatrix()
38
39
    : rowFeSpace(NULL),
      colFeSpace(NULL),
40
41
42
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
43
      nnzPerRow(0),
44
45
      inserter(NULL)
  {}
46

Thomas Witkowski's avatar
Thomas Witkowski committed
47

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

60
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
61

62
63
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
64

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

68
    boundaryManager = new BoundaryManager(rowFeSpace);
69

70
71
    nRow = rowFeSpace->getBasisFcts()->getNumber();
    nCol = colFeSpace->getBasisFcts()->getNumber();
72
73
74
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
75
76
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
77

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
91

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
104

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
113

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

118
    double tol = 1e-5;
119

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
139

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

144
    int non_symmetric = !symmetric();
145

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

Thomas Witkowski's avatar
Thomas Witkowski committed
152

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

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
177

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

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

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

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

210
211
    using namespace mtl;

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

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

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

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

253

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

259

260
  void DOFMatrix::freeDOFContent(int index)
261
  {}
262

263

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

270
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
271

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
    if (operators.size())
      addElementMatrix(elementMatrix, bound, elInfo, NULL); 
284
285
  }

286

287
288
289
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
290
291
			   Operator *op)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
294
295
296
297
298
299
    FUNCNAME("DOFMatrix::assemble()");
    
    TEST_EXIT_DBG(op)("No operator!\n");
    
    set_to_zero(elementMatrix);
    op->getElementMatrix(elInfo, elementMatrix, factor);
    
    if (factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
300
	elementMatrix *= factor;
Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
    
    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
    if (factor != 1.0)
      elementMatrix *= factor;

Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
    if (op || operators.size())
      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
  }

341

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
383

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

Thomas Witkowski's avatar
Thomas Witkowski committed
392

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

401
402
403
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
404

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

410

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

Thomas Witkowski's avatar
Thomas Witkowski committed
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
420
421
422
423
  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
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
482
  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;
    }    
  }

483

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

489

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

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

497
    inserter_type &ins = *inserter;  
Thomas Witkowski's avatar
Thomas Witkowski committed
498
499
    for (std::set<int>::iterator it = dirichletDofs.begin(); 
	 it != dirichletDofs.end(); ++it)
500
      ins[*it][*it] = 1.0;    
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
  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
519
520

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

Thomas Witkowski's avatar
Thomas Witkowski committed
523

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

531
}