ProblemVec.cc 57.5 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"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "TraverseParallel.h"
38
39
#include "io/VtkWriter.h"
#include "io/ValueReader.h"
40
#include "ProblemVecDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
41
#include "Debug.h"
42
43
44

namespace AMDiS {

45
46
  using boost::lexical_cast;

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

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

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

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

81
82
83
84

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

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

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

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

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

115
116
117
      }
    }

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

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

    // === create solver ===
Thomas Witkowski's avatar
Thomas Witkowski committed
132
    if (solver) {
133
134
135
136
137
138
      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
139
140
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
141
142
143
      }
    }

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

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

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

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

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

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

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

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

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

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

209
	// If AMDiS is compiled for parallel computations, the global refinements are
210
211
212
	// 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
213
	GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d", 
214
		      &globalRefinements);
215
#endif
216

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

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

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

	// === do global refinements ===
234

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

    doOtherStuff();
  }

244

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

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

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

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

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

293

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

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

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

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

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

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

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

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

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

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

344

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

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

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

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

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

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

372

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

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

388

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

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

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

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


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

424

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

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

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

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

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

446

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

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

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

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

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

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

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

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

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

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

498

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

503

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

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

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

    clock_t first = clock();
518
    
519
    solver->solveSystem(solverMatrix, *solution, *rhs);
520

521
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
522
523
    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);
524
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
525
526
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
527
528
#endif

529
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
530
531
532
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
533
534
  }

535

536
537
538
539
  void ProblemVec::estimate(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::estimate()");

540
541
542
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#else
543
#ifdef _OPENMP
544
545
546
547
    double first = omp_get_wtime();
#else
    clock_t first = clock();
#endif
548
549
#endif

550
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
551
      computeError(adaptInfo);
552
553
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
554
	Estimator *scalEstimator = estimator[i];
555
556
	
	if (scalEstimator) {
557
558
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
559
	  scalEstimator->estimate(adaptInfo->getTimestep());
560

561
562
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
563
564
565
566
567

	  if (adaptInfo->getRosenbrockMode() == false) {
	    adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	    adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	  }
568
	}
569
570
571
      }
    }

572
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
573
    MPI::COMM_WORLD.Barrier();
574
575
576
    INFO(info, 8)("estimation of the error needed %.5f seconds\n", 
		  MPI::Wtime() - first);
#else
577
#ifdef _OPENMP
578
579
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  omp_get_wtime() - first);
580
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
581
582
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
583
#endif
584
#endif
585
586
  }

587

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

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

595
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
596
    for (int i = 0; i < nComponents; i++) {
597
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
598
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
599
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
600
601
	WARNING("No marker for component %d\n", i);              
    }
602
    
603
604
605
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
606

607
608
609
610
  Flag ProblemVec::refineMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::refineMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
611
    int nMeshes = static_cast<int>(meshes.size());
612
    Flag refineFlag = 0;
613
614
615
616
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

617
618
619
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
620

621
622
623
624
  Flag ProblemVec::coarsenMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::coarsenMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
625
    int nMeshes = static_cast<int>(meshes.size());
626
    Flag coarsenFlag = 0;
627
628
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
629
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
630
631
632
633

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
634

635
636
637
638
  Flag ProblemVec::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ProblemVec::oneIteration()");

639
640
641
642
643
    for (int i = 0; i < nComponents; i++)
      if (adaptInfo->spaceToleranceReached(i))
	adaptInfo->allowRefinement(false, i);
      else
	adaptInfo->allowRefinement(true, i);	
644
645
646
647

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

Thomas Witkowski's avatar
Thomas Witkowski committed
648

649
650
  void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				     bool asmMatrix, bool asmVector)
651
652
653
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

654
655
656
//     buildAfterCoarsen_sebastianMode(adaptInfo, flag);
//     return;

657
    if (dualMeshTraverseRequired()) {
658
659
660
661
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      ERROR_EXIT("Dual mesh assemble does not work in parallel code!\n");
#endif

662
      dualAssemble(adaptInfo, flag, asmMatrix, asmVector);
663
664
      return;
    }
665

666
667
668
669
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#endif

670

Thomas Witkowski's avatar
Thomas Witkowski committed
671
672
673
674
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


675
676
677
678
679
680
681
#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
682
#ifdef _OPENMP
683
684
685
686
    double first = omp_get_wtime();
#else
    clock_t first = clock();
#endif
687
688
#endif

689
690
691

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
694
695
696
697
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
698
699
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
700
    if (useGetBound)
701
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
702
703

  
Thomas Witkowski's avatar
Thomas Witkowski committed
704
    traverseInfo.updateStatus();
705
706
707
    // Used to calculate the overall number of non zero entries.
    int nnz = 0;

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

710
      MSG("%d DOFs for %s\n", 
711
	  componentSpaces[i]->getAdmin()->getUsedDofs(), 
712
713
714
715
	  componentSpaces[i]->getName().c_str());

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

716
      for (int j = 0; j < nComponents; j++) {
717

718
719
	if (writeAsmInfo)
	  MSG("--------- %d %d -------------\n", i, j);
720

721
722
723
724
	// 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).
725
	DOFMatrix *matrix = (asmMatrix ? (*systemMatrix)[i][j] : NULL);
726

727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
	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;
	  }
	}

742
743
744
	if (matrix) 
	  matrix->calculateNnz();
	
745
746
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
747
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
748
	  assembleMatrix = false;
749
750
751
752
753
754
755
	} 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
756
757

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

761
762
763
	// 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.
764
765
766
767
	if (!assembleMatrix && i != j) {
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();

768
	  continue;
769
	}
770

771
772
773
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

774
775
776
777
778
779
780
781
	// 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
782

783
784
	if (assembleMatrix)
	  matrix->finishInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
785

786
787
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	
788

789
	if (matrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
790
	  nnz += matrix->getBaseMatrix().nnz();
791
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
792

793
794
      // And now assemble boundary conditions on the vectors
      assembleBoundaryConditions(rhs->getDOFVector(i),
Thomas Witkowski's avatar
Thomas Witkowski committed
795
796
797
 				 solution->getDOFVector(i),
 				 componentMeshes[i],
 				 assembleFlag);     
798
    }  
799

800
    if (asmMatrix) {
801
      solverMatrix.setMatrix(*systemMatrix);
802
      createPrecon();
Thomas Witkowski's avatar
Thomas Witkowski committed
803

804
      INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
805
    }
806

807
808
809
810
811
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    MPI::COMM_WORLD.Barrier();
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  MPI::Wtime() - first);
#else
812
#ifdef _OPENMP
813
814
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
		  omp_get_wtime() - wtime);
815
816
817
#else
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));    
818
819
#endif
#endif
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
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
  void ProblemVec::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

    clock_t first = clock();
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

    for (int i = 0; i < static_cast<int>(meshes.size()); i++)
      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;

    traverseInfo.updateStatus();

    // Used to calculate the overall number of non zero entries.
    int nnz = 0;  


    /// === INITIALIZE ===

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

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

      for (int j = 0; j < nComponents; j++) {
	// 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).
	DOFMatrix *matrix = (*systemMatrix)[i][j];

	if (matrix) 
	  matrix->calculateNnz();
	
	// 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;
	} else if (matrix) {
	  matrix->getBaseMatrix().
	    change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
		       componentSpaces[j]->getAdmin()->getUsedSize());

	  set_to_zero(matrix->getBaseMatrix());	  
	}

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

	// 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) {
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();

	  continue;
	}

	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

	if (matrix && assembleMatrix) 
	  matrix->startInsertion(matrix->getNnz());
      }
    }

    // === TRAVERSE ===

    Mesh *mesh = componentMeshes[0];
    const FiniteElemSpace *feSpace = componentSpaces[0];
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    ElementMatrix elMat(basisFcts->getNumber(), basisFcts->getNumber());
    ElementMatrix tmpElMat(elMat);
    ElementVector elVec(basisFcts->getNumber());
    ElementVector tmpElVec(elVec);
    TraverseStack stack;
    BoundaryType *bound = 
      useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
    while (elInfo) {
      if (useGetBound)
	basisFcts->getBound(elInfo, bound);

      for (std::map<Operator*, std::vector<OperatorPos> >::iterator opIt = operators.begin();
	   opIt != operators.end(); ++opIt) {
	if (opIt->first->getNeedDualTraverse() == true)
	  continue;

	if (opFlags[opIt->first].isSet(Operator::MATRIX_OPERATOR)) {
	  set_to_zero(elMat);
	  opIt->first->getElementMatrix(elInfo, elMat, 1.0);
	}
	if (opFlags[opIt->first].isSet(Operator::VECTOR_OPERATOR)) {
	  set_to_zero(elVec);
	  opIt->first->getElementVector(elInfo, elVec, 1.0);
	}
	
	for (std::vector<OperatorPos>::iterator posIt = opIt->second.begin();
	     posIt != opIt->second.end(); ++posIt) {

	  if (posIt->operatorType.isSet(Operator::MATRIX_OPERATOR)) {
	    if (!posIt->factor || *(posIt->factor) == 1.0) {
	      (*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(elMat, bound, elInfo, NULL);
	    } else {
	      tmpElMat = *(posIt->factor) * elMat;
	      (*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(tmpElMat, bound, elInfo, NULL);
	    }
	  }
	  
	  if (posIt->operatorType.isSet(Operator::VECTOR_OPERATOR)) {
	    if (!posIt->factor || *(posIt->factor) == 1.0) {
	      rhs->getDOFVector(posIt->row)->addElementVector(1.0, elVec, bound, elInfo);
	    } else {
	      tmpElVec = *(posIt->factor) * elVec;
	      rhs->getDOFVector(posIt->row)->addElementVector(1.0, tmpElVec, bound, elInfo);
	    }
	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    } 
      
    if (useGetBound)
      delete [] bound;

    // === FINELIZE ===

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	bool assembleMatrix = true;
	DOFMatrix *matrix = (*systemMatrix)[i][j];

	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j])
	  assembleMatrix = false;
	if (!matrix || !assembleMatrix)
	  assembleMatrix = false;
	if (!assembleMatrix && i != j)
	  continue;

	assembledMatrix[i][j] = true;

	if (assembleMatrix) {
	  matrix->removeRowsWithDBC(matrix->getApplyDBCs());
	  matrix->finishInsertion();
	}
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);
	if (matrix)
	  nnz += matrix->getBaseMatrix().nnz();
      }

      assembleBoundaryConditions(rhs->getDOFVector(i),