DOFMatrix.cc 12.9 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
35
36
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
37
38
    : rowFeSpace(NULL),
      colFeSpace(NULL),
39
40
41
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      nnzPerRow(0),
43
44
      inserter(NULL)
  {}
45

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

58
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
59
  
60
61
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
62

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

66
    boundaryManager = new BoundaryManager(rowFeSpace);
67

68
69
    nRow = rowFeSpace->getBasisFcts()->getNumber();
    nCol = colFeSpace->getBasisFcts()->getNumber();
70
71
72
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
73
74

    applyDBCs.clear();
75
76
77
  }

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

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

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

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

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

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

    if (inserter) 
      inserter->print();
108
109
110
111
112
113
  }

  bool DOFMatrix::symmetric()
  {
    FUNCNAME("DOFMatrix::symmetric()");

114
    double tol = 1e-5;
115

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

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

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

  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

139
    int non_symmetric = !symmetric();
140

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

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

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

    return *this;
  }

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

178
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
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;

203
204
205
206
207
208
209
210
211
212
213
214
215
#if 0
    std::cout << "----- PRINT MAT --------" << 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;
#endif
         
Thomas Witkowski's avatar
Thomas Witkowski committed
216
    for (int i = 0; i < nRow; i++)  {
217
      DegreeOfFreedom row = rowIndices[i];
218

219
220
221
222
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
223
224
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
225
226
	  if ((*rankDofs)[rowIndices[i]]) 
	    applyDBCs.insert(static_cast<int>(row));
227
#else
228
	  applyDBCs.insert(static_cast<int>(row));
229
230
#endif
	}
231
      } else {
232
233
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
234
	  ins[row][col] += elMat[i][j];
235
	}
236
      }
237
    }
238
  }
239

240

241
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
242
  {
243
    return matrix[a][b];
244
245
  }

246

247
  void DOFMatrix::freeDOFContent(int index)
248
  {}
249

250

251
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
252
  {
253
    FUNCNAME("DOFMatrix::assemble()");
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)
Thomas Witkowski's avatar
Thomas Witkowski committed
260
      if ((*it)->getNeedDualTraverse() == false)
261
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);      
262

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

266
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
267
268
  }

269

270
271
272
273
274
275
276
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

277
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
278
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
281
282
283

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

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

286

287
288
289
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
290
291
292
293
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

294
    if (!op && operators.size() == 0)
295
296
      return;

297
298
    set_to_zero(elementMatrix);

299
    if (op) {
300
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
301
			   false, elementMatrix);
302
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
305
      for (; it != operators.end(); ++it, ++factorIt)
306
307
308
309
310
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
311
				elementMatrix, 
312
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
313
314
    }

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

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

321

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

Thomas Witkowski's avatar
Thomas Witkowski committed
329
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
      return;

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

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

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
361
362
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
363

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

Thomas Witkowski's avatar
Thomas Witkowski committed
372

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

381
382
383
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
384

385
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
386
  {
387
    matrix+= a * x.matrix + y.matrix;
388
389
  }

390

391
392
  void DOFMatrix::scal(double b) 
  {
393
    matrix*= b;
394
395
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
396

Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
399
400
401
402
403
  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
404

405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
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
  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;
    }    
  }

463

464
465
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
466
    matrix = rhs.matrix;
467
468
  }

469

Thomas Witkowski's avatar
Thomas Witkowski committed
470
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
471
  {      
472
473
474
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
475
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
476
      ins[*it][*it] = 1.0;    
477
478

    rows->clear();
479
480
  }

481

482
483
  int DOFMatrix::memsize() 
  {   
484
485
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
486
487
  }

488

Thomas Witkowski's avatar
Thomas Witkowski committed
489
490
491
492
493
494
495
496
497
498
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

499
}