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

Thomas Witkowski's avatar
Thomas Witkowski committed
44

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

57
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
58

59
60
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
61

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

65
    boundaryManager = new BoundaryManager(rowFeSpace);
66

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

Thomas Witkowski's avatar
Thomas Witkowski committed
74

75
  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
76
    : name(rhs.name + "copy")
77
  {
78
79
    FUNCNAME("DOFMatrix::DOFMatrix()");

80
    *this = rhs;
81
82
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFeSpace->getAdmin()))->addDOFIndexed(this);
83

84
85
    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
    inserter = 0;
86
87
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
88

89
90
  DOFMatrix::~DOFMatrix()
  {
91
92
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
93
94
95
96
    if (boundaryManager) 
      delete boundaryManager;
    if (inserter) 
      delete inserter;
97
98
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
99

100
101
  void DOFMatrix::print() const
  {
102
103
    if (inserter) 
      inserter->print();
104
105
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
106

107
108
  bool DOFMatrix::symmetric()
  {
109
    double tol = 1e-5;
110

111
112
113
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
114

115
116
117
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
118

119
120
121
122
123
124
125
    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)
126
127
128
129
	  return false;
    return true;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
130

131
132
133
134
  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

135
    int non_symmetric = !symmetric();
136

Thomas Witkowski's avatar
Thomas Witkowski committed
137
    if (non_symmetric)
138
      MSG("Matrix `%s' not symmetric.\n", name.data());
Thomas Witkowski's avatar
Thomas Witkowski committed
139
    else
140
      MSG("Matrix `%s' is symmetric.\n", name.data());
141
142
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
143

144
145
  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
146
147
    rowFeSpace = rhs.rowFeSpace;
    colFeSpace = rhs.colFeSpace;
148
149
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
150
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
154
155

    /// 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
156
    if (rhs.boundaryManager)
Thomas Witkowski's avatar
Thomas Witkowski committed
157
      boundaryManager = new BoundaryManager(*rhs.boundaryManager);
Thomas Witkowski's avatar
Thomas Witkowski committed
158
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
159
      boundaryManager = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
160
    
161
162
163
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
164
165
166
167

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
168

Thomas Witkowski's avatar
Thomas Witkowski committed
169
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
170
				   const BoundaryType *bound,
171
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
172
				   ElInfo* colElInfo)
173
  {
174
    FUNCNAME_DBG("DOFMatrix::addElementMatrix()");
175

176
177
178
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode\n");
    TEST_EXIT_DBG(rowFeSpace)("Have now rowFeSpace!\n");

179
180
    inserter_type &ins= *inserter;
 
181
    // === Get indices mapping from local to global matrix indices. ===
182

183
184
    rowFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFeSpace->getAdmin(),
185
						rowIndices);
186
    if (rowFeSpace == colFeSpace) {
187
      colIndices = rowIndices;
188
189
    } else {
      if (colElInfo) {
190
191
	colFeSpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
192
						    colIndices);
193
194
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
195
196
	colFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
197
						    colIndices);
198
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
199
    }
200

201
202
    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
203
    for (int i = 0; i < nRow; i++)  {
204
      DegreeOfFreedom row = rowIndices[i];
205

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

209
      if (condition && condition->isDirichlet()) {	
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
210
	if (condition->applyBoundaryCondition())
Thomas Witkowski's avatar
Thomas Witkowski committed
211
	  dirichletDofs.insert(row);
212
      } else {
213
214
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
215
	  ins[row][col] += elMat[i][j];
216
	}
217
      }
218
    }
219
  }
220

221

222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
  void DOFMatrix::assembleOperator(Operator &op)
  {
    FUNCNAME("DOFMatrix::assembleOperator()");

    TEST_EXIT(rowFeSpace->getMesh() == colFeSpace->getMesh())
      ("This function does not support for multi mesh procedure!\n");
    TEST_EXIT(op.getRowFeSpace() == rowFeSpace)
      ("Row FE spaces do not fit together!\n");
    TEST_EXIT(op.getColFeSpace() == colFeSpace)
      ("Column FE spaces do not fit together!\n");

    clearOperators();
    addOperator(&op);

    matrix.change_dim(rowFeSpace->getAdmin()->getUsedSize(),
		      colFeSpace->getAdmin()->getUsedSize());

    Mesh *mesh = rowFeSpace->getMesh();
    mesh->dofCompress();
    const BasisFunction *basisFcts = rowFeSpace->getBasisFcts();

    Flag assembleFlag = getAssembleFlag() | 
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH |
      Mesh::FILL_BOUND;

    BoundaryType *bound = new BoundaryType[basisFcts->getNumber()];

    calculateNnz();
    if (getBoundaryManager())
      getBoundaryManager()->initMatrix(this);

    startInsertion(getNnz());

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
    while (elInfo) {
      basisFcts->getBound(elInfo, bound);

      assemble(1.0, elInfo, bound);

      if (getBoundaryManager())
	getBoundaryManager()->fillBoundaryConditions(elInfo, this);

      elInfo = stack.traverseNext(elInfo);
    }

    clearDirichletRows();
    finishAssembling();
    finishInsertion();
    getBoundaryManager()->exitMatrix(this);

    delete [] bound;
  }


281
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
282
  {
283
    return matrix[a][b];
284
285
  }

286

287
  void DOFMatrix::freeDOFContent(int index)
288
  {}
289

290

291
292
293
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
294
  {
295
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
296

297
298
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
299
    for (; it != operators.end(); ++it, ++factorIt)
300
301
      if ((*it)->getNeedDualTraverse() == false && 
	  (*factorIt == NULL || **factorIt != 0.0))
302
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
303

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

Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
    if (operators.size())
      addElementMatrix(elementMatrix, bound, elInfo, NULL); 
309
310
  }

311

312
313
314
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
315
316
			   Operator *op)
  {
317
    FUNCNAME_DBG("DOFMatrix::assemble()");
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
320
321
322
323
324
    
    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
325
	elementMatrix *= factor;
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
    
    addElementMatrix(elementMatrix, bound, elInfo, NULL);
328
329
  }

330

331
332
333
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
334
335
			   const BoundaryType *bound, Operator *op)
  {
336
    if (!op && operators.size() == 0)
337
338
      return;

339
340
    set_to_zero(elementMatrix);

341
    if (op) {
342
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
343
			   false, elementMatrix);
344
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
347
      for (; it != operators.end(); ++it, ++factorIt)
348
349
350
351
352
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
353
				elementMatrix, 
354
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
359
    if (factor != 1.0)
      elementMatrix *= factor;

Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
    if (op || operators.size())
      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
  }

364

Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
367
368
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
371
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
372
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
      return;

375
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
376
    
377
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
380
381
382
383
384
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
385
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
389
390
391
392
393
	  (*it)->getElementMatrix(mainElInfo, 
				  auxElInfo,
				  smallElInfo, 
				  largeElInfo,
				  rowFeSpace == colFeSpace,
394
395
396
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
397
398
      }      
    }
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
404
405
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
406

Thomas Witkowski's avatar
Thomas Witkowski committed
407
408
409
410
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
411
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
414
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
415

416
  // Should work as before
417
418
419
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
420
421
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
422
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
423

424
425
426
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
427

428
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
429
  {
430
    matrix += a * x.matrix + y.matrix;
431
432
  }

433

434
435
  void DOFMatrix::scal(double b) 
  {
436
    matrix *= b;
437
438
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
439

Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
442
443
444
445
446
  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
447

448
449
450
451
452
453
454
455
  void DOFMatrix::clearOperators()
  {
    operators.clear();
    operatorFactor.clear();
    operatorEstFactor.clear();
  }


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
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
  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;
    }    
  }

514

515
516
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
517
    matrix = rhs.matrix;
518
519
  }

520

Thomas Witkowski's avatar
Thomas Witkowski committed
521
  void DOFMatrix::clearDirichletRows()
522
  {      
523
524
    // Do the following only in sequential code. In parallel mode, the specific
    // solver method must care about dirichlet boundary conditions.
525
    inserter_type &ins = *inserter;  
Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
    for (std::set<int>::iterator it = dirichletDofs.begin(); 
	 it != dirichletDofs.end(); ++it)
528
      ins[*it][*it] = 1.0;
529
530
  }

531

532
533
  int DOFMatrix::memsize() 
  {   
534
535
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
536
537
  }

538

Thomas Witkowski's avatar
Thomas Witkowski committed
539
540
541
542
543
544
545
546
  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
547
548

    dirichletDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
549
550
  }

551
}