DOFMatrix.cc 15.1 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 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
  void DOFMatrix::assembleOperator(Operator &op)
  {
    FUNCNAME("DOFMatrix::assembleOperator()");

    TEST_EXIT(rowFeSpace->getMesh() == colFeSpace->getMesh())
      ("This function does not support for multi mesh procedure!\n");
    TEST_EXIT(op.getRowFeSpace() == rowFeSpace)
      ("Row FE spaces do not fit together!\n");
    TEST_EXIT(op.getColFeSpace() == colFeSpace)
      ("Column FE spaces do not fit together!\n");

    clearOperators();
    addOperator(&op);

    matrix.change_dim(rowFeSpace->getAdmin()->getUsedSize(),
		      colFeSpace->getAdmin()->getUsedSize());

    Mesh *mesh = rowFeSpace->getMesh();
    mesh->dofCompress();
    const BasisFunction *basisFcts = rowFeSpace->getBasisFcts();

    Flag assembleFlag = getAssembleFlag() | 
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH |
      Mesh::FILL_BOUND;

    BoundaryType *bound = new BoundaryType[basisFcts->getNumber()];

    calculateNnz();
    if (getBoundaryManager())
      getBoundaryManager()->initMatrix(this);

    startInsertion(getNnz());

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
    while (elInfo) {
      basisFcts->getBound(elInfo, bound);

      assemble(1.0, elInfo, bound);

      if (getBoundaryManager())
	getBoundaryManager()->fillBoundaryConditions(elInfo, this);

      elInfo = stack.traverseNext(elInfo);
    }

    clearDirichletRows();
    finishAssembling();
    finishInsertion();
    getBoundaryManager()->exitMatrix(this);

    delete [] bound;
  }


305
  double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
306
  {
307
    return matrix[a][b];
308 309
  }

310

311
  void DOFMatrix::freeDOFContent(int index)
312
  {}
313

314

315 316 317
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound)
318
  {
319
    FUNCNAME("DOFMatrix::assemble()");
320

321
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
322

323 324
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
325
    for (; it != operators.end(); ++it, ++factorIt)
326 327
      if ((*it)->getNeedDualTraverse() == false && 
	  (*factorIt == NULL || **factorIt != 0.0))
328
	(*it)->getElementMatrix(elInfo,	elementMatrix, *factorIt ? **factorIt : 1.0);
329

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

Thomas Witkowski's avatar
Thomas Witkowski committed
333 334
    if (operators.size())
      addElementMatrix(elementMatrix, bound, elInfo, NULL); 
335 336
  }

337

338 339 340
  void DOFMatrix::assemble(double factor, 
			   ElInfo *elInfo, 
			   const BoundaryType *bound,
341 342
			   Operator *op)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
343 344 345 346 347 348 349 350
    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
351
	elementMatrix *= factor;
Thomas Witkowski's avatar
Thomas Witkowski committed
352 353
    
    addElementMatrix(elementMatrix, bound, elInfo, NULL);
354 355
  }

356

357 358 359
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
360 361 362 363
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

364
    if (!op && operators.size() == 0)
365 366
      return;

367 368
    set_to_zero(elementMatrix);

369
    if (op) {
370
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
371
			   false, elementMatrix);
372
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
373 374
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
375
      for (; it != operators.end(); ++it, ++factorIt)
376 377 378 379 380
	(*it)->getElementMatrix(rowElInfo, 
				colElInfo,
				smallElInfo, 
				largeElInfo,
				false,
Thomas Witkowski's avatar
Thomas Witkowski committed
381
				elementMatrix, 
382
				*factorIt ? **factorIt : 1.0);	     
Thomas Witkowski's avatar
Thomas Witkowski committed
383 384
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
385 386 387
    if (factor != 1.0)
      elementMatrix *= factor;

Thomas Witkowski's avatar
Thomas Witkowski committed
388 389
    if (op || operators.size())
      addElementMatrix(elementMatrix, bound, rowElInfo, colElInfo);       
Thomas Witkowski's avatar
Thomas Witkowski committed
390 391
  }

392

Thomas Witkowski's avatar
Thomas Witkowski committed
393 394
  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
395 396
			    ElInfo *smallElInfo, ElInfo *largeElInfo,
			    const BoundaryType *bound, Operator *op)
Thomas Witkowski's avatar
Thomas Witkowski committed
397 398 399
  {
    FUNCNAME("DOFMatrix::assemble2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
400
    if (!op && operators.size() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
401 402
      return;

403
    set_to_zero(elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
404
    
405
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
406 407 408 409 410 411 412
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
413
      for(it = operators.begin(), factorIt = operatorFactor.begin();	
Thomas Witkowski's avatar
Thomas Witkowski committed
414 415 416
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
417 418 419 420 421
	  (*it)->getElementMatrix(mainElInfo, 
				  auxElInfo,
				  smallElInfo, 
				  largeElInfo,
				  rowFeSpace == colFeSpace,
422 423 424
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
425 426
      }      
    }
427

Thomas Witkowski's avatar
Thomas Witkowski committed
428 429 430 431
    if (factor != 1.0)
      elementMatrix *= factor;

    addElementMatrix(elementMatrix, bound, mainElInfo, NULL);       
432 433
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
434

Thomas Witkowski's avatar
Thomas Witkowski committed
435 436 437 438
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
439
	 it != operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
440 441 442
      (*it)->finishAssembling();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
443

444
  // Should work as before
445 446 447
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
448 449
    for (std::vector<Operator*>::iterator op = operators.begin(); 
	 op != operators.end(); ++op)
450
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
451

452 453 454
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
455

456
  void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
457
  {
458
    matrix += a * x.matrix + y.matrix;
459 460
  }

461

462 463
  void DOFMatrix::scal(double b) 
  {
464
    matrix *= b;
465 466
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
467

Thomas Witkowski's avatar
Thomas Witkowski committed
468 469 470 471 472 473 474
  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
475

476 477 478 479 480 481 482 483
  void DOFMatrix::clearOperators()
  {
    operators.clear();
    operatorFactor.clear();
    operatorEstFactor.clear();
  }


484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541
  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;
    }    
  }

542

543 544
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
545
    matrix = rhs.matrix;
546 547
  }

548

Thomas Witkowski's avatar
Thomas Witkowski committed
549
  void DOFMatrix::clearDirichletRows()
550
  {      
Thomas Witkowski's avatar
Thomas Witkowski committed
551
    FUNCNAME("DOFMatrix::clearDirichletRows()");
552

553 554 555
    // Do the following only in sequential code. In parallel mode, the specific
    // solver method must care about dirichlet boundary conditions.

556
    inserter_type &ins = *inserter;  
Thomas Witkowski's avatar
Thomas Witkowski committed
557 558
    for (std::set<int>::iterator it = dirichletDofs.begin(); 
	 it != dirichletDofs.end(); ++it)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
559 560 561 562
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      if (dofMap->isRankDof(*it))
#endif
	ins[*it][*it] = 1.0;    
563 564
  }

565

566 567
  int DOFMatrix::memsize() 
  {   
568 569
    return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type)
      + matrix.nnz() * sizeof(base_matrix_type::value_type);
570 571
  }

572

Thomas Witkowski's avatar
Thomas Witkowski committed
573 574 575 576 577 578 579 580
  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
581 582

    dirichletDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
583 584
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
585

586
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
587 588 589 590
  void DOFMatrix::setDofMap(FeSpaceDofMap& m)
  {
    dofMap = &m; 
  }
591
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
592

593
}