ProblemScal.cc 19.7 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
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include "ProblemScal.h"
#include "AbstractFunction.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "FixVec.h"
#include "Flag.h"
#include "Serializer.h"
#include "RecoveryEstimator.h"
#include "Operator.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "OEMSolver.h"
#include "DOFVector.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "ElInfo.h"
30
#include "io/FileWriter.h"
31
32
33
34
#include "RefinementManager.h"
#include "CoarseningManager.h"
#include "Lagrange.h"
#include "PeriodicBC.h"
35
36
#include "io/ValueReader.h"
#include "io/ElementFileWriter.h"
37
#include "ProblemVec.h"
38
#include "Debug.h"
39
40
41

namespace AMDiS {

42
43
44
45
  ProblemScal::~ProblemScal()
  {
    FUNCNAME("ProblemScal::~ProblemScal()");

46
47
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++)
      delete fileWriters[i];
48
49

    if (systemMatrix)
50
      delete systemMatrix;
51
    if (rhs) 
52
      delete rhs;
53
    if (solution)
54
      delete solution;
55
    if (estimator)
56
      delete estimator;
57
    if (marker)
58
      delete marker;
59
    if (solver)
60
      delete solver;
61
    if (mesh)
62
      delete mesh;
63
64
65
66
67

    FiniteElemSpace::clear();
    Lagrange::clear();
  }

68
69
70
  void ProblemScal::writeFiles(AdaptInfo *adaptInfo, bool force) 
  {
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++)
71
      fileWriters[i]->writeFiles(adaptInfo, force);
Thomas Witkowski's avatar
Thomas Witkowski committed
72
73
  }

74
75
76
77
78
  void ProblemScal::writeFiles(AdaptInfo &adaptInfo, bool force) 
  {
    writeFiles(&adaptInfo, force);
  }

79
80
  void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
81
    solution->interpol(fct);
82
83
  }

84
  void ProblemScal::addMatrixOperator(Operator *op, double *factor, double *estFactor) 
85
  {
86
    systemMatrix->addOperator(op, factor, estFactor);
87
88
  }

89
90
91
92
93
94
  void ProblemScal::addMatrixOperator(Operator &op, double *factor, double *estFactor) 
  {
    systemMatrix->addOperator(&op, factor, estFactor);
  }

  void ProblemScal::addVectorOperator(Operator *op, double *factor, double *estFactor) 
95
  {
96
    rhs->addOperator(op, factor, estFactor);
97
  }
98
99
100
101
102

  void ProblemScal::addVectorOperator(Operator &op, double *factor, double *estFactor) 
  {
    rhs->addOperator(&op, factor, estFactor);
  }
103
104
105
106
  
  void ProblemScal::addDirichletBC(BoundaryType type, 
				   AbstractFunction<double, WorldVector<double> >* b) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
107
    DirichletBC *dirichlet = new DirichletBC(type, b, feSpace);
108
109
110
111
112
113
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
114
115
  }

116
  void ProblemScal::addDirichletBC(BoundaryType type, DOFVector<double> *vec) 
117
118
  {
    DirichletBC *dirichlet = new DirichletBC(type, vec);
119
120
121
122
123
124
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
125
126
127
128
129
130
  }

  void ProblemScal::addRobinBC(BoundaryType type, 
			       AbstractFunction<double, WorldVector<double> > *n,
			       AbstractFunction<double, WorldVector<double> > *r)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
131
    RobinBC *robin = new RobinBC(type, n, r, feSpace);
132
133
134
135
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
136
137
138
139
140
  }

  void ProblemScal::addNeumannBC(BoundaryType type, 
				 AbstractFunction<double, WorldVector<double> > *n)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
141
    NeumannBC *neumann = new NeumannBC(type, n, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
142
    if (rhs) 
143
      rhs->getBoundaryManager()->addBoundaryCondition(neumann);
144
145
146
  }

  void ProblemScal::addRobinBC(BoundaryType type, 
147
			       DOFVector<double> *n, DOFVector<double> *r)
148
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    RobinBC *robin = new RobinBC(type, n, r, feSpace);
150
151
152
153
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
154
155
156
157
  }

  void ProblemScal::addPeriodicBC(BoundaryType type) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
158
    PeriodicBC *periodic = new PeriodicBC(type, feSpace);
159

160
161
162
163
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(periodic);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(periodic);
164
165
  }

166

167
168
169
170
171
  void ProblemScal::createMesh()
  {
    TEST_EXIT(Parameters::initialized())("parameters not initialized\n");

    // === create problems mesh ===
Thomas Witkowski's avatar
Thomas Witkowski committed
172
    std::string meshName("");
173

Thomas Witkowski's avatar
Thomas Witkowski committed
174
    GET_PARAMETER(0, name + "->info", "%d", &info);
Thomas Witkowski's avatar
Thomas Witkowski committed
175
    GET_PARAMETER(0, name + "->mesh", &meshName);
176

177
    TEST_EXIT(meshName != "")("no mesh name specified\n");
178
179
180

    // get problem dimension
    int dim = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
182
    TEST_EXIT(dim)("no problem dimension specified!\n");
183
184

    // create the mesh
185
    mesh = new Mesh(meshName, dim);
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
    switch (dim) {
188
    case 1:
189
190
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
191
192
      break;
    case 2:
193
194
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
195
196
      break;
    case 3:
197
198
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
199
200
201
202
203
204
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

205
206
  void ProblemScal::setMeshFromProblemVec(ProblemVec* pv, int i) 
  { 
207
    mesh = pv->getMesh(i);
208
209
    coarseningManager = pv->getCoarseningManager(i);
    refinementManager = pv->getRefinementManager(i);
210
211
  }

212
213
  Flag ProblemScal::markElements(AdaptInfo *adaptInfo) 
  { 
214
215
    if (marker)
      return marker->markMesh(adaptInfo, mesh);
216
217
    else
      WARNING("no marker\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
218

219
220
221
222
223
    return 0;
  }

  Flag ProblemScal::refineMesh(AdaptInfo *adaptInfo) 
  { 
224
    return refinementManager->refineMesh(mesh); 
225
226
227
228
  }

  Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    if (adaptInfo->isCoarseningAllowed(0)) {
230
      return coarseningManager->coarsenMesh(mesh);
Thomas Witkowski's avatar
Thomas Witkowski committed
231
    } else {
232
233
234
235
236
      WARNING("coarsening not allowed\n");
      return 0;
    }
  }

237
  void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
238
239
  {
    FUNCNAME("Problem::solve()");
240

241
    if (!solver) {
242
243
244
245
      WARNING("no solver\n");
      return;
    }

246
247
    clock_t first = clock();    

248
249
    SolverMatrix<DOFMatrix> solverMatrix;
    solverMatrix.setMatrix(*systemMatrix);
250
    solver->solveSystem(solverMatrix, *solution, *rhs);
251

Thomas Witkowski's avatar
Thomas Witkowski committed
252
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
253
254
		   TIME_USED(first, clock()));

255
    adaptInfo->setSolverIterations(solver->getIterations());
256
257
258
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
259
260
261
262
263
264
265
266
  }

  void ProblemScal::initialize(Flag initFlag,
			       ProblemScal *adoptProblem,
			       Flag adoptFlag) 
  {
    FUNCNAME("Problem::initialize()");
    // === create mesh ===
267
    if (mesh) { 
268
269
      WARNING("mesh already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
272
      if (initFlag.isSet(CREATE_MESH) || 
	  ((!adoptFlag.isSet(INIT_MESH))&&
	   (initFlag.isSet(INIT_SYSTEM)||initFlag.isSet(INIT_FE_SPACE)))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
275
	createMesh();
      } 
      
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
279
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_MESH) || 
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
280
281
	TEST_EXIT(!mesh)("mesh already created\n");
	mesh = adoptProblem->getMesh();
282
283
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
Thomas Witkowski's avatar
Thomas Witkowski committed
284
      }
285
286
    }

287
    if (!mesh) 
288
      WARNING("no mesh created\n");
289

290
291
292
    if (refinementManager == NULL)
	WARNING("refinement manager not set\n");

293
    // === create fespace ===
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    if (feSpace) {
295
296
      WARNING("feSpace already created\n");
    } else {
297
298
      if (initFlag.isSet(INIT_FE_SPACE) || 
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
299
	createFeSpace();
300
      } 
301

302
303
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
304
	TEST_EXIT(!feSpace)("feSpace already created");
305
	feSpace = dynamic_cast<ProblemScal*>(adoptProblem)->getFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
306
      }
307
308
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
309
    if (!feSpace) 
310
      WARNING("no feSpace created\n");
311
312
313
314
315
316

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM)) {
      createMatricesAndVectors();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
317
318
319
320
      TEST_EXIT(!solution)("solution already created\n");
      TEST_EXIT(!rhs)("rhs already created\n");
      TEST_EXIT(!systemMatrix)("systemMatrix already created\n");
      solution = adoptProblem->getSolution();
321
      rhs = adoptProblem->getRhs();
322
      systemMatrix  = adoptProblem->getSystemMatrix();
323
324
325
    }

    // === create solver ===
326
    if (solver) {
327
328
329
330
331
332
      WARNING("solver already created\n");
    } else {
      if (initFlag.isSet(INIT_SOLVER)) {
	createSolver();
      } 
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
333
334
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
335
336
337
      }
    }

338
    if (!solver) 
339
340
341
      WARNING("no solver created\n");

    // === create estimator ===
342
    if (estimator) {
343
344
      WARNING("estimator already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
345
      if (initFlag.isSet(INIT_ESTIMATOR)) {
346
347
	createEstimator();
      } 
Thomas Witkowski's avatar
Thomas Witkowski committed
348
      if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) {
349
350
	TEST_EXIT(!estimator)("estimator already created\n");
	estimator = adoptProblem->getEstimator();
351
352
353
      } 
    }

354
    if (!estimator) 
355
      WARNING("no estimator created\n");
356
357

    // === create marker ===
358
    if (marker) {
359
360
      WARNING("marker already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
361
      if (initFlag.isSet(INIT_MARKER))
362
	createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
363
364

      if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) {
365
366
	TEST_EXIT(!marker)("marker already created\n");
	marker = adoptProblem->getMarker();
367
368
369
      } 
    }

370
    if (!marker) 
371
      WARNING("no marker created\n");
372
373

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
374
    if (initFlag.isSet(INIT_FILEWRITER))
375
      createFileWriter();
Thomas Witkowski's avatar
Thomas Witkowski committed
376
 
377
378
379
380
381
382
383
384
385
386

    // === read serialization and init mesh ===

    // There are two possiblities where the user can define a serialization
    // to be read from disk. Either by providing the parameter -rs when executing
    // the program or in the init file. The -rs parameter is always checked first,
    // because it can be added automatically when  rescheduling the program
    // before timeout of the runqueue.

    int readSerialization = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
387
    std::string serializationFilename = "";
388
389
390
391
392
    GET_PARAMETER(0, "argv->rs", &serializationFilename);

    // If the parameter -rs is set, we do nothing here, because the problem will be
    // deserialized in the constructor of a following AdaptInstationary initialization.
    if (!serializationFilename.compare("")) {
Thomas Witkowski's avatar
Thomas Witkowski committed
393
      GET_PARAMETER(0, name + "->input->read serialization", "%d", 
394
395
		    &readSerialization);
      if (readSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
396
	GET_PARAMETER(0, name + "->input->serialization filename", 
397
398
399
400
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
401
	std::ifstream in(serializationFilename.c_str());
402
403
404
	deserialize(in);
	in.close();
      } else {
405
406
	if (initFlag.isSet(INIT_MESH) && mesh && !mesh->isInitialized()) {
	  mesh->initialize();
407
408
409
	}
    
	// === read value file and use it for the mesh values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
410
	std::string valueFilename("");
411
	GET_PARAMETER(0, mesh->getName() + "->value file name", &valueFilename);     
412
	if (valueFilename.length())
413
	  ValueReader::readValue(valueFilename,
414
415
416
				 mesh,
				 solution,
				 mesh->getMacroFileInfo());
417
	
418
	// === do global refinements ===
419
420
	if (initFlag.isSet(INIT_GLOBAL_REFINES)) {
	  int globalRefinements = 0;
421
422

#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
423
	  GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinements);
424
425
#endif

426
	  refinementManager->globalRefine(mesh, globalRefinements);	
427
	}
428
429
      }
    }
430
	
431
432
  }

433
  void ProblemScal::createFeSpace()
434
435
436
  {
    // create finite element space
    int degree = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
437
    GET_PARAMETER(0, name + "->polynomial degree" ,"%d", &degree);
438
    feSpace = FiniteElemSpace::provideFeSpace(NULL,
439
440
					       Lagrange::getLagrange(mesh->getDim(), degree),
					       mesh,
Thomas Witkowski's avatar
Thomas Witkowski committed
441
					       name + "->feSpace");  
442
443

    // create dof admin for vertex dofs if neccessary
444
    if (mesh->getNumberOfDofs(VERTEX) == 0) {
445
      DimVec<int> ln_dof(mesh->getDim(), DEFAULT_VALUE, 0);
446
      ln_dof[VERTEX]= 1;
447
      mesh->createDOFAdmin("vertex dofs", ln_dof);
448
449
450
451
452
453
    }
  }

  void ProblemScal::createMatricesAndVectors()
  {
    // === create vectors and system matrix ===
454
455
456
    systemMatrix = new DOFMatrix(feSpace, feSpace, "system matrix");
    rhs = new DOFVector<double>(feSpace, "rhs");
    solution = new DOFVector<double>(feSpace, "solution");
457

458
459
    solution->setCoarsenOperation(COARSE_INTERPOL);
    solution->set(0.0);      /*  initialize u_h  !                      */
460
461
462
463
464
  }

  void ProblemScal::createSolver()
  {
    // === create solver ===
465
    std::string solverType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
466
    GET_PARAMETER(0, name + "->solver", &solverType);
467
468
    OEMSolverCreator *solverCreator = 
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
469
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
470
    solverCreator->setName(name + "->solver");
471
472
    solver = solverCreator->create();
    solver->initParameters();
473
474
475
476
477
  }

  void ProblemScal::createEstimator()
  {
    // create and set leaf data prototype
478
    mesh->setElementDataPrototype(new LeafDataEstimatable(new LeafDataCoarsenable));
479
480

    // === create estimator ===
481
    std::string estimatorType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
482
    GET_PARAMETER(0, name + "->estimator", &estimatorType);
483
484
485
    EstimatorCreator *estimatorCreator = 
      dynamic_cast<EstimatorCreator*>(
				      CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    if (estimatorCreator) {
Thomas Witkowski's avatar
Thomas Witkowski committed
487
      estimatorCreator->setName(name + "->estimator");
Thomas Witkowski's avatar
Thomas Witkowski committed
488
      if (estimatorType == "recovery") {
489
	dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->setSolution(solution);
490
      }
491
      estimator = estimatorCreator->create();
492
493

      // init estimator
494
      estimator->addSystem(systemMatrix, solution, rhs);
495
496
497
498
499
    }
  }

  void ProblemScal::createMarker()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
500
    marker = dynamic_cast<Marker*>(Marker::createMarker(name + "->marker", -1));
501
502
503
504
  }

  void ProblemScal::createFileWriter()
  {
505
    fileWriters.push_back(new FileWriter(name + "->output", mesh, solution));
506
    int writeSerialization = 0;
507
508
509
510
511
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);

    // Serialization is not allowed to be done by the problem, if its part of a parallel
    // problem definition. Than, the parallel problem object must be serialized.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
512
513
    if (writeSerialization)
      fileWriters.push_back(new Serializer<ProblemScal>(this));
514
#endif
515
516
517
518
519
520
  }

  void ProblemScal::estimate(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("Problem::estimate()");

521
    if (estimator) {
522
      clock_t first = clock();
523
      estimator->estimate(adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
524
      INFO(info,8)("estimation of the error needed %.5f seconds\n",
525
526
		    TIME_USED(first,clock()));

527
      adaptInfo->setEstSum(estimator->getErrorSum(), 0);
528
529
530
      adaptInfo->setTimeEstSum(estimator->getTimeEst(), 0);
      adaptInfo->setEstMax(estimator->getErrorMax(), 0);
      adaptInfo->setTimeEstMax(estimator->getTimeEstMax(), 0);
531
532
533
534
535
    } else {
      WARNING("no estimator\n");
    }
  }

536
537
  void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				      bool assembleMatrix, bool assembleVector)
538
  {
539
    FUNCNAME("ProblemScal::buildAfterCoarsen()");
540
541
542

    clock_t first = clock();

543
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
544
    mesh->dofCompress();
Thomas Witkowski's avatar
Thomas Witkowski committed
545
#endif
546
547

    MSG("%d DOFs for %s\n", 
548
	feSpace->getAdmin()->getUsedDofs(), feSpace->getName().c_str());
549

550
    rhs->set(0.0);
551

552
553
554
555
556
557
    DOFMatrix::base_matrix_type& base_matrix= systemMatrix->getBaseMatrix(); 

    int nnz_per_row= 0;
    if (num_rows(base_matrix) != 0)
      nnz_per_row= int(double(base_matrix.nnz()) / num_rows(base_matrix) * 1.2); 
    if (nnz_per_row < 5) 
Thomas Witkowski's avatar
Thomas Witkowski committed
558
      nnz_per_row= 5;
559
560
561
562

    // Correct dimensionality of matrix
    base_matrix.change_dim(feSpace->getAdmin()->getUsedSize(), 
			   feSpace->getAdmin()->getUsedSize());
563
    set_to_zero(base_matrix);
564
565
566
567
   
    // Reuse old sparsity information (if available) or default
    systemMatrix->startInsertion(nnz_per_row);

568
    // fill boundary conditions
569
570
571
572
573
574
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->initMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->initVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->initVector(solution);
575

576
    BoundaryType *bound;
577
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
578
579
580
581
582
583
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | 
					           Mesh::FILL_BOUND |
					           Mesh::FILL_COORDS |
					           Mesh::FILL_DET |
					           Mesh::FILL_GRD_LAMBDA |
					           Mesh::FILL_NEIGH);
584

585
    while (elInfo) {
586
      if (useGetBound) {
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	bound = new BoundaryType[feSpace->getBasisFcts()->getNumber()];
Thomas Witkowski's avatar
Thomas Witkowski committed
588
	feSpace->getBasisFcts()->getBound(elInfo, bound);
589
590
591
      } else {
	bound = NULL;
      }
592

593
594
595
      systemMatrix->assemble(1.0, elInfo, bound);
      rhs->assemble(1.0, elInfo, bound);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
596
      if (useGetBound)
Thomas Witkowski's avatar
Thomas Witkowski committed
597
	delete [] bound;
598
      
599
600
601
602
603
604
      if (systemMatrix->getBoundaryManager()) 
	systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
      if (rhs->getBoundaryManager())
	rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
      if (solution->getBoundaryManager())
	solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
605
606
607
      elInfo = stack.traverseNext(elInfo);
    }

608
    systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
609
    systemMatrix->finishInsertion();
610

611
    // TODO: ExitMatrix should be called after finishInsertion!
612
613
614
615
616
617
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->exitVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->exitVector(solution);
618

619
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first,clock()));
620
 
621
    createPrecon();
622
623
  }

624

Thomas Witkowski's avatar
Thomas Witkowski committed
625
  void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, std::string name)
626
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
627
    FUNCNAME("ProblemScal::writeResidualMesh()");
628

Thomas Witkowski's avatar
Thomas Witkowski committed
629
    std::map<int, double> vec;    
630
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
631
632
    ElInfo *elInfo = stack.traverseFirst(this->getMesh(),  -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
633
634
    
    while (elInfo) {		  
Thomas Witkowski's avatar
Thomas Witkowski committed
635
      vec[elInfo->getElement()->getIndex()] = elInfo->getElement()->getEstimation(0);
636
637
638
      elInfo = stack.traverseNext(elInfo);
    }
    
639
    ElementFileWriter fw(name, this->getMesh(), vec);
640
641
642
    fw.writeFiles(adaptInfo, true);    
  }

643

644
645
  void ProblemScal::createPrecon()
  {
646
    std::string preconType("no");
647
648
649
650
651
    GET_PARAMETER(0, name + "->solver->left precon", &preconType);

    CreatorInterface<ITL_BasePreconditioner> *preconCreator = 
      CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);

652
    solver->setLeftPrecon(preconCreator->create(systemMatrix->getBaseMatrix()));
653
654
655
656
657

    preconType= "no";
    GET_PARAMETER(0, name + "->solver->right precon", &preconType);

    preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
658
    solver->setRightPrecon(preconCreator->create(systemMatrix->getBaseMatrix()));
659
660
  }

661

Thomas Witkowski's avatar
Thomas Witkowski committed
662
  void ProblemScal::serialize(std::ostream &out) 
663
664
665
  {
    FUNCNAME("ProblemScal::serialize()");

666
667
    mesh->serialize(out);
    solution->serialize(out);
668
669
  }

670

Thomas Witkowski's avatar
Thomas Witkowski committed
671
  void ProblemScal::deserialize(std::istream &in) 
672
673
674
  {
    FUNCNAME("ProblemScal::deserialize()");

675
676
    mesh->deserialize(in);    
    solution->deserialize(in);
677
678
679
  }

}