DOFMatrix.cc 14.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"
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
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
92

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

Thomas Witkowski's avatar
Thomas Witkowski committed
101

102
103
104
  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
105
106
107

    if (inserter) 
      inserter->print();
108
109
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
110

111
112
113
114
  bool DOFMatrix::symmetric()
  {
    FUNCNAME("DOFMatrix::symmetric()");

115
    double tol = 1e-5;
116

117
118
119
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
120

121
122
123
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
124

125
126
127
128
129
130
131
    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)
132
133
134
135
	  return false;
    return true;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
136

137
138
139
140
  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

141
    int non_symmetric = !symmetric();
142

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

Thomas Witkowski's avatar
Thomas Witkowski committed
149

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

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
174

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

182
183
184
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode\n");
    TEST_EXIT_DBG(rowFeSpace)("Have now rowFeSpace!\n");

185
186
    inserter_type &ins= *inserter;
 
187
    // === Get indices mapping from local to global matrix indices. ===
188

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

207
208
    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
209
    for (int i = 0; i < nRow; i++)  {
210
      DegreeOfFreedom row = rowIndices[i];
211

212
213
214
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

215
      if (condition && condition->isDirichlet()) {	
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
216
	if (condition->applyBoundaryCondition())
Thomas Witkowski's avatar
Thomas Witkowski committed
217
	  dirichletDofs.insert(row);
218
      } else {
219
220
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
221
	  ins[row][col] += elMat[i][j];
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
281
282
283
284
285
286
  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;
  }


287
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
288
  {
289
    return matrix[a][b];
290
291
  }

292

293
  void DOFMatrix::freeDOFContent(int index)
294
  {}
295

296

297
298
299
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
300
  {
301
    FUNCNAME("DOFMatrix::assemble()");
302

303
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
304

305
306
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
307
    for (; it != operators.end(); ++it, ++factorIt)
308
309
      if ((*it)->getNeedDualTraverse() == false && 
	  (*factorIt == NULL || **factorIt != 0.0))
310
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
311

Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
    if (factor != 1.0)
      elementMatrix *= factor;
314

Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
    if (operators.size())
      addElementMatrix(elementMatrix, bound, elInfo, NULL); 
317
318
  }

319

320
321
322
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
323
324
			   Operator *op)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
328
329
330
331
332
    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
333
	elementMatrix *= factor;
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
    
    addElementMatrix(elementMatrix, bound, elInfo, NULL);
336
337
  }

338

339
340
341
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
342
343
344
345
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

346
    if (!op && operators.size() == 0)
347
348
      return;

349
350
    set_to_zero(elementMatrix);

351
    if (op) {
352
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
353
			   false, elementMatrix);
354
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
357
      for (; it != operators.end(); ++it, ++factorIt)
358
359
360
361
362
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
363
				elementMatrix, 
364
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
    if (op || operators.size())
      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
  }

374

Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
377
378
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
      return;

385
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
386
    
387
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
390
391
392
393
394
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
395
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
398
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
399
400
401
402
403
	  (*it)->getElementMatrix(mainElInfo, 
				  auxElInfo,
				  smallElInfo, 
				  largeElInfo,
				  rowFeSpace == colFeSpace,
404
405
406
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
407
408
      }      
    }
409

Thomas Witkowski's avatar
Thomas Witkowski committed
410
411
412
413
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
414
415
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
420
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
421
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
424
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
425

426
  // Should work as before
427
428
429
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
430
431
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
432
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
433

434
435
436
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
437

438
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
439
  {
440
    matrix += a * x.matrix + y.matrix;
441
442
  }

443

444
445
  void DOFMatrix::scal(double b) 
  {
446
    matrix *= b;
447
448
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
449

Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
452
453
454
455
456
  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
457

458
459
460
461
462
463
464
465
  void DOFMatrix::clearOperators()
  {
    operators.clear();
    operatorFactor.clear();
    operatorEstFactor.clear();
  }


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
514
515
516
517
518
519
520
521
522
523
  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;
    }    
  }

524

525
526
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
527
    matrix = rhs.matrix;
528
529
  }

530

Thomas Witkowski's avatar
Thomas Witkowski committed
531
  void DOFMatrix::clearDirichletRows()
532
  {      
Thomas Witkowski's avatar
Thomas Witkowski committed
533
    FUNCNAME("DOFMatrix::clearDirichletRows()");
534

535
536
    // Do the following only in sequential code. In parallel mode, the specific
    // solver method must care about dirichlet boundary conditions.
537
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
538
    inserter_type &ins = *inserter;  
Thomas Witkowski's avatar
Thomas Witkowski committed
539
540
    for (std::set<int>::iterator it = dirichletDofs.begin(); 
	 it != dirichletDofs.end(); ++it)
541
	ins[*it][*it] = 1.0;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
542
#endif
543
544
  }

545

546
547
  int DOFMatrix::memsize() 
  {   
548
549
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
550
551
  }

552

Thomas Witkowski's avatar
Thomas Witkowski committed
553
554
555
556
557
558
559
560
  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
561
562

    dirichletDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
563
564
  }

565
}