ProblemStat.cc 49.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <sstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include <boost/lexical_cast.hpp>
15
16

#include "ProblemStat.h"
17
18
19
20
21
22
23
24
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Marker.h"
#include "AdaptInfo.h"
25
#include "io/FileWriter.h"
26
27
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "DualTraverse.h"
29
30
31
32
33
34
#include "Mesh.h"
#include "OEMSolver.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
35
#include "Flag.h"
36
#include "est/Estimator.h"
37
38
#include "io/VtkWriter.h"
#include "io/ValueReader.h"
39
#include "ProblemStatDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
40
#include "Debug.h"
41
42
43

namespace AMDiS {

44
  using namespace std;
45
46
  using boost::lexical_cast;

47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
  ProblemStatSeq::ProblemStatSeq(string nameStr,
				 ProblemIterationInterface *problemIteration)
    : StandardProblemIteration(this),
      name(nameStr),
      nComponents(-1),
      nMeshes(0),
      traverseInfo(0),
      solver(NULL),
      solution(NULL),
      rhs(NULL),
      systemMatrix(NULL),
      useGetBound(true),
      refinementManager(NULL),
      coarseningManager(NULL),
      info(10),
      deserialized(false),
      computeExactError(false),
      boundaryConditionSet(false),
      writeAsmInfo(false)
  {
    Parameters::get(name + "->components", nComponents);
    TEST_EXIT(nComponents > 0)("components not set!\n");    
    estimator.resize(nComponents, NULL);
    marker.resize(nComponents, NULL);
    
    assembleMatrixOnlyOnce.resize(nComponents);
    assembledMatrix.resize(nComponents);
    for (int i = 0; i < nComponents; i++) {
      assembleMatrixOnlyOnce[i].resize(nComponents);
      assembledMatrix[i].resize(nComponents);
      for (int j = 0; j < nComponents; j++) {
	assembleMatrixOnlyOnce[i][j] = false;
	assembledMatrix[i][j] = false;
      }
    }
    
    exactSolutionFcts.resize(nComponents);


    // === Initialize name of components. ===

    componentNames.resize(nComponents, "");
    for (int i = 0; i < nComponents; i++)
      componentNames[i] = "solution[" + lexical_cast<string>(i) + "]";

    Parameters::get(name + "->name", componentNames);
93
94
    componentNames.resize(nComponents);

95
96
97
    for (int i = 0; i < nComponents; i++)
      Parameters::get(name + "->name[" + lexical_cast<string>(i) + "]",
		      componentNames[i]);
Thomas Witkowski's avatar
Thomas Witkowski committed
98
99

    Parameters::get("debug->write asm info", writeAsmInfo);
100
101
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
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

  ProblemStatSeq::~ProblemStatSeq()
  {
    delete rhs;
    rhs = NULL;
    delete solution;
    solution = NULL;
    
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*systemMatrix)[i][j]) {
	  delete (*systemMatrix)[i][j];
	  (*systemMatrix)[i][j] = NULL;
	}

    delete systemMatrix;
    systemMatrix = NULL;

    for (unsigned int i = 0; i < meshes.size(); i++)
      delete meshes[i];
    
    for (unsigned int i = 0; i < estimator.size(); i++)
      delete estimator[i];

    for (unsigned int i = 0; i < marker.size(); i++)
      delete marker[i];
  }


131
132
133
  void ProblemStatSeq::initialize(Flag initFlag,
				  ProblemStatSeq *adoptProblem,
				  Flag adoptFlag)
134
  {
135
    FUNCNAME("ProblemStat::initialize()");
136

137
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
138
    if (meshes.size() != 0) { 
139
140
141
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
Thomas Witkowski's avatar
Thomas Witkowski committed
142
	  (!adoptFlag.isSet(INIT_MESH) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
143
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
144
	createMesh();      
Thomas Witkowski's avatar
Thomas Witkowski committed
145

146
      if (adoptProblem && 
Thomas Witkowski's avatar
Thomas Witkowski committed
147
	  (adoptFlag.isSet(INIT_MESH) ||
148
149
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
150
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
151
	componentMeshes = adoptProblem->componentMeshes;
152
153
154

	// If the adopt problem has fewer components than this problem, but only one
	// mesh for all component, than scal up the componentMeshes array.
Thomas Witkowski's avatar
Thomas Witkowski committed
155
	if (adoptProblem->getNumComponents() < nComponents) {
Thomas Witkowski's avatar
Thomas Witkowski committed
156
	  TEST_EXIT(meshes.size() == 1)("Daran muss ich noch arbeiten!\n");
157
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
158
	  componentMeshes.resize(nComponents);
159
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
160
	    componentMeshes[i] = componentMeshes[0];
161
	}
162
163
164
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    if (meshes.size() == 0) 
166
167
      WARNING("no mesh created\n");

168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
    // === create refinement/corasening-manager ===
    if (refinementManager != NULL && coarseningManager != NULL) { 
      WARNING("refinement-/coarseningmanager already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
          (!adoptFlag.isSet(INIT_MESH) &&
           (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
        createRefCoarseManager();

      if (adoptProblem && 
          (adoptFlag.isSet(INIT_MESH) ||
           adoptFlag.isSet(INIT_SYSTEM) ||
           adoptFlag.isSet(INIT_FE_SPACE))) {
        refinementManager = adoptProblem->refinementManager;
        coarseningManager = adoptProblem->coarseningManager;
      }
    }

    if (refinementManager == NULL || coarseningManager == NULL) 
      WARNING("no refinement-/coarseningmanager created\n");

189
    // === create fespace ===
190
    if (feSpaces.size() != 0) {
191
192
193
      WARNING("feSpaces already created\n");
    } else {
      if (initFlag.isSet(INIT_FE_SPACE) || 
194
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
195
	createFeSpace(NULL);
196

197
198
      if (adoptProblem &&
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
199
	feSpaces = adoptProblem->getFeSpaces();
200
	componentSpaces = adoptProblem->componentSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
201
	traverseInfo = adoptProblem->traverseInfo;
202
203
204

	// If the adopt problem has fewer components than this problem, but only one
	// fe space for all component, than scal up the componentSpaces array.
Thomas Witkowski's avatar
Thomas Witkowski committed
205
206
	if (adoptProblem->getNumComponents() < nComponents) {
	  TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n");
207
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
208
	  componentSpaces.resize(nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
209
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
210
	    componentSpaces[i] = componentSpaces[0];
211
212
	}

213
214
215
      }
    }

216
    if (feSpaces.size() == 0) 
217
218
219
      WARNING("no feSpace created\n");

    // === create system ===
220
    if (initFlag.isSet(INIT_SYSTEM)) 
221
      createMatricesAndVectors();
222
    
223
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
224
      solution = adoptProblem->getSolution();
225
      rhs = adoptProblem->getRhs();
Thomas Witkowski's avatar
Thomas Witkowski committed
226
      systemMatrix = adoptProblem->getSystemMatrix();
227
228
229
    }

    // === create solver ===
Thomas Witkowski's avatar
Thomas Witkowski committed
230
    if (solver) {
231
232
      WARNING("solver already created\n");
    } else {
233
      if (initFlag.isSet(INIT_SOLVER))
234
	createSolver();
235
      
236
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
239
240
241
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
242
    if (!solver) 
243
244
245
      WARNING("no solver created\n");

    // === create estimator ===
Thomas Witkowski's avatar
Thomas Witkowski committed
246
    if (initFlag.isSet(INIT_ESTIMATOR))
247
      createEstimator();
Thomas Witkowski's avatar
Thomas Witkowski committed
248
249

    if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
250
      estimator = adoptProblem->getEstimators();
251
252

    // === create marker ===
Thomas Witkowski's avatar
Thomas Witkowski committed
253
    if (initFlag.isSet(INIT_MARKER))
254
      createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
255

Thomas Witkowski's avatar
Thomas Witkowski committed
256
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
257
      marker = adoptProblem->getMarkers();
258
259

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
260
    if (initFlag.isSet(INIT_FILEWRITER))
261
262
263
264
265
      createFileWriter();
    
    // === read serialization and init mesh ===
    
    // There are two possiblities where the user can define a serialization
266
267
268
269
    // 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.
270
271

    int readSerialization = 0;
272
273
    string serializationFilename = "";
    Parameters::get("argv->rs", serializationFilename);
274
275
276
277
278
279

    // 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;

280
281
282
      Parameters::get(name + "->input->read serialization", readSerialization);
      Parameters::get(name + "->input->serialization with adaptinfo",
		      readSerializationWithAdaptInfo);
283
284
285
286
287
288

      // 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) {
289
290
	Parameters::get(name + "->input->serialization filename", 
			serializationFilename);
291
292
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

293
294
	// If AMDiS is compiled for parallel computations, the deserialization is
	// controlled by the parallel problem object.
295
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
296
	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
297
	ifstream in(serializationFilename.c_str());
298
299
300
301
302
303

	// Read the revision number of the AMDiS version which was used to create 
	// the serialization file.
	int revNumber = -1;
	SerUtil::deserialize(in, revNumber);

304
305
	deserialize(in);
	in.close();
306
#endif
307
308

	deserialized = true;
309
      } else {
310
	int globalRefinements = 0;
311

312
	// If AMDiS is compiled for parallel computations, the global refinements are
313
314
315
	// ignored here. Later, each rank will add the global refinements to its 
	// private mesh.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
316
317
	Parameters::get(meshes[0]->getName() + "->global refinements", 
			globalRefinements);
318
#endif
319

320
	bool initMesh = initFlag.isSet(INIT_MESH);
321

Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
	// Initialize the meshes if there is no serialization file.
	for (int i = 0; i < static_cast<int>(meshes.size()); i++)
324
	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
Thomas Witkowski's avatar
Thomas Witkowski committed
325
	    meshes[i]->initialize();	    	
326
327

	// === read value file and use it for the mesh values ===
328
329
	string valueFilename("");
	Parameters::get(meshes[0]->getName() + "->value file name", valueFilename); 
330
	if (valueFilename.length())
331
332
333
334
335
336
	  ValueReader::readValue(valueFilename,
				 meshes[0],
				 solution->getDOFVector(0),
				 meshes[0]->getMacroFileInfo());

	// === do global refinements ===
337

Thomas Witkowski's avatar
Thomas Witkowski committed
338
	for (unsigned int i = 0; i < meshes.size(); i++)
339
	  if (initMesh && meshes[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
340
	    refinementManager->globalRefine(meshes[i], globalRefinements);	
341
342
343
344
345
346
      }
    }

    doOtherStuff();
  }

347

348
  void ProblemStatSeq::createMesh() 
349
  {
350
    FUNCNAME("ProblemStat::createMesh()");
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
    componentMeshes.resize(nComponents);
353
354
355
356
357
358
    map<int, Mesh*> meshForRefinementSet;

    string meshName("");
    Parameters::get(name + "->mesh", meshName);
    TEST_EXIT(meshName != "")("No mesh name specified for \"%s->mesh\"!\n", 
			      name.c_str());
359
360

    int dim = 0;
361
362
363
    Parameters::get(name + "->dim", dim);
    TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
		   name.c_str());
364

365
    for (int i = 0; i < nComponents; i++) {
366
      int refSet = -1;
367
368
      Parameters::get(name + "->refinement set[" + lexical_cast<string>(i) + "]", 
		      refSet);
369
      if (refSet < 0)
370
	refSet = 0;
371
      
372
      if (meshForRefinementSet[refSet] == NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
373
	Mesh *newMesh = new Mesh(meshName, dim);
374
	meshForRefinementSet[refSet] = newMesh;
Thomas Witkowski's avatar
Thomas Witkowski committed
375
376
	meshes.push_back(newMesh);
	nMeshes++;
377
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
378
      componentMeshes[i] = meshForRefinementSet[refSet];
379
    }
380
381
382
383
384
385
386
387
388
389
390
  }


  void ProblemStatSeq::createRefCoarseManager() 
  {
    FUNCNAME("ProblemStat::createRefCoarseManager()");

    int dim = 0;
    Parameters::get(name + "->dim", dim);
    TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
                   name.c_str());
391
392

    switch (dim) {
393
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
394
      refinementManager = new RefinementManager1d();
395
      coarseningManager = new CoarseningManager1d();
396
397
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
398
      refinementManager = new RefinementManager2d();
399
      coarseningManager = new CoarseningManager2d();
400
401
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
402
      refinementManager = new RefinementManager3d();
403
      coarseningManager = new CoarseningManager3d();
404
405
406
407
408
409
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

410

411
  void ProblemStatSeq::createFeSpace(DOFAdmin *admin)
412
  {
413
    FUNCNAME("ProblemStat::createFeSpace()");
414

415
    map<pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
416
    int dim = -1;
417
    Parameters::get(name + "->dim", dim);
418
    TEST_EXIT(dim != -1)("no problem dimension specified!\n");
419

420
    componentSpaces.resize(nComponents, NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
421
    traverseInfo.resize(nComponents);
422

423
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
424
      int degree = 1;
425
426
      Parameters::get(name + "->polynomial degree[" + 
		      boost::lexical_cast<string>(i) + "]", degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
      TEST_EXIT(degree > 0)
	("Poynomial degree in component %d must be larger than zero!\n", i);
429
      TEST_EXIT(componentSpaces[i] == NULL)("feSpace already created\n");
430

431
      if (feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
434
	stringstream s;
	s << name << "->feSpace[" << i << "]";

435
436
	FiniteElemSpace *newFeSpace = 
	  FiniteElemSpace::provideFeSpace(admin, Lagrange::getLagrange(dim, degree),
Thomas Witkowski's avatar
Thomas Witkowski committed
437
					  componentMeshes[i], s.str());
438
	feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)] = newFeSpace;
439
	feSpaces.push_back(newFeSpace);
440
      }
441

442
      componentSpaces[i] = feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)];
Thomas Witkowski's avatar
Thomas Witkowski committed
443
444
445
    }

    for (int i = 0; i < nComponents; i++) {
446
      for (int j = 0; j < nComponents; j++)
447
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
448
      
449
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
450
451
452
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
453
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
454
      if (meshes[i]->getNumberOfDofs(VERTEX) == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
455
	DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
456
	ln_dof[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
457
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
458
459
460
461
      }
    }
  }

462

463
  void ProblemStatSeq::createMatricesAndVectors()
464
  {
465
    FUNCNAME("ProblemStat::createMatricesAndVectors()");
466
467
468

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

Thomas Witkowski's avatar
Thomas Witkowski committed
469
    systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
470
    systemMatrix->set(NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
471
472
    rhs = new SystemVector("rhs", componentSpaces, nComponents);
    solution = new SystemVector("solution", componentSpaces, nComponents);
473

Thomas Witkowski's avatar
Thomas Witkowski committed
474
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
475
      (*systemMatrix)[i][i] = new DOFMatrix(componentSpaces[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
476
					    componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
477
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
Thomas Witkowski's avatar
Thomas Witkowski committed
478

479
480
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], 
						 "rhs[" + lexical_cast<string>(i) + "]"));
481

Thomas Witkowski's avatar
Thomas Witkowski committed
482
      solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
483
      						      componentNames[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
486
487
488
    }
  }

489

490
  void ProblemStatSeq::createSolver()
491
  {
492
    FUNCNAME("ProblemStat::createSolver()");
493
494

    // === create solver ===
495
    string solverType("0");
496
497
    string initFileStr = name + "->solver";
    Parameters::get(initFileStr, solverType);
498
    OEMSolverCreator *solverCreator = 
499
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType, initFileStr));
500
    TEST_EXIT(solverCreator)("no solver type\n");
501
    solverCreator->setName(initFileStr);
Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
    solver = solverCreator->create();
    solver->initParameters();
504
505
  }

506

507
  void ProblemStatSeq::createEstimator()
508
  {
509
    FUNCNAME("ProblemStat::createEstimator()");
510
511

    // create and set leaf data prototype
Thomas Witkowski's avatar
Thomas Witkowski committed
512
    for (unsigned int i = 0; i < meshes.size(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
513
      meshes[i]->setElementDataPrototype
Thomas Witkowski's avatar
Thomas Witkowski committed
514
	(new LeafDataEstimatableVec(new LeafDataCoarsenableVec));
515

Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
    for (int i = 0; i < nComponents; i++) {
      TEST_EXIT(estimator[i] == NULL)("estimator already created\n");
518
519
      string estName = 
	name + "->estimator[" + boost::lexical_cast<string>(i) + "]";
520
521

      // === create estimator ===
522
523
      string estimatorType("0");
      Parameters::get(estName, estimatorType);
524
      EstimatorCreator *estimatorCreator = 
525
	dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType, estName));
Thomas Witkowski's avatar
Thomas Witkowski committed
526
      if (estimatorCreator) {
527
528
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
529
	estimatorCreator->setSolution(solution->getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
530
	estimator[i] = estimatorCreator->create();
531
532
533
      }


Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
536
537
538
      if (estimator[i])
	for (int j = 0; j < nComponents; j++)
	  estimator[i]->addSystem((*systemMatrix)[i][j], 
				  solution->getDOFVector(j), 
				  rhs->getDOFVector(j));      
539
540
541
    }
  }

542

543
  void ProblemStatSeq::createMarker()
544
  {
545
    FUNCNAME("ProblemStat::createMarker()");
546

Thomas Witkowski's avatar
Thomas Witkowski committed
547
    int nMarkersCreated = 0;
548

549
    for (int i = 0; i < nComponents; i++) {
550
      marker[i] = Marker::createMarker
551
	(name + "->marker[" + boost::lexical_cast<string>(i) + "]", i);
552

553
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
554
555
556
557
	nMarkersCreated++;

	// If there is more than one marker, and all components are defined
	// on the same mesh, the maximum marking has to be enabled.
558
 	if (nMarkersCreated > 1 && nMeshes == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
559
 	  marker[i]->setMaximumMarking(true);
560
561
562
563
      }
    }
  }

564

565
  void ProblemStatSeq::createFileWriter()
566
  {
567
    FUNCNAME("ProblemStat::createFileWriter()");
568
569
  
    // Create one filewriter for all components of the problem
570
571
572
    string numberedName  = name + "->output";
    string filename = "";
    Parameters::get(numberedName + "->filename", filename);
573
574

    if (filename != "") {
575
      vector< DOFVector<double>* > solutionList(nComponents);
576

577
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
578
	TEST_EXIT(componentMeshes[0] == componentMeshes[i])
579
580
	  ("All Meshes have to be equal to write a vector file.\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
581
	solutionList[i] = solution->getDOFVector(i);
582
583
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
584
      fileWriters.push_back(new FileWriter(numberedName,
585
586
					   componentMeshes[0],
					   solutionList));
587
588
589
    }

    // Create own filewriters for each components of the problem
590
    for (int i = 0; i < nComponents; i++) {
591
      numberedName = name + "->output[" + boost::lexical_cast<string>(i) + "]";
592
      filename = "";
593
      Parameters::get(numberedName + "->filename", filename);
594

595
      if (filename != "")
Thomas Witkowski's avatar
Thomas Witkowski committed
596
	fileWriters.push_back(new FileWriter(numberedName, 
597
					     componentMeshes[i], 
598
					     solution->getDOFVector(i)));      
599
600
601
    }

    int writeSerialization = 0;
602
    Parameters::get(name + "->output->write serialization", writeSerialization);
603
    if (writeSerialization)
604
      fileWriters.push_back(new Serializer<ProblemStatSeq>(this));
605
606
  }

607

608
  void ProblemStatSeq::doOtherStuff()
609
  {}
610

611

612
613
614
  void ProblemStatSeq::solve(AdaptInfo *adaptInfo, 	       
			     bool createMatrixData,
			     bool storeMatrixData)
615
616
617
  {
    FUNCNAME("Problem::solve()");

Thomas Witkowski's avatar
Thomas Witkowski committed
618
    if (!solver) {
619
620
621
622
623
      WARNING("no solver\n");
      return;
    }

    clock_t first = clock();
624
625
626

    solver->solveSystem(solverMatrix, *solution, *rhs, 
			createMatrixData, storeMatrixData);
627

Thomas Witkowski's avatar
Thomas Witkowski committed
628
629
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
630

631
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
635
636
  }

637

638
  void ProblemStatSeq::estimate(AdaptInfo *adaptInfo) 
639
  {
640
    FUNCNAME("ProblemStat::estimate()");
641

642
643
644
645
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#else
    clock_t first = clock();
646
647
#endif

648
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
649
      computeError(adaptInfo);
650
651
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
652
	Estimator *scalEstimator = estimator[i];
653
654
	
	if (scalEstimator) {
655
656
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
657
	  scalEstimator->estimate(adaptInfo->getTimestep());
658

659
660
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
661
662
663
664
665

	  if (adaptInfo->getRosenbrockMode() == false) {
	    adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	    adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	  }
666
	}
667
668
669
      }
    }

670
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
671
    MPI::COMM_WORLD.Barrier();
672
673
    INFO(info, 8)("estimation of the error needed %.5f seconds\n", 
		  MPI::Wtime() - first);
674
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
675
676
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
677
#endif
678
679
  }

680

681
  Flag ProblemStatSeq::markElements(AdaptInfo *adaptInfo) 
682
  {
683
    FUNCNAME("ProblemStat::markElements()");
684

685
    TEST_EXIT_DBG(static_cast<unsigned int>(nComponents) == marker.size())
Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
      ("Wrong number of markers!\n");

688
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
689
    for (int i = 0; i < nComponents; i++)
690
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
691
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
692
    
693
694
695
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
696

697
  Flag ProblemStatSeq::refineMesh(AdaptInfo *adaptInfo) 
698
  {
699
    FUNCNAME("ProblemStat::refineMesh()");
700

Thomas Witkowski's avatar
Thomas Witkowski committed
701
    int nMeshes = static_cast<int>(meshes.size());
702
    Flag refineFlag = 0;
703
704
705
706
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

707
708
709
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
710

711
  Flag ProblemStatSeq::coarsenMesh(AdaptInfo *adaptInfo) 
712
  {
713
    FUNCNAME("ProblemStat::coarsenMesh()");
714

Thomas Witkowski's avatar
Thomas Witkowski committed
715
    int nMeshes = static_cast<int>(meshes.size());
716
    Flag coarsenFlag = 0;
717
718
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
719
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
720
721
722
723

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
724

725
  Flag ProblemStatSeq::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
726
  {
727
    FUNCNAME("ProblemStat::oneIteration()");
728

729
730
    for (int i = 0; i < nComponents; i++)
      if (adaptInfo->spaceToleranceReached(i))
731
 	adaptInfo->allowRefinement(false, i);
732
      else
733
 	adaptInfo->allowRefinement(true, i);	
734
735
736
737

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

Thomas Witkowski's avatar
Thomas Witkowski committed
738

739
740
  void ProblemStatSeq::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
					 bool asmMatrix, bool asmVector)
741
  {
742
    FUNCNAME("ProblemStat::buildAfterCoarsen()");
743

744
    if (dualMeshTraverseRequired()) {
745
746
747
748
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      ERROR_EXIT("Dual mesh assemble does not work in parallel code!\n");
#endif

749
      dualAssemble(adaptInfo, flag, asmMatrix, asmVector);
750
751
      return;
    }
752

753
754
755
756
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#endif

757

Thomas Witkowski's avatar
Thomas Witkowski committed
758
759
760
761
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


762
763
764
765
766
767
768
769
770
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    MPI::COMM_WORLD.Barrier();
    INFO(info, 8)("dof compression needed %.5f seconds\n", 
		  MPI::Wtime() - first);

    first = MPI::Wtime();
#else
    clock_t first = clock();
#endif
771

772
773
774

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
777
778
779
780
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
781
782
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
783
    if (useGetBound)
784
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786

  
Thomas Witkowski's avatar
Thomas Witkowski committed
787
    traverseInfo.updateStatus();
788
789
790
    // Used to calculate the overall number of non zero entries.
    int nnz = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
791
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
792

793
      MSG("%d DOFs for %s\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
794
795
	  componentSpaces[rowComponent]->getAdmin()->getUsedDofs(), 
	  componentSpaces[rowComponent]->getName().c_str());
796

Thomas Witkowski's avatar
Thomas Witkowski committed
797
      rhs->getDOFVector(rowComponent)->set(0.0);
798

Thomas Witkowski's avatar
Thomas Witkowski committed
799
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
800

801
	if (writeAsmInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
802
	  MSG("--------- %d %d -------------\n", rowComponent, colComponent);
803