ProblemVec.cc 37.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include "ProblemVec.h"
#include "RecoveryEstimator.h"
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "FileWriter.h"
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
15
#include "DualTraverse.h"
16
17
18
19
20
21
22
23
#include "Mesh.h"
#include "OEMSolver.h"
#include "Preconditioner.h"
#include "MatVecMultiplier.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
#include "Flag.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "TraverseParallel.h"
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

namespace AMDiS {

  ProblemVec *ProblemVec::traversePtr_ = NULL;

  void ProblemVec::initialize(Flag initFlag,
			      ProblemVec *adoptProblem,
			      Flag adoptFlag)
  {
    FUNCNAME("ProblemVec::initialize()");
    
    // === create meshes ===
    if (meshes_.size() != 0) { 
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
	  ((!adoptFlag.isSet(INIT_MESH))&&
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
	createMesh();
      } 
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_MESH) || 
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
	meshes_ = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
51
	componentMeshes = adoptProblem->componentMeshes;
52
53
	refinementManager_ = adoptProblem->refinementManager_;
	coarseningManager_ = adoptProblem->coarseningManager_;
54
55
56
57
58
59
60
61
62
63
64
65

	// If the adopt problem has fewer components than this problem, but only one
	// mesh for all component, than scal up the componentMeshes array.
	if (adoptProblem->getNumComponents() < numComponents_) {
	  TEST_EXIT(meshes_.size() == 1)("Daran muss ich noch arbeiten!\n");
	  
	  componentMeshes_.resize(numComponents_);
	  for (int i = adoptProblem->getNumComponents(); i < numComponents_; i++) {
	    componentMeshes_[i] = componentMeshes_[0];
	  }
	}

66
67
68
69
70
71
72
      }
    }

    if (meshes_.size() == 0) 
      WARNING("no mesh created\n");

    // === create fespace ===
73
    if (feSpaces.size() != 0) {
74
75
76
77
78
79
80
81
      WARNING("feSpaces already created\n");
    } else {
      if (initFlag.isSet(INIT_FE_SPACE) || 
	  (initFlag.isSet(INIT_SYSTEM)&&!adoptFlag.isSet(INIT_FE_SPACE))) {
	createFESpace();
      } 
      if (adoptProblem &&
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
82
83
	feSpaces = adoptProblem->getFESpaces();
	componentSpaces = adoptProblem->componentSpaces;
84
85
86
87
88
89
90
91
92
93
94
95

	// If the adopt problem has fewer components than this problem, but only one
	// fe space for all component, than scal up the componentSpaces array.
	if (adoptProblem->getNumComponents() < numComponents_) {
	  TEST_EXIT(feSpaces_.size() == 1)("Daran muss ich noch arbeiten!\n");
	  
	  componentSpaces_.resize(numComponents_);
	  for (int i = adoptProblem->getNumComponents(); i < numComponents_; i++) {
	    componentSpaces_[i] = componentSpaces_[0];
	  }
	}

96
97
98
      }
    }

99
    if (feSpaces.size() == 0) 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
      WARNING("no feSpace created\n");

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM)) {
      createMatricesAndVectors();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
      solution_ = adoptProblem->getSolution();
      rhs_ = adoptProblem->getRHS();
      systemMatrix_ = adoptProblem->getSystemMatrix();
    }

    // === create solver ===
    if (solver_) {
      WARNING("solver already created\n");
    } else {
      if (initFlag.isSet(INIT_SOLVER)) {
	createSolver();
      } 
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
	TEST_EXIT(!solver_)("solver already created\n");
	solver_ = adoptProblem->getSolver();
      }
    }

    if (!solver_) 
      WARNING("no solver created\n");

    // === create estimator ===
    if (initFlag.isSet(INIT_ESTIMATOR)) {
      createEstimator();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) {
      estimator_ = adoptProblem->getEstimator();
    } 

    // === create marker ===
    if (initFlag.isSet(INIT_MARKER)) {
      createMarker();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) {
141
      marker = adoptProblem->getMarker();
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
    } 


    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER)) {
      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;
160
    std::string serializationFilename = "";
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
    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("")) {
      int readSerializationWithAdaptInfo = 0;

      GET_PARAMETER(0, name_ + "->input->read serialization", "%d", 
		    &readSerialization);
      GET_PARAMETER(0, name_ + "->input->serialization with adaptinfo", "%d",
		    &readSerializationWithAdaptInfo);

      // The serialization file is only read, if the adaptInfo part should not be used.
      // If the adaptInfo part should be also read, the serialization file will be read
      // in the constructor of the AdaptInstationary problem, because we do not have here
      // the adaptInfo object.
      if (readSerialization && !readSerializationWithAdaptInfo) {
	GET_PARAMETER(0, name_ + "->input->serialization filename", 
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
183
	std::ifstream in(serializationFilename.c_str());
184
185
186
	deserialize(in);
	in.close();
      } else {
187
188
189
190
	int globalRefinements = 0;
	GET_PARAMETER(0, meshes_[0]->getName() + "->global refinements", "%d", 
		      &globalRefinements);

191
192
193
194
195
	// Initialize the meshes if there is no serialization file.
	for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
	  if (initFlag.isSet(INIT_MESH) && 
	      meshes_[i] && 
	      !(meshes_[i]->isInitialized())) {
196
197
	    meshes_[i]->initialize();	    
	    refinementManager_->globalRefine(meshes_[i], globalRefinements);
198
199
200
201
202
203
204
205
206
207
208
209
	  }
	}	
      }
    }

    doOtherStuff();
  }

  void ProblemVec::createMesh() 
  {
    FUNCNAME("ProblemVec::createMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
210
    componentMeshes.resize(nComponents);
211
    std::map<int, Mesh*> meshForRefinementSet;
212
213
    char number[3];

214
    std::string meshName("");
215
    GET_PARAMETER(0, name_ + "->mesh", &meshName);
216
    TEST_EXIT(meshName != "")("no mesh name specified\n");
217
218
    int dim = 0;
    GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
219
    TEST_EXIT(dim)("no problem dimension specified!\n");
220

221
    for (int i = 0; i < nComponents; i++) {
222
      sprintf(number, "%d", i);
223
224
225
      int refSet = -1;
      GET_PARAMETER(0, name_ + "->refinement set[" + number + "]", "%d", &refSet);
      if (refSet < 0) {
226
227
	refSet = 0;
      }
228
      if (meshForRefinementSet[refSet] == NULL) {
229
230
231
232
	Mesh *newMesh = NEW Mesh(meshName, dim);
	meshForRefinementSet[refSet] = newMesh;
	meshes_.push_back(newMesh);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
233
      componentMeshes[i] = meshForRefinementSet[refSet];
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
    }
    switch(dim) {
    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");
    }
  }

  void ProblemVec::createFESpace()
  {
    FUNCNAME("ProblemVec::createFESpace()");

    int degree = 1;
    char number[3];

260
    std::map< std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
261
    int dim = -1;
262
    GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
263
    TEST_EXIT(dim != -1)("no problem dimension specified!\n");
264

265
    componentSpaces.resize(nComponents, NULL);
266

267
    for (int i = 0; i < nComponents; i++) {
268
269
270
      sprintf(number, "%d", i);
      GET_PARAMETER(0, name_ + "->polynomial degree[" + number + "]","%d", &degree);

271
      TEST_EXIT(componentSpaces[i] == NULL)("feSpace already created\n");
272

Thomas Witkowski's avatar
Thomas Witkowski committed
273
      if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
274
275
276
	FiniteElemSpace *newFESpace = 
	  FiniteElemSpace::provideFESpace(NULL,
					  Lagrange::getLagrange(dim, degree),
Thomas Witkowski's avatar
Thomas Witkowski committed
277
					  componentMeshes[i],
278
					  name_ + "->feSpace");
Thomas Witkowski's avatar
Thomas Witkowski committed
279
	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFESpace;
280
	feSpaces.push_back(newFESpace);
281
      }
282
      componentSpaces[i] = 
Thomas Witkowski's avatar
Thomas Witkowski committed
283
	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
284
285
286
    }

    // create dof admin for vertex dofs if neccessary
287
    for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
288
289
290
291
292
293
294
295
296
297
298
299
300
301
      if (meshes_[i]->getNumberOfDOFs(VERTEX) == 0) {
	DimVec<int> ln_dof(meshes_[i]->getDim(), DEFAULT_VALUE, 0);
	ln_dof[VERTEX]= 1;
	meshes_[i]->createDOFAdmin("vertex dofs", ln_dof);      
      }
    }
  }

  void ProblemVec::createMatricesAndVectors()
  {
    FUNCNAME("ProblemVec::createMatricesAndVectors()");

    // === create vectors and system matrix ===

302
    systemMatrix_ = NEW Matrix<DOFMatrix*>(nComponents, nComponents);
303
    systemMatrix_->set(NULL);
304
305
    rhs_ = NEW SystemVector("rhs", componentSpaces, nComponents);
    solution_ = NEW SystemVector("solution", componentSpaces, nComponents);
306
307

    char number[10];
308
    std::string numberedName;
Thomas Witkowski's avatar
Thomas Witkowski committed
309
    for (int i = 0; i < nComponents; i++) {
310
311
      (*systemMatrix_)[i][i] = NEW DOFMatrix(componentSpaces[i], 
					     componentSpaces[i], "A_ii");
312
313
      (*systemMatrix_)[i][i]->setCoupleMatrix(false);
      sprintf(number, "[%d]", i);
314
      numberedName = "rhs" + std::string(number);
315
      rhs_->setDOFVector(i, NEW DOFVector<double>(componentSpaces[i], numberedName));
316
      numberedName = name_ + std::string(number);
317
      solution_->setDOFVector(i, NEW DOFVector<double>(componentSpaces[i], 
318
319
320
321
322
323
324
325
326
327
328
329
330
331
						       numberedName));
      solution_->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution_->getDOFVector(i)->set(0.0);
    }

    // === create matVec ===
    matVec_ = NEW StandardMatVec<Matrix<DOFMatrix*>, SystemVector>(systemMatrix_);
  }

  void ProblemVec::createSolver()
  {
    FUNCNAME("ProblemVec::createSolver()");

    // === create solver ===
332
    std::string solverType("no");
333
334
335
336
337
338
339
340
341
342
343
344
    GET_PARAMETER(0, name_ + "->solver", &solverType);
    OEMSolverCreator<SystemVector> *solverCreator = 
      dynamic_cast<OEMSolverCreator<SystemVector>*>(
						    CreatorMap<OEMSolver<SystemVector> >
						    ::getCreator(solverType)
						    );
    TEST_EXIT(solverCreator)("no solver type\n");
    solverCreator->setName(name_ + "->solver");
    solver_ = solverCreator->create();
    solver_->initParameters();

    // === create preconditioners ===
345
    std::string preconType("no");
346
347

    PreconditionerScal *scalPrecon;
348
    PreconditionerVec *vecPrecon = NEW PreconditionerVec(nComponents);
349
350
351
352
353
354
355
356
357
358
359

    GET_PARAMETER(0, name_ + "->solver->left precon", &preconType);
    CreatorInterface<PreconditionerScal> *preconCreator =
      CreatorMap<PreconditionerScal>::getCreator(preconType);

    int i, j;

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

360
      for(i = 0; i < nComponents; i++) {
361
	dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
362
	  setSizeAndRow(nComponents, i);
363
364
    
	scalPrecon = preconCreator->create();
365
	for(j = 0; j < nComponents; j++) {
366
367
368
369
370
371
372
373
	  scalPrecon->setMatrix(&(*systemMatrix_)[i][j], j);
	}
	vecPrecon->setScalarPrecon(i, scalPrecon);
      }
      leftPrecon_ = vecPrecon;
    }


374
    vecPrecon = NEW PreconditionerVec(nComponents);
375
376
377
378
379
380
381
382
383
384

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

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


385
      for(i = 0; i < nComponents; i++) {
386
	dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
387
	  setSizeAndRow(nComponents, i);
388
389
    
	scalPrecon = preconCreator->create();
390
	for(j = 0; j < nComponents; j++) {
391
392
393
394
395
396
397
398
399
400
	  scalPrecon->setMatrix(&(*systemMatrix_)[i][j], j);
	}
	vecPrecon->setScalarPrecon(i, scalPrecon);
      }
      rightPrecon_ = vecPrecon;
    }


    // === create vector creator ===
    solver_->setVectorCreator(NEW SystemVector::Creator("temp",
401
							componentSpaces, 
402
							nComponents));
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
  }

  void ProblemVec::createEstimator()
  {
    FUNCNAME("ProblemVec::createEstimator()");

    int i, j;

    // create and set leaf data prototype
    for(i = 0; i < static_cast<int>(meshes_.size()); i++) {
      meshes_[i]->setElementDataPrototype
	(NEW LeafDataEstimatableVec(NEW LeafDataCoarsenableVec));
    }  

    char number[3];
418
    std::string estName;
419

420
    for(i = 0; i < nComponents; i++) {
421
422
      TEST_EXIT(estimator_[i] == NULL)("estimator already created\n");
      sprintf(number, "%d", i);
423
      estName = name_ + "->estimator[" + std::string(number) + "]";
424
425

      // === create estimator ===
426
      std::string estimatorType("no");
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
      GET_PARAMETER(0, estName, &estimatorType);
      EstimatorCreator *estimatorCreator = 
	dynamic_cast<EstimatorCreator*>(
					CreatorMap<Estimator>::getCreator(estimatorType));
      if(estimatorCreator) {
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
	if(estimatorType == "recovery") {
	  dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->
	    setSolution(solution_->getDOFVector(i));
	}
	estimator_[i] = estimatorCreator->create();
      }


      if(estimator_[i]) {
443
	for(j=0; j < nComponents; j++) {
444
445
446
447
448
449
450
451
452
453
454
455
	  estimator_[i]->addSystem((*systemMatrix_)[i][j], 
				   solution_->getDOFVector(j), 
				   rhs_->getDOFVector(j));
	}
      }
    }
  }

  void ProblemVec::createMarker()
  {
    FUNCNAME("ProblemVec::createMarker()");

456
    std::string numberedName;
457
458
    char number[10];
    int numMarkersCreated = 0;
459

460
    for (int i = 0; i < nComponents; i++) {
461
      sprintf(number, "[%d]", i);
462
      numberedName = name_ + "->marker" + std::string(number);
463
464
      marker[i] = Marker::createMarker(numberedName, i);
      if (marker[i]) {
465
466
	numMarkersCreated++;
	if (numMarkersCreated > 1)
467
	  marker[i]->setMaximumMarking(true);
468
469
470
471
472
473
474
475
476
477
      }
    }
  }

  void ProblemVec::createFileWriter()
  {
    FUNCNAME("ProblemVec::createFileWriter()");
  

    // Create one filewriter for all components of the problem
478
479
    std::string numberedName  = name_ + "->output";
    std::string filename = "";
480
481
482
    GET_PARAMETER(0, numberedName + "->filename", &filename);

    if (filename != "") {
483
      std::vector< DOFVector<double>* > solutionList(nComponents);
484

485
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
486
	TEST_EXIT(componentMeshes[0] == componentMeshes[i])
487
488
489
490
491
492
	  ("All Meshes have to be equal to write a vector file.\n");

	solutionList[i] = solution_->getDOFVector(i);
      }

      fileWriters_.push_back(NEW FileWriter(numberedName,
Thomas Witkowski's avatar
Thomas Witkowski committed
493
					    componentMeshes[0],
494
495
496
497
498
499
					    solutionList));
    }


    // Create own filewriters for each components of the problem
    char number[10];
500
    for (int i = 0; i < nComponents; i++) {
501
      sprintf(number, "[%d]", i);
502
      numberedName  = name_ + "->output" + std::string(number);
503
504
505
506
507
      filename = "";
      GET_PARAMETER(0, numberedName + "->filename", &filename);

      if (filename != "") {
	fileWriters_.push_back(NEW FileWriter(numberedName, 
Thomas Witkowski's avatar
Thomas Witkowski committed
508
					      componentMeshes[i], 
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
					      solution_->getDOFVector(i)));
      }
    }


    // Check for serializer
    int writeSerialization = 0;
    GET_PARAMETER(0, name_ + "->write serialization", "%d", &writeSerialization);
    if (writeSerialization) {
      MSG("Use are using the obsolete parameter: %s->write serialization\n", name_.c_str());
      MSG("Please use instead the following parameter: %s->output->write serialization\n", name_.c_str());
      ERROR_EXIT("Usage of an obsolete parameter (see message above)!\n");
    }

    GET_PARAMETER(0, name_ + "->output->write serialization", "%d", &writeSerialization);
    if (writeSerialization) {
      fileWriters_.push_back(NEW Serializer<ProblemVec>(this));
    }
  }

  void ProblemVec::doOtherStuff()
  {
  }

533
  void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
534
535
536
537
538
539
540
541
542
543
544
545
546
  {
    FUNCNAME("Problem::solve()");

    if (!solver_) {
      WARNING("no solver\n");
      return;
    }

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

    clock_t first = clock();
547
548
    int iter = solver_->solve(matVec_, solution_, rhs_, 
			      leftPrecon_, rightPrecon_, fixedMatrix);
549
    
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
#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);
    adaptInfo->setMaxSolverIterations(solver_->getMaxIterations());
    adaptInfo->setSolverTolerance(solver_->getTolerance());
    adaptInfo->setSolverResidual(solver_->getResidual());
  }

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

    clock_t first = clock();

572
573
574
575
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

576
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
577
      computeError(adaptInfo);
578
579
580
581
582
583
584
585
586
587
588
589
590
    } else {
      for (int i = 0; i < nComponents; i++) {
	Estimator *scalEstimator = estimator_[i];
	
	if (scalEstimator) {
	  scalEstimator->estimate(adaptInfo->getTimestep());
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
	  adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	  adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	} else {
	  WARNING("no estimator for component %d\n" , i);
	}
591
592
593
      }
    }

594
595
596
597
598
599
600
601
602
603
#ifdef _OPENMP
    INFO(info_, 8)("estimation of the error needed %.5f seconds system time / %.5f seconds wallclock time\n",
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
    INFO(info_, 8)("estimation of the error needed %.5f seconds\n",
		   TIME_USED(first, clock()));

#endif

604
605
606
607
608
609
610
611
612
613
614
  }

  Flag ProblemVec::markElements(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::markElements()");

    // to enforce albert-like behavior: refinement even if space tolerance
    // here is reached already because of time adaption
    allowFirstRefinement();

    Flag markFlag = 0;
615
    for (int i = 0; i < nComponents; i++) {
616
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
617
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
618
619
620
621
      } else {
	WARNING("no marker for component %d\n", i);
      }
    }
622
    
623
624
625
626
627
628
629
    return markFlag;
  }

  Flag ProblemVec::refineMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::refineMesh()");

630
    int nMeshes = static_cast<int>(meshes_.size());
631
    Flag refineFlag = 0;
632
    for (int i = 0; i < nMeshes; i++) {
633
634
635
636
637
638
639
640
641
      refineFlag |= refinementManager_->refineMesh(meshes_[i]);
    }
    return refineFlag;
  }

  Flag ProblemVec::coarsenMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::coarsenMesh()");

642
    int nMeshes = static_cast<int>(meshes_.size());
643
    Flag coarsenFlag = 0;
644
645
    for (int i = 0; i < nMeshes; i++) {
      if (adaptInfo->isCoarseningAllowed(i)) {
646
647
648
649
650
651
652
653
654
655
656
657
658
	coarsenFlag |= coarseningManager_->coarsenMesh(meshes_[i]);

	WARNING("coarsening for component %d no allowed\n", i);
      }
    }
    return coarsenFlag;
  }

  Flag ProblemVec::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ProblemVec::oneIteration()");

    if (allowFirstRef_) {
659
      for (int i = 0; i < nComponents; i++) {
660
661
662
663
	adaptInfo->allowRefinement(true, i);
      }
      allowFirstRef_ = false;
    } else {
664
      for (int i = 0; i < nComponents; i++) {
665
666
667
668
669
670
671
672
673
674
675
	if (adaptInfo->spaceToleranceReached(i)) {
	  adaptInfo->allowRefinement(false, i);
	} else {
	  adaptInfo->allowRefinement(true, i);	
	}
      }
    }

    return StandardProblemIteration::oneIteration(adaptInfo, toDo);
  }

676
677
  void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				     bool asmMatrix, bool asmVector)
678
679
680
681
682
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

    clock_t first = clock();

683
684
685
686
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
687
    for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
      meshes_[i]->dofCompress();
    }

    Flag assembleFlag = 
      flag | 
      (*systemMatrix_)[0][0]->getAssembleFlag() | 
      rhs_->getDOFVector(0)->getAssembleFlag()   |
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
705

706
    for (int i = 0; i < nComponents; i++) {
707
      MSG("%d DOFs for %s\n", 
708
709
	  componentSpaces[i]->getAdmin()->getUsedSize(), 
	  componentSpaces[i]->getName().c_str());
710
711

      rhs_->getDOFVector(i)->set(0.0);
712
      for (int j = 0; j < nComponents; j++) {
713
714
715
	if ((*systemMatrix_)[i][j]) {
	  // The matrix should not be deleted, if it was assembled before
	  // and it is marked to be assembled only once.
716
	  if (!(assembleMatrixOnlyOnce_[i][j] && assembledMatrix_[i][j]) && asmMatrix) {
717
718
719
	    (*systemMatrix_)[i][j]->clear();
	  }
	}
720
721
722
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
723
    for (int i = 0; i < nComponents; i++) {
724
      for (int j = 0; j < nComponents; j++) {
725
726
727
728
	// Only if this variable is true, the current matrix will be assembled.	
	bool assembleMatrix = true;
	// The DOFMatrix which should be assembled (or not, if assembleMatrix
	// will be set to false).
729
730
	DOFMatrix *matrix = (*systemMatrix_)[i][j];

731
732
733
734
735
736
737
738
739
740
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
	if (assembleMatrixOnlyOnce_[i][j] && assembledMatrix_[i][j]) {
	  assembleMatrix = false;
	}
	// If there is no DOFMatrix (e.g. if it is completly 0), do not assemble.
	if (!matrix) {
	  assembleMatrix = false;
	}

741
742
743
744
	if (!asmMatrix) {
	  assembleMatrix = false;
	}

745
746
747
748
749
750
751
	// If the matrix should not be assembled, the rhs vector has to be considered.
	// This will be only done, if i == j. So, if both is not true, we can jump
	// to the next matrix.
	if (!assembleMatrix && i != j) {
	  continue;
	}

752
753
754
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

Thomas Witkowski's avatar
Thomas Witkowski committed
755
756
757
758
	if (componentSpaces[i] == componentSpaces[j]) {
	  assembleOnOneMesh(componentSpaces[i],
			    assembleFlag,
			    assembleMatrix ? matrix : NULL,
759
			    ((i == j) && asmVector) ? rhs_->getDOFVector(i) : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
760
761
762
763
	} else {
	  assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
			      assembleFlag,
			      assembleMatrix ? matrix : NULL,
764
			      ((i == j) && asmVector) ? rhs_->getDOFVector(i) : NULL);	  
765
766
767
768
	  TEST_EXIT(matrix->getUsedSize() == componentSpaces[i]->getAdmin()->getUsedSize())
	    ("Assembled matrix has wrong dimension!\n");
	  TEST_EXIT(matrix->getNumCols() == componentSpaces[j]->getAdmin()->getUsedSize())
	    ("Assembled matrix has wrong dimension!\n");
769
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
770

771
772
773
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->exitMatrix(matrix);	  
	
774
	assembledMatrix_[i][j] = true;
775
776
      }

777
778
779
780
781
      // And now assemble boundary conditions on the vectors
      assembleBoundaryConditions(rhs_->getDOFVector(i),
				 solution_->getDOFVector(i),
				 componentMeshes[i],
				 assembleFlag);
782
    }
783
784
785
786
787
788

#ifdef _OPENMP
    INFO(info_, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
789
    INFO(info_, 8)("buildAfterCoarsen needed %.5f seconds\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
790
		   TIME_USED(first, clock()));
791
#endif
792
793
794
795
796
797
  }

  void ProblemVec::writeFiles(AdaptInfo *adaptInfo, bool force) 
  {
    FUNCNAME("ProblemVec::writeFiles()");

798
799
800
801
802
803
804
805
806
807
808
809
810
    clock_t first = clock();

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

    int i;
    int size = static_cast<int>(fileWriters_.size());
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 1)
#endif
    for (i = 0; i < size; i++) {
      fileWriters_[i]->writeFiles(adaptInfo, force);
811
    }
812
813
814
815
816
817
818
819
820
    
#ifdef _OPENMP
    INFO(info_, 8)("writeFiles needed %.5f seconds system time / %.5f seconds wallclock time\n",
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
    INFO(info_, 8)("writeFiles needed %.5f seconds\n",
		   TIME_USED(first, clock()));
#endif
821
822
  }

823
  void ProblemVec::interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct) 
824
825
826
827
828
829
830
831
832
833
834
835
836
  {
    FUNCNAME("ProblemVec::interpolInitialSolution()");

    solution_->interpol(fct);
  }

  void ProblemVec::addMatrixOperator(Operator *op, 
				     int i, int j,
				     double *factor,
				     double *estFactor)
  {
    FUNCNAME("ProblemVec::addMatrixOperator()");

837
    if (!(*systemMatrix_)[i][j]) {
838
      TEST_EXIT(i != j)("should have been created already\n");
839
840
      (*systemMatrix_)[i][j] = NEW DOFMatrix(componentSpaces[i],
					     componentSpaces[j],
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
					     "");
      (*systemMatrix_)[i][j]->setCoupleMatrix(true);

      (*systemMatrix_)[i][j]->getBoundaryManager()->
	setBoundaryConditionMap((*systemMatrix_)[i][i]->getBoundaryManager()->
				getBoundaryConditionMap());
    }    
    (*systemMatrix_)[i][j]->addOperator(op, factor, estFactor);
  }

  void ProblemVec::addVectorOperator(Operator *op, int i,
				     double *factor,
				     double *estFactor)
  {
    FUNCNAME("ProblemVec::addVectorOperator()");

    rhs_->getDOFVector(i)->addOperator(op, factor, estFactor);
  }

  void ProblemVec::addDirichletBC(BoundaryType type, int system,
				  AbstractFunction<double, WorldVector<double> >* b)
  {
    FUNCNAME("ProblemVec::addDirichletBC()");

    DirichletBC *dirichlet = new DirichletBC(type, 
					     b, 
867
					     componentSpaces[system]);
868
    for (int i = 0; i < nComponents; i++) {
869
870
871
872
      if (systemMatrix_ && (*systemMatrix_)[system][i]) {
	(*systemMatrix_)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet);
      }
    }
873

874
875
    if (rhs_)
      rhs_->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
876

877
878
879
880
881
882
883
884
885
886
887
    if (solution_)
      solution_->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
  }

  void ProblemVec::addNeumannBC(BoundaryType type, int row, int col, 
				AbstractFunction<double, WorldVector<double> > *n)
  {
    FUNCNAME("ProblemVec::addNeumannBC()");

    NeumannBC *neumann = 
      new NeumannBC(type, n, 
888
889
		    componentSpaces[row], 
		    componentSpaces[col]);
890
    if (rhs_)
891
892
893
894
895
896
897
898
899
900
901
      rhs_->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
  }

  void ProblemVec::addRobinBC(BoundaryType type, int row, int col, 
			      AbstractFunction<double, WorldVector<double> > *n,
			      AbstractFunction<double, WorldVector<double> > *r)
  {
    FUNCNAME("ProblemVec::addRobinBC()");

    RobinBC *robin = 
      new RobinBC(type, n, r, 
902
903
		  componentSpaces[row], 
		  componentSpaces[col]);
904
    if (rhs_)
905
      rhs_->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
906
907

    if (systemMatrix_ && (*systemMatrix_)[row][col]) {
908
909
910
911
912
913
914
915
      (*systemMatrix_)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
    }
  }

  void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col) 
  {
    FUNCNAME("ProblemVec::addPeriodicBC()");

916
    FiniteElemSpace *feSpace = componentSpaces[row];
917
918
919

    PeriodicBC *periodic = new PeriodicBC(type, feSpace);

920
    if (systemMatrix_ && (*systemMatrix_)[row][col]) 
921
      (*systemMatrix_)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);
922
923

    if (rhs_) 
924
925
926
927
      rhs_->getDOFVector(row)->getBoundaryManager()->
	addBoundaryCondition(periodic);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
928
929
930
931
932
  void ProblemVec::assembleOnOneMesh(FiniteElemSpace *feSpace, Flag assembleFlag,
				     DOFMatrix *matrix, DOFVector<double> *vector)
  {
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
933
934

#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
935
    TraverseParallelStack stack;
936
937
938
#else
    TraverseStack stack;
#endif   
Thomas Witkowski's avatar
Thomas Witkowski committed
939

940
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
941
#pragma omp parallel
942
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
943
944
945
946
947
948
949
    {
      BoundaryType *bound = useGetBound_ ? GET_MEMORY(BoundaryType, basisFcts->getNumber()) : NULL;

      // Create for every thread its private matrix and vector, on that
      // the thread will assemble its part of the mesh.
      DOFMatrix *tmpMatrix = NULL;
      DOFVector<double> *tmpVector = NULL; 
Thomas Witkowski's avatar
Thomas Witkowski committed
950
951

      if (matrix) {
Thomas Witkowski's avatar
Thomas Witkowski committed
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
	tmpMatrix = NEW DOFMatrix(matrix->getRowFESpace(),
				  matrix->getColFESpace(),
				  "tmp");

	// Copy the global matrix to the private matrix, because we need the
	// operators defined on the global matrix in the private one. Only the
	// values have to be set to zero.
	*tmpMatrix = *matrix;
	tmpMatrix->clear();
      }

      if (vector) {
	tmpVector = NEW DOFVector<double>(vector->getFESpace(), "tmp");

	// Copy the global vector to the private vector, because we need the
	// operatirs defined on the global vector in the private one. But set
	// the values to zero of the private vector after copying.
	*tmpVector = *vector;
	tmpVector->set(0.0);
      }

973
974
975
976
      // Because we are using the parallel traverse stack, each thread will
      // traverse only a part of the mesh.
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);

977
978
979
980
981
982
983
      // After creating privat copies of the DOFMatrix and the DOFVector, all threads
      // have to wait at this barrier. Especially for small problems this is required,
      // because otherwise one thread may be finished with assembling, before another
      // has made his private copy.
#ifdef _OPENMP
#pragma omp barrier
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
      while (elInfo) {
	if (useGetBound_) {
	  basisFcts->getBound(elInfo, bound);
	}
	
	if (matrix) {
	  tmpMatrix->assemble(1.0, elInfo, bound);
	  
	  // Take the matrix boundary manager from the public matrix,
	  // but assemble the boundary conditions on the thread private matrix.
	  if (matrix->getBoundaryManager()) {
	    matrix->getBoundaryManager()->
	      fillBoundaryConditions(elInfo, tmpMatrix);
	  }		      
	}
	
	if (vector) {
	  tmpVector->assemble(1.0, elInfo, bound);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
1003
	
Thomas Witkowski's avatar
Thomas Witkowski committed
1004
	elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
1005
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1006
1007
1008
1009
1010
1011
1012

      // After mesh traverse, all thread have to added their private matrices and
      // vectors to the global public matrix and public vector. Therefore, this is 
      // a critical section, which is allowed to be executed by on thread only at 
      // the same time.

      if (matrix) {
1013
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
1014
#pragma omp critical
1015
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
	{
	  addDOFMatrix(matrix, tmpMatrix);

	  // Remove rows corresponding to DOFs on a Dirichlet boundary.
	  matrix->removeRowsWithDBC(tmpMatrix->getApplyDBCs());
	}

	DELETE tmpMatrix;
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
1026
      if (vector) {
1027
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
1028
#pragma omp critical
1029
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
1030
1031
1032
	*vector += *tmpVector;

	DELETE tmpVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
1033
1034
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
1035
1036
1037
1038
1039
1040
      if (useGetBound_) {
	FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
      }	      

    } // pragma omp parallel

Thomas Witkowski's avatar
Thomas Witkowski committed
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
  }

  void ProblemVec::assembleOnDifMeshes(FiniteElemSpace *rowFeSpace, FiniteElemSpace *colFeSpace,
				       Flag assembleFlag,
				       DOFMatrix *matrix, DOFVector<double> *vector)
  {
    Mesh *rowMesh = rowFeSpace->getMesh();
    Mesh *colMesh = colFeSpace->getMesh();

    const BasisFunction *basisFcts = rowFeSpace->getBasisFcts();
    BoundaryType *bound = NULL;
    if (useGetBound_) {
      bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
    }

    DualTraverse dualTraverse;
    ElInfo *rowElInfo, *colElInfo;
    ElInfo *largeElInfo, *smallElInfo;

1060
    dualTraverse.setFillSubElemMat(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
1061
1062
1063
1064
1065
1066
1067
1068
    bool cont = dualTraverse.traverseFirst(rowMesh, colMesh, -1, -1,
					   assembleFlag, assembleFlag,
					   &rowElInfo, &colElInfo,
					   &smallElInfo, &largeElInfo);
    while (cont) {
      Element *rowElem = rowElInfo->getElement();
      Element *colElem = colElInfo->getElement();
      
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
      if (useGetBound_) {
	basisFcts->getBound(rowElInfo, bound);
      }
      
      if (matrix) {
	matrix->assemble(1.0, rowElInfo, colElInfo, smallElInfo, largeElInfo, bound);
	
	if (matrix->getBoundaryManager()) {
	  matrix->getBoundaryManager()->
	    fillBoundaryConditions(rowElInfo, matrix);
	}		      
Thomas Witkowski's avatar
Thomas Witkowski committed
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
      }
      
      if (vector) {
	vector->assemble(1.0, rowElInfo, bound);
      }

      cont = dualTraverse.traverseNext(&rowElInfo, &colElInfo,
				       &smallElInfo, &largeElInfo);
    }

    if (useGetBound_) {
      FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
    }
  }

1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
  void ProblemVec::assembleBoundaryConditions(DOFVector<double> *rhs,
					      DOFVector<double> *solution,
					      Mesh *mesh,
					      Flag assembleFlag)
  {
    /* ================ Initialization of vectors ==================== */

    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->initVector(rhs);      
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->initVector(solution);
    
#ifdef _OPENMP
    TraverseParallelStack stack;
#else
    TraverseStack stack;
#endif

    /* ================= Parallel Boundary Assemblage ================= */
#ifdef _OPENMP
#pragma omp parallel
#endif
    {
      // Each thread assembles on its own dof-vectors.
      DOFVector<double> *tmpRhsVec = NEW DOFVector<double>(rhs->getFESpace(), "tmpRhs");
      DOFVector<double> *tmpSolVec = NEW DOFVector<double>(solution->getFESpace(), "tmpSol");
      tmpRhsVec->set(0.0);
      tmpSolVec->set(0.0);


      // (Parallel) traverse of mesh.
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
      while (elInfo) {
	if (rhs->getBoundaryManager())
	  rhs->getBoundaryManager()-> fillBoundaryConditions(elInfo, tmpRhsVec);
	
	if (solution->getBoundaryManager())
	  solution->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpSolVec);
      
	elInfo = stack.traverseNext(elInfo);
      }


      // After (parallel) mesh traverse, the result is applied to the final
      // vectors. This section is not allowed to be executed by more than one
      // thread at the same time.
#ifdef _OPENMP
#pragma omp critical
#endif
      {
	DOFVector<double>::Iterator rhsIt(rhs, USED_DOFS);
	DOFVector<double>::Iterator solIt(solution, USED_DOFS);
	DOFVector<double>::Iterator tmpRhsIt(tmpRhsVec, USED_DOFS);
	DOFVector<double>::Iterator tmpSolIt(tmpSolVec, USED_DOFS);
	for (rhsIt.reset(), solIt.reset(), tmpRhsIt.reset(), tmpSolIt.reset();
	     !rhsIt.end();
	     ++rhsIt, ++solIt, ++tmpRhsIt, ++tmpSolIt) {	     
	  if (*tmpRhsIt != 0.0)
	    *rhsIt = *tmpRhsIt;
	  
	  if (*tmpSolIt != 0.0)
	    *solIt = *tmpSolIt;	  
	}
      } // pragma omp critical


      DELETE tmpRhsVec;
      DELETE tmpSolVec;
    } // pragma omp parallel
     

    /* ======================= Finalize vectors ================== */

    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->exitVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->exitVector(solution);
  }


1175
  void ProblemVec::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
  {
    FUNCNAME("ProblemVec::writeResidualMesh()");

    Mesh *mesh = this->getMesh(0);
    FiniteElemSpace *fe = this->getFESpace(0);
    
    std::map<int, double> vec;
    
    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);    
  }

1202
  void ProblemVec::serialize(std::ostream &out) 
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
  {
    FUNCNAME("ProblemVec::serialize()");

    SerializerUtil::serializeBool(out, &allowFirstRef_);
    
    for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
      meshes_[i]->serialize(out);
    }

    solution_->serialize(out);
  }

1215
  void ProblemVec::deserialize(std::istream &in) 
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
  {
    FUNCNAME("ProblemVec::deserialize()");

    SerializerUtil::deserializeBool(in, &allowFirstRef_);

    for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
      meshes_[i]->deserialize(in);
    }

    solution_->deserialize(in);
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282

  void ProblemVec::computeError(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::computeError()");

    for (int i = 0; i < nComponents; i++) {		
      TEST_EXIT(exactSolutionFcts[i])("No solution function given!\n");

      // Compute the difference between exact and computed solution
      DOFVector<double> *tmp = NEW DOFVector<double>(componentSpaces[i], "tmp");
      tmp->interpol(exactSolutionFcts[i]);
      double solMax = tmp->absMax();
      *tmp -= *(solution_->getDOFVector(i));
      
      MSG("L2    error = %.8e\n", tmp->L2Norm());
      MSG("L-inf error = %.8e\n", tmp->absMax() / solMax);
      
      adaptInfo->setEstSum(tmp->absMax() / solMax, i);
      adaptInfo->setEstMax(tmp->absMax() / solMax, i);
      
      // To set element estimates, compute a vector with the difference
      // between exact and computed solution for each DOF.
      DOFVector<double> *sol = NEW DOFVector<double>(componentSpaces[i], "tmp");
      sol->interpol(exactSolutionFcts[i]);
      DOFVector<double>::Iterator it1(sol, USED_DOFS);
      DOFVector<double>::Iterator it2(tmp, USED_DOFS);
      for (it1.reset(), it2.reset(); !it1.end(); ++it1, ++it2) {
	if ((abs(*it1) <= DBL_TOL) || (abs(*it2) <= DBL_TOL)) {
	  *it2 = 0.0;
	} else {
	  *it2 = abs(*it2 / *it1);
	}
      }

      // Compute estimate for every mesh element
      Vector<DegreeOfFreedom> locInd(componentSpaces[i]->getBasisFcts()->getNumber());
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	componentSpaces[i]->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
							       componentSpaces[i]->getAdmin(),
							       &locInd);
	double estimate = 0.0;
	for (int j = 0; j < componentSpaces[i]->getBasisFcts()->getNumber(); j++) {
	  estimate += (*tmp)[locInd[j]];
	}
	elInfo->getElement()->setEstimation(estimate, i);
	elInfo->getElement()->setMark(0);
								
	elInfo = stack.traverseNext(elInfo);
      }  
      
      DELETE tmp;	
      DELETE sol;
    }						           
  }
1283
1284
}