ProblemScal.cc 20.2 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
24
25
#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 "Preconditioner.h"
#include "MatVecMultiplier.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"
26
#include "ElementFileWriter.h"
27
28
29
30
31

namespace AMDiS {

  ProblemScal *ProblemScal::traversePtr_ = NULL;

32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
  ProblemScal::~ProblemScal()
  {
    FUNCNAME("ProblemScal::~ProblemScal()");

    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++) {
      DELETE fileWriters[i];
    }

    if (systemMatrix)
      DELETE systemMatrix;
    if (rhs) 
      DELETE rhs;
    if (solution)
      DELETE solution;
    if (estimator)
      DELETE estimator;
    if (marker)
      DELETE marker;
    if (solver)
      DELETE solver;
    if (leftPrecon)
      DELETE leftPrecon;
    if (rightPrecon)
      DELETE rightPrecon;
    if (mesh)
      DELETE mesh;

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

63
  void ProblemScal::writeFiles(AdaptInfo *adaptInfo, bool force) {
64
65
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++) {
      fileWriters[i]->writeFiles(adaptInfo, force);
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
    }
  }

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

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

  void ProblemScal::addVectorOperator(Operator *op, 
				      double *factor,
				      double *estFactor) 
  {
85
    rhs->addOperator(op, factor, estFactor);
86
87
88
89
90
91
  }
  
  void ProblemScal::addDirichletBC(BoundaryType type, 
				   AbstractFunction<double, WorldVector<double> >* b) 
  {
    DirichletBC *dirichlet = new DirichletBC(type, b, feSpace_);
92
93
94
95
96
97
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
98
99
100
101
102
103
  }

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

  void ProblemScal::addRobinBC(BoundaryType type, 
			       AbstractFunction<double, WorldVector<double> > *n,
			       AbstractFunction<double, WorldVector<double> > *r)
  {
    RobinBC *robin = new RobinBC(type, n, r, feSpace_);
117
118
119
120
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
121
122
123
124
125
126
  }

  void ProblemScal::addNeumannBC(BoundaryType type, 
				 AbstractFunction<double, WorldVector<double> > *n)
  {
    NeumannBC *neumann = new NeumannBC(type, n, feSpace_);
127
128
    if(rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(neumann);
129
130
131
132
133
134
135
  }

  void ProblemScal::addRobinBC(BoundaryType type, 
			       DOFVector<double> *n,
			       DOFVector<double> *r)
  {
    RobinBC *robin = new RobinBC(type, n, r, feSpace_);
136
137
138
139
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
140
141
142
143
144
145
  }

  void ProblemScal::addPeriodicBC(BoundaryType type) 
  {
    PeriodicBC *periodic = new PeriodicBC(type, feSpace_);

146
147
148
149
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(periodic);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(periodic);
150
151
152
153
154
155
156
  }

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

    // === create problems mesh ===
Thomas Witkowski's avatar
Thomas Witkowski committed
157
    std::string meshName("");
158
159
160
161

    GET_PARAMETER(0, name_ + "->info", "%d", &info_);
    GET_PARAMETER(0, name_ + "->mesh", &meshName);

162
    TEST_EXIT(meshName != "")("no mesh name specified\n");
163
164
165
166

    // get problem dimension
    int dim = 0;
    GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
167
    TEST_EXIT(dim)("no problem dimension specified!\n");
168
169

    // create the mesh
170
    mesh = NEW Mesh(meshName, dim);
171

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

190
191
192

  void ProblemScal::setMeshFromProblemVec(ProblemVec* pv, int i) 
  { 
193
    mesh = pv->getMesh(i);
194
195
196
197
198
    coarseningManager_ = pv->getCoarseningManager(i);
    refinementManager_ = pv->getRefinementManager(i);
  }


199
200
  Flag ProblemScal::markElements(AdaptInfo *adaptInfo) 
  { 
201
202
    if (marker)
      return marker->markMesh(adaptInfo, mesh);
203
204
    else
      WARNING("no marker\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
205

206
207
208
    return 0;
  }

209

210
211
  Flag ProblemScal::refineMesh(AdaptInfo *adaptInfo) 
  { 
212
    return refinementManager_->refineMesh(mesh); 
213
214
  }

215

216
217
  Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
218
    if (adaptInfo->isCoarseningAllowed(0)) {
219
      return coarseningManager_->coarsenMesh(mesh);
Thomas Witkowski's avatar
Thomas Witkowski committed
220
    } else {
221
222
223
224
225
      WARNING("coarsening not allowed\n");
      return 0;
    }
  }

226
  void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
227
228
  {
    FUNCNAME("Problem::solve()");
229
    if (!solver) {
230
231
232
233
234
235
236
237
238
      WARNING("no solver\n");
      return;
    }

#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

    clock_t first = clock();
239
    int iter = solver->solve(matVec_, solution, rhs, leftPrecon, rightPrecon); 
240
241
242
243
244
245
246
247
248
249
250

#ifdef _OPENMP
    INFO(info_, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
    INFO(info_, 8)("solution of discrete system needed %.5f seconds\n",
		   TIME_USED(first, clock()));
#endif

    adaptInfo->setSolverIterations(iter);
251
252
253
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
254
255
256
257
258
259
260
261
262
  }

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

    // === create mesh ===
263
    if (mesh) { 
264
265
      WARNING("mesh already created\n");
    } else {
266
      if (initFlag.isSet(CREATE_MESH) || 
267
	 ((!adoptFlag.isSet(INIT_MESH))&&
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
270
271
	  (initFlag.isSet(INIT_SYSTEM)||initFlag.isSet(INIT_FE_SPACE)))) {
	createMesh();
      } 
      
272
      if (adoptProblem && 
273
274
	 (adoptFlag.isSet(INIT_MESH) || 
	  adoptFlag.isSet(INIT_SYSTEM) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
275
	  adoptFlag.isSet(INIT_FE_SPACE))) {
276
277
	TEST_EXIT(!mesh)("mesh already created\n");
	mesh = adoptProblem->getMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
280
	refinementManager_ = adoptProblem->refinementManager_;
	coarseningManager_ = adoptProblem->coarseningManager_;
      }
281
282
    }

283
    if (!mesh) 
284
      WARNING("no mesh created\n");
285
286

    // === create fespace ===
287
    if (feSpace_) {
288
289
      WARNING("feSpace already created\n");
    } else {
290
      if (initFlag.isSet(INIT_FE_SPACE) || (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
291
292
	createFESpace();
      } 
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
295
296
      if (adoptProblem && (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
	TEST_EXIT(!feSpace_)("feSpace already created");
	feSpace_ = dynamic_cast<ProblemScal*>(adoptProblem)->getFESpace();
      }
297
298
    }

299
300
    if (!feSpace_) 
      WARNING("no feSpace created\n");
301
302
303
304
305
306

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM)) {
      createMatricesAndVectors();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
307
308
309
310
311
312
      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();
313
314
315
    }

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

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

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

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

    // === create marker ===
348
    if (marker) {
349
350
      WARNING("marker already created\n");
    } else {
351
      if (initFlag.isSet(INIT_MARKER)) {
352
353
	createMarker();
      } 
354
      if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) {
355
356
	TEST_EXIT(!marker)("marker already created\n");
	marker = adoptProblem->getMarker();
357
358
359
      } 
    }

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

    // === create file writer ===
364
    if (initFlag.isSet(INIT_FILEWRITER)) {
365
366
367
368
369
370
371
372
373
374
375
376
377
378
      createFileWriter();
    }

  

    // === 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
379
    std::string serializationFilename = "";
380
381
382
383
384
385
386
387
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("")) {
      GET_PARAMETER(0, name_ + "->input->read serialization", "%d", 
		    &readSerialization);
      if (readSerialization) {
	GET_PARAMETER(0, name_ + "->input->serialization filename", 
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

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

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

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

  void ProblemScal::createMatricesAndVectors()
  {
    // === create vectors and system matrix ===
446
447
448
    systemMatrix = NEW DOFMatrix(feSpace_, feSpace_, "system matrix");
    rhs = NEW DOFVector<double>(feSpace_, "rhs");
    solution = NEW DOFVector<double>(feSpace_, "solution");
449

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

    // === create matVec ===
454
    matVec_ = NEW StandardMatVec<DOFMatrix, DOFVector<double> >(systemMatrix);
455
456
457
458
459
  }

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

    // === create preconditioners ===
Thomas Witkowski's avatar
Thomas Witkowski committed
470
    std::string preconType("no");
471
472
473
474
475
476
477
478
479
480
481
482
483
    Preconditioner<DOFVector<double> > *precon;
    GET_PARAMETER(0, name_ + "->solver->left precon", &preconType);
    CreatorInterface<PreconditionerScal> *preconCreator = 
      CreatorMap<PreconditionerScal>::getCreator(preconType);

    if (!preconCreator->isNullCreator()) {
      dynamic_cast<PreconditionerScalCreator*>(preconCreator)->setSizeAndRow(1, 0);
      dynamic_cast<PreconditionerScalCreator*>(preconCreator)->setName(name_ + "->solver->left precon");
    }
  
    precon = preconCreator->create();

    if (precon) {
484
485
      dynamic_cast<PreconditionerScal*>(precon)->setMatrix(&systemMatrix);
      leftPrecon = precon;
486
487
488
489
490
491
492
493
494
495
496
497
498
    }

    preconType.assign("no");
    GET_PARAMETER(0, name_ + "->solver->right precon", &preconType);
    preconCreator = CreatorMap<PreconditionerScal>::getCreator(preconType);

    if (!preconCreator->isNullCreator()) {
      dynamic_cast<PreconditionerScalCreator*>(preconCreator)->setSizeAndRow(1, 0);
      dynamic_cast<PreconditionerScalCreator*>(preconCreator)->setName(name_ + "->solver->left precon");
    }

    precon = preconCreator->create();
    if (precon) {
499
500
      dynamic_cast<PreconditionerScal*>(precon)->setMatrix(&systemMatrix);
      rightPrecon = precon;
501
502
503
    }

    // === create vector creator ===
504
    solver->setVectorCreator(new DOFVector<double>::Creator(feSpace_));
505
506
507
508
509
510

  }

  void ProblemScal::createEstimator()
  {
    // create and set leaf data prototype
511
    mesh->setElementDataPrototype(NEW LeafDataEstimatable(NEW LeafDataCoarsenable));
512
513

    // === create estimator ===
Thomas Witkowski's avatar
Thomas Witkowski committed
514
    std::string estimatorType("no");
515
516
517
518
    GET_PARAMETER(0, name_ + "->estimator", &estimatorType);
    EstimatorCreator *estimatorCreator = 
      dynamic_cast<EstimatorCreator*>(
				      CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
519
    if (estimatorCreator) {
520
      estimatorCreator->setName(name_ + "->estimator");
Thomas Witkowski's avatar
Thomas Witkowski committed
521
      if (estimatorType == "recovery") {
522
	dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->setSolution(solution);
523
      }
524
      estimator = estimatorCreator->create();
525
526

      // init estimator
527
      estimator->addSystem(systemMatrix, solution, rhs);
528
529
530
531
532
    }
  }

  void ProblemScal::createMarker()
  {
533
    marker = dynamic_cast<Marker*>(Marker::createMarker(name_ + "->marker", -1));
534
535
536
537
  }

  void ProblemScal::createFileWriter()
  {
538
    fileWriters.push_back(NEW FileWriter(name_ + "->output", mesh, solution));
539
540
541
    int writeSerialization = 0;
    GET_PARAMETER(0, name_ + "->output->write serialization", "%d", 
		  &writeSerialization);
Thomas Witkowski's avatar
Thomas Witkowski committed
542
    if (writeSerialization) {
543
      fileWriters.push_back(NEW Serializer<ProblemScal>(this));
544
545
546
547
548
549
550
    }
  }

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

551
    if (estimator) {
552
      clock_t first = clock();
553
      estimator->estimate(adaptInfo->getTimestep());
554
555
556
      INFO(info_,8)("estimation of the error needed %.5f seconds\n",
		    TIME_USED(first,clock()));

557
      adaptInfo->setEstSum(estimator->getErrorSum(), 0);
558
559

      adaptInfo->
560
	setTimeEstSum(estimator->getTimeEst(), 0);
561
562

      adaptInfo->
563
	setEstMax(estimator->getErrorMax(), 0);
564
565

      adaptInfo->
566
	setTimeEstMax(estimator->getTimeEstMax(), 0);
567
568
569
570
571
572

    } else {
      WARNING("no estimator\n");
    }
  }

573
574
  void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				      bool assembleMatrix, bool assembleVector)
575
  {
576
    FUNCNAME("ProblemScal::buildAfterCoarsen()");
577
578
579

    clock_t first = clock();

580
    mesh->dofCompress();
581
582
583
584
585
586
587

    MSG("%d DOFs for %s\n", 
	feSpace_->getAdmin()->getUsedSize(), 
	feSpace_->getName().c_str());

    Flag assembleFlag = 
      flag | 
588
589
      systemMatrix->getAssembleFlag() | 
      rhs->getAssembleFlag();
590
591
592
593

    if (useGetBound_)
      assembleFlag |= Mesh::FILL_BOUND;

594
595
    systemMatrix->clear();
    rhs->set(0.0);
596
597

    traversePtr_ = this;
598
    mesh->traverse(-1, assembleFlag, &buildAfterCoarsenFct);
599
600

    // fill boundary conditions
601
602
603
604
605
606
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->initMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->initVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->initVector(solution);
607
608

    TraverseStack stack;
609
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
610
611
612
613
614
615
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_BOUND |
					 Mesh::FILL_COORDS |
					 Mesh::FILL_DET |
					 Mesh::FILL_GRD_LAMBDA |
					 Mesh::FILL_NEIGH);
Thomas Witkowski's avatar
Thomas Witkowski committed
616

617
618
    // for all elements ...
    while (elInfo) {
619
620
621
622
623
624
      if (systemMatrix->getBoundaryManager()) 
	systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
      if (rhs->getBoundaryManager())
	rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
      if (solution->getBoundaryManager())
	solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
625
626
627
      elInfo = stack.traverseNext(elInfo);
    }

628
    systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
Thomas Witkowski's avatar
Thomas Witkowski committed
629

630
631
632
633
634
635
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->exitVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->exitVector(solution);
636
637
638
639
640
641
642
643
644

    INFO(info_, 8)("buildAfterCoarsen needed %.5f seconds\n",
		   TIME_USED(first,clock()));

    return;
  }

  int ProblemScal::buildAfterCoarsenFct(ElInfo *elInfo) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
645
    BoundaryType *bound;
646

647
    if (traversePtr_->getBoundUsed()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
648
649
      bound = GET_MEMORY(BoundaryType, traversePtr_->getFESpace()->getBasisFcts()->getNumber());
      traversePtr_->getFESpace()->getBasisFcts()->getBound(elInfo, bound);
650
    } else {
651
      bound = NULL;
652
    }
653
654
655
656
  
    traversePtr_->getSystemMatrix()->assemble(1.0, elInfo, bound);
    traversePtr_->getRHS()->assemble(1.0, elInfo, bound);

Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
659
660
    if (traversePtr_->getBoundUsed()) {
      FREE_MEMORY(bound, BoundaryTyppe, traversePtr_->getFESpace()->getBasisFcts()->getNumber());
    }

661
662
663
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
664
  void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
665
666
667
668
  {
    FUNCNAME("ProblemVec::writeResidualMesh()");

    Mesh *mesh = this->getMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
    FiniteElemSpace *fe = this->getFESpace();   
    std::map<int, double> vec;    
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh,
					 -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS);
    
    while (elInfo) {		  
      Element *el = elInfo->getElement();
      double lError = el->getEstimation(0);
      
      vec[elInfo->getElement()->getIndex()] = lError;
      elInfo = stack.traverseNext(elInfo);
    }
    
    ElementFileWriter fw(name, mesh, fe, vec);
    fw.writeFiles(adaptInfo, true);    
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
689
  void ProblemScal::serialize(std::ostream &out) 
690
691
692
  {
    FUNCNAME("ProblemScal::serialize()");

693
694
    mesh->serialize(out);
    solution->serialize(out);
695
696
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
697
  void ProblemScal::deserialize(std::istream &in) 
698
699
700
  {
    FUNCNAME("ProblemScal::deserialize()");

701
702
    mesh->deserialize(in);    
    solution->deserialize(in);
703
704
705
  }

}