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

18
19

//#include <iostream>
20
21
namespace AMDiS {

22
23
  using namespace mtl;

24
25
26
  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
27
28
    : rowFESpace(NULL),
      colFESpace(NULL),
29
30
31
      elementMatrix(3, 3),
      nRow(0),
      nCol(0),
32
33
      inserter(NULL)
  {}
34
35
36

  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowFESpace_,
		       const FiniteElemSpace* colFESpace_,
37
		       std::string name_)
38
39
40
    : rowFESpace(rowFESpace_),
      colFESpace(colFESpace_),
      name(name_), 
41
42
      coupleMatrix(false),
      inserter(NULL)
43
44
45
  {
    TEST_EXIT(rowFESpace)("no rowFESpace\n");
  
Thomas Witkowski's avatar
Thomas Witkowski committed
46
    if (!colFESpace)
47
48
49
      colFESpace = rowFESpace;

    if (rowFESpace && rowFESpace->getAdmin())
Thomas Witkowski's avatar
Thomas Witkowski committed
50
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
51

Thomas Witkowski's avatar
Thomas Witkowski committed
52
    boundaryManager = new BoundaryManager(rowFESpace_);
53
54
55
56
57
58

    nRow = rowFESpace->getBasisFcts()->getNumber();
    nCol = colFESpace->getBasisFcts()->getNumber();
    elementMatrix.change_dim(nRow, nCol);
    rowIndices.resize(nRow);
    colIndices.resize(nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
59
60

    applyDBCs.clear();
61
62
63
  }

  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
64
    : name(rhs.name + "copy")
65
  {
66
    *this = rhs;
67
68
    if (rowFESpace && rowFESpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
69
70
71

    TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion");
    inserter= 0;
72
73
74
75
  }

  DOFMatrix::~DOFMatrix()
  {
76
    FUNCNAME("DOFMatrix::~DOFMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
77

78
    if (rowFESpace && rowFESpace->getAdmin())
79
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
80
81
82
83
    if (boundaryManager) 
      delete boundaryManager;
    if (inserter) 
      delete inserter;
84
85
86
87
88
  }

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");
89
90
  if(inserter) inserter->print();
  /*  using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
91
92
    namespace traits= mtl::traits;
    typedef base_matrix_type Matrix;
93

94
95
96
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
97

98
99
    typedef traits::range_generator<major, Matrix>::type      cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
100
101

    std::cout.precision(10);
102
103
    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)
104
	if (value(*icursor) != 0.0)
Thomas Witkowski's avatar
nix.  
Thomas Witkowski committed
105
106
	    std::cout << "(" << row(*icursor) << "," << col(*icursor) << "," << value(*icursor) << ") ";
      
107
      std::cout << "\n";
108
    }*/
109
110
111
112
113
114
  }

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

115
    double tol = 1e-5;
116

117
118
119
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef base_matrix_type   Matrix;
120

121
122
123
    traits::row<Matrix>::type                                 row(matrix);
    traits::col<Matrix>::type                                 col(matrix);
    traits::const_value<Matrix>::type                         value(matrix);
124

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

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

140
    int non_symmetric = !symmetric();
141

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

  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
154
    coupleMatrix = rhs.coupleMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
155
156
157
158
159

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

    return *this;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
172
  void DOFMatrix::addElementMatrix(const ElementMatrix& elMat, 
173
				   const BoundaryType *bound,
174
				   ElInfo* rowElInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
175
				   ElInfo* colElInfo)
176
  {
177
    FUNCNAME("DOFMatrix::addElementMatrix()");
178

179
    TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode");
180

181
    // === Get indices mapping from local to global matrix indices. ===
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
185
    rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						rowFESpace->getAdmin(),
						rowIndices);
186
    if (rowFESpace == colFESpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
187
      *inserter << element_matrix(elMat, rowIndices, rowIndices);
188
189
    } else {
      if (colElInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
190
191
192
	colFESpace->getBasisFcts()->getLocalIndices(colElInfo->getElement(),
						    colFESpace->getAdmin(),
						    colIndices);
193
194
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
197
	colFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
						    colFESpace->getAdmin(),
						    colIndices);
198
      }
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
      *inserter << element_matrix(elMat, rowIndices, colIndices);
    }
202

203

Thomas Witkowski's avatar
Thomas Witkowski committed
204
    for (int i = 0; i < nRow; i++)  {
205
206
207
208
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
209
210
	if (condition->applyBoundaryCondition()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
211
212
	  if (rankDofs[rowIndices[i]]) 
	    applyDBCs.insert(static_cast<int>(rowIndices[i]));
213
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
214
	  applyDBCs.insert(static_cast<int>(rowIndices[i]));
215
216
#endif
	}
217
      }
218
    }
219
 }
220

221
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
222
  {
223
    return matrix[a][b];
224
225
  }

226
  void DOFMatrix::freeDOFContent(int index)
227
  {}
228

229
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
230
  {
231
    FUNCNAME("DOFMatrix::assemble()");
232

233
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
234

235
236
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
237
238
239
    for (; it != operators.end(); ++it, ++factorIt)
      if ((*it)->getNeedDualTraverse() == false)
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
240

Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
    if (factor != 1.0)
      elementMatrix *= factor;
243
244
/*MSG("elementMatrix\n");
std::cout<<elementMatrix<<"\n";*/
Thomas Witkowski's avatar
Thomas Witkowski committed
245
    addElementMatrix(elementMatrix, bound, elInfo, NULL);   
246
247
  }

248
249
250
251
252
253
254
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

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

255
      set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
256
      op->getElementMatrix(elInfo, elementMatrix, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
257
258
259
260
261

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

      addElementMatrix(elementMatrix, bound, elInfo, NULL);
262
263
  }

264
265
266
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
267
268
269
270
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

271
    if (!op && operators.size() == 0)
272
273
      return;

274
275
    set_to_zero(elementMatrix);

276
    if (op) {
277
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
278
			   elementMatrix);
279
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
      for (; it != operators.end(); ++it, ++factorIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
285
286
287
288
289
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
292
293
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
296
297
298
299
300
301
302
  }

  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
			    ElInfo *smallElInfo, ElInfo *largeElInfo,			    
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
303
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
304
305
      return;

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

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

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
332
333
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
337
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
338
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
341
      (*it)->finishAssembling();
  }

342
  // Should work as before
343
344
345
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
346
347
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
348
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
349

350
351
352
    return fillFlag;
  }

353
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
354
  {
355
    matrix+= a * x.matrix + y.matrix;
356
357
  }

358

359
360
  void DOFMatrix::scal(double b) 
  {
361
    matrix*= b;
362
363
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
367
368
369
370
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

371
372
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
373
    matrix= rhs.matrix;
374
375
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
376
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
377
  {      
378
379
380
    inserter_type &ins= *inserter;
   
    for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it)
381
382
383
384
385
386
387
388
389
390
391
392
393
      {
	//ins.make_empty(*it);
	//clear row
	for(int i=0;i<(*it);i++)
	      ins[*it][i]=0;
	//set diag
	ins[*it][*it] = 1.0;
	//clear rest of row
	
	for(int i=(*it)+1;i<matrix.num_cols();i++)
	      ins[*it][i]=0;
      }
     rows->clear();
394
395
  }

396
397
  int DOFMatrix::memsize() 
  {   
398
399
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
400
401
402
  }

}