ProblemStat.cc 51.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

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

namespace AMDiS {

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

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

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

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

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

82
83
84
85

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

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

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

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

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

116
117
118
      }
    }

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

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

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

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

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

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

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

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

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
163
    if (initFlag.isSet(INIT_FILEWRITER))
164
165
166
167
168
      createFileWriter();
    
    // === read serialization and init mesh ===
    
    // There are two possiblities where the user can define a serialization
169
170
171
172
    // 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.
173
174

    int readSerialization = 0;
175
176
    string serializationFilename = "";
    Parameters::get("argv->rs", serializationFilename);
177
178
179
180
181
182

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

183
184
185
      Parameters::get(name + "->input->read serialization", readSerialization);
      Parameters::get(name + "->input->serialization with adaptinfo",
		      readSerializationWithAdaptInfo);
186
187
188
189
190
191

      // 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) {
192
193
	Parameters::get(name + "->input->serialization filename", 
			serializationFilename);
194
195
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

196
197
	// If AMDiS is compiled for parallel computations, the deserialization is
	// controlled by the parallel problem object.
198
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
199
	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
200
	ifstream in(serializationFilename.c_str());
201
202
	deserialize(in);
	in.close();
203
204

#endif
205
206

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

210
	// If AMDiS is compiled for parallel computations, the global refinements are
211
212
213
	// ignored here. Later, each rank will add the global refinements to its 
	// private mesh.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
214
215
	Parameters::get(meshes[0]->getName() + "->global refinements", 
			globalRefinements);
216
#endif
217

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

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

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

	// === do global refinements ===
235

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

    doOtherStuff();
  }

245

246
  void ProblemStatSeq::createMesh() 
247
  {
248
    FUNCNAME("ProblemStat::createMesh()");
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    componentMeshes.resize(nComponents);
251
252
253
254
255
256
    map<int, Mesh*> meshForRefinementSet;

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

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

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

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

297

298
  void ProblemStatSeq::createFeSpace(DOFAdmin *admin)
299
  {
300
    FUNCNAME("ProblemStat::createFeSpace()");
301

302
    map<pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
303
    int dim = -1;
304
    Parameters::get(name + "->dim", dim);
305
    TEST_EXIT(dim != -1)("no problem dimension specified!\n");
306

307
    componentSpaces.resize(nComponents, NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
308
    traverseInfo.resize(nComponents);
309

310
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
311
      int degree = 1;
312
313
      Parameters::get(name + "->polynomial degree[" + 
		      boost::lexical_cast<string>(i) + "]", degree);
314

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

317
      if (feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
320
	stringstream s;
	s << name << "->feSpace[" << i << "]";

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

328
      componentSpaces[i] = feSpaceMap[pair<Mesh*, int>(componentMeshes[i], degree)];
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
    }

    for (int i = 0; i < nComponents; i++) {
332
      for (int j = 0; j < nComponents; j++)
333
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
334
      
335
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
336
337
338
    }

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

348

349
  void ProblemStatSeq::createMatricesAndVectors()
350
  {
351
    FUNCNAME("ProblemStat::createMatricesAndVectors()");
352
353
354

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

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
Thomas Witkowski's avatar
Thomas Witkowski committed
356
    systemMatrix->set(NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
    rhs = new SystemVector("rhs", componentSpaces, nComponents);
    solution = new SystemVector("solution", componentSpaces, nComponents);
359

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

365
      string numberedName = "rhs[" + boost::lexical_cast<string>(i) + "]";
Thomas Witkowski's avatar
Thomas Witkowski committed
366
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName));
367

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

376

377
  void ProblemStatSeq::createSolver()
378
  {
379
    FUNCNAME("ProblemStat::createSolver()");
380
381

    // === create solver ===
382
383
    string solverType("0");
    Parameters::get(name + "->solver", solverType);
384
385
    OEMSolverCreator *solverCreator = 
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
386
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
389
    solverCreator->setName(name + "->solver");
    solver = solverCreator->create();
    solver->initParameters();
390
391
  }

392

393
  void ProblemStatSeq::createEstimator()
394
  {
395
    FUNCNAME("ProblemStat::createEstimator()");
396
397

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

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

      // === create estimator ===
408
409
      string estimatorType("0");
      Parameters::get(estName, estimatorType);
410
      EstimatorCreator *estimatorCreator = 
Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
411
	dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
412
      if (estimatorCreator) {
413
414
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
415
	estimatorCreator->setSolution(solution->getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
416
	estimator[i] = estimatorCreator->create();
417
418
419
      }


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

428

429
  void ProblemStatSeq::createMarker()
430
  {
431
    FUNCNAME("ProblemStat::createMarker()");
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    int nMarkersCreated = 0;
434

435
    for (int i = 0; i < nComponents; i++) {
436
      marker[i] = Marker::createMarker
437
	(name + "->marker[" + boost::lexical_cast<string>(i) + "]", i);
438

439
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
442
443
	nMarkersCreated++;

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

450

451
  void ProblemStatSeq::createFileWriter()
452
  {
453
    FUNCNAME("ProblemStat::createFileWriter()");
454
455
  
    // Create one filewriter for all components of the problem
456
457
458
    string numberedName  = name + "->output";
    string filename = "";
    Parameters::get(numberedName + "->filename", filename);
459
460

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
467
	solutionList[i] = solution->getDOFVector(i);
468
469
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
470
      fileWriters.push_back(new FileWriter(numberedName,
471
472
					   componentMeshes[0],
					   solutionList));
473
474
475
    }

    // Create own filewriters for each components of the problem
476
    for (int i = 0; i < nComponents; i++) {
477
      numberedName = name + "->output[" + boost::lexical_cast<string>(i) + "]";
478
      filename = "";
479
      Parameters::get(numberedName + "->filename", filename);
480

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

    // Check for serializer
    int writeSerialization = 0;
489
    Parameters::get(name + "->write serialization", writeSerialization);
490
    if (writeSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
      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());
493
494
495
      ERROR_EXIT("Usage of an obsolete parameter (see message above)!\n");
    }

496
    Parameters::get(name + "->output->write serialization", writeSerialization);
497

498
    if (writeSerialization)
499
      fileWriters.push_back(new Serializer<ProblemStatSeq>(this));
500
501
  }

502

503
  void ProblemStatSeq::doOtherStuff()
504
505
506
  {
  }

507

508
  void ProblemStatSeq::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
509
510
511
  {
    FUNCNAME("Problem::solve()");

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

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

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

524
    adaptInfo->setSolverIterations(solver->getIterations());
Thomas Witkowski's avatar
Thomas Witkowski committed
525
526
527
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
528
529
  }

530

531
  void ProblemStatSeq::estimate(AdaptInfo *adaptInfo) 
532
  {
533
    FUNCNAME("ProblemStat::estimate()");
534

535
536
537
538
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#else
    clock_t first = clock();
539
540
#endif

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

552
553
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
554
555
556
557
558

	  if (adaptInfo->getRosenbrockMode() == false) {
	    adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	    adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	  }
559
	}
560
561
562
      }
    }

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

573

574
  Flag ProblemStatSeq::markElements(AdaptInfo *adaptInfo) 
575
  {
576
    FUNCNAME("ProblemStat::markElements()");
577

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

581
    Flag markFlag = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
582
    for (int i = 0; i < nComponents; i++)
583
      if (marker[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
584
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
585
    
586
587
588
    return markFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
589

590
  Flag ProblemStatSeq::refineMesh(AdaptInfo *adaptInfo) 
591
  {
592
    FUNCNAME("ProblemStat::refineMesh()");
593

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

600
601
602
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
603

604
  Flag ProblemStatSeq::coarsenMesh(AdaptInfo *adaptInfo) 
605
  {
606
    FUNCNAME("ProblemStat::coarsenMesh()");
607

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

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
617

618
  Flag ProblemStatSeq::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
619
  {
620
    FUNCNAME("ProblemStat::oneIteration()");
621

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
631

632
633
  void ProblemStatSeq::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
					 bool asmMatrix, bool asmVector)
634
  {
635
    FUNCNAME("ProblemStat::buildAfterCoarsen()");
636

637
638
639
//     buildAfterCoarsen_sebastianMode(adaptInfo, flag);
//     return;

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

645
      dualAssemble(adaptInfo, flag, asmMatrix, asmVector);
646
647
      return;
    }
648

649
650
651
652
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double first = MPI::Wtime();
#endif

653

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


658
659
660
661
662
663
664
665
666
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    MPI::COMM_WORLD.Barrier();
    INFO(info, 8)("dof compression needed %.5f seconds\n", 
		  MPI::Wtime() - first);

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

668
669
670

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

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

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

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

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

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

695
      for (int j = 0; j < nComponents; j++) {
696

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

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

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

721
722
723
	if (matrix) 
	  matrix->calculateNnz();
	
724
725
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
726
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
727
	  assembleMatrix = false;
728
729
730
731
732
733
734
	} 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
735
736

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

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

747
	  continue;
748
	}
749

750
751
752
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

753
754
755
756
757
758
759
760
	// 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
761

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

765
766
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	
767

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

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

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

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

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

795
796

  void ProblemStatSeq::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
797
  {
798
    FUNCNAME("ProblemStat::buildAfterCoarsen()");
799
800
801
802
803
804
805
806
807
808

    clock_t first = clock();

    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()  |
809
810
811
812
      Mesh::CALL_LEAF_EL                       | 
      Mesh::FILL_COORDS                        |
      Mesh::FILL_DET                           |
      Mesh::FILL_GRD_LAMBDA                    |
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
<