DOFMatrix.cc 13 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
    std::cout << "----- PRINT MAT " << rowElInfo->getElement()->getIndex() << "--------" << std::endl;
215
216
217
218
219
220
221
222
223
224
225
    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
245
	  if (fabs(elMat[i][j]) > 1e-10)
	    ins[row][col] += elMat[i][j];
246
	}
247
      }
248
    }
249
  }
250

251

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

257

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

261

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

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

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

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

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

280

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

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

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

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

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

297

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

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

308
309
    set_to_zero(elementMatrix);

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

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

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

332

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
374

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

Thomas Witkowski's avatar
Thomas Witkowski committed
383

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

392
393
394
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
395

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

401

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

Thomas Witkowski's avatar
Thomas Witkowski committed
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
410
411
412
413
414
  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
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
473
  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;
    }    
  }

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::removeRowsWithDBC(std::set<int> *rows)
482
  {      
483
484
485
    FUNCNAME("DOFMatrix::removeRowsWithDBC()");

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

    rows->clear();
490
491
  }

492

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

499

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

    inserter = new inserter_type(matrix, nnz_per_row);
  }

510
}