DOFMatrix.cc 15.2 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
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
560
561
      //      if (dofMap->isRankDof(*it))
      if (dofMap->isSet(*it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
562
563
#endif
	ins[*it][*it] = 1.0;    
564
565
  }

566

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

573

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
586

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

594
}