DOFMatrix.cc 13.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
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

Thomas Witkowski's avatar
Thomas Witkowski committed
46

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

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

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

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

67
    boundaryManager = new BoundaryManager(rowFeSpace);
68

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

    applyDBCs.clear();
76
77
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
78

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
92

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
105

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
114

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

119
    double tol = 1e-5;
120

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
140

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

145
    int non_symmetric = !symmetric();
146

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

Thomas Witkowski's avatar
Thomas Witkowski committed
153

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

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
178

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

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

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

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

211
212
    using namespace mtl;

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

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

234
      if (condition && condition->isDirichlet()) {	
235
236
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
237
238
	  if ((*rankDofs)[rowIndices[i]]) 
	    applyDBCs.insert(static_cast<int>(row));
239
#else
240
	  applyDBCs.insert(static_cast<int>(row));
241
242
#endif
	}
243
      } else {
244
245
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
246
	  if (fabs(elMat[i][j]) > 1e-10) {
247
	    ins[row][col] += elMat[i][j];
248
	  }
249
	}
250
      }
251
    }
252
  }
253

254

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

260

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

264

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

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

273
274
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
275
    for (; it != operators.end(); ++it, ++factorIt)
276
277
      if ((*it)->getNeedDualTraverse() == false) 
	(*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

282
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
283
284
  }

285

286
287
288
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
289
290
291
292
293
294
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

295
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
296
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
299
300
301

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
302
303
  }

304

305
306
307
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
308
309
310
311
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

312
    if (!op && operators.size() == 0)
313
314
      return;

315
316
    set_to_zero(elementMatrix);

317
    if (op) {
318
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
319
			   false, elementMatrix);
320
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
323
      for (; it != operators.end(); ++it, ++factorIt)
324
325
326
327
328
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
329
				elementMatrix, 
330
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
335
336
    if (factor != 1.0)
      elementMatrix *= factor;

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

339

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

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

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

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

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
379
380
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
381

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

Thomas Witkowski's avatar
Thomas Witkowski committed
390

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

399
400
401
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
402

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

408

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

Thomas Witkowski's avatar
Thomas Witkowski committed
414

Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
421
  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
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
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
  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;
    }    
  }

481

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

487

Thomas Witkowski's avatar
Thomas Witkowski committed
488
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
489
  {      
490
491
492
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
493
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
494
      ins[*it][*it] = 1.0;    
495
496

    rows->clear();
497
498
  }

499

500
501
  int DOFMatrix::memsize() 
  {   
502
503
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
504
505
  }

506

Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
509
510
511
512
513
514
515
516
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

517
}