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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
  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 (mesh)
      DELETE mesh;

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

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

62
63
  void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
64
    solution->interpol(fct);
65
66
67
68
69
70
  }

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

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

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

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

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

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

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

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

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


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

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

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

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

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

    // create the mesh
173
    mesh = NEW Mesh(meshName, dim);
174

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

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

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

207
208
209
210
211
    return 0;
  }

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

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

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

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

    clock_t first = clock();
238
239
240
241

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  void ProblemScal::createMatricesAndVectors()
  {
    // === create vectors and system matrix ===
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
450
    systemMatrix = NEW DOFMatrix(feSpace, feSpace, "system matrix");
    rhs = NEW DOFVector<double>(feSpace, "rhs");
    solution = NEW DOFVector<double>(feSpace, "solution");
451

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

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

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

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

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

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

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

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

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

518
      adaptInfo->setEstSum(estimator->getErrorSum(), 0);
519
520

      adaptInfo->
521
	setTimeEstSum(estimator->getTimeEst(), 0);
522
523

      adaptInfo->
524
	setEstMax(estimator->getErrorMax(), 0);
525
526

      adaptInfo->
527
	setTimeEstMax(estimator->getTimeEstMax(), 0);
528
529
530
531
532
533

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

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

    clock_t first = clock();

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

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

548
    rhs->set(0.0);
549

550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
    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) 
	nnz_per_row= 5;

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

565
    // fill boundary conditions
566
567
568
569
570
571
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->initMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->initVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->initVector(solution);
572

573
    BoundaryType *bound;
574
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
575
576
577
578
579
580
    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);
581

582
    while (elInfo) {
583
      if (useGetBound) {
Thomas Witkowski's avatar
Thomas Witkowski committed
584
585
	bound = GET_MEMORY(BoundaryType, feSpace->getBasisFcts()->getNumber());
	feSpace->getBasisFcts()->getBound(elInfo, bound);
586
587
588
589
590
591
592
593
      } else {
	bound = NULL;
      }
      
      systemMatrix->assemble(1.0, elInfo, bound);
      rhs->assemble(1.0, elInfo, bound);
      
      if (useGetBound) {
Thomas Witkowski's avatar
Thomas Witkowski committed
594
	FREE_MEMORY(bound, BoundaryType, feSpace->getBasisFcts()->getNumber());	    
595
596
      }
      
597
598
599
600
601
602
      if (systemMatrix->getBoundaryManager()) 
	systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
      if (rhs->getBoundaryManager())
	rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
      if (solution->getBoundaryManager())
	solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
603
604
605
      elInfo = stack.traverseNext(elInfo);
    }

606
    systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
Thomas Witkowski's avatar
Thomas Witkowski committed
607

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

616
    systemMatrix->finishInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
617

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

620
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
623
624
625
626
//     PetscErrorCode ierr;
//     Mat A;

//     ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr);
#endif
    
627
    createPrecon();
628
629
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
630
  void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
631
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
632
    FUNCNAME("ProblemScal::writeResidualMesh()");
633

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

648
649
  void ProblemScal::createPrecon()
  {
650
    std::string preconType("no");
651
652
653
654
655
656
657
658
659
660
661
662
663
664
    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
665
  void ProblemScal::serialize(std::ostream &out) 
666
667
668
  {
    FUNCNAME("ProblemScal::serialize()");

669
670
    mesh->serialize(out);
    solution->serialize(out);
671
672
  }

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

677
678
    mesh->deserialize(in);    
    solution->deserialize(in);
679
680
681
  }

}