ProblemScal.cc 19.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#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"
#include "FileWriter.h"
#include "RefinementManager.h"
#include "CoarseningManager.h"
#include "Lagrange.h"
#include "PeriodicBC.h"
#include "ValueReader.h"
24
#include "ElementFileWriter.h"
25
#include "ProblemVec.h"
26
27
28

namespace AMDiS {

29
30
31
32
  ProblemScal::~ProblemScal()
  {
    FUNCNAME("ProblemScal::~ProblemScal()");

33
34
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++)
      delete fileWriters[i];
35
36

    if (systemMatrix)
37
      delete systemMatrix;
38
    if (rhs) 
39
      delete rhs;
40
    if (solution)
41
      delete solution;
42
    if (estimator)
43
      delete estimator;
44
    if (marker)
45
      delete marker;
46
    if (solver)
47
      delete solver;
48
    if (mesh)
49
      delete mesh;
50
51
52
53
54

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

55
56
57
  void ProblemScal::writeFiles(AdaptInfo *adaptInfo, bool force) 
  {
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++)
58
      fileWriters[i]->writeFiles(adaptInfo, force);
Thomas Witkowski's avatar
Thomas Witkowski committed
59
60
  }

61
62
63
64
65
  void ProblemScal::writeFiles(AdaptInfo &adaptInfo, bool force) 
  {
    writeFiles(&adaptInfo, force);
  }

66
67
  void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
68
    solution->interpol(fct);
69
70
71
72
73
74
  }

  void ProblemScal::addMatrixOperator(Operator *op, 
				      double *factor,
				      double *estFactor) 
  {
75
    systemMatrix->addOperator(op, factor, estFactor);
76
77
78
79
80
81
  }

  void ProblemScal::addVectorOperator(Operator *op, 
				      double *factor,
				      double *estFactor) 
  {
82
    rhs->addOperator(op, factor, estFactor);
83
84
85
86
87
  }
  
  void ProblemScal::addDirichletBC(BoundaryType type, 
				   AbstractFunction<double, WorldVector<double> >* b) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
88
    DirichletBC *dirichlet = new DirichletBC(type, b, feSpace);
89
90
91
92
93
94
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
95
96
97
98
99
100
  }

  void ProblemScal::addDirichletBC(BoundaryType type, 
				   DOFVector<double> *vec) 
  {
    DirichletBC *dirichlet = new DirichletBC(type, vec);
101
102
103
104
105
106
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
107
108
109
110
111
112
  }

  void ProblemScal::addRobinBC(BoundaryType type, 
			       AbstractFunction<double, WorldVector<double> > *n,
			       AbstractFunction<double, WorldVector<double> > *r)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
113
    RobinBC *robin = new RobinBC(type, n, r, feSpace);
114
115
116
117
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
118
119
120
121
122
  }

  void ProblemScal::addNeumannBC(BoundaryType type, 
				 AbstractFunction<double, WorldVector<double> > *n)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
123
    NeumannBC *neumann = new NeumannBC(type, n, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
124
    if (rhs) 
125
      rhs->getBoundaryManager()->addBoundaryCondition(neumann);
126
127
128
129
130
131
  }

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

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

143
144
145
146
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(periodic);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(periodic);
147
148
149
150
151
152
153
  }

  void ProblemScal::createMesh()
  {
    TEST_EXIT(Parameters::initialized())("parameters not initialized\n");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
156
    GET_PARAMETER(0, name + "->info", "%d", &info);
Thomas Witkowski's avatar
Thomas Witkowski committed
157
    GET_PARAMETER(0, name + "->mesh", &meshName);
158

159
    TEST_EXIT(meshName != "")("no mesh name specified\n");
160
161
162

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

    // create the mesh
167
    mesh = new Mesh(meshName, dim);
168

Thomas Witkowski's avatar
Thomas Witkowski committed
169
    switch (dim) {
170
    case 1:
171
172
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
173
174
      break;
    case 2:
175
176
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
177
178
      break;
    case 3:
179
180
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
181
182
183
184
185
186
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

187
188
  void ProblemScal::setMeshFromProblemVec(ProblemVec* pv, int i) 
  { 
189
    mesh = pv->getMesh(i);
190
191
    coarseningManager = pv->getCoarseningManager(i);
    refinementManager = pv->getRefinementManager(i);
192
193
  }

194
195
  Flag ProblemScal::markElements(AdaptInfo *adaptInfo) 
  { 
196
197
    if (marker)
      return marker->markMesh(adaptInfo, mesh);
198
199
    else
      WARNING("no marker\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
200

201
202
203
204
205
    return 0;
  }

  Flag ProblemScal::refineMesh(AdaptInfo *adaptInfo) 
  { 
206
    return refinementManager->refineMesh(mesh); 
207
208
209
210
  }

  Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    if (adaptInfo->isCoarseningAllowed(0)) {
212
      return coarseningManager->coarsenMesh(mesh);
Thomas Witkowski's avatar
Thomas Witkowski committed
213
    } else {
214
215
216
217
218
      WARNING("coarsening not allowed\n");
      return 0;
    }
  }

219
  void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
220
221
  {
    FUNCNAME("Problem::solve()");
222

223
    if (!solver) {
224
225
226
227
228
229
230
231
      WARNING("no solver\n");
      return;
    }

#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif
    clock_t first = clock();
232
233
234

    SolverMatrix<DOFMatrix> solverMatrix;
    solverMatrix.setMatrix(*systemMatrix);
235
    solver->solveSystem(solverMatrix, *solution, *rhs);
236
237

#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
238
    INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
239
240
241
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
243
244
245
		   TIME_USED(first, clock()));
#endif

246
    adaptInfo->setSolverIterations(solver->getIterations());
247
248
249
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
250
251
252
253
254
255
256
257
258
  }

  void ProblemScal::initialize(Flag initFlag,
			       ProblemScal *adoptProblem,
			       Flag adoptFlag) 
  {
    FUNCNAME("Problem::initialize()");

    // === create mesh ===
259
    if (mesh) { 
260
261
      WARNING("mesh already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
262
263
264
      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
265
266
267
	createMesh();
      } 
      
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
270
271
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_MESH) || 
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
272
273
	TEST_EXIT(!mesh)("mesh already created\n");
	mesh = adoptProblem->getMesh();
274
275
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
Thomas Witkowski's avatar
Thomas Witkowski committed
276
      }
277
278
    }

279
    if (!mesh) 
280
      WARNING("no mesh created\n");
281
282

    // === create fespace ===
Thomas Witkowski's avatar
Thomas Witkowski committed
283
    if (feSpace) {
284
285
      WARNING("feSpace already created\n");
    } else {
286
287
      if (initFlag.isSet(INIT_FE_SPACE) || 
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
288
289
	createFESpace();
      } 
290

291
292
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
	TEST_EXIT(!feSpace)("feSpace already created");
	feSpace = dynamic_cast<ProblemScal*>(adoptProblem)->getFESpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
295
      }
296
297
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
298
    if (!feSpace) 
299
      WARNING("no feSpace created\n");
300
301
302
303
304
305

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM)) {
      createMatricesAndVectors();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
306
307
308
309
310
311
      TEST_EXIT(!solution)("solution already created\n");
      TEST_EXIT(!rhs)("rhs already created\n");
      TEST_EXIT(!systemMatrix)("systemMatrix already created\n");
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRHS();
      systemMatrix  = adoptProblem->getSystemMatrix();
312
313
314
    }

    // === create solver ===
315
    if (solver) {
316
317
318
319
320
321
      WARNING("solver already created\n");
    } else {
      if (initFlag.isSet(INIT_SOLVER)) {
	createSolver();
      } 
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
322
323
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
324
325
326
      }
    }

327
    if (!solver) 
328
329
330
      WARNING("no solver created\n");

    // === create estimator ===
331
    if (estimator) {
332
333
      WARNING("estimator already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
334
      if (initFlag.isSet(INIT_ESTIMATOR)) {
335
336
	createEstimator();
      } 
Thomas Witkowski's avatar
Thomas Witkowski committed
337
      if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) {
338
339
	TEST_EXIT(!estimator)("estimator already created\n");
	estimator = adoptProblem->getEstimator();
340
341
342
      } 
    }

343
    if (!estimator) 
344
      WARNING("no estimator created\n");
345
346

    // === create marker ===
347
    if (marker) {
348
349
      WARNING("marker already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
350
      if (initFlag.isSet(INIT_MARKER))
351
	createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353

      if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) {
354
355
	TEST_EXIT(!marker)("marker already created\n");
	marker = adoptProblem->getMarker();
356
357
358
      } 
    }

359
    if (!marker) 
360
      WARNING("no marker created\n");
361
362

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
363
    if (initFlag.isSet(INIT_FILEWRITER))
364
      createFileWriter();
Thomas Witkowski's avatar
Thomas Witkowski committed
365
 
366
367
368
369
370
371
372
373
374
375

    // === 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
376
    std::string serializationFilename = "";
377
378
379
380
381
    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
382
      GET_PARAMETER(0, name + "->input->read serialization", "%d", 
383
384
		    &readSerialization);
      if (readSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
385
	GET_PARAMETER(0, name + "->input->serialization filename", 
386
387
388
389
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	std::ifstream in(serializationFilename.c_str());
391
392
393
	deserialize(in);
	in.close();
      } else {
394
395
	if (initFlag.isSet(INIT_MESH) && mesh && !mesh->isInitialized()) {
	  mesh->initialize();
396
397
398
	}
    
	// === read value file and use it for the mesh values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
399
	std::string valueFilename("");
400
	GET_PARAMETER(0, mesh->getName() + "->value file name", &valueFilename);     
401
	if (valueFilename.length()) {     
402
	  ValueReader::readValue(valueFilename,
403
404
405
406
				 mesh,
				 solution,
				 mesh->getMacroFileInfo());
	  mesh->clearMacroFileInfo();
407
	}
408
	
409
	// === do global refinements ===
410
411
	if (initFlag.isSet(INIT_GLOBAL_REFINES)) {
	  int globalRefinements = 0;
412
413

#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
414
	  GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinements);
415
416
#endif

417
	  refinementManager->globalRefine(mesh, globalRefinements);	
418
	}
419
420
421
422
423
424
425
426
      }
    }
  }

  void ProblemScal::createFESpace()
  {
    // create finite element space
    int degree = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
427
    GET_PARAMETER(0, name + "->polynomial degree" ,"%d", &degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
    feSpace = FiniteElemSpace::provideFESpace(NULL,
429
430
					       Lagrange::getLagrange(mesh->getDim(), degree),
					       mesh,
Thomas Witkowski's avatar
Thomas Witkowski committed
431
					       name + "->feSpace");  
432
433

    // create dof admin for vertex dofs if neccessary
434
435
    if (mesh->getNumberOfDOFs(VERTEX) == 0) {
      DimVec<int> ln_dof(mesh->getDim(), DEFAULT_VALUE, 0);
436
      ln_dof[VERTEX]= 1;
437
      mesh->createDOFAdmin("vertex dofs", ln_dof);
438
439
440
441
442
443
    }
  }

  void ProblemScal::createMatricesAndVectors()
  {
    // === create vectors and system matrix ===
444
445
446
    systemMatrix = new DOFMatrix(feSpace, feSpace, "system matrix");
    rhs = new DOFVector<double>(feSpace, "rhs");
    solution = new DOFVector<double>(feSpace, "solution");
447

448
449
    solution->setCoarsenOperation(COARSE_INTERPOL);
    solution->set(0.0);      /*  initialize u_h  !                      */
450
451
452
453
454
  }

  void ProblemScal::createSolver()
  {
    // === create solver ===
455
    std::string solverType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
456
    GET_PARAMETER(0, name + "->solver", &solverType);
457
458
    OEMSolverCreator *solverCreator = 
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
459
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
460
    solverCreator->setName(name + "->solver");
461
462
    solver = solverCreator->create();
    solver->initParameters();
463
464
465
466
467
  }

  void ProblemScal::createEstimator()
  {
    // create and set leaf data prototype
468
    mesh->setElementDataPrototype(new LeafDataEstimatable(new LeafDataCoarsenable));
469
470

    // === create estimator ===
471
    std::string estimatorType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
472
    GET_PARAMETER(0, name + "->estimator", &estimatorType);
473
474
475
    EstimatorCreator *estimatorCreator = 
      dynamic_cast<EstimatorCreator*>(
				      CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
476
    if (estimatorCreator) {
Thomas Witkowski's avatar
Thomas Witkowski committed
477
      estimatorCreator->setName(name + "->estimator");
Thomas Witkowski's avatar
Thomas Witkowski committed
478
      if (estimatorType == "recovery") {
479
	dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->setSolution(solution);
480
      }
481
      estimator = estimatorCreator->create();
482
483

      // init estimator
484
      estimator->addSystem(systemMatrix, solution, rhs);
485
486
487
488
489
    }
  }

  void ProblemScal::createMarker()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
490
    marker = dynamic_cast<Marker*>(Marker::createMarker(name + "->marker", -1));
491
492
493
494
  }

  void ProblemScal::createFileWriter()
  {
495
    fileWriters.push_back(new FileWriter(name + "->output", mesh, solution));
496
    int writeSerialization = 0;
497
498
499
500
501
    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
502
503
    if (writeSerialization)
      fileWriters.push_back(new Serializer<ProblemScal>(this));
504
#endif
505
506
507
508
509
510
  }

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

511
    if (estimator) {
512
      clock_t first = clock();
513
      estimator->estimate(adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
514
      INFO(info,8)("estimation of the error needed %.5f seconds\n",
515
516
		    TIME_USED(first,clock()));

517
      adaptInfo->setEstSum(estimator->getErrorSum(), 0);
518
519
520
      adaptInfo->setTimeEstSum(estimator->getTimeEst(), 0);
      adaptInfo->setEstMax(estimator->getErrorMax(), 0);
      adaptInfo->setTimeEstMax(estimator->getTimeEstMax(), 0);
521
522
523
524
525
    } else {
      WARNING("no estimator\n");
    }
  }

526
527
  void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				      bool assembleMatrix, bool assembleVector)
528
  {
529
    FUNCNAME("ProblemScal::buildAfterCoarsen()");
530
531
532

    clock_t first = clock();

533
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
534
    mesh->dofCompress();
Thomas Witkowski's avatar
Thomas Witkowski committed
535
#endif
536
537

    MSG("%d DOFs for %s\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
538
	feSpace->getAdmin()->getUsedDOFs(), feSpace->getName().c_str());
539

540
    rhs->set(0.0);
541

542
543
544
545
546
547
    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
548
      nnz_per_row= 5;
549
550
551
552

    // Correct dimensionality of matrix
    base_matrix.change_dim(feSpace->getAdmin()->getUsedSize(), 
			   feSpace->getAdmin()->getUsedSize());
553
    set_to_zero(base_matrix);
554
555
556
557
   
    // Reuse old sparsity information (if available) or default
    systemMatrix->startInsertion(nnz_per_row);

558
    // fill boundary conditions
559
560
561
562
563
564
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->initMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->initVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->initVector(solution);
565

566
    BoundaryType *bound;
567
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
568
569
570
571
572
573
    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);
574

575
    while (elInfo) {
576
      if (useGetBound) {
Thomas Witkowski's avatar
Thomas Witkowski committed
577
	bound = new BoundaryType[feSpace->getBasisFcts()->getNumber()];
Thomas Witkowski's avatar
Thomas Witkowski committed
578
	feSpace->getBasisFcts()->getBound(elInfo, bound);
579
580
581
      } else {
	bound = NULL;
      }
582

583
584
585
      systemMatrix->assemble(1.0, elInfo, bound);
      rhs->assemble(1.0, elInfo, bound);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
586
      if (useGetBound)
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	delete [] bound;
588
      
589
590
591
592
593
594
      if (systemMatrix->getBoundaryManager()) 
	systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
      if (rhs->getBoundaryManager())
	rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
      if (solution->getBoundaryManager())
	solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
595
596
597
      elInfo = stack.traverseNext(elInfo);
    }

598
    systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
599
    systemMatrix->finishInsertion();
600

601
    // TODO: ExitMatrix should be called after finishInsertion!
602
603
604
605
606
607
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->exitVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->exitVector(solution);
608

609
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first,clock()));
610
 
611
    createPrecon();
612
613
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
614
  void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, std::string name)
615
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
616
    FUNCNAME("ProblemScal::writeResidualMesh()");
617

Thomas Witkowski's avatar
Thomas Witkowski committed
618
    std::map<int, double> vec;    
619
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
    ElInfo *elInfo = stack.traverseFirst(this->getMesh(),  -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
622
623
    
    while (elInfo) {		  
Thomas Witkowski's avatar
Thomas Witkowski committed
624
      vec[elInfo->getElement()->getIndex()] = elInfo->getElement()->getEstimation(0);
625
626
627
      elInfo = stack.traverseNext(elInfo);
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
628
    ElementFileWriter fw(name, this->getFESpace(), vec);
629
630
631
    fw.writeFiles(adaptInfo, true);    
  }

632
633
  void ProblemScal::createPrecon()
  {
634
    std::string preconType("no");
635
636
637
638
639
640
641
642
643
644
645
646
647
648
    GET_PARAMETER(0, name + "->solver->left precon", &preconType);

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

    solver->setLeftPrecon( preconCreator->create(systemMatrix->getBaseMatrix()) );

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
649
  void ProblemScal::serialize(std::ostream &out) 
650
651
652
  {
    FUNCNAME("ProblemScal::serialize()");

653
654
    mesh->serialize(out);
    solution->serialize(out);
655
656
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
657
  void ProblemScal::deserialize(std::istream &in) 
658
659
660
  {
    FUNCNAME("ProblemScal::deserialize()");

661
662
    mesh->deserialize(in);    
    solution->deserialize(in);
663
664
665
  }

}