ProblemVec.cc 48.8 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <sstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
2
#include <boost/lexical_cast.hpp>
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "ProblemVec.h"
#include "RecoveryEstimator.h"
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "FileWriter.h"
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
17
#include "DualTraverse.h"
18
19
20
21
22
23
#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
24
#include "Flag.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "TraverseParallel.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include "VtkWriter.h"
27
#include "ValueReader.h"
28
#include "ProblemVecDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
29
#include "Debug.h"
30
31
32

namespace AMDiS {

33
34
  using boost::lexical_cast;

35
36
37
38
39
  void ProblemVec::initialize(Flag initFlag,
			      ProblemVec *adoptProblem,
			      Flag adoptFlag)
  {
    FUNCNAME("ProblemVec::initialize()");
40

41
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
42
    if (meshes.size() != 0) { 
43
44
45
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
Thomas Witkowski's avatar
Thomas Witkowski committed
46
	  (!adoptFlag.isSet(INIT_MESH) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
47
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
48
	createMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
49

50
      if (adoptProblem && 
Thomas Witkowski's avatar
Thomas Witkowski committed
51
	  (adoptFlag.isSet(INIT_MESH) ||
52
53
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
54
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
55
	componentMeshes = adoptProblem->componentMeshes;
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
58
59
60

	// 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
61
	if (adoptProblem->getNumComponents() < nComponents) {
Thomas Witkowski's avatar
Thomas Witkowski committed
62
	  TEST_EXIT(meshes.size() == 1)("Daran muss ich noch arbeiten!\n");
63
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
64
	  componentMeshes.resize(nComponents);
65
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
66
	    componentMeshes[i] = componentMeshes[0];
67
68
	}

69
70
71
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    if (meshes.size() == 0) 
73
74
75
      WARNING("no mesh created\n");

    // === create fespace ===
76
    if (feSpaces.size() != 0) {
77
78
79
      WARNING("feSpaces already created\n");
    } else {
      if (initFlag.isSet(INIT_FE_SPACE) || 
80
81
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
	createFESpace(NULL);
82

83
84
      if (adoptProblem &&
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
85
86
	feSpaces = adoptProblem->getFESpaces();
	componentSpaces = adoptProblem->componentSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
87
	traverseInfo = adoptProblem->traverseInfo;
88
89
90

	// 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
91
92
	if (adoptProblem->getNumComponents() < nComponents) {
	  TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n");
93
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
94
	  componentSpaces.resize(nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
95
	  for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
96
	    componentSpaces[i] = componentSpaces[0];
97
98
	}

99
100
101
      }
    }

102
    if (feSpaces.size() == 0) 
103
104
105
      WARNING("no feSpace created\n");

    // === create system ===
106
    if (initFlag.isSet(INIT_SYSTEM)) 
107
      createMatricesAndVectors();
108
    
109
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
112
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRHS();
      systemMatrix = adoptProblem->getSystemMatrix();
113
114
115
    }

    // === create solver ===
Thomas Witkowski's avatar
Thomas Witkowski committed
116
    if (solver) {
117
118
119
120
121
122
      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
123
124
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
125
126
127
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
128
    if (!solver) 
129
130
131
      WARNING("no solver created\n");

    // === create estimator ===
Thomas Witkowski's avatar
Thomas Witkowski committed
132
    if (initFlag.isSet(INIT_ESTIMATOR))
133
      createEstimator();
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135

    if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
Thomas Witkowski's avatar
Thomas Witkowski committed
136
      estimator = adoptProblem->getEstimator();
137
138

    // === create marker ===
Thomas Witkowski's avatar
Thomas Witkowski committed
139
    if (initFlag.isSet(INIT_MARKER))
140
      createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
141
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
142
      marker = adoptProblem->getMarker();
143
144

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
145
    if (initFlag.isSet(INIT_FILEWRITER))
146
147
148
149
150
151
152
153
154
155
156
      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;
157
    std::string serializationFilename = "";
158
159
160
161
162
163
164
    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;

165
      GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
Thomas Witkowski's avatar
Thomas Witkowski committed
166
      GET_PARAMETER(0, name + "->input->serialization with adaptinfo", "%d",
167
168
169
170
171
172
173
		    &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
174
	GET_PARAMETER(0, name + "->input->serialization filename", 
175
176
177
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

178
179
	// If AMDiS is compiled for parallel computations, the deserialization is
	// controlled by the parallel problem object.
180
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
181
	std::ifstream in(serializationFilename.c_str());
182
183
	deserialize(in);
	in.close();
184
185
186

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
#endif
187
      } else {
188
	int globalRefinements = 0;
189
190
191
192
193

	// If AMDiS is compied for parallel computations, the global refinements are
	// 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
194
	GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d", 
195
		      &globalRefinements);
196
#endif
197

198
	// Initialize the meshes if there is no serialization file.
199
200
201
202
203
	for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
	  bool initMesh = initFlag.isSet(INIT_MESH) ||
	    (adoptProblem && adoptFlag.isSet(INIT_MESH));

	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
Thomas Witkowski's avatar
Thomas Witkowski committed
204
	    meshes[i]->initialize();	    
205
	}
206
207
208
209

	// === read value file and use it for the mesh values ===
	std::string valueFilename("");
	GET_PARAMETER(0, meshes[0]->getName() + "->value file name", &valueFilename); 
Thomas Witkowski's avatar
nix.  
Thomas Witkowski committed
210
	if (valueFilename.length()) {
211
212
213
214
215
216
217
218
	  ValueReader::readValue(valueFilename,
				 meshes[0],
				 solution->getDOFVector(0),
				 meshes[0]->getMacroFileInfo());
	  meshes[0]->clearMacroFileInfo();
	}

	// === do global refinements ===
219
220

	for (unsigned int i = 0; i < meshes.size(); i++)
221
	  if (initFlag.isSet(INIT_MESH) && meshes[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
222
	    refinementManager->globalRefine(meshes[i], globalRefinements);
223
224
225
226
227
228
229
230
231
232
      }
    }

    doOtherStuff();
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
233
    componentMeshes.resize(nComponents);
234
    std::map<int, Mesh*> meshForRefinementSet;
235

236
    std::string meshName("");
Thomas Witkowski's avatar
Thomas Witkowski committed
237
    GET_PARAMETER(0, name + "->mesh", &meshName);
238
    TEST_EXIT(meshName != "")("no mesh name specified\n");
239
    int dim = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
240
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
241
    TEST_EXIT(dim)("no problem dimension specified!\n");
242

243
    for (int i = 0; i < nComponents; i++) {
244
      int refSet = -1;
245
246
      GET_PARAMETER(0, name + "->refinement set[" + 
		       lexical_cast<std::string>(i) + "]", "%d", &refSet);
247
      if (refSet < 0)
248
	refSet = 0;
249
      
250
      if (meshForRefinementSet[refSet] == NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
251
	Mesh *newMesh = new Mesh(meshName, dim);
252
	meshForRefinementSet[refSet] = newMesh;
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254
	meshes.push_back(newMesh);
	nMeshes++;
255
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
256
      componentMeshes[i] = meshForRefinementSet[refSet];
257
258
259
    }
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
262
263
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
266
267
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
270
271
272
273
274
275
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

276
  void ProblemVec::createFESpace(DOFAdmin *admin)
277
278
279
  {
    FUNCNAME("ProblemVec::createFESpace()");

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

285
    componentSpaces.resize(nComponents, NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
286
    traverseInfo.resize(nComponents);
287

288
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
289
      int degree = 1;
290
291
      GET_PARAMETER(0, name + "->polynomial degree[" + 
		    boost::lexical_cast<std::string>(i) + "]","%d", &degree);
292

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

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

299
	FiniteElemSpace *newFESpace = 
300
	  FiniteElemSpace::provideFESpace(admin, Lagrange::getLagrange(dim, degree),
Thomas Witkowski's avatar
Thomas Witkowski committed
301
					  componentMeshes[i], s.str());
Thomas Witkowski's avatar
Thomas Witkowski committed
302
	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFESpace;
303
	feSpaces.push_back(newFESpace);
304
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
      componentSpaces[i] = feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
    }

    for (int i = 0; i < nComponents; i++) {
309
      for (int j = 0; j < nComponents; j++)
310
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
311
      
312
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
313
314
315
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
318
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
      if (meshes[i]->getNumberOfDOFs(VERTEX) == 0) {
	DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
319
	ln_dof[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
320
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
321
322
323
324
325
326
327
328
329
330
      }
    }
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
331
    systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    systemMatrix->set(NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
    rhs = new SystemVector("rhs", componentSpaces, nComponents);
    solution = new SystemVector("solution", componentSpaces, nComponents);
335

Thomas Witkowski's avatar
Thomas Witkowski committed
336
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
337
      (*systemMatrix)[i][i] = new DOFMatrix(componentSpaces[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
338
					    componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
339
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
Thomas Witkowski's avatar
Thomas Witkowski committed
340

341
      std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]";
Thomas Witkowski's avatar
Thomas Witkowski committed
342
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName));
343
344

      numberedName = "solution[" + boost::lexical_cast<std::string>(i) + "]";
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
      solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
      						      numberedName));
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
349
350
351
352
353
354
355
356
    }
  }

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

    // === create solver ===
357
    std::string solverType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    GET_PARAMETER(0, name + "->solver", &solverType);
359
360
    OEMSolverCreator *solverCreator = 
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
361
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
364
    solverCreator->setName(name + "->solver");
    solver = solverCreator->create();
    solver->initParameters();
365
366
  }

367

368
369
370
371
372
  void ProblemVec::createEstimator()
  {
    FUNCNAME("ProblemVec::createEstimator()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
    for (int i = 0; i < nComponents; i++) {
      TEST_EXIT(estimator[i] == NULL)("estimator already created\n");
379
380
      std::string estName = 
	name + "->estimator[" + boost::lexical_cast<std::string>(i) + "]";
381
382

      // === create estimator ===
383
      std::string estimatorType("0");
384
385
      GET_PARAMETER(0, estName, &estimatorType);
      EstimatorCreator *estimatorCreator = 
Thomas Witkowski's avatar
nix.  
Thomas Witkowski committed
386
	dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
387
      if (estimatorCreator) {
388
389
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
390
	if (estimatorType == "recovery")
391
	  dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->
392
	    setSolution(solution->getDOFVector(i));       
Thomas Witkowski's avatar
Thomas Witkowski committed
393
	estimator[i] = estimatorCreator->create();
394
395
396
      }


Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
399
400
401
      if (estimator[i])
	for (int j = 0; j < nComponents; j++)
	  estimator[i]->addSystem((*systemMatrix)[i][j], 
				  solution->getDOFVector(j), 
				  rhs->getDOFVector(j));      
402
403
404
    }
  }

405

406
407
408
409
  void ProblemVec::createMarker()
  {
    FUNCNAME("ProblemVec::createMarker()");

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    int nMarkersCreated = 0;
411

412
    for (int i = 0; i < nComponents; i++) {
413
414
415
      marker[i] = Marker::createMarker
	(name + "->marker[" + boost::lexical_cast<std::string>(i) + "]", i);

416
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
420
	nMarkersCreated++;

	// If there is more than one marker, and all components are defined
	// on the same mesh, the maximum marking has to be enabled.
421
 	if (nMarkersCreated > 1 && nMeshes == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
422
 	  marker[i]->setMaximumMarking(true);
423
424
425
426
427
428
429
430
431
      }
    }
  }

  void ProblemVec::createFileWriter()
  {
    FUNCNAME("ProblemVec::createFileWriter()");
  
    // Create one filewriter for all components of the problem
Thomas Witkowski's avatar
Thomas Witkowski committed
432
    std::string numberedName  = name + "->output";
433
    std::string filename = "";
434
435
436
    GET_PARAMETER(0, numberedName + "->filename", &filename);

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
443
	solutionList[i] = solution->getDOFVector(i);
444
445
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
446
      fileWriters.push_back(new FileWriter(numberedName,
Thomas Witkowski's avatar
Thomas Witkowski committed
447
					    componentMeshes[0],
448
449
450
451
					    solutionList));
    }

    // Create own filewriters for each components of the problem
452
    for (int i = 0; i < nComponents; i++) {
453
      numberedName = name + "->output[" + boost::lexical_cast<std::string>(i) + "]";
454
455
456
      filename = "";
      GET_PARAMETER(0, numberedName + "->filename", &filename);

457
      if (filename != "")
Thomas Witkowski's avatar
Thomas Witkowski committed
458
	fileWriters.push_back(new FileWriter(numberedName, 
459
460
					     componentMeshes[i], 
					     solution->getDOFVector(i)));
461
462
463
464
    }

    // Check for serializer
    int writeSerialization = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
465
    GET_PARAMETER(0, name + "->write serialization", "%d", &writeSerialization);
466
    if (writeSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
      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());
469
470
471
      ERROR_EXIT("Usage of an obsolete parameter (see message above)!\n");
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
472
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
473
474
475
476

    // Serialization is not allowed to be done by the problem, if its part of a parallel
    // problem definition. Than, the parallel problem object must be serialized.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
477
    if (writeSerialization)
Thomas Witkowski's avatar
Thomas Witkowski committed
478
      fileWriters.push_back(new Serializer<ProblemVec>(this));
479
#endif
480
481
482
483
484
485
  }

  void ProblemVec::doOtherStuff()
  {
  }

486
  void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
487
488
489
  {
    FUNCNAME("Problem::solve()");

Thomas Witkowski's avatar
Thomas Witkowski committed
490
    if (!solver) {
491
492
493
494
495
496
497
498
499
      WARNING("no solver\n");
      return;
    }

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

    clock_t first = clock();
500
    
501
    solver->solveSystem(solverMatrix, *solution, *rhs);
502

503
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    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);
506
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
509
510
#endif

511
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
515
516
  }

517

518
519
520
521
522
523
  void ProblemVec::estimate(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::estimate()");

    clock_t first = clock();

524
525
526
527
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

528
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
529
      computeError(adaptInfo);
530
531
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
532
	Estimator *scalEstimator = estimator[i];
533
534
	
	if (scalEstimator) {
535
536
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
537
	  scalEstimator->estimate(adaptInfo->getTimestep());
538

539
540
541
542
543
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
	  adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	  adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	} else {
544
	  WARNING("No estimator for component %d\n" , i);
545
	}
546
547
548
      }
    }

549
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
550
551
    INFO(info, 8)("estimation of the error needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
552
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
553
554
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
555
#endif
556
557
  }

558

559
560
561
562
  Flag ProblemVec::markElements(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::markElements()");

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

566
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
567
    for (int i = 0; i < nComponents; i++) {
568
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
569
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
570
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
571
572
	WARNING("No marker for component %d\n", i);              
    }
573
    
574
575
576
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
577

578
579
580
581
  Flag ProblemVec::refineMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::refineMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
582
    int nMeshes = static_cast<int>(meshes.size());
583
    Flag refineFlag = 0;
584
585
586
587
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

588
589
590
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
591

592
593
594
595
  Flag ProblemVec::coarsenMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::coarsenMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
596
    int nMeshes = static_cast<int>(meshes.size());
597
    Flag coarsenFlag = 0;
598
599
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
600
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
601
602
603
604

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
605

606
607
608
609
  Flag ProblemVec::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ProblemVec::oneIteration()");

610
611
612
613
614
    for (int i = 0; i < nComponents; i++)
      if (adaptInfo->spaceToleranceReached(i))
	adaptInfo->allowRefinement(false, i);
      else
	adaptInfo->allowRefinement(true, i);	
615
616
617
618

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

Thomas Witkowski's avatar
Thomas Witkowski committed
619

620
621
  void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				     bool asmMatrix, bool asmVector)
622
623
624
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

625
626
627
628
    if (dualMeshTraverseRequired()) {
      dualAssemble(adaptInfo, flag);
      return;
    }
629

630
    //    printOpenmpTraverseInfo(this, true);
631
632
    //    std::cout << "ElInfo = " << ElInfo::subElemMatrices.size() << std::endl;

Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
635
636
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


637
    clock_t first = clock();
638
639
640
641
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

642
643
644

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
647
648
649
650
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
651
652
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
653
    if (useGetBound)
654
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656

  
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
    traverseInfo.updateStatus();

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

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

664
665
666
667
668
669
      MSG("%d DOFs for %s\n", 
	  componentSpaces[i]->getAdmin()->getUsedSize(), 
	  componentSpaces[i]->getName().c_str());

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

670
      for (int j = 0; j < nComponents; j++) {
671

672
673
	if (writeAsmInfo)
	  std::cout << "-------" << i << " " << j << "----------------" << std::endl;
674

675
676
677
678
	// 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).
679
	DOFMatrix *matrix = (asmMatrix ? (*systemMatrix)[i][j] : NULL);
680

681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
	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;
	  }
	}

696
697
698
	if (matrix) 
	  matrix->calculateNnz();
	
699
700
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
701
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
702
	  assembleMatrix = false;
703
704
705
706
707
708
709
	} 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
710
711

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

715
716
717
	// 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.
718
719
720
721
	if (!assembleMatrix && i != j) {
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();

722
	  continue;
723
	}
724

725
726
727
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

728
729
730
731
732
733
734
735
	// 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
736

Thomas Witkowski's avatar
Thomas Witkowski committed
737
	
738
739
	if (assembleMatrix)
	  matrix->finishInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
740

741
742
743
744
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	
	
	if (matrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
745
	  nnz += matrix->getBaseMatrix().nnz();
746
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
747

748
749
      // And now assemble boundary conditions on the vectors
      assembleBoundaryConditions(rhs->getDOFVector(i),
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
752
 				 solution->getDOFVector(i),
 				 componentMeshes[i],
 				 assembleFlag);     
753
    }
754

755
756
757
    if (asmMatrix) {
      solverMatrix.setMatrix(*systemMatrix);      
      createPrecon();
Thomas Witkowski's avatar
Thomas Witkowski committed
758

759
760
      INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
    }
761

762
763
764
765
766
767
768
769
#ifdef _OPENMP
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
#else
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));    
#endif    
  }
770

771
772
773
774
775
776
777
778
779
780
781
782
783
784
785

  bool ProblemVec::dualMeshTraverseRequired()
  {
    FUNCNAME("ProblemVec::dualMeshTraverseRequired()");

    TEST_EXIT(meshes.size() <= 2)("More than two mneshes are not supported yet!\n");    

    return (meshes.size() == 2);
  }
  

  void ProblemVec::dualAssemble(AdaptInfo *adaptInfo, Flag flag)
  {
    FUNCNAME("ProblemVec::dualAssemble()");

Thomas Witkowski's avatar
Thomas Witkowski committed
786
787
788
789
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


790
791
792
793
794
    clock_t first = clock();
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
795
    
796
797
798
799
    Flag assembleFlag = 
      flag | 
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
800
801
802
803
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
804
805
806
807
      Mesh::FILL_NEIGH;

    if (useGetBound)
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
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
    traverseInfo.updateStatus();

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

    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++) {

	if (writeAsmInfo)
	  std::cout << "-------" << i << " " << j << "----------------" << std::endl;

	// 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 (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;
Thomas Witkowski's avatar
Thomas Witkowski committed
845
	  }
846
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
847

848
849
850
851
852
853
854
855
856
857
858
	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());
859

860
861
	  set_to_zero(matrix->getBaseMatrix());
	  matrix->startInsertion(matrix->getNnz());
862

863
864
865
	  if (matrix->getBoundaryManager())
	    matrix->getBoundaryManager()->initMatrix(matrix);
	}
866

867
868
869
	// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
	if (!matrix || !assembleMatrix)
	  assembleMatrix = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
870

871
872
873
874
875
876
877
878
	// 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();
      }
    }
879

880
    TEST_EXIT(meshes.size() == 2)("There is something wrong!\n");
881

882
883
884
885
886
887
888
889
890
891
892
893
    const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
    BoundaryType *bound = 
      useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;

    DualTraverse dualTraverse;
    DualElInfo dualElInfo;
    int oldElIndex0 = -1;
    int oldElIndex1 = -1;
    dualTraverse.setFillSubElemMat(true, basisFcts);
    bool cont = dualTraverse.traverseFirst(meshes[0], meshes[1], -1, -1, 
					   assembleFlag, assembleFlag, dualElInfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
894
    while (cont) {     
895
896
897
      bool newEl0 = (dualElInfo.rowElInfo->getElement()->getIndex() != oldElIndex0);
      bool newEl1 = (dualElInfo.colElInfo->getElement()->getIndex() != oldElIndex1);
      oldElIndex0 = dualElInfo.rowElInfo->getElement()->getIndex();