ProblemVec.cc 51.8 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
17
18
19
20
21
22
23
24
25
#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"
26
#include "io/FileWriter.h"
27
28
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
29
#include "DualTraverse.h"
30
31
32
33
34
35
#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
36
#include "Flag.h"
37
38
#include "io/VtkWriter.h"
#include "io/ValueReader.h"
39
#include "ProblemVecDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
40
#include "Debug.h"
41
42
43

namespace AMDiS {

44
45
  using boost::lexical_cast;

46
47
48
49
50
  void ProblemVec::initialize(Flag initFlag,
			      ProblemVec *adoptProblem,
			      Flag adoptFlag)
  {
    FUNCNAME("ProblemVec::initialize()");
51

52
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
53
    if (meshes.size() != 0) { 
54
55
56
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
Thomas Witkowski's avatar
Thomas Witkowski committed
57
	  (!adoptFlag.isSet(INIT_MESH) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
58
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
59
	createMesh();      
Thomas Witkowski's avatar
Thomas Witkowski committed
60

61
      if (adoptProblem && 
Thomas Witkowski's avatar
Thomas Witkowski committed
62
	  (adoptFlag.isSet(INIT_MESH) ||
63
64
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
65
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
66
	componentMeshes = adoptProblem->componentMeshes;
Thomas Witkowski's avatar
Thomas Witkowski committed
67
68
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
69
70
71

	// 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
72
	if (adoptProblem->getNumComponents() < nComponents) {
Thomas Witkowski's avatar
Thomas Witkowski committed
73
	  TEST_EXIT(meshes.size() == 1)("Daran muss ich noch arbeiten!\n");
74
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
75
	  componentMeshes.resize(nComponents);
76
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
77
	    componentMeshes[i] = componentMeshes[0];
78
79
	}

80
81
82
83

	// If the problem adopts the mesh but creates an own FE space.
	if (initFlag.isSet(INIT_FE_SPACE)) {
	}
84
85
86
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
87
    if (meshes.size() == 0) 
88
89
90
      WARNING("no mesh created\n");

    // === create fespace ===
91
    if (feSpaces.size() != 0) {
92
93
94
      WARNING("feSpaces already created\n");
    } else {
      if (initFlag.isSet(INIT_FE_SPACE) || 
95
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
96
	createFeSpace(NULL);
97

98
99
      if (adoptProblem &&
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
100
	feSpaces = adoptProblem->getFeSpaces();
101
	componentSpaces = adoptProblem->componentSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
102
	traverseInfo = adoptProblem->traverseInfo;
103
104
105

	// 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
106
107
	if (adoptProblem->getNumComponents() < nComponents) {
	  TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n");
108
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
109
	  componentSpaces.resize(nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
110
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
111
	    componentSpaces[i] = componentSpaces[0];
112
113
	}

114
115
116
      }
    }

117
    if (feSpaces.size() == 0) 
118
119
120
      WARNING("no feSpace created\n");

    // === create system ===
121
    if (initFlag.isSet(INIT_SYSTEM)) 
122
      createMatricesAndVectors();
123
    
124
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
125
      solution = adoptProblem->getSolution();
126
      rhs = adoptProblem->getRhs();
Thomas Witkowski's avatar
Thomas Witkowski committed
127
      systemMatrix = adoptProblem->getSystemMatrix();
128
129
130
    }

    // === create solver ===
Thomas Witkowski's avatar
Thomas Witkowski committed
131
    if (solver) {
132
133
134
135
136
137
      WARNING("solver already created\n");
    } else {
      if (initFlag.isSet(INIT_SOLVER)) {
	createSolver();
      } 
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
140
141
142
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
143
    if (!solver) 
144
145
146
      WARNING("no solver created\n");

    // === create estimator ===
Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (initFlag.isSet(INIT_ESTIMATOR))
148
      createEstimator();
Thomas Witkowski's avatar
Thomas Witkowski committed
149
150

    if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
151
      estimator = adoptProblem->getEstimators();
152
153

    // === create marker ===
Thomas Witkowski's avatar
Thomas Witkowski committed
154
    if (initFlag.isSet(INIT_MARKER))
155
      createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
156

Thomas Witkowski's avatar
Thomas Witkowski committed
157
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
158
      marker = adoptProblem->getMarkers();
159
160

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
161
    if (initFlag.isSet(INIT_FILEWRITER))
162
163
164
165
166
167
168
169
170
171
172
      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;
173
    std::string serializationFilename = "";
174
175
176
177
178
179
180
    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;

181
      GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
Thomas Witkowski's avatar
Thomas Witkowski committed
182
      GET_PARAMETER(0, name + "->input->serialization with adaptinfo", "%d",
183
184
185
186
187
188
189
		    &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) {
Thomas Witkowski's avatar
Thomas Witkowski committed
190
	GET_PARAMETER(0, name + "->input->serialization filename", 
191
192
193
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

194
195
	// If AMDiS is compiled for parallel computations, the deserialization is
	// controlled by the parallel problem object.
196
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
197
	std::ifstream in(serializationFilename.c_str());
198
199
	deserialize(in);
	in.close();
200
201
202

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
#endif
203
204

	deserialized = true;
205
      } else {
206
	int globalRefinements = 0;
207

208
	// If AMDiS is compiled for parallel computations, the global refinements are
209
210
211
	// ignored here. Later, each rank will add the global refinements to its 
	// private mesh.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
212
	GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d", 
213
		      &globalRefinements);
214
#endif
215

216
	bool initMesh = initFlag.isSet(INIT_MESH);
217

Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
	// Initialize the meshes if there is no serialization file.
	for (int i = 0; i < static_cast<int>(meshes.size()); i++)
220
	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
Thomas Witkowski's avatar
Thomas Witkowski committed
221
	    meshes[i]->initialize();	    	
222
223
224
225

	// === read value file and use it for the mesh values ===
	std::string valueFilename("");
	GET_PARAMETER(0, meshes[0]->getName() + "->value file name", &valueFilename); 
226
	if (valueFilename.length())
227
228
229
230
231
232
	  ValueReader::readValue(valueFilename,
				 meshes[0],
				 solution->getDOFVector(0),
				 meshes[0]->getMacroFileInfo());

	// === do global refinements ===
233

Thomas Witkowski's avatar
Thomas Witkowski committed
234
	for (unsigned int i = 0; i < meshes.size(); i++)
235
	  if (initMesh && meshes[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
236
	    refinementManager->globalRefine(meshes[i], globalRefinements);	
237
238
239
240
241
242
      }
    }

    doOtherStuff();
  }

243

244
245
246
247
  void ProblemVec::createMesh() 
  {
    FUNCNAME("ProblemVec::createMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
248
    componentMeshes.resize(nComponents);
249
    std::map<int, Mesh*> meshForRefinementSet;
250

251
    std::string meshName("");
Thomas Witkowski's avatar
Thomas Witkowski committed
252
    GET_PARAMETER(0, name + "->mesh", &meshName);
253
    TEST_EXIT(meshName != "")("no mesh name specified\n");
254
    int dim = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
255
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
256
    TEST_EXIT(dim)("no problem dimension specified!\n");
257

258
    for (int i = 0; i < nComponents; i++) {
259
      int refSet = -1;
260
261
      GET_PARAMETER(0, name + "->refinement set[" + 
		       lexical_cast<std::string>(i) + "]", "%d", &refSet);
262
      if (refSet < 0)
263
	refSet = 0;
264
      
265
      if (meshForRefinementSet[refSet] == NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
266
	Mesh *newMesh = new Mesh(meshName, dim);
267
	meshForRefinementSet[refSet] = newMesh;
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
	meshes.push_back(newMesh);
	nMeshes++;
270
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
271
      componentMeshes[i] = meshForRefinementSet[refSet];
272
    }
273
274

    switch (dim) {
275
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
278
279
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
282
283
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
284
285
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
286
287
288
289
290
291
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

292

293
  void ProblemVec::createFeSpace(DOFAdmin *admin)
294
  {
295
    FUNCNAME("ProblemVec::createFeSpace()");
296

297
    std::map<std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
298
    int dim = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
299
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
300
    TEST_EXIT(dim != -1)("no problem dimension specified!\n");
301

302
    componentSpaces.resize(nComponents, NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
303
    traverseInfo.resize(nComponents);
304

305
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
306
      int degree = 1;
307
308
      GET_PARAMETER(0, name + "->polynomial degree[" + 
		    boost::lexical_cast<std::string>(i) + "]","%d", &degree);
309

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

Thomas Witkowski's avatar
Thomas Witkowski committed
312
      if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
313
314
315
	stringstream s;
	s << name << "->feSpace[" << i << "]";

316
317
	FiniteElemSpace *newFeSpace = 
	  FiniteElemSpace::provideFeSpace(admin, Lagrange::getLagrange(dim, degree),
Thomas Witkowski's avatar
Thomas Witkowski committed
318
					  componentMeshes[i], s.str());
319
320
	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFeSpace;
	feSpaces.push_back(newFeSpace);
321
      }
322

Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
325
326
      componentSpaces[i] = feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
    }

    for (int i = 0; i < nComponents; i++) {
327
      for (int j = 0; j < nComponents; j++)
328
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
329
      
330
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
331
332
333
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
335
      if (meshes[i]->getNumberOfDofs(VERTEX) == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
336
	DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
337
	ln_dof[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
338
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
339
340
341
342
      }
    }
  }

343

344
345
346
347
348
349
  void ProblemVec::createMatricesAndVectors()
  {
    FUNCNAME("ProblemVec::createMatricesAndVectors()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
350
    systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
351
    systemMatrix->set(NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
    rhs = new SystemVector("rhs", componentSpaces, nComponents);
    solution = new SystemVector("solution", componentSpaces, nComponents);
354

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      (*systemMatrix)[i][i] = new DOFMatrix(componentSpaces[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
357
					    componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
358
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
Thomas Witkowski's avatar
Thomas Witkowski committed
359

360
      std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]";
Thomas Witkowski's avatar
Thomas Witkowski committed
361
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName));
362
363

      numberedName = "solution[" + boost::lexical_cast<std::string>(i) + "]";
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
      solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
      						      numberedName));
Thomas Witkowski's avatar
Thomas Witkowski committed
366
367
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
368
369
370
    }
  }

371

372
373
374
375
376
  void ProblemVec::createSolver()
  {
    FUNCNAME("ProblemVec::createSolver()");

    // === create solver ===
377
    std::string solverType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
378
    GET_PARAMETER(0, name + "->solver", &solverType);
379
380
    OEMSolverCreator *solverCreator = 
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
381
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
384
    solverCreator->setName(name + "->solver");
    solver = solverCreator->create();
    solver->initParameters();
385
386
  }

387

388
389
390
391
392
  void ProblemVec::createEstimator()
  {
    FUNCNAME("ProblemVec::createEstimator()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
    for (int i = 0; i < nComponents; i++) {
      TEST_EXIT(estimator[i] == NULL)("estimator already created\n");
399
400
      std::string estName = 
	name + "->estimator[" + boost::lexical_cast<std::string>(i) + "]";
401
402

      // === create estimator ===
403
      std::string estimatorType("0");
404
405
      GET_PARAMETER(0, estName, &estimatorType);
      EstimatorCreator *estimatorCreator = 
Thomas Witkowski's avatar
nix.  
Thomas Witkowski committed
406
	dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
407
      if (estimatorCreator) {
408
409
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
410
	estimatorCreator->setSolution(solution->getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
411
	estimator[i] = estimatorCreator->create();
412
413
414
      }


Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
      if (estimator[i])
	for (int j = 0; j < nComponents; j++)
	  estimator[i]->addSystem((*systemMatrix)[i][j], 
				  solution->getDOFVector(j), 
				  rhs->getDOFVector(j));      
420
421
422
    }
  }

423

424
425
426
427
  void ProblemVec::createMarker()
  {
    FUNCNAME("ProblemVec::createMarker()");

Thomas Witkowski's avatar
Thomas Witkowski committed
428
    int nMarkersCreated = 0;
429

430
    for (int i = 0; i < nComponents; i++) {
431
432
433
      marker[i] = Marker::createMarker
	(name + "->marker[" + boost::lexical_cast<std::string>(i) + "]", i);

434
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
437
438
	nMarkersCreated++;

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

445

446
447
448
449
450
  void ProblemVec::createFileWriter()
  {
    FUNCNAME("ProblemVec::createFileWriter()");
  
    // Create one filewriter for all components of the problem
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    std::string numberedName  = name + "->output";
452
    std::string filename = "";
453
454
455
    GET_PARAMETER(0, numberedName + "->filename", &filename);

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
462
	solutionList[i] = solution->getDOFVector(i);
463
464
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
465
      fileWriters.push_back(new FileWriter(numberedName,
466
467
					   componentMeshes[0],
					   solutionList));
468
469
470
    }

    // Create own filewriters for each components of the problem
471
    for (int i = 0; i < nComponents; i++) {
472
      numberedName = name + "->output[" + boost::lexical_cast<std::string>(i) + "]";
473
474
475
      filename = "";
      GET_PARAMETER(0, numberedName + "->filename", &filename);

476
      if (filename != "")
Thomas Witkowski's avatar
Thomas Witkowski committed
477
	fileWriters.push_back(new FileWriter(numberedName, 
478
					     componentMeshes[i], 
479
					     solution->getDOFVector(i)));      
480
481
482
483
    }

    // Check for serializer
    int writeSerialization = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
484
    GET_PARAMETER(0, name + "->write serialization", "%d", &writeSerialization);
485
    if (writeSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
      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());
488
489
490
      ERROR_EXIT("Usage of an obsolete parameter (see message above)!\n");
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
491
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
492

493
    if (writeSerialization)
Thomas Witkowski's avatar
Thomas Witkowski committed
494
      fileWriters.push_back(new Serializer<ProblemVec>(this));
495
496
  }

497

498
499
500
501
  void ProblemVec::doOtherStuff()
  {
  }

502

503
  void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
504
505
506
  {
    FUNCNAME("Problem::solve()");

Thomas Witkowski's avatar
Thomas Witkowski committed
507
    if (!solver) {
508
509
510
511
512
      WARNING("no solver\n");
      return;
    }

    clock_t first = clock();
513
    
514
    solver->solveSystem(solverMatrix, *solution, *rhs);
515

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

519
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521
522
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
523
524
  }

525

526
527
528
529
  void ProblemVec::estimate(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::estimate()");

530
531
532
533
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#else
    clock_t first = clock();
534
535
#endif

536
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
537
      computeError(adaptInfo);
538
539
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
540
	Estimator *scalEstimator = estimator[i];
541
542
	
	if (scalEstimator) {
543
544
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
545
	  scalEstimator->estimate(adaptInfo->getTimestep());
546

547
548
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
549
550
551
552
553

	  if (adaptInfo->getRosenbrockMode() == false) {
	    adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	    adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	  }
554
	}
555
556
557
      }
    }

558
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
559
    MPI::COMM_WORLD.Barrier();
560
561
    INFO(info, 8)("estimation of the error needed %.5f seconds\n", 
		  MPI::Wtime() - first);
562
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
563
564
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
565
#endif
566
567
  }

568

569
570
571
572
  Flag ProblemVec::markElements(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::markElements()");

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

576
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
577
    for (int i = 0; i < nComponents; i++) {
578
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
579
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
580
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
581
582
	WARNING("No marker for component %d\n", i);              
    }
583
    
584
585
586
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
587

588
589
590
591
  Flag ProblemVec::refineMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::refineMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
592
    int nMeshes = static_cast<int>(meshes.size());
593
    Flag refineFlag = 0;
594
595
596
597
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

598
599
600
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
601

602
603
604
605
  Flag ProblemVec::coarsenMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::coarsenMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
606
    int nMeshes = static_cast<int>(meshes.size());
607
    Flag coarsenFlag = 0;
608
609
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
610
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
611
612
613
614

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
615

616
617
618
619
  Flag ProblemVec::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ProblemVec::oneIteration()");

620
621
622
623
624
    for (int i = 0; i < nComponents; i++)
      if (adaptInfo->spaceToleranceReached(i))
	adaptInfo->allowRefinement(false, i);
      else
	adaptInfo->allowRefinement(true, i);	
625
626
627
628

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

Thomas Witkowski's avatar
Thomas Witkowski committed
629

630
631
  void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				     bool asmMatrix, bool asmVector)
632
633
634
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

635
636
637
//     buildAfterCoarsen_sebastianMode(adaptInfo, flag);
//     return;

638
    if (dualMeshTraverseRequired()) {
639
640
641
642
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      ERROR_EXIT("Dual mesh assemble does not work in parallel code!\n");
#endif

643
      dualAssemble(adaptInfo, flag, asmMatrix, asmVector);
644
645
      return;
    }
646

647
648
649
650
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#endif

651

Thomas Witkowski's avatar
Thomas Witkowski committed
652
653
654
655
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


656
657
658
659
660
661
662
663
664
#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
665

666
667
668

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
671
672
673
674
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
675
676
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
677
    if (useGetBound)
678
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680

  
Thomas Witkowski's avatar
Thomas Witkowski committed
681
    traverseInfo.updateStatus();
682
683
684
    // Used to calculate the overall number of non zero entries.
    int nnz = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
685
    for (int i = 0; i < nComponents; i++) {
686

687
      MSG("%d DOFs for %s\n", 
688
	  componentSpaces[i]->getAdmin()->getUsedDofs(), 
689
690
691
692
	  componentSpaces[i]->getName().c_str());

      rhs->getDOFVector(i)->set(0.0);

693
      for (int j = 0; j < nComponents; j++) {
694

695
696
	if (writeAsmInfo)
	  MSG("--------- %d %d -------------\n", i, j);
697

698
699
700
701
	// 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).
702
	DOFMatrix *matrix = (asmMatrix ? (*systemMatrix)[i][j] : NULL);
703

704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
	if (writeAsmInfo && matrix) {
	  for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
	       it != matrix->getOperatorsEnd(); ++it) {
	    Assembler *assembler = (*it)->getAssembler();
	    if (assembler->getZeroOrderAssembler()) 
	      std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
	    if (assembler->getFirstOrderAssembler(GRD_PSI)) 
	      std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
	    if (assembler->getFirstOrderAssembler(GRD_PHI)) 
	      std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
	    if (assembler->getSecondOrderAssembler()) 
	      std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
	  }
	}

719
720
721
	if (matrix) 
	  matrix->calculateNnz();
	
722
723
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
724
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
725
	  assembleMatrix = false;
726
727
728
729
730
731
732
	} else if (matrix) {
	  matrix->getBaseMatrix().
	    change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
		       componentSpaces[j]->getAdmin()->getUsedSize());

	  set_to_zero(matrix->getBaseMatrix());
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
733
734

	// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
735
	if (!matrix || !assembleMatrix)
736
737
	  assembleMatrix = false;

738
739
740
	// 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.
741
742
743
744
	if (!assembleMatrix && i != j) {
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();

745
	  continue;
746
	}
747

748
749
750
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

751
752
753
754
755
756
757
758
	// The simplest case: either the right hand side has no operaters, no aux
	// fe spaces, or all aux fe spaces are equal to the row and col fe space.
	assembleOnOneMesh(componentSpaces[i],
			  assembleFlag,
			  assembleMatrix ? matrix : NULL,
			  ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
      
	assembledMatrix[i][j] = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
759

760
761
	if (assembleMatrix)
	  matrix->finishInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
762

763
764
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	
765

766
	if (matrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
767
	  nnz += matrix->getBaseMatrix().nnz();
768
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
769

770
771
      // And now assemble boundary conditions on the vectors
      assembleBoundaryConditions(rhs->getDOFVector(i),
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
774
 				 solution->getDOFVector(i),
 				 componentMeshes[i],
 				 assembleFlag);     
775
    }  
776

777
    if (asmMatrix) {
778
      solverMatrix.setMatrix(*systemMatrix);
779
      createPrecon();
Thomas Witkowski's avatar
Thomas Witkowski committed
780

781
      INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
782
    }
783

784
785
786
787
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    MPI::COMM_WORLD.Barrier();
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  MPI::Wtime() - first);
788
789
790
#else
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));    
791
#endif
792
  }
793

794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
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
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935