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

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
214
215
216
217
218
219
220
221
222
223
224
225
#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
226
    for (int i = 0; i < nRow; i++)  {
227
      DegreeOfFreedom row = rowIndices[i];
228

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

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

250

251
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
252
  {
253
    return matrix[a][b];
254
255
  }

256

257
  void DOFMatrix::freeDOFContent(int index)
258
  {}
259

260

261
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
262
  {
263
    FUNCNAME("DOFMatrix::assemble()");
264

265
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
266

267
268
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
269
    for (; it != operators.end(); ++it, ++factorIt)
270
271
      if ((*it)->getNeedDualTraverse() == false) 
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
272

Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
    if (factor != 1.0)
      elementMatrix *= factor;
275

276
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
277
278
  }

279

280
281
282
283
284
285
286
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

287
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
288
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
289
290
291
292
293

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
294
295
  }

296

297
298
299
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
300
301
302
303
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

304
    if (!op && operators.size() == 0)
305
306
      return;

307
308
    set_to_zero(elementMatrix);

309
    if (op) {
310
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
311
			   false, elementMatrix);
312
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
313
314
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
315
      for (; it != operators.end(); ++it, ++factorIt)
316
317
318
319
320
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
321
				elementMatrix, 
322
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
328
    if (factor != 1.0)
      elementMatrix *= factor;

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

331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
334
335
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
338
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
339
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
      return;

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

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

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
371
372
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
376
377
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
378
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
382

383
  // Should work as before
384
385
386
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
387
388
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
389
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
390

391
392
393
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
394

395
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
396
  {
397
    matrix+= a * x.matrix + y.matrix;
398
399
  }

400

401
402
  void DOFMatrix::scal(double b) 
  {
403
    matrix*= b;
404
405
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
406

Thomas Witkowski's avatar
Thomas Witkowski committed
407
408
409
410
411
412
413
  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
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
463
464
465
466
467
468
469
470
471
472
  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;
    }    
  }

473

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

479

Thomas Witkowski's avatar
Thomas Witkowski committed
480
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
481
  {      
482
483
484
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
485
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
486
      ins[*it][*it] = 1.0;    
487
488

    rows->clear();
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
503
504
505
506
507
508
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

509
}