DOFMatrix.cc 13.4 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
29
30
31
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/ParallelDofMapping.h"
#endif
32

33
34
namespace AMDiS {

35
36
  using namespace mtl;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
47

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

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

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

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

68
    boundaryManager = new BoundaryManager(rowFeSpace);
69

70
71
    nRow = rowFeSpace->getBasisFcts()->getNumber();
    nCol = colFeSpace->getBasisFcts()->getNumber();
72
73
74
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
75
76
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
77

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
91

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
104

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
113

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

118
    double tol = 1e-5;
119

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
139

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

144
    int non_symmetric = !symmetric();
145

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

Thomas Witkowski's avatar
Thomas Witkowski committed
152

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

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
177

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

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

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

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

210
211
    using namespace mtl;

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

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

233
      if (condition && condition->isDirichlet()) {	
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
234
	if (condition->applyBoundaryCondition())
Thomas Witkowski's avatar
Thomas Witkowski committed
235
	  dirichletDofs.insert(row);
236
      } else {
237
238
	for (int j = 0; j < nCol; j++) {
	  DegreeOfFreedom col = colIndices[j];
239
	  ins[row][col] += elMat[i][j];
240
	}
241
      }
242
    }
243
  }
244

245

246
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
247
  {
248
    return matrix[a][b];
249
250
  }

251

252
  void DOFMatrix::freeDOFContent(int index)
253
  {}
254

255

256
257
258
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
259
  {
260
    FUNCNAME("DOFMatrix::assemble()");
261

262
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
263

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
    if (operators.size())
      addElementMatrix(elementMatrix, bound, elInfo, NULL); 
276
277
  }

278

279
280
281
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
282
283
			   Operator *op)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
284
285
286
287
288
289
290
291
    FUNCNAME("DOFMatrix::assemble()");
    
    TEST_EXIT_DBG(op)("No operator!\n");
    
    set_to_zero(elementMatrix);
    op->getElementMatrix(elInfo, elementMatrix, factor);
    
    if (factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
292
	elementMatrix *= factor;
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
    
    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
    if (factor != 1.0)
      elementMatrix *= factor;

Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
    if (op || operators.size())
      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
  }

333

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
375

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

Thomas Witkowski's avatar
Thomas Witkowski committed
384

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

393
394
395
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
396

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

402

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

Thomas Witkowski's avatar
Thomas Witkowski committed
408

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

475

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

481

Thomas Witkowski's avatar
Thomas Witkowski committed
482
  void DOFMatrix::clearDirichletRows()
483
  {      
Thomas Witkowski's avatar
Thomas Witkowski committed
484
    FUNCNAME("DOFMatrix::clearDirichletRows()");
485

486
487
488
    // Do the following only in sequential code. In parallel mode, the specific
    // solver method must care about dirichlet boundary conditions.

489
    inserter_type &ins = *inserter;  
Thomas Witkowski's avatar
Thomas Witkowski committed
490
491
    for (std::set<int>::iterator it = dirichletDofs.begin(); 
	 it != dirichletDofs.end(); ++it)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
492
493
494
495
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      if (dofMap->isRankDof(*it))
#endif
	ins[*it][*it] = 1.0;    
496
497
  }

498

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

505

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

    inserter = new inserter_type(matrix, nnz_per_row);
Thomas Witkowski's avatar
Thomas Witkowski committed
514
515

    dirichletDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
518

519
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
522
523
  void DOFMatrix::setDofMap(FeSpaceDofMap& m)
  {
    dofMap = &m; 
  }
524
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
525

526
}