Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

ProblemScal.cc 19.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
#include "ProblemScal.h"
#include "AbstractFunction.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "FixVec.h"
#include "Flag.h"
#include "Serializer.h"
#include "RecoveryEstimator.h"
#include "Operator.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "OEMSolver.h"
#include "DOFVector.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "ElInfo.h"
#include "FileWriter.h"
#include "RefinementManager.h"
#include "CoarseningManager.h"
#include "Lagrange.h"
#include "PeriodicBC.h"
#include "ValueReader.h"
24
#include "ElementFileWriter.h"
25
#include "ProblemVec.h"
26 27 28

namespace AMDiS {

29 30 31 32
  ProblemScal::~ProblemScal()
  {
    FUNCNAME("ProblemScal::~ProblemScal()");

33 34
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++)
      delete fileWriters[i];
35 36

    if (systemMatrix)
37
      delete systemMatrix;
38
    if (rhs) 
39
      delete rhs;
40
    if (solution)
41
      delete solution;
42
    if (estimator)
43
      delete estimator;
44
    if (marker)
45
      delete marker;
46
    if (solver)
47
      delete solver;
48
    if (mesh)
49
      delete mesh;
50 51 52 53 54

    FiniteElemSpace::clear();
    Lagrange::clear();
  }

55 56 57
  void ProblemScal::writeFiles(AdaptInfo *adaptInfo, bool force) 
  {
    for (int i = 0; i < static_cast<int>(fileWriters.size()); i++)
58
      fileWriters[i]->writeFiles(adaptInfo, force);
Thomas Witkowski's avatar
Thomas Witkowski committed
59 60
  }

61 62
  void ProblemScal::interpolInitialSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
63
    solution->interpol(fct);
64 65 66 67 68 69
  }

  void ProblemScal::addMatrixOperator(Operator *op, 
				      double *factor,
				      double *estFactor) 
  {
70
    systemMatrix->addOperator(op, factor, estFactor);
71 72 73 74 75 76
  }

  void ProblemScal::addVectorOperator(Operator *op, 
				      double *factor,
				      double *estFactor) 
  {
77
    rhs->addOperator(op, factor, estFactor);
78 79 80 81 82
  }
  
  void ProblemScal::addDirichletBC(BoundaryType type, 
				   AbstractFunction<double, WorldVector<double> >* b) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    DirichletBC *dirichlet = new DirichletBC(type, b, feSpace);
84 85 86 87 88 89
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
90 91 92 93 94 95
  }

  void ProblemScal::addDirichletBC(BoundaryType type, 
				   DOFVector<double> *vec) 
  {
    DirichletBC *dirichlet = new DirichletBC(type, vec);
96 97 98 99 100 101
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getBoundaryManager()->addBoundaryCondition(dirichlet);
102 103 104 105 106 107
  }

  void ProblemScal::addRobinBC(BoundaryType type, 
			       AbstractFunction<double, WorldVector<double> > *n,
			       AbstractFunction<double, WorldVector<double> > *r)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    RobinBC *robin = new RobinBC(type, n, r, feSpace);
109 110 111 112
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
113 114 115 116 117
  }

  void ProblemScal::addNeumannBC(BoundaryType type, 
				 AbstractFunction<double, WorldVector<double> > *n)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
118
    NeumannBC *neumann = new NeumannBC(type, n, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
119
    if (rhs) 
120
      rhs->getBoundaryManager()->addBoundaryCondition(neumann);
121 122 123 124 125 126
  }

  void ProblemScal::addRobinBC(BoundaryType type, 
			       DOFVector<double> *n,
			       DOFVector<double> *r)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
127
    RobinBC *robin = new RobinBC(type, n, r, feSpace);
128 129 130 131
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(robin);
    if (systemMatrix)
      systemMatrix->getBoundaryManager()->addBoundaryCondition(robin);
132 133 134 135
  }

  void ProblemScal::addPeriodicBC(BoundaryType type) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
136
    PeriodicBC *periodic = new PeriodicBC(type, feSpace);
137

138 139 140 141
    if (systemMatrix) 
      systemMatrix->getBoundaryManager()->addBoundaryCondition(periodic);
    if (rhs) 
      rhs->getBoundaryManager()->addBoundaryCondition(periodic);
142 143 144 145 146 147 148
  }

  void ProblemScal::createMesh()
  {
    TEST_EXIT(Parameters::initialized())("parameters not initialized\n");

    // === create problems mesh ===
Thomas Witkowski's avatar
Thomas Witkowski committed
149
    std::string meshName("");
150

Thomas Witkowski's avatar
Thomas Witkowski committed
151
    GET_PARAMETER(0, name + "->info", "%d", &info);
Thomas Witkowski's avatar
Thomas Witkowski committed
152
    GET_PARAMETER(0, name + "->mesh", &meshName);
153

154
    TEST_EXIT(meshName != "")("no mesh name specified\n");
155 156 157

    // get problem dimension
    int dim = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
158
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
159
    TEST_EXIT(dim)("no problem dimension specified!\n");
160 161

    // create the mesh
162
    mesh = new Mesh(meshName, dim);
163

Thomas Witkowski's avatar
Thomas Witkowski committed
164
    switch (dim) {
165
    case 1:
166 167
      coarseningManager = new CoarseningManager1d();
      refinementManager = new RefinementManager1d();
168 169
      break;
    case 2:
170 171
      coarseningManager = new CoarseningManager2d();
      refinementManager = new RefinementManager2d();
172 173
      break;
    case 3:
174 175
      coarseningManager = new CoarseningManager3d();
      refinementManager = new RefinementManager3d();
176 177 178 179 180 181
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }

182 183
  void ProblemScal::setMeshFromProblemVec(ProblemVec* pv, int i) 
  { 
184
    mesh = pv->getMesh(i);
185 186
    coarseningManager = pv->getCoarseningManager(i);
    refinementManager = pv->getRefinementManager(i);
187 188
  }

189 190
  Flag ProblemScal::markElements(AdaptInfo *adaptInfo) 
  { 
191 192
    if (marker)
      return marker->markMesh(adaptInfo, mesh);
193 194
    else
      WARNING("no marker\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
195

196 197 198 199 200
    return 0;
  }

  Flag ProblemScal::refineMesh(AdaptInfo *adaptInfo) 
  { 
201
    return refinementManager->refineMesh(mesh); 
202 203 204 205
  }

  Flag ProblemScal::coarsenMesh(AdaptInfo *adaptInfo) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
206
    if (adaptInfo->isCoarseningAllowed(0)) {
207
      return coarseningManager->coarsenMesh(mesh);
Thomas Witkowski's avatar
Thomas Witkowski committed
208
    } else {
209 210 211 212 213
      WARNING("coarsening not allowed\n");
      return 0;
    }
  }

214
  void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
215 216
  {
    FUNCNAME("Problem::solve()");
217

218
    if (!solver) {
219 220 221 222 223 224 225 226
      WARNING("no solver\n");
      return;
    }

#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif
    clock_t first = clock();
227 228 229 230

    SolverMatrix<DOFMatrix> solverMatrix;
    solverMatrix.setMatrix(*systemMatrix);
    int iter = solver->solveSystem(solverMatrix, *solution, *rhs);
231 232

#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
233
    INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
234 235 236
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
237
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
238 239 240 241
		   TIME_USED(first, clock()));
#endif

    adaptInfo->setSolverIterations(iter);
242 243 244
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
245 246 247 248 249 250 251 252 253
  }

  void ProblemScal::initialize(Flag initFlag,
			       ProblemScal *adoptProblem,
			       Flag adoptFlag) 
  {
    FUNCNAME("Problem::initialize()");

    // === create mesh ===
254
    if (mesh) { 
255 256
      WARNING("mesh already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
257 258 259
      if (initFlag.isSet(CREATE_MESH) || 
	  ((!adoptFlag.isSet(INIT_MESH))&&
	   (initFlag.isSet(INIT_SYSTEM)||initFlag.isSet(INIT_FE_SPACE)))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
260 261 262
	createMesh();
      } 
      
Thomas Witkowski's avatar
Thomas Witkowski committed
263 264 265 266
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_MESH) || 
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
267 268
	TEST_EXIT(!mesh)("mesh already created\n");
	mesh = adoptProblem->getMesh();
269 270
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
Thomas Witkowski's avatar
Thomas Witkowski committed
271
      }
272 273
    }

274
    if (!mesh) 
275
      WARNING("no mesh created\n");
276 277

    // === create fespace ===
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    if (feSpace) {
279 280
      WARNING("feSpace already created\n");
    } else {
281 282
      if (initFlag.isSet(INIT_FE_SPACE) || 
	  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
283 284
	createFESpace();
      } 
285

286 287
      if (adoptProblem && 
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
288 289
	TEST_EXIT(!feSpace)("feSpace already created");
	feSpace = dynamic_cast<ProblemScal*>(adoptProblem)->getFESpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
290
      }
291 292
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
293
    if (!feSpace) 
294
      WARNING("no feSpace created\n");
295 296 297 298 299 300

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM)) {
      createMatricesAndVectors();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
301 302 303 304 305 306
      TEST_EXIT(!solution)("solution already created\n");
      TEST_EXIT(!rhs)("rhs already created\n");
      TEST_EXIT(!systemMatrix)("systemMatrix already created\n");
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRHS();
      systemMatrix  = adoptProblem->getSystemMatrix();
307 308 309
    }

    // === create solver ===
310
    if (solver) {
311 312 313 314 315 316
      WARNING("solver already created\n");
    } else {
      if (initFlag.isSet(INIT_SOLVER)) {
	createSolver();
      } 
      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
317 318
	TEST_EXIT(!solver)("solver already created\n");
	solver = adoptProblem->getSolver();
319 320 321
      }
    }

322
    if (!solver) 
323 324 325
      WARNING("no solver created\n");

    // === create estimator ===
326
    if (estimator) {
327 328
      WARNING("estimator already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
329
      if (initFlag.isSet(INIT_ESTIMATOR)) {
330 331
	createEstimator();
      } 
Thomas Witkowski's avatar
Thomas Witkowski committed
332
      if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) {
333 334
	TEST_EXIT(!estimator)("estimator already created\n");
	estimator = adoptProblem->getEstimator();
335 336 337
      } 
    }

338
    if (!estimator) 
339
      WARNING("no estimator created\n");
340 341

    // === create marker ===
342
    if (marker) {
343 344
      WARNING("marker already created\n");
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
345
      if (initFlag.isSet(INIT_MARKER))
346
	createMarker();
Thomas Witkowski's avatar
Thomas Witkowski committed
347 348

      if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) {
349 350
	TEST_EXIT(!marker)("marker already created\n");
	marker = adoptProblem->getMarker();
351 352 353
      } 
    }

354
    if (!marker) 
355
      WARNING("no marker created\n");
356 357

    // === create file writer ===
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    if (initFlag.isSet(INIT_FILEWRITER))
359
      createFileWriter();
Thomas Witkowski's avatar
Thomas Witkowski committed
360
 
361 362 363 364 365 366 367 368 369 370

    // === 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
371
    std::string serializationFilename = "";
372 373 374 375 376
    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("")) {
Thomas Witkowski's avatar
Thomas Witkowski committed
377
      GET_PARAMETER(0, name + "->input->read serialization", "%d", 
378 379
		    &readSerialization);
      if (readSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
380
	GET_PARAMETER(0, name + "->input->serialization filename", 
381 382 383 384
		      &serializationFilename);
	TEST_EXIT(serializationFilename != "")("no serialization file\n");

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
Thomas Witkowski's avatar
Thomas Witkowski committed
385
	std::ifstream in(serializationFilename.c_str());
386 387 388
	deserialize(in);
	in.close();
      } else {
389 390
	if (initFlag.isSet(INIT_MESH) && mesh && !mesh->isInitialized()) {
	  mesh->initialize();
391 392 393
	}
    
	// === read value file and use it for the mesh values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
394
	std::string valueFilename("");
395
	GET_PARAMETER(0, mesh->getName() + "->value file name", &valueFilename);     
396
	if (valueFilename.length()) {     
397
	  ValueReader::readValue(valueFilename,
398 399 400 401
				 mesh,
				 solution,
				 mesh->getMacroFileInfo());
	  mesh->clearMacroFileInfo();
402
	}
403
	
404
	// === do global refinements ===
405 406
	if (initFlag.isSet(INIT_GLOBAL_REFINES)) {
	  int globalRefinements = 0;
407
	  GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinements);
408
	  refinementManager->globalRefine(mesh, globalRefinements);	
409
	}
410 411 412 413 414 415 416 417
      }
    }
  }

  void ProblemScal::createFESpace()
  {
    // create finite element space
    int degree = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
418
    GET_PARAMETER(0, name + "->polynomial degree" ,"%d", &degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
419
    feSpace = FiniteElemSpace::provideFESpace(NULL,
420 421
					       Lagrange::getLagrange(mesh->getDim(), degree),
					       mesh,
Thomas Witkowski's avatar
Thomas Witkowski committed
422
					       name + "->feSpace");  
423 424

    // create dof admin for vertex dofs if neccessary
425 426
    if (mesh->getNumberOfDOFs(VERTEX) == 0) {
      DimVec<int> ln_dof(mesh->getDim(), DEFAULT_VALUE, 0);
427
      ln_dof[VERTEX]= 1;
428
      mesh->createDOFAdmin("vertex dofs", ln_dof);
429 430 431 432 433 434
    }
  }

  void ProblemScal::createMatricesAndVectors()
  {
    // === create vectors and system matrix ===
435 436 437
    systemMatrix = new DOFMatrix(feSpace, feSpace, "system matrix");
    rhs = new DOFVector<double>(feSpace, "rhs");
    solution = new DOFVector<double>(feSpace, "solution");
438

439 440
    solution->setCoarsenOperation(COARSE_INTERPOL);
    solution->set(0.0);      /*  initialize u_h  !                      */
441 442 443 444 445
  }

  void ProblemScal::createSolver()
  {
    // === create solver ===
446
    std::string solverType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
447
    GET_PARAMETER(0, name + "->solver", &solverType);
448 449
    OEMSolverCreator *solverCreator = 
      dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
450
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    solverCreator->setName(name + "->solver");
452 453
    solver = solverCreator->create();
    solver->initParameters();
454 455 456 457 458
  }

  void ProblemScal::createEstimator()
  {
    // create and set leaf data prototype
459
    mesh->setElementDataPrototype(new LeafDataEstimatable(new LeafDataCoarsenable));
460 461

    // === create estimator ===
462
    std::string estimatorType("0");
Thomas Witkowski's avatar
Thomas Witkowski committed
463
    GET_PARAMETER(0, name + "->estimator", &estimatorType);
464 465 466
    EstimatorCreator *estimatorCreator = 
      dynamic_cast<EstimatorCreator*>(
				      CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
467
    if (estimatorCreator) {
Thomas Witkowski's avatar
Thomas Witkowski committed
468
      estimatorCreator->setName(name + "->estimator");
Thomas Witkowski's avatar
Thomas Witkowski committed
469
      if (estimatorType == "recovery") {
470
	dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->setSolution(solution);
471
      }
472
      estimator = estimatorCreator->create();
473 474

      // init estimator
475
      estimator->addSystem(systemMatrix, solution, rhs);
476 477 478 479 480
    }
  }

  void ProblemScal::createMarker()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
481
    marker = dynamic_cast<Marker*>(Marker::createMarker(name + "->marker", -1));
482 483 484 485
  }

  void ProblemScal::createFileWriter()
  {
486
    fileWriters.push_back(new FileWriter(name + "->output", mesh, solution));
487
    int writeSerialization = 0;
488 489 490 491 492
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);

    // 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
493 494
    if (writeSerialization)
      fileWriters.push_back(new Serializer<ProblemScal>(this));
495
#endif
496 497 498 499 500 501
  }

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

502
    if (estimator) {
503
      clock_t first = clock();
504
      estimator->estimate(adaptInfo->getTimestep());
Thomas Witkowski's avatar
Thomas Witkowski committed
505
      INFO(info,8)("estimation of the error needed %.5f seconds\n",
506 507
		    TIME_USED(first,clock()));

508
      adaptInfo->setEstSum(estimator->getErrorSum(), 0);
509 510 511
      adaptInfo->setTimeEstSum(estimator->getTimeEst(), 0);
      adaptInfo->setEstMax(estimator->getErrorMax(), 0);
      adaptInfo->setTimeEstMax(estimator->getTimeEstMax(), 0);
512 513 514 515 516
    } else {
      WARNING("no estimator\n");
    }
  }

517 518
  void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				      bool assembleMatrix, bool assembleVector)
519
  {
520
    FUNCNAME("ProblemScal::buildAfterCoarsen()");
521 522 523

    clock_t first = clock();

524
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
525
    mesh->dofCompress();
Thomas Witkowski's avatar
Thomas Witkowski committed
526
#endif
527 528

    MSG("%d DOFs for %s\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
529
	feSpace->getAdmin()->getUsedDOFs(), feSpace->getName().c_str());
530

531
    rhs->set(0.0);
532

533 534 535 536 537 538
    DOFMatrix::base_matrix_type& base_matrix= systemMatrix->getBaseMatrix(); 

    int nnz_per_row= 0;
    if (num_rows(base_matrix) != 0)
      nnz_per_row= int(double(base_matrix.nnz()) / num_rows(base_matrix) * 1.2); 
    if (nnz_per_row < 5) 
Thomas Witkowski's avatar
Thomas Witkowski committed
539
      nnz_per_row= 5;
540 541 542 543

    // Correct dimensionality of matrix
    base_matrix.change_dim(feSpace->getAdmin()->getUsedSize(), 
			   feSpace->getAdmin()->getUsedSize());
544
    set_to_zero(base_matrix);
545 546 547 548
   
    // Reuse old sparsity information (if available) or default
    systemMatrix->startInsertion(nnz_per_row);

549
    // fill boundary conditions
550 551 552 553 554 555
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->initMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->initVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->initVector(solution);
556

557
    BoundaryType *bound;
558
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
559 560 561 562 563 564
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | 
					           Mesh::FILL_BOUND |
					           Mesh::FILL_COORDS |
					           Mesh::FILL_DET |
					           Mesh::FILL_GRD_LAMBDA |
					           Mesh::FILL_NEIGH);
565

566
    while (elInfo) {
567
      if (useGetBound) {
Thomas Witkowski's avatar
Thomas Witkowski committed
568
	bound = new BoundaryType[feSpace->getBasisFcts()->getNumber()];
Thomas Witkowski's avatar
Thomas Witkowski committed
569
	feSpace->getBasisFcts()->getBound(elInfo, bound);
570 571 572
      } else {
	bound = NULL;
      }
573

574 575 576
      systemMatrix->assemble(1.0, elInfo, bound);
      rhs->assemble(1.0, elInfo, bound);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
577
      if (useGetBound)
Thomas Witkowski's avatar
Thomas Witkowski committed
578
	delete [] bound;
579
      
580 581 582 583 584 585
      if (systemMatrix->getBoundaryManager()) 
	systemMatrix->getBoundaryManager()->fillBoundaryConditions(elInfo, systemMatrix);
      if (rhs->getBoundaryManager())
	rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
      if (solution->getBoundaryManager())
	solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
586 587 588
      elInfo = stack.traverseNext(elInfo);
    }

589
    systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
Thomas Witkowski's avatar
Thomas Witkowski committed
590

591
    systemMatrix->finishInsertion();
592
  
593
    // TODO: ExitMatrix should be called after finishInsertion!
594 595 596 597 598 599
    if (systemMatrix->getBoundaryManager())
      systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix);
    if (rhs->getBoundaryManager())
      rhs->getBoundaryManager()->exitVector(rhs);
    if (solution->getBoundaryManager())
      solution->getBoundaryManager()->exitVector(solution);
600

601 602
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first,clock()));

603
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
604 605 606 607 608
//     PetscErrorCode ierr;
//     Mat A;

//     ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr);
#endif
609
  
610
    createPrecon();
611 612
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
613
  void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, std::string name)
614
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
615
    FUNCNAME("ProblemScal::writeResidualMesh()");
616

Thomas Witkowski's avatar
Thomas Witkowski committed
617
    std::map<int, double> vec;    
618
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
619 620
    ElInfo *elInfo = stack.traverseFirst(this->getMesh(),  -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
621 622
    
    while (elInfo) {		  
Thomas Witkowski's avatar
Thomas Witkowski committed
623
      vec[elInfo->getElement()->getIndex()] = elInfo->getElement()->getEstimation(0);
624 625 626
      elInfo = stack.traverseNext(elInfo);
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
627
    ElementFileWriter fw(name, this->getFESpace(), vec);
628 629 630
    fw.writeFiles(adaptInfo, true);    
  }

631 632
  void ProblemScal::createPrecon()
  {
633
    std::string preconType("no");
634 635 636 637 638 639 640 641 642 643 644 645 646 647
    GET_PARAMETER(0, name + "->solver->left precon", &preconType);

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

    solver->setLeftPrecon( preconCreator->create(systemMatrix->getBaseMatrix()) );

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

    preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
    solver->setRightPrecon( preconCreator->create(systemMatrix->getBaseMatrix()) );
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
648
  void ProblemScal::serialize(std::ostream &out) 
649 650 651
  {
    FUNCNAME("ProblemScal::serialize()");

652 653
    mesh->serialize(out);
    solution->serialize(out);
654 655
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
656
  void ProblemScal::deserialize(std::istream &in) 
657 658 659
  {
    FUNCNAME("ProblemScal::deserialize()");

660 661
    mesh->deserialize(in);    
    solution->deserialize(in);
662 663 664
  }

}