ProblemVec.cc 48.5 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <sstream>
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"
29
30
31

namespace AMDiS {

32
33
  using boost::lexical_cast;

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

40
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
41
    if (meshes.size() != 0) { 
42
43
44
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
Thomas Witkowski's avatar
Thomas Witkowski committed
45
	  (!adoptFlag.isSet(INIT_MESH) &&
46
47
48
49
50
51
52
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
	createMesh();
      } 
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_MESH) || 
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
53
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
54
	componentMeshes = adoptProblem->componentMeshes;
Thomas Witkowski's avatar
Thomas Witkowski committed
55
56
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
57
58
59

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

68
69
70
      }
    }

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

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

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

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

98
99
100
      }
    }

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

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

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

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

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

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

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

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

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

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

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

	// 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
193
	GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d", 
194
		      &globalRefinements);
195
#endif
196

197
	// Initialize the meshes if there is no serialization file.
198
199
200
201
202
	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
203
	    meshes[i]->initialize();	    
204
	}
205
206
207
208

	// === 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
209
	if (valueFilename.length()) {
210
211
212
213
214
215
216
217
218
219
	  ValueReader::readValue(valueFilename,
				 meshes[0],
				 solution->getDOFVector(0),
				 meshes[0]->getMacroFileInfo());
	  meshes[0]->clearMacroFileInfo();
	}

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

    doOtherStuff();
  }

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

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

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

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

274
  void ProblemVec::createFESpace(DOFAdmin *admin)
275
276
277
  {
    FUNCNAME("ProblemVec::createFESpace()");

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

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

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

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

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

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

    for (int i = 0; i < nComponents; i++) {
307
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
308
	traverseInfo.getMatrix(i, j).setFESpace(componentSpaces[i], componentSpaces[j]);
309
      
Thomas Witkowski's avatar
Thomas Witkowski committed
310
      traverseInfo.getVector(i).setFESpace(componentSpaces[i]);
311
312
313
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
    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);
317
	ln_dof[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
318
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
319
320
321
322
323
324
325
326
327
328
      }
    }
  }

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

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

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

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

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

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

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

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

    // create and set leaf data prototype
369
    for (int i = 0; i < static_cast<int>(meshes.size()); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
370
      meshes[i]->setElementDataPrototype
Thomas Witkowski's avatar
Thomas Witkowski committed
371
	(new LeafDataEstimatableVec(new LeafDataCoarsenableVec));
372

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

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


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
407
    int nMarkersCreated = 0;
408

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

413
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
417
	nMarkersCreated++;

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
440
	solutionList[i] = solution->getDOFVector(i);
441
442
      }

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

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

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

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

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

    // 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
474
    if (writeSerialization)
Thomas Witkowski's avatar
Thomas Witkowski committed
475
      fileWriters.push_back(new Serializer<ProblemVec>(this));
476
#endif
477
478
479
480
481
482
  }

  void ProblemVec::doOtherStuff()
  {
  }

483
  void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
484
485
486
  {
    FUNCNAME("Problem::solve()");

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

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

    clock_t first = clock();
497
    solver->solveSystem(solverMatrix, *solution, *rhs);
498

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

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

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

    clock_t first = clock();

519
520
521
522
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

523
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
524
      computeError(adaptInfo);
525
526
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
527
	Estimator *scalEstimator = estimator[i];
528
529
530
531
532
533
534
535
536
537
	
	if (scalEstimator) {
	  scalEstimator->estimate(adaptInfo->getTimestep());
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
	  adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	  adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	} else {
	  WARNING("no estimator for component %d\n" , i);
	}
538
539
540
      }
    }

541
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
542
543
    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);
544
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
545
546
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
547
548
549

#endif

550
551
552
553
554
555
556
557
558
559
560
  }

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
576
    int nMeshes = static_cast<int>(meshes.size());
577
    Flag refineFlag = 0;
578
579
580
581
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

582
583
584
585
586
587
588
    return refineFlag;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
589
    int nMeshes = static_cast<int>(meshes.size());
590
    Flag coarsenFlag = 0;
591
592
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isCoarseningAllowed(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
593
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
594
595
596
597
598
599
600
601

    return coarsenFlag;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
602
    if (allowFirstRef) {
Thomas Witkowski's avatar
Thomas Witkowski committed
603
      for (int i = 0; i < nComponents; i++)
604
	adaptInfo->allowRefinement(true, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
605

Thomas Witkowski's avatar
Thomas Witkowski committed
606
      allowFirstRef = false;
607
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
608
609
      for (int i = 0; i < nComponents; i++)
	if (adaptInfo->spaceToleranceReached(i))
610
	  adaptInfo->allowRefinement(false, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
611
	else
612
613
614
615
616
617
	  adaptInfo->allowRefinement(true, i);	
    }

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

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

623
    //    buildAfterCoarsen_sebastianMode(adaptInfo, flag);
624

625
    clock_t first = clock();
626
627
628
629
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

630
    for (int i = 0; i < static_cast<int>(meshes.size()); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
631
      meshes[i]->dofCompress();
632
633
634

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
637
638
639
640
641
642
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
643
    if (useGetBound)
644
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
645

Thomas Witkowski's avatar
Thomas Witkowski committed
646
647
    traverseInfo.updateStatus();

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

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

653
654
655
656
657
658
      MSG("%d DOFs for %s\n", 
	  componentSpaces[i]->getAdmin()->getUsedSize(), 
	  componentSpaces[i]->getName().c_str());

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

659
      for (int j = 0; j < nComponents; j++) {
660

661
662
663
664
	// 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).
Thomas Witkowski's avatar
Thomas Witkowski committed
665
	DOFMatrix *matrix = (*systemMatrix)[i][j];
666

667
668
669
	if (matrix) 
	  matrix->calculateNnz();
	
670
671
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
672
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
673
	  assembleMatrix = false;
674
675
676
677
678
679
680
	} 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
681
682

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

686
687
688
	// 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.
689
690
691
692
	if (!assembleMatrix && i != j) {
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();

693
	  continue;
694
	}
695

696
697
698
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
	if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_NO_AUX ||
	    traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_AUX) {

	  // Row fe space and col fe space are both equal

	  if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_NO_AUX ||
	      traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_AUX) {

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

	  } else if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {

	    // Row fe space and col fe space are both equal, but right hand side has at
718
	    // least one another aux fe space. 
Thomas Witkowski's avatar
Thomas Witkowski committed
719
720
721
722
723
724

	    assembleOnOneMesh(componentSpaces[i],
			      assembleFlag,
			      assembleMatrix ? matrix : NULL,
			      ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);

725
726
	    assembleOnDifMeshes2(componentSpaces[i], 
				 traverseInfo.getAuxFESpace(i, j),
Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
				 assembleFlag,
				 NULL,
729
				 ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
730
731
732
733
734
735
736

	  } else {
	    ERROR_EXIT("Possible? If yes, not yet implemented!\n");
	  }

	} else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
739
	  assembleOnOneMesh(componentSpaces[i],
			    assembleFlag,
			    assembleMatrix ? matrix : NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
740
			    ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
741
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
742
743
744
745
	  assembleOnDifMeshes2(componentSpaces[i],
			       traverseInfo.getAuxFESpace(i, j),
			       assembleFlag,
			       assembleMatrix ? matrix : NULL,
746
			       ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
750

	} else if (traverseInfo.getStatus(i, j) ==  SingleComponentInfo::DIF_SPACES_NO_AUX ||
		   traverseInfo.getStatus(i, j) ==  SingleComponentInfo::DIF_SPACES_WITH_AUX) {

Thomas Witkowski's avatar
Thomas Witkowski committed
751
752
753
	  assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
			      assembleFlag,
			      assembleMatrix ? matrix : NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
756
			      ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);	  
	} else {
	  ERROR_EXIT("Not yet implemented!\n");
757
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
758

Thomas Witkowski's avatar
Thomas Witkowski committed
759
	assembledMatrix[i][j] = true;
760
761
762

	if (assembleMatrix)
	  matrix->finishInsertion();
763
764

 	if (assembleMatrix && matrix->getBoundaryManager())
765
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	
766
767
	
	if (matrix)
768
	  nnz += matrix->getBaseMatrix().nnz();	  
769
770
      }

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

778
    solverMatrix.setMatrix(*systemMatrix);
779

Thomas Witkowski's avatar
Thomas Witkowski committed
780
    createPrecon();
781

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

784
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
787
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
790
#endif    
791
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
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
  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) == 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) == 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),
				 solution->getDOFVector(i),
				 componentMeshes[i],
				 assembleFlag);     
    }

    solverMatrix.setMatrix(*systemMatrix);

    createPrecon();

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

#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     

    exit(0);

  }

988
989
  void ProblemVec::createPrecon()
  {
990
    std::string preconType("no");
991
992
993
994
995
    GET_PARAMETER(0, name + "->solver->left precon", &preconType);

    CreatorInterface<ITL_BasePreconditioner> *preconCreator = 
      CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);

Thomas Witkowski's avatar
Thomas Witkowski committed
996
    solver->setLeftPrecon(preconCreator->create(solverMatrix.getMatrix()));
997
998
999
1000
1001

    preconType= "no";
    GET_PARAMETER(0, name + "->solver->right precon", &preconType);

    preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
Thomas Witkowski's avatar
Thomas Witkowski committed
1002
    solver->setRightPrecon(preconCreator->create(solverMatrix.getMatrix()));
1003
1004
1005
  }