Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

ProblemVec.cc 42 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <sstream>
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#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
16
#include "DualTraverse.h"
17
18
19
20
21
22
#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
23
#include "Flag.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
#include "TraverseParallel.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "VtkWriter.h"
26
#include "ValueReader.h"
27
28
29
30
31
32
33
34
35
36

namespace AMDiS {

  void ProblemVec::initialize(Flag initFlag,
			      ProblemVec *adoptProblem,
			      Flag adoptFlag)
  {
    FUNCNAME("ProblemVec::initialize()");
    
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
37
    if (meshes.size() != 0) { 
38
39
40
41
42
43
44
45
46
47
48
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
	  ((!adoptFlag.isSet(INIT_MESH))&&
	   (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
49
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
50
	componentMeshes = adoptProblem->componentMeshes;
Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
53
54
55

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

65
66
67
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
68
    if (meshes.size() == 0) 
69
70
71
      WARNING("no mesh created\n");

    // === create fespace ===
72
    if (feSpaces.size() != 0) {
73
74
75
76
77
78
79
80
      WARNING("feSpaces already created\n");
    } else {
      if (initFlag.isSet(INIT_FE_SPACE) || 
	  (initFlag.isSet(INIT_SYSTEM)&&!adoptFlag.isSet(INIT_FE_SPACE))) {
	createFESpace();
      } 
      if (adoptProblem &&
	  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
81
82
	feSpaces = adoptProblem->getFESpaces();
	componentSpaces = adoptProblem->componentSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
83
	traverseInfo = adoptProblem->traverseInfo;
84
85
86

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

96
97
98
      }
    }

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

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

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

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

    // === create estimator ===
    if (initFlag.isSet(INIT_ESTIMATOR)) {
      createEstimator();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
133
      estimator = adoptProblem->getEstimator();
134
135
136
137
138
139
140
    } 

    // === create marker ===
    if (initFlag.isSet(INIT_MARKER)) {
      createMarker();
    } 
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) {
141
      marker = adoptProblem->getMarker();
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
    } 


    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER)) {
      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;
160
    std::string serializationFilename = "";
161
162
163
164
165
166
167
    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;

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

	MSG("Deserialization from file: %s\n", serializationFilename.c_str());
183
	std::ifstream in(serializationFilename.c_str());
184
185
186
	deserialize(in);
	in.close();
      } else {
187
	int globalRefinements = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
188
	GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d", 
189
190
		      &globalRefinements);

191
	// Initialize the meshes if there is no serialization file.
192
193
	for (int i = 0; i < static_cast<int>(meshes.size()); i++)
	  if (initFlag.isSet(INIT_MESH) && meshes[i] && !(meshes[i]->isInitialized()))
Thomas Witkowski's avatar
Thomas Witkowski committed
194
	    meshes[i]->initialize();	    
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

	// === read value file and use it for the mesh values ===
	std::string valueFilename("");
	GET_PARAMETER(0, meshes[0]->getName() + "->value file name", &valueFilename); 
	if (valueFilename.length()) {     
	  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
210
	    refinementManager->globalRefine(meshes[i], globalRefinements);
211
212
213
214
215
216
217
218
219
220
      }
    }

    doOtherStuff();
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
221
    componentMeshes.resize(nComponents);
222
    std::map<int, Mesh*> meshForRefinementSet;
223
224
    char number[3];

225
    std::string meshName("");
Thomas Witkowski's avatar
Thomas Witkowski committed
226
    GET_PARAMETER(0, name + "->mesh", &meshName);
227
    TEST_EXIT(meshName != "")("no mesh name specified\n");
228
    int dim = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
230
    TEST_EXIT(dim)("no problem dimension specified!\n");
231

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

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

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

274
    componentSpaces.resize(nComponents, NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
275
    traverseInfo.resize(nComponents);
276

277
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
278
      char number[3];
279
      sprintf(number, "%d", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
      int degree = 1;
      GET_PARAMETER(0, name + "->polynomial degree[" + number + "]","%d", &degree);
282

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

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

289
	FiniteElemSpace *newFESpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
	  FiniteElemSpace::provideFESpace(NULL, Lagrange::getLagrange(dim, degree),
					  componentMeshes[i], s.str());
Thomas Witkowski's avatar
Thomas Witkowski committed
292
	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFESpace;
293
	feSpaces.push_back(newFESpace);
294
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
297
298
299
300
301
302
      componentSpaces[i] = feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
    }

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	traverseInfo.getMatrix(i, j).setFESpace(componentSpaces[i], componentSpaces[j]);
      }
      traverseInfo.getVector(i).setFESpace(componentSpaces[i]);
303
304
305
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
308
    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);
309
	ln_dof[VERTEX]= 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
310
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
311
312
313
314
315
316
317
318
319
320
      }
    }
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
323
324
    systemMatrix = NEW Matrix<DOFMatrix*>(nComponents, nComponents);
    systemMatrix->set(NULL);
    rhs = NEW SystemVector("rhs", componentSpaces, nComponents);
    solution = NEW SystemVector("solution", componentSpaces, nComponents);
325
326

    char number[10];
327
    std::string numberedName;
Thomas Witkowski's avatar
Thomas Witkowski committed
328
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
329
      (*systemMatrix)[i][i] = NEW DOFMatrix(componentSpaces[i], 
330
					     componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
331
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
332
      sprintf(number, "[%d]", i);
333
      numberedName = "rhs" + std::string(number);
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
      rhs->setDOFVector(i, NEW DOFVector<double>(componentSpaces[i], numberedName));
      numberedName = name + std::string(number);
      solution->setDOFVector(i, NEW DOFVector<double>(componentSpaces[i], 
337
						       numberedName));
Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
340
341
342
343
344
345
346
347
    }
  }

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

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

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

    // create and set leaf data prototype
Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
      meshes[i]->setElementDataPrototype
365
366
367
368
	(NEW LeafDataEstimatableVec(NEW LeafDataCoarsenableVec));
    }  

    char number[3];
369
    std::string estName;
370

Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
    for (int i = 0; i < nComponents; i++) {
      TEST_EXIT(estimator[i] == NULL)("estimator already created\n");
373
      sprintf(number, "%d", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
374
      estName = name + "->estimator[" + std::string(number) + "]";
375
376

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


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

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

407
    std::string numberedName;
408
    char number[10];
Thomas Witkowski's avatar
Thomas Witkowski committed
409
    int nMarkersCreated = 0;
410

411
    for (int i = 0; i < nComponents; i++) {
412
      sprintf(number, "[%d]", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
413
      numberedName = name + "->marker" + std::string(number);
414
415
      marker[i] = Marker::createMarker(numberedName, i);
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
418
419
420
421
	nMarkersCreated++;

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

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

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

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

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

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

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


    // Create own filewriters for each components of the problem
    char number[10];
454
    for (int i = 0; i < nComponents; i++) {
455
      sprintf(number, "[%d]", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
456
      numberedName  = name + "->output" + std::string(number);
457
458
459
460
      filename = "";
      GET_PARAMETER(0, numberedName + "->filename", &filename);

      if (filename != "") {
Thomas Witkowski's avatar
Thomas Witkowski committed
461
	fileWriters.push_back(NEW FileWriter(numberedName, 
Thomas Witkowski's avatar
Thomas Witkowski committed
462
					      componentMeshes[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
463
					      solution->getDOFVector(i)));
464
465
466
467
468
469
      }
    }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
477
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
478
    if (writeSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
479
      fileWriters.push_back(NEW Serializer<ProblemVec>(this));
480
481
482
483
484
485
486
    }
  }

  void ProblemVec::doOtherStuff()
  {
  }

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

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

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

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

503
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
506
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
509
510
511
512
#endif


    adaptInfo->setSolverIterations(iter);
Thomas Witkowski's avatar
Thomas Witkowski committed
513
514
515
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
516
517
518
519
520
521
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
522
523
    VtkWriter::writeFile(solution->getDOFVector(0), "test.vtu");

524
525
    clock_t first = clock();

526
527
528
529
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

530
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
531
      computeError(adaptInfo);
532
533
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
534
	Estimator *scalEstimator = estimator[i];
535
536
537
538
539
540
541
542
543
544
	
	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);
	}
545
546
547
      }
    }

548
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
549
550
    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);
551
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
554
555
556

#endif

557
558
559
560
561
562
563
564
565
566
567
  }

  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;
568
    for (int i = 0; i < nComponents; i++) {
569
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
570
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
571
572
573
574
      } else {
	WARNING("no marker for component %d\n", i);
      }
    }
575
    
576
577
578
579
580
581
582
    return markFlag;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
583
    int nMeshes = static_cast<int>(meshes.size());
584
    Flag refineFlag = 0;
585
    for (int i = 0; i < nMeshes; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
586
      refineFlag |= refinementManager->refineMesh(meshes[i]);
587
588
589
590
591
592
593
594
    }
    return refineFlag;
  }

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

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

	WARNING("coarsening for component %d no allowed\n", i);
      }
    }
    return coarsenFlag;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
611
    if (allowFirstRef) {
612
      for (int i = 0; i < nComponents; i++) {
613
614
	adaptInfo->allowRefinement(true, i);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
615
      allowFirstRef = false;
616
    } else {
617
      for (int i = 0; i < nComponents; i++) {
618
619
620
621
622
623
624
625
626
627
628
	if (adaptInfo->spaceToleranceReached(i)) {
	  adaptInfo->allowRefinement(false, i);
	} else {
	  adaptInfo->allowRefinement(true, i);	
	}
      }
    }

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

629
630
  void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				     bool asmMatrix, bool asmVector)
631
632
633
634
635
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

    clock_t first = clock();

636
637
638
639
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
640
641
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
      meshes[i]->dofCompress();
642
643
644
645
    }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
658

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

Thomas Witkowski's avatar
Thomas Witkowski committed
664
      rhs->getDOFVector(i)->set(0.0);
665
      for (int j = 0; j < nComponents; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
666
	if ((*systemMatrix)[i][j]) {
667
668
669
670
671
672
673
674
675
676
677
678
679
680
	  DOFMatrix*                   dof_matrix= (*systemMatrix)[i][j];
	  DOFMatrix::base_matrix_type& base_matrix= dof_matrix->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) 
	      nnz_per_row= 5;

	  // Correct dimensionality of matrix
	  base_matrix.change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
				 componentSpaces[j]->getAdmin()->getUsedSize());
	  
	  // Reuse old sparsity information (if available) or default
681
	  //	  dof_matrix->startInsertion(nnz_per_row);
682
	}
683
684
685
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
    traverseInfo.updateStatus();

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

691
692
693
694
	// 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
695
	DOFMatrix *matrix = (*systemMatrix)[i][j];
696

697
698
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
699
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j])
700
	  assembleMatrix = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
701
702

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

706
707
708
709
710
711
712
	// 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) {
	  continue;
	}

713
714
715
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
	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
735
	    // least one another aux fe space. 
Thomas Witkowski's avatar
Thomas Witkowski committed
736
737
738
739
740
741

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

742
743
	    assembleOnDifMeshes2(componentSpaces[i], 
				 traverseInfo.getAuxFESpace(i, j),
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
				 assembleFlag,
				 NULL,
746
				 ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
750
751
752
753

	  } 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
754
755
756
	  assembleOnOneMesh(componentSpaces[i],
			    assembleFlag,
			    assembleMatrix ? matrix : NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
757
			    ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
758
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
761
762
	  assembleOnDifMeshes2(componentSpaces[i],
			       traverseInfo.getAuxFESpace(i, j),
			       assembleFlag,
			       assembleMatrix ? matrix : NULL,
763
			       ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
764
765
766
767

	} 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
768
769
770
	  assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
			      assembleFlag,
			      assembleMatrix ? matrix : NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
771
772
773
			      ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);	  
	} else {
	  ERROR_EXIT("Not yet implemented!\n");
774
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
775

776
777
778
	// TODO: ExitMatrix should be called after finishInsertion!
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	  
779
	
Thomas Witkowski's avatar
Thomas Witkowski committed
780
	assembledMatrix[i][j] = true;
781
782
      }

783
      // And now assemble boundary conditions on the vectors
Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
      assembleBoundaryConditions(rhs->getDOFVector(i),
				 solution->getDOFVector(i),
786
787
				 componentMeshes[i],
				 assembleFlag);
788
    }
789

790
791
    int nnz = 0;

792
793
794
    // Finish insertion
    for (int i = 0; i < nComponents; i++) 
      for (int j = 0; j < nComponents; j++) 
795
 	if ((*systemMatrix)[i][j]) {
796
 	  (*systemMatrix)[i][j]->finishInsertion(); 
797
798
	  nnz += (*systemMatrix)[i][j]->getBaseMatrix().nnz();
	}
799

800
    //    clock_t first1 = clock();
801
    solverMatrix.setMatrix(*systemMatrix);
802
    //    clock_t first2 = clock();
803
    createPrecon();
804
805
806
807
    //    clock_t first3 = clock();

    //    std::cout << "T1 = " << TIME_USED(first1, first2) << std::endl;
    //    std::cout << "T2 = " << TIME_USED(first2, first3) << std::endl;
808

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

811
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
812
813
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
814
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
817
#endif
818
819
  }

820
821
  void ProblemVec::createPrecon()
  {
822
    std::string preconType("no");
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
    GET_PARAMETER(0, name + "->solver->left precon", &preconType);

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

    solver->setLeftPrecon( preconCreator->create(solverMatrix.getMatrix()) );

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

    preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
    solver->setRightPrecon( preconCreator->create(solverMatrix.getMatrix()) );
  }


838
839
840
841
  void ProblemVec::writeFiles(AdaptInfo *adaptInfo, bool force) 
  {
    FUNCNAME("ProblemVec::writeFiles()");

842
843
844
845
846
847
848
    clock_t first = clock();

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

    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
849
    int size = static_cast<int>(fileWriters.size());
850
851
852
853
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 1)
#endif
    for (i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
854
      fileWriters[i]->writeFiles(adaptInfo, force);
855
    }
856
857
    
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
858
859
860
    INFO(info, 8)("writeFiles needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()),
		  omp_get_wtime() - wtime);
861
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
862
863
    INFO(info, 8)("writeFiles needed %.5f seconds\n",
		  TIME_USED(first, clock()));
864
#endif
865
866
  }

867
  void ProblemVec::interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct) 
868
869
870
  {
    FUNCNAME("ProblemVec::interpolInitialSolution()");

Thomas Witkowski's avatar
Thomas Witkowski committed
871
    solution->interpol(fct);
872
873
874
875
876
  }

  void ProblemVec::addMatrixOperator(Operator *op, 
				     int i, int j,
				     double *factor,
877
878
				     double *estFactor,
				     bool fake)
879
880
  {
    FUNCNAME("ProblemVec::addMatrixOperator()");
Thomas Witkowski's avatar
Thomas Witkowski committed
881
   
Thomas Witkowski's avatar
Thomas Witkowski committed
882
    if (!(*systemMatrix)[i][j]) {
883
      TEST_EXIT(i != j)("should have been created already\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
884
885
886
887
      (*systemMatrix)[i][j] = NEW DOFMatrix(componentSpaces[i], componentSpaces[j], "");
      (*systemMatrix)[i][j]->setCoupleMatrix(true);
      (*systemMatrix)[i][j]->getBoundaryManager()->
	setBoundaryConditionMap((*systemMatrix)[i][i]->getBoundaryManager()->
888
889
				getBoundaryConditionMap());
    }    
Thomas Witkowski's avatar
Thomas Witkowski committed
890

891
    (*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
Thomas Witkowski's avatar
Thomas Witkowski committed
892

893
894
895
896
897
898
899
900
901
902
903
    if (!fake) {    
      traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces()); 
      
      for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
	if ((op->getAuxFESpaces())[k]->getMesh() != componentSpaces[i]->getMesh() ||
	    (op->getAuxFESpaces())[k]->getMesh() != componentSpaces[j]->getMesh()) {
	  op->setNeedDualTraverse(true);
	  break;
	}
      }    
    } 
904
905
906
907
  }

  void ProblemVec::addVectorOperator(Operator *op, int i,
				     double *factor,
908
909
				     double *estFactor,
				     bool fake)
910
911
912
  {
    FUNCNAME("ProblemVec::addVectorOperator()");

Thomas Witkowski's avatar
Thomas Witkowski committed
913
914
    rhs->getDOFVector(i)->addOperator(op, factor, estFactor);

915
916
917
918
919
920
921
922
923
924
    if (!fake) {
      traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces()); 
      
      for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
	if ((op->getAuxFESpaces())[j]->getMesh() != componentSpaces[i]->getMesh()) {
	  op->setNeedDualTraverse(true);
	  break;
	}
      }    
    }
925
926
927
928
929
930
931
932
933
  }

  void ProblemVec::addDirichletBC(BoundaryType type, int system,
				  AbstractFunction<double, WorldVector<double> >* b)
  {
    FUNCNAME("ProblemVec::addDirichletBC()");

    DirichletBC *dirichlet = new DirichletBC(type, 
					     b, 
934
					     componentSpaces[system]);
935
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
936
937
      if (systemMatrix && (*systemMatrix)[system][i]) {
	(*systemMatrix)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet);
938
939
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
940
941
942
943
    if (rhs)
      rhs->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
    if (solution)
      solution->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
944
945
946
947
948
949
950
951
952
  }

  void ProblemVec::addNeumannBC(BoundaryType type, int row, int col, 
				AbstractFunction<double, WorldVector<double> > *n)
  {
    FUNCNAME("ProblemVec::addNeumannBC()");

    NeumannBC *neumann = 
      new NeumannBC(type, n, 
953
954
		    componentSpaces[row], 
		    componentSpaces[col]);
Thomas Witkowski's avatar
Thomas Witkowski committed
955
956
    if (rhs)
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
957
958
959
960
961
962
963
964
965
966
  }

  void ProblemVec::addRobinBC(BoundaryType type, int row, int col, 
			      AbstractFunction<double, WorldVector<double> > *n,
			      AbstractFunction<double, WorldVector<double> > *r)
  {
    FUNCNAME("ProblemVec::addRobinBC()");

    RobinBC *robin = 
      new RobinBC(type, n, r, 
967
968
		  componentSpaces[row], 
		  componentSpaces[col]);
Thomas Witkowski's avatar
Thomas Witkowski committed
969
970
    if (rhs)
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
971

Thomas Witkowski's avatar
Thomas Witkowski committed
972
973
    if (systemMatrix && (*systemMatrix)[row][col]) {
      (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
974
975
976
977
978
979
980
    }
  }

  void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col) 
  {
    FUNCNAME("ProblemVec::addPeriodicBC()");

981
    FiniteElemSpace *feSpace = componentSpaces[row];
982
983
984

    PeriodicBC *periodic = new PeriodicBC(type, feSpace);

Thomas Witkowski's avatar
Thomas Witkowski committed
985
986
    if (systemMatrix && (*systemMatrix)[row][col]) 
      (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);
987

Thomas Witkowski's avatar
Thomas Witkowski committed
988
    if (rhs) 
Thomas Witkowski's avatar
Thomas Witkowski committed
989
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(periodic);
990
991
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
992
993
  void ProblemVec::assembleOnOneMesh(FiniteElemSpace *feSpace, 
				     Flag assembleFlag,
Thomas Witkowski's avatar
Thomas Witkowski committed
994
995
996
997
				     DOFMatrix *matrix, DOFVector<double> *vector)
  {
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
998
999

#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
1000
    TraverseParallelStack stack;
1001
1002
1003
#else
    TraverseStack stack;
#endif   
Thomas Witkowski's avatar
Thomas Witkowski committed
1004

1005
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
1006
#pragma omp parallel
1007
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
1008
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
1009
1010
      BoundaryType *bound = 
	useGetBound ? GET_MEMORY(BoundaryType, basisFcts->getNumber()) : NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
1011
1012
1013
1014
1015

      // Create for every thread its private matrix and vector, on that
      // the thread will assemble its part of the mesh.
      DOFMatrix *tmpMatrix = NULL;
      DOFVector<double> *tmpVector = NULL; 
Thomas Witkowski's avatar
Thomas Witkowski committed
1016
1017

      if (matrix) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1018
	tmpMatrix = NEW DOFMatrix(matrix->getRowFESpace(), matrix->getColFESpace(), "tmp");
Thomas Witkowski's avatar
Thomas Witkowski committed
1019
1020
1021
1022
1023
1024

	// Copy the global matrix to the private matrix, because we need the
	// operators defined on the global matrix in the private one. Only the
	// values have to be set to zero.
	*tmpMatrix = *matrix;
	tmpMatrix->clear();
1025
1026
1027
1028

	tmpMatrix->getBaseMatrix().change_dim(matrix->getRowFESpace()->getAdmin()->getUsedSize(),
					      matrix->getColFESpace()->getAdmin()->getUsedSize());
	tmpMatrix->startInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
      }

      if (vector) {
	tmpVector = NEW DOFVector<double>(vector->