DOFMatrix.cc 13.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include <algorithm>
23
#include <boost/numeric/mtl/mtl.hpp>
24
#include "DOFMatrix.h"
25
26
27
28
29
30
31
32
33
34
35
36
#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"
37
#include "Serializer.h"
38

39
40
namespace AMDiS {

41
42
  using namespace mtl;

43
  DOFMatrix::DOFMatrix()
44
45
    : rowFeSpace(nullptr),
      colFeSpace(nullptr),
46
47
48
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      nnzPerRow(0),
50
      inserter(nullptr)
51
  {}
52

Thomas Witkowski's avatar
Thomas Witkowski committed
53

54
55
56
  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowSpace,
		       const FiniteElemSpace* colSpace,
		       std::string n)
57
58
    : rowFeSpace(rowSpace),
      colFeSpace(colSpace),
59
      name(n), 
60
      coupleMatrix(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
61
      nnzPerRow(0),
62
      inserter(nullptr)
63
  {
64
65
    FUNCNAME("DOFMatrix::DOFMatrix()");

66
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
67

68
69
70
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
    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
    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
    inserter = 0;
88
89
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
90

91
92
  DOFMatrix::~DOFMatrix()
  {
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
  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
109
110
    rowFeSpace = rhs.rowFeSpace;
    colFeSpace = rhs.colFeSpace;
111
112
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
113
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
116
117
118

    /// 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
119
    if (rhs.boundaryManager)
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      boundaryManager = new BoundaryManager(*rhs.boundaryManager);
Thomas Witkowski's avatar
Thomas Witkowski committed
121
    else
122
      boundaryManager = nullptr;
Thomas Witkowski's avatar
Thomas Witkowski committed
123
    
124
125
126
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
127
128
129
130

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
133
				   const BoundaryType *bound,
134
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
135
				   ElInfo* colElInfo)
136
  {
137
    FUNCNAME_DBG("DOFMatrix::addElementMatrix()");
138

139
140
141
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode\n");
    TEST_EXIT_DBG(rowFeSpace)("Have now rowFeSpace!\n");

142
143
    inserter_type &ins= *inserter;
 
144
    // === Get indices mapping from local to global matrix indices. ===
145

146
147
    rowFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFeSpace->getAdmin(),
148
						rowIndices);
149
    if (rowFeSpace == colFeSpace) {
150
      colIndices = rowIndices;
151
152
    } else {
      if (colElInfo) {
153
154
	colFeSpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
155
						    colIndices);
156
157
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
158
159
	colFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
160
						    colIndices);
161
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    }
163

164
165
    using namespace mtl;

Thomas Witkowski's avatar
Thomas Witkowski committed
166
    for (int i = 0; i < nRow; i++)  {
167
      DegreeOfFreedom row = rowIndices[i];
168

169
      BoundaryCondition *condition = 
170
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : nullptr;
171

172
      if (condition && condition->isDirichlet()) {	
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
173
	if (condition->applyBoundaryCondition())
Thomas Witkowski's avatar
Thomas Witkowski committed
174
	  dirichletDofs.insert(row);
175
      } else {
176
177
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
178
	  ins[row][col] += elMat[i][j];
179
	}
180
      }
181
    }
182
  }
183

184

185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
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
  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;
  }


244
#if 0
245
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
246
  {
247
    return matrix[a][b];
248
  }
249
#endif
250

251
252
253
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
254
  {
255
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
256

257
258
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
259
    for (; it != operators.end(); ++it, ++factorIt)
260
      if ((*it)->getNeedDualTraverse() == false && 
261
	  (*factorIt == nullptr || **factorIt != 0.0))
262
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
    if (factor != 1.0)
      elementMatrix *= factor;
266

Thomas Witkowski's avatar
Thomas Witkowski committed
267
    if (operators.size())
268
      addElementMatrix(elementMatrix, bound, elInfo, nullptr); 
269
270
  }

271

272
273
274
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
275
276
			   Operator *op)
  {
277
    FUNCNAME_DBG("DOFMatrix::assemble()");
Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
280
281
282
283
284
    
    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
285
	elementMatrix *= factor;
Thomas Witkowski's avatar
Thomas Witkowski committed
286
    
287
    addElementMatrix(elementMatrix, bound, elInfo, nullptr);
288
289
  }

290

291
292
293
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
294
295
			   const BoundaryType *bound, Operator *op)
  {
296
    if (!op && operators.size() == 0)
297
298
      return;

299
300
    set_to_zero(elementMatrix);

301
    if (op) {
302
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
303
			   false, elementMatrix);
304
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
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
310
311
312
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
313
				elementMatrix, 
314
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
319
    if (factor != 1.0)
      elementMatrix *= factor;

Thomas Witkowski's avatar
Thomas Witkowski committed
320
321
    if (op || operators.size())
      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
  }

324

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
327
328
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
332
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
      return;

335
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    
337
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
340
341
342
343
344
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
345
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
348
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
349
350
351
352
353
	  (*it)->getElementMatrix(mainElInfo, 
				  auxElInfo,
				  smallElInfo, 
				  largeElInfo,
				  rowFeSpace == colFeSpace,
354
355
356
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
357
358
      }      
    }
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362
    if (factor != 1.0)
      elementMatrix *= factor;

363
    addElementMatrix(elementMatrix, bound, mainElInfo, nullptr);       
364
365
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
  void DOFMatrix::finishAssembling()
  {
369
    // call the operators cleanup procedures
Thomas Witkowski's avatar
Thomas Witkowski committed
370
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
371
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
375

376
  // Should work as before
377
378
379
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
380
381
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
382
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
383

384
385
386
    return fillFlag;
  }

387
#if 0
388
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
389
  {
390
    matrix += a * x.matrix + y.matrix;
391
392
  }

393

394
395
  void DOFMatrix::scal(double b) 
  {
396
    matrix *= b;
397
  }
398
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
404
405
406
  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
407

408
409
410
411
412
413
414
415
  void DOFMatrix::clearOperators()
  {
    operators.clear();
    operatorFactor.clear();
    operatorEstFactor.clear();
  }


416
417
418
419
  void DOFMatrix::serialize(std::ostream &out)
  {
    using namespace mtl; 
    
420
421
    typedef mtl::traits::range_generator<mtl::tag::major, base_matrix_type>::type c_type;
    typedef mtl::traits::range_generator<mtl::tag::nz, c_type>::type ic_type;
422
423
424
425
    
    typedef Collection<base_matrix_type>::size_type size_type;
    typedef Collection<base_matrix_type>::value_type value_type;
    
426
427
428
    mtl::traits::row<base_matrix_type>::type row(matrix); 
    mtl::traits::col<base_matrix_type>::type col(matrix);
    mtl::traits::const_value<base_matrix_type>::type value(matrix); 
429
430
431
432
433
434
    
    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);
    
435
436
437
438
    for (c_type cursor(mtl::begin<mtl::tag::major>(matrix)), 
	   cend(mtl::end<mtl::tag::major>(matrix)); cursor != cend; ++cursor)
      for (ic_type icursor(mtl::begin<mtl::tag::nz>(cursor)), 
	     icend(mtl::end<mtl::tag::nz>(cursor)); icursor != icend; ++icursor) {
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
	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;
    }    
  }

474

475
476
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
477
    matrix = rhs.matrix;
478
479
  }

480

Thomas Witkowski's avatar
Thomas Witkowski committed
481
  void DOFMatrix::clearDirichletRows()
482
  {      
483
484
    // Do the following only in sequential code. In parallel mode, the specific
    // solver method must care about dirichlet boundary conditions.
485
    inserter_type &ins = *inserter;  
486
    for (std::set<DegreeOfFreedom>::iterator it = dirichletDofs.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
487
	 it != dirichletDofs.end(); ++it)
488
      ins[(*it)][(*it)] = 1.0;
489
490
  }

491

492
493
  int DOFMatrix::memsize() 
  {   
494
495
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
496
497
  }

498

Thomas Witkowski's avatar
Thomas Witkowski committed
499
500
501
502
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
503
      inserter = nullptr; 
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
506
    }

    inserter = new inserter_type(matrix, nnz_per_row);
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508

    dirichletDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
  }

511
}