Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

ProblemScal.cc 19.5 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
  void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
63
    solution->interpol(fct);
64
65
66
67
68
69
  }

  void ProblemScal::addMatrixOperator(Operator *op, 
				      double *factor,
				      double *estFactor) 
  {
70
    systemMatrix->addOperator(op, factor, estFactor);
71
72
73
74
75
76
  }

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

  void ProblemScal::addDirichletBC(BoundaryType type, 
				   DOFVector<double> *vec) 
  {
    DirichletBC *dirichlet = new DirichletBC(type, vec);
96
97
98
99
100
101
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
102
103
104
105
106
107
  }

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

  void ProblemScal::addNeumannBC(BoundaryType type, 
				 AbstractFunction<double, WorldVector<double> > *n)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
118
    NeumannBC *neumann = new NeumannBC(type, n, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
119
    if (rhs) 
120
      rhs->getBoundaryManager()->addBoundaryCondition(neumann);
121
122
123
124
125
126
  }

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

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

138
139
140
141
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(periodic);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(periodic);
142
143
  }

144
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
145
  void ProblemScal::useApplicationOrdering(AO *ao) 
146
147
148
149
150
151
152
153
  {
    applicationOrdering = ao;
    systemMatrix->useApplicationOrdering(ao);
    rhs->useApplicationOrdering(ao);
  }
#endif


154
155
156
157
158
  void ProblemScal::createMesh()
  {
    TEST_EXIT(Parameters::initialized())("parameters not initialized\n");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
161
    GET_PARAMETER(0, name + "->info", "%d", &info);
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    GET_PARAMETER(0, name + "->mesh", &meshName);
163

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

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

    // create the mesh
172
    mesh = new Mesh(meshName, dim);
173

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

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

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
209
210
    return 0;
  }

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

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

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

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

    clock_t first = clock();
237
238
239
240

    SolverMatrix<DOFMatrix> solverMatrix;
    solverMatrix.setMatrix(*systemMatrix);
    int iter = solver->solveSystem(solverMatrix, *solution, *rhs);
241
242

#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
243
    INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
244
245
246
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
247
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
248
249
250
251
		   TIME_USED(first, clock()));
#endif

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

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

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

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

    // === create fespace ===
Thomas Witkowski's avatar
Thomas Witkowski committed
288
    if (feSpace) {
289
290
      WARNING("feSpace already created\n");
    } else {
291
292
      if (initFlag.isSet(INIT_FE_SPACE) || 
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
293
294
	createFESpace();
      } 
295

296
297
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
	TEST_EXIT(!feSpace)("feSpace already created");
	feSpace = dynamic_cast<ProblemScal*>(adoptProblem)->getFESpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
300
      }
301
302
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
303
    if (!feSpace) 
304
      WARNING("no feSpace created\n");
305
306
307
308
309
310

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM)) {
      createMatricesAndVectors();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
311
312
313
314
315
316
      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();
317
318
319
    }

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

332
    if (!solver) 
333
334
335
      WARNING("no solver created\n");

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

348
    if (!estimator) 
349
      WARNING("no estimator created\n");
350
351

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

364
    if (!marker) 
365
      WARNING("no marker created\n");
366
367

    // === create file writer ===
368
    if(initFlag.isSet(INIT_FILEWRITER)) {
369
370
371
372
373
374
375
376
377
378
379
380
381
382
      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
383
    std::string serializationFilename = "";
384
385
386
387
388
    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
389
      GET_PARAMETER(0, name + "->input->read serialization", "%d", 
390
391
		    &readSerialization);
      if (readSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
392
	GET_PARAMETER(0, name + "->input->serialization filename", 
393
394
395
396
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

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

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

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

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

451
452
    solution->setCoarsenOperation(COARSE_INTERPOL);
    solution->set(0.0);      /*  initialize u_h  !                      */
453
454
455
456
457
  }

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

  void ProblemScal::createEstimator()
  {
    // create and set leaf data prototype
471
    mesh->setElementDataPrototype(new LeafDataEstimatable(new LeafDataCoarsenable));
472
473

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

      // init estimator
487
      estimator->addSystem(systemMatrix, solution, rhs);
488
489
490
491
492
    }
  }

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

  void ProblemScal::createFileWriter()
  {
498
    fileWriters.push_back(new FileWriter(name + "->output", mesh, solution));
499
    int writeSerialization = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
500
    GET_PARAMETER(0, name + "->output->write serialization", "%d", 
501
		  &writeSerialization);
502
503
    if (writeSerialization)
      fileWriters.push_back(new Serializer<ProblemScal>(this));
504
505
506
507
508
509
  }

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

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

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

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

    clock_t first = clock();

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

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

539
    rhs->set(0.0);
540

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

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

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

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

574
    while (elInfo) {
575
      if (useGetBound) {
Thomas Witkowski's avatar
Thomas Witkowski committed
576
577
	bound = GET_MEMORY(BoundaryType, feSpace->getBasisFcts()->getNumber());
	feSpace->getBasisFcts()->getBound(elInfo, bound);
578
579
580
581
582
583
584
      } else {
	bound = NULL;
      }
      
      systemMatrix->assemble(1.0, elInfo, bound);
      rhs->assemble(1.0, elInfo, bound);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
585
      if (useGetBound)
Thomas Witkowski's avatar
Thomas Witkowski committed
586
	FREE_MEMORY(bound, BoundaryType, feSpace->getBasisFcts()->getNumber());	    
587
      
588
589
590
591
592
593
      if (systemMatrix->getBoundaryManager()) 
	systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
      if (rhs->getBoundaryManager())
	rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
      if (solution->getBoundaryManager())
	solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
594
595
596
      elInfo = stack.traverseNext(elInfo);
    }

597
    systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
Thomas Witkowski's avatar
Thomas Witkowski committed
598

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

607
    systemMatrix->finishInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
608

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

611
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
612
613
614
615
616
617
//     PetscErrorCode ierr;
//     Mat A;

//     ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr);
#endif
    
618
    createPrecon();
619
620
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
621
  void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
622
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
623
    FUNCNAME("ProblemScal::writeResidualMesh()");
624

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

639
640
  void ProblemScal::createPrecon()
  {
641
    std::string preconType("no");
642
643
644
645
646
647
648
649
650
651
652
653
654
655
    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
656
  void ProblemScal::serialize(std::ostream &out) 
657
658
659
  {
    FUNCNAME("ProblemScal::serialize()");

660
661
    mesh->serialize(out);
    solution->serialize(out);
662
663
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
664
  void ProblemScal::deserialize(std::istream &in) 
665
666
667
  {
    FUNCNAME("ProblemScal::deserialize()");

668
669
    mesh->deserialize(in);    
    solution->deserialize(in);
670
671
672
  }

}