DOFMatrix.cc 12.7 KB
Newer Older
1
#include <algorithm>
2
#include <boost/numeric/mtl/mtl.hpp>
3
#include "DOFMatrix.h"
4
5
6
7
8
9
10
11
12
13
14
15
#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"
16
#include "Serializer.h"
17

18
19
namespace AMDiS {

20
21
  using namespace mtl;

22
23
24
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
25
26
    : rowFeSpace(NULL),
      colFeSpace(NULL),
27
28
29
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
30
      nnzPerRow(0),
31
32
      inserter(NULL)
  {}
33

34
35
36
  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowSpace,
		       const FiniteElemSpace* colSpace,
		       std::string n)
37
38
    : rowFeSpace(rowSpace),
      colFeSpace(colSpace),
39
      name(n), 
40
      coupleMatrix(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
41
      nnzPerRow(0),
42
      inserter(NULL)
43
  {
44
45
    FUNCNAME("DOFMatrix::DOFMatrix()");

46
    TEST_EXIT(rowFeSpace)("No fe space for row!\n");
47
  
48
49
    if (!colFeSpace)
      colFeSpace = rowFeSpace;
50

51
52
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->addDOFIndexed(this);
53

54
    boundaryManager = new BoundaryManager(rowFeSpace);
55

56
57
    nRow = rowFeSpace->getBasisFcts()->getNumber();
    nCol = colFeSpace->getBasisFcts()->getNumber();
58
59
60
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
61
62

    applyDBCs.clear();
63
64
65
  }

  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
66
    : name(rhs.name + "copy")
67
  {
68
69
    FUNCNAME("DOFMatrix::DOFMatrix()");

70
    *this = rhs;
71
72
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFeSpace->getAdmin()))->addDOFIndexed(this);
73

74
75
    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
    inserter = 0;
76
77
78
79
  }

  DOFMatrix::~DOFMatrix()
  {
80
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
81

82
83
    if (rowFeSpace && rowFeSpace->getAdmin())
      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
86
87
    if (boundaryManager) 
      delete boundaryManager;
    if (inserter) 
      delete inserter;
88
89
90
91
92
  }

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
93
94
95

    if (inserter) 
      inserter->print();
96
97
98
99
100
101
  }

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

102
    double tol = 1e-5;
103

104
105
106
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
107

108
109
110
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
111

112
113
114
115
116
117
118
    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)
119
120
121
122
123
124
125
126
	  return false;
    return true;
  }

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

127
    int non_symmetric = !symmetric();
128

Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (non_symmetric)
130
      MSG("Matrix `%s' not symmetric.\n", name.data());
Thomas Witkowski's avatar
Thomas Witkowski committed
131
    else
132
      MSG("Matrix `%s' is symmetric.\n", name.data());
133
134
135
136
  }

  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
137
138
    rowFeSpace = rhs.rowFeSpace;
    colFeSpace = rhs.colFeSpace;
139
140
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
141
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
144
145
146

    /// 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
147
    if (rhs.boundaryManager)
Thomas Witkowski's avatar
Thomas Witkowski committed
148
      boundaryManager = new BoundaryManager(*rhs.boundaryManager);
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
150
      boundaryManager = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
151
    
152
153
154
    nRow = rhs.nRow;
    nCol = rhs.nCol;
    elementMatrix.change_dim(nRow, nCol);
155
156
157
158

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
159
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
160
				   const BoundaryType *bound,
161
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
162
				   ElInfo* colElInfo)
163
  {
164
    FUNCNAME("DOFMatrix::addElementMatrix()");
165

166
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
167
168
    inserter_type &ins= *inserter;
 
169
    // === Get indices mapping from local to global matrix indices. ===
170

171
172
    rowFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFeSpace->getAdmin(),
173
						rowIndices);
174
    if (rowFeSpace == colFeSpace) {
175
      colIndices = rowIndices;
176
177
    } else {
      if (colElInfo) {
178
179
	colFeSpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
180
						    colIndices);
181
182
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
183
184
	colFeSpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFeSpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
185
						    colIndices);
186
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    }
188

189
190
    using namespace mtl;

191
192
193
194
195
196
197
198
199
200
201
202
203
#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
204
    for (int i = 0; i < nRow; i++)  {
205
      DegreeOfFreedom row = rowIndices[i];
206

207
208
209
210
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
211
212
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
213
 	  if ((*rankDofs)[rowIndices[i]]) 
214
 	    applyDBCs.insert(static_cast<int>(row));
215
#else
216
	  applyDBCs.insert(static_cast<int>(row));
217
218
#endif
	}
219
      } else {
220
221
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
222
223
224
225
226

// 	  if (MPI::COMM_WORLD.Get_rank() == 0  && row <= 10 && col <= 10) {
// 	    MSG("%d/%d entry: %e\n", row, col, elMat[i][j]);
// 	  }

227
	  ins[row][col] += elMat[i][j];
228
	}
229
      }
230
    }
231
  }
232

233

234
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
235
  {
236
    return matrix[a][b];
237
238
  }

239

240
  void DOFMatrix::freeDOFContent(int index)
241
  {}
242

243

244
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
245
  {
246
    FUNCNAME("DOFMatrix::assemble()");
247

248
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
249

250
251
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
252
253
    for (; it != operators.end(); ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse() == false)
254
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);      
255

Thomas Witkowski's avatar
Thomas Witkowski committed
256
257
    if (factor != 1.0)
      elementMatrix *= factor;
258

259
    addElementMatrix(elementMatrix, bound, elInfo, NULL); 
260
261
  }

262

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

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

270
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
271
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
274
275
276

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

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

279

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

287
    if (!op && operators.size() == 0)
288
289
      return;

290
291
    set_to_zero(elementMatrix);

292
    if (op) {
293
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
294
			   elementMatrix);
295
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
298
      for (; it != operators.end(); ++it, ++factorIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
302
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
    if (factor != 1.0)
      elementMatrix *= factor;

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

311

Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
314
315
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
318
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
319
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
320
321
      return;

322
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
323
    
324
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
328
329
330
331
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
332
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
335
336
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
337
338
339
340
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
341
342
      }      
    }
343

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
347
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
348
349
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
350

Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
354
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
355
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
358
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
359

360
  // Should work as before
361
362
363
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
364
365
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
366
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
367

368
369
370
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
371

372
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
373
  {
374
    matrix+= a * x.matrix + y.matrix;
375
376
  }

377

378
379
  void DOFMatrix::scal(double b) 
  {
380
    matrix*= b;
381
382
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
383

Thomas Witkowski's avatar
Thomas Witkowski committed
384
385
386
387
388
389
390
  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
391

392
393
394
395
396
397
398
399
400
401
402
403
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
  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;
    }    
  }

450

451
452
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
453
    matrix = rhs.matrix;
454
455
  }

456

Thomas Witkowski's avatar
Thomas Witkowski committed
457
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
458
  {      
459
460
461
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

    inserter_type &ins = *inserter;  
462
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
463
      ins[*it][*it] = 1.0;    
464
465

    rows->clear();
466
467
  }

468

469
470
  int DOFMatrix::memsize() 
  {   
471
472
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
473
474
  }

475

Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
478
479
480
481
482
483
484
485
  void DOFMatrix::startInsertion(int nnz_per_row)
  {
    if (inserter) {
      delete inserter;
      inserter = NULL; 
    }

    inserter = new inserter_type(matrix, nnz_per_row);
  }

486
}