ProblemVec.cc 41.7 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
#include <sstream>

3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "ProblemVec.h"
#include "RecoveryEstimator.h"
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "FileWriter.h"
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
17
#include "DualTraverse.h"
18
19
20
21
22
23
24
25
#include "Mesh.h"
#include "OEMSolver.h"
#include "Preconditioner.h"
#include "MatVecMultiplier.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include "Flag.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include "TraverseParallel.h"
28
29
30
31
32
33
34
35
36
37

namespace AMDiS {

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

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

66
67
68
      }
    }

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

    // === create fespace ===
73
    if (feSpaces.size() != 0) {
74
75
76
77
78
79
80
81
      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))) {
82
83
	feSpaces = adoptProblem->getFESpaces();
	componentSpaces = adoptProblem->componentSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
84
	traverseInfo = adoptProblem->traverseInfo;
85
86
87

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

97
98
99
      }
    }

100
    if (feSpaces.size() == 0) 
101
102
103
104
105
106
107
      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
108
109
110
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRHS();
      systemMatrix = adoptProblem->getSystemMatrix();
111
112
113
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
126
    if (!solver) 
127
128
129
130
131
132
133
      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
134
      estimator = adoptProblem->getEstimator();
135
136
137
138
139
140
141
    } 

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


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

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

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

    doOtherStuff();
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
211
    componentMeshes.resize(nComponents);
212
    std::map<int, Mesh*> meshForRefinementSet;
213
214
    char number[3];

215
    std::string meshName("");
Thomas Witkowski's avatar
Thomas Witkowski committed
216
    GET_PARAMETER(0, name + "->mesh", &meshName);
217
    TEST_EXIT(meshName != "")("no mesh name specified\n");
218
    int dim = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    GET_PARAMETER(0, name + "->dim", "%d", &dim);
220
    TEST_EXIT(dim)("no problem dimension specified!\n");
221

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

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

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

264
    componentSpaces.resize(nComponents, NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
265
    traverseInfo.resize(nComponents);
266

267
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      char number[3];
269
      sprintf(number, "%d", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
      int degree = 1;
      GET_PARAMETER(0, name + "->polynomial degree[" + number + "]","%d", &degree);
272

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

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

279
	FiniteElemSpace *newFESpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
	  FiniteElemSpace::provideFESpace(NULL, Lagrange::getLagrange(dim, degree),
					  componentMeshes[i], s.str());
Thomas Witkowski's avatar
Thomas Witkowski committed
282
	feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFESpace;
283
	feSpaces.push_back(newFESpace);
284
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
287
288
289
290
291
292
      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]);
293
294
295
    }

    // create dof admin for vertex dofs if neccessary
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
    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);
299
	ln_dof[VERTEX]= 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
300
	meshes[i]->createDOFAdmin("vertex dofs", ln_dof);      
301
302
303
304
305
306
307
308
309
310
      }
    }
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
311
312
313
314
    systemMatrix = NEW Matrix<DOFMatrix*>(nComponents, nComponents);
    systemMatrix->set(NULL);
    rhs = NEW SystemVector("rhs", componentSpaces, nComponents);
    solution = NEW SystemVector("solution", componentSpaces, nComponents);
315
316

    char number[10];
317
    std::string numberedName;
Thomas Witkowski's avatar
Thomas Witkowski committed
318
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
319
      (*systemMatrix)[i][i] = NEW DOFMatrix(componentSpaces[i], 
320
					     componentSpaces[i], "A_ii");
Thomas Witkowski's avatar
Thomas Witkowski committed
321
      (*systemMatrix)[i][i]->setCoupleMatrix(false);
322
      sprintf(number, "[%d]", i);
323
      numberedName = "rhs" + std::string(number);
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
326
      rhs->setDOFVector(i, NEW DOFVector<double>(componentSpaces[i], numberedName));
      numberedName = name + std::string(number);
      solution->setDOFVector(i, NEW DOFVector<double>(componentSpaces[i], 
327
						       numberedName));
Thomas Witkowski's avatar
Thomas Witkowski committed
328
329
      solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
      solution->getDOFVector(i)->set(0.0);
330
331
332
    }

    // === create matVec ===
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    matVec = NEW StandardMatVec<Matrix<DOFMatrix*>, SystemVector>(systemMatrix);
334
335
336
337
338
339
340
  }

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

    // === create solver ===
341
    std::string solverType("no");
Thomas Witkowski's avatar
Thomas Witkowski committed
342
    GET_PARAMETER(0, name + "->solver", &solverType);
343
344
345
346
347
348
    OEMSolverCreator<SystemVector> *solverCreator = 
      dynamic_cast<OEMSolverCreator<SystemVector>*>(
						    CreatorMap<OEMSolver<SystemVector> >
						    ::getCreator(solverType)
						    );
    TEST_EXIT(solverCreator)("no solver type\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
351
    solverCreator->setName(name + "->solver");
    solver = solverCreator->create();
    solver->initParameters();
352
353

    // === create preconditioners ===
354
    std::string preconType("no");
355
356

    PreconditionerScal *scalPrecon;
357
    PreconditionerVec *vecPrecon = NEW PreconditionerVec(nComponents);
358

Thomas Witkowski's avatar
Thomas Witkowski committed
359
    GET_PARAMETER(0, name + "->solver->left precon", &preconType);
360
361
362
363
364
    CreatorInterface<PreconditionerScal> *preconCreator =
      CreatorMap<PreconditionerScal>::getCreator(preconType);

    if (!preconCreator->isNullCreator()) {
      dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
Thomas Witkowski's avatar
Thomas Witkowski committed
365
	setName(name + "->solver->left precon");
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
      for (int i = 0; i < nComponents; i++) {
368
	dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
369
	  setSizeAndRow(nComponents, i);
370
371
    
	scalPrecon = preconCreator->create();
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
	for (int j = 0; j < nComponents; j++) {
	  scalPrecon->setMatrix(&(*systemMatrix)[i][j], j);
374
375
376
	}
	vecPrecon->setScalarPrecon(i, scalPrecon);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
377
      leftPrecon = vecPrecon;
378
379
380
    }


381
    vecPrecon = NEW PreconditionerVec(nComponents);
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    GET_PARAMETER(0, name + "->solver->right precon", &preconType);
384
385
386
    preconCreator = 
      CreatorMap<PreconditionerScal>::getCreator(preconType);

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    if (!preconCreator->isNullCreator()) {
388
      dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
Thomas Witkowski's avatar
Thomas Witkowski committed
389
	setName(name + "->solver->left precon");
390
391


Thomas Witkowski's avatar
Thomas Witkowski committed
392
      for (int i = 0; i < nComponents; i++) {
393
	dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
394
	  setSizeAndRow(nComponents, i);
395
396
    
	scalPrecon = preconCreator->create();
Thomas Witkowski's avatar
Thomas Witkowski committed
397
398
	for (int j = 0; j < nComponents; j++) {
	  scalPrecon->setMatrix(&(*systemMatrix)[i][j], j);
399
400
401
	}
	vecPrecon->setScalarPrecon(i, scalPrecon);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
402
      rightPrecon = vecPrecon;
403
404
405
406
    }


    // === create vector creator ===
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    solver->setVectorCreator(NEW SystemVector::Creator("temp",
408
							componentSpaces, 
409
							nComponents));
410
411
412
413
414
415
416
  }

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

    // create and set leaf data prototype
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
      meshes[i]->setElementDataPrototype
419
420
421
422
	(NEW LeafDataEstimatableVec(NEW LeafDataCoarsenableVec));
    }  

    char number[3];
423
    std::string estName;
424

Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
    for (int i = 0; i < nComponents; i++) {
      TEST_EXIT(estimator[i] == NULL)("estimator already created\n");
427
      sprintf(number, "%d", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
      estName = name + "->estimator[" + std::string(number) + "]";
429
430

      // === create estimator ===
431
      std::string estimatorType("no");
432
433
434
435
      GET_PARAMETER(0, estName, &estimatorType);
      EstimatorCreator *estimatorCreator = 
	dynamic_cast<EstimatorCreator*>(
					CreatorMap<Estimator>::getCreator(estimatorType));
Thomas Witkowski's avatar
Thomas Witkowski committed
436
      if (estimatorCreator) {
437
438
	estimatorCreator->setName(estName);
	estimatorCreator->setRow(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
439
	if (estimatorType == "recovery") {
440
	  dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->
Thomas Witkowski's avatar
Thomas Witkowski committed
441
	    setSolution(solution->getDOFVector(i));
442
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
443
	estimator[i] = estimatorCreator->create();
444
445
446
      }


Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449
450
451
      if (estimator[i]) {
	for (int j = 0; j < nComponents; j++) {
	  estimator[i]->addSystem((*systemMatrix)[i][j], 
				   solution->getDOFVector(j), 
				   rhs->getDOFVector(j));
452
453
454
455
456
457
458
459
460
	}
      }
    }
  }

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

461
    std::string numberedName;
462
    char number[10];
Thomas Witkowski's avatar
Thomas Witkowski committed
463
    int nMarkersCreated = 0;
464

465
    for (int i = 0; i < nComponents; i++) {
466
      sprintf(number, "[%d]", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
467
      numberedName = name + "->marker" + std::string(number);
468
469
      marker[i] = Marker::createMarker(numberedName, i);
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
472
473
474
475
	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);
476
477
478
479
480
481
482
483
484
485
      }
    }
  }

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

    // Create one filewriter for all components of the problem
Thomas Witkowski's avatar
Thomas Witkowski committed
486
    std::string numberedName  = name + "->output";
487
    std::string filename = "";
488
489
490
    GET_PARAMETER(0, numberedName + "->filename", &filename);

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
497
	solutionList[i] = solution->getDOFVector(i);
498
499
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
500
      fileWriters.push_back(NEW FileWriter(numberedName,
Thomas Witkowski's avatar
Thomas Witkowski committed
501
					    componentMeshes[0],
502
503
504
505
506
507
					    solutionList));
    }


    // Create own filewriters for each components of the problem
    char number[10];
508
    for (int i = 0; i < nComponents; i++) {
509
      sprintf(number, "[%d]", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
510
      numberedName  = name + "->output" + std::string(number);
511
512
513
514
      filename = "";
      GET_PARAMETER(0, numberedName + "->filename", &filename);

      if (filename != "") {
Thomas Witkowski's avatar
Thomas Witkowski committed
515
	fileWriters.push_back(NEW FileWriter(numberedName, 
Thomas Witkowski's avatar
Thomas Witkowski committed
516
					      componentMeshes[i], 
Thomas Witkowski's avatar
Thomas Witkowski committed
517
					      solution->getDOFVector(i)));
518
519
520
521
522
523
      }
    }


    // Check for serializer
    int writeSerialization = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
524
    GET_PARAMETER(0, name + "->write serialization", "%d", &writeSerialization);
525
    if (writeSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
      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());
528
529
530
      ERROR_EXIT("Usage of an obsolete parameter (see message above)!\n");
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
531
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
532
    if (writeSerialization) {
Thomas Witkowski's avatar
Thomas Witkowski committed
533
      fileWriters.push_back(NEW Serializer<ProblemVec>(this));
534
535
536
537
538
539
540
    }
  }

  void ProblemVec::doOtherStuff()
  {
  }

541
  void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
542
543
544
  {
    FUNCNAME("Problem::solve()");

Thomas Witkowski's avatar
Thomas Witkowski committed
545
    if (!solver) {
546
547
548
549
550
551
552
553
554
      WARNING("no solver\n");
      return;
    }

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

    clock_t first = clock();
Thomas Witkowski's avatar
Thomas Witkowski committed
555
    int iter = solver->solve(matVec, solution, rhs, leftPrecon, rightPrecon, fixedMatrix);
556
    
557
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
    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);
560
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
561
562
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
563
564
565
566
#endif


    adaptInfo->setSolverIterations(iter);
Thomas Witkowski's avatar
Thomas Witkowski committed
567
568
569
    adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
    adaptInfo->setSolverTolerance(solver->getTolerance());
    adaptInfo->setSolverResidual(solver->getResidual());
570
571
572
573
574
575
576
577
  }

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

    clock_t first = clock();

578
579
580
581
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

582
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
583
      computeError(adaptInfo);
584
585
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
586
	Estimator *scalEstimator = estimator[i];
587
588
589
590
591
592
593
594
595
596
	
	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);
	}
597
598
599
      }
    }

600
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
    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);
603
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
604
605
    INFO(info, 8)("estimation of the error needed %.5f seconds\n",
		  TIME_USED(first, clock()));
606
607
608

#endif

609
610
611
612
613
614
615
616
617
618
619
  }

  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;
620
    for (int i = 0; i < nComponents; i++) {
621
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
622
	markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
623
624
625
626
      } else {
	WARNING("no marker for component %d\n", i);
      }
    }
627
    
628
629
630
631
632
633
634
    return markFlag;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
635
    int nMeshes = static_cast<int>(meshes.size());
636
    Flag refineFlag = 0;
637
    for (int i = 0; i < nMeshes; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
638
      refineFlag |= refinementManager->refineMesh(meshes[i]);
639
640
641
642
643
644
645
646
    }
    return refineFlag;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
647
    int nMeshes = static_cast<int>(meshes.size());
648
    Flag coarsenFlag = 0;
649
650
    for (int i = 0; i < nMeshes; i++) {
      if (adaptInfo->isCoarseningAllowed(i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
651
	coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
652
653
654
655
656
657
658
659
660
661
662

	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
663
    if (allowFirstRef) {
664
      for (int i = 0; i < nComponents; i++) {
665
666
	adaptInfo->allowRefinement(true, i);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
667
      allowFirstRef = false;
668
    } else {
669
      for (int i = 0; i < nComponents; i++) {
670
671
672
673
674
675
676
677
678
679
680
	if (adaptInfo->spaceToleranceReached(i)) {
	  adaptInfo->allowRefinement(false, i);
	} else {
	  adaptInfo->allowRefinement(true, i);	
	}
      }
    }

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

681
682
  void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
				     bool asmMatrix, bool asmVector)
683
684
685
686
687
  {
    FUNCNAME("ProblemVec::buildAfterCoarsen()");

    clock_t first = clock();

688
689
690
691
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
    for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
      meshes[i]->dofCompress();
694
695
696
697
    }

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
700
701
702
703
704
705
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
706
    if (useGetBound) {
707
708
709
      assembleFlag |= Mesh::FILL_BOUND;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
710

711
    for (int i = 0; i < nComponents; i++) {
712
      MSG("%d DOFs for %s\n", 
713
714
	  componentSpaces[i]->getAdmin()->getUsedSize(), 
	  componentSpaces[i]->getName().c_str());
715

Thomas Witkowski's avatar
Thomas Witkowski committed
716
      rhs->getDOFVector(i)->set(0.0);
717
      for (int j = 0; j < nComponents; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
718
	if ((*systemMatrix)[i][j]) {
719
720
	  // The matrix should not be deleted, if it was assembled before
	  // and it is marked to be assembled only once.
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
	  if (!(assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) && asmMatrix) {
	    (*systemMatrix)[i][j]->clear();
723
724
	  }
	}
725
726
727
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
728
729
    traverseInfo.updateStatus();

Thomas Witkowski's avatar
Thomas Witkowski committed
730
    for (int i = 0; i < nComponents; i++) {
731
      for (int j = 0; j < nComponents; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
732

733
734
735
736
	// 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
737
	DOFMatrix *matrix = (*systemMatrix)[i][j];
738

739
740
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
Thomas Witkowski's avatar
Thomas Witkowski committed
741
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
742
743
	  assembleMatrix = false;
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745

	// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
746
747
748
749
	if (!matrix) {
	  assembleMatrix = false;
	}

750
751
752
753
	if (!asmMatrix) {
	  assembleMatrix = false;
	}

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

761
762
763
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

Thomas Witkowski's avatar
Thomas Witkowski committed
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
	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
	    // least one another aux fe space. In this case, do not assemble the rhs,
	    // this will be done afterwards.

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

	    assembleOnDifMeshes2(componentSpaces[1], componentSpaces[0],
				 assembleFlag,
				 NULL,
				 rhs->getDOFVector(1));

	  } 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
802
803
804
	  assembleOnOneMesh(componentSpaces[i],
			    assembleFlag,
			    assembleMatrix ? matrix : NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
805
806
807
808
809
810
811
812
813
814
815
			    ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);

	  assembleOnDifMeshes2(componentSpaces[i],
			       traverseInfo.getAuxFESpace(i, j),
			       assembleFlag,
			       assembleMatrix ? matrix : NULL,
			       NULL);

	} 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
816
817
818
	  assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
			      assembleFlag,
			      assembleMatrix ? matrix : NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
819
820
821
			      ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);	  

	  TEST_EXIT_DBG(matrix->getUsedSize() == componentSpaces[i]->getAdmin()->getUsedSize())
822
	    ("Assembled matrix has wrong dimension!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
825
826
827
	  TEST_EXIT_DBG(matrix->getNumCols() == componentSpaces[j]->getAdmin()->getUsedSize())
	    ("Assembled matrix has wrong dimension!\n");	  

	} else {
	  ERROR_EXIT("Not yet implemented!\n");
828
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
829

830
831
832
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->exitMatrix(matrix);	  
	
Thomas Witkowski's avatar
Thomas Witkowski committed
833
	assembledMatrix[i][j] = true;
834
835
      }

836
      // And now assemble boundary conditions on the vectors
Thomas Witkowski's avatar
Thomas Witkowski committed
837
838
      assembleBoundaryConditions(rhs->getDOFVector(i),
				 solution->getDOFVector(i),
839
840
				 componentMeshes[i],
				 assembleFlag);
841
    }
842
843

#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
844
845
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
846
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
847
848
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));
849
#endif
850
851
852
853
854
855
  }

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

856
857
858
859
860
861
862
    clock_t first = clock();

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

    int i;
Thomas Witkowski's avatar
Thomas Witkowski committed
863
    int size = static_cast<int>(fileWriters.size());
864
865
866
867
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 1)
#endif
    for (i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
868
      fileWriters[i]->writeFiles(adaptInfo, force);
869
    }
870
871
    
#ifdef _OPENMP
Thomas Witkowski's avatar
Thomas Witkowski committed
872
873
874
    INFO(info, 8)("writeFiles needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()),
		  omp_get_wtime() - wtime);
875
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
876
877
    INFO(info, 8)("writeFiles needed %.5f seconds\n",
		  TIME_USED(first, clock()));
878
#endif
879
880
  }

881
  void ProblemVec::interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct) 
882
883
884
  {
    FUNCNAME("ProblemVec::interpolInitialSolution()");

Thomas Witkowski's avatar
Thomas Witkowski committed
885
    solution->interpol(fct);
886
887
888
889
890
891
892
893
  }

  void ProblemVec::addMatrixOperator(Operator *op, 
				     int i, int j,
				     double *factor,
				     double *estFactor)
  {
    FUNCNAME("ProblemVec::addMatrixOperator()");
Thomas Witkowski's avatar
Thomas Witkowski committed
894
   
Thomas Witkowski's avatar
Thomas Witkowski committed
895
    if (!(*systemMatrix)[i][j]) {
896
      TEST_EXIT(i != j)("should have been created already\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
897
898
899
900
      (*systemMatrix)[i][j] = NEW DOFMatrix(componentSpaces[i], componentSpaces[j], "");
      (*systemMatrix)[i][j]->setCoupleMatrix(true);
      (*systemMatrix)[i][j]->getBoundaryManager()->
	setBoundaryConditionMap((*systemMatrix)[i][i]->getBoundaryManager()->
901
902
				getBoundaryConditionMap());
    }    
Thomas Witkowski's avatar
Thomas Witkowski committed
903
904
905
906
907
908
909
910
911
912
913
    (*systemMatrix)[i][j]->addOperator(op, factor, estFactor);

    traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces()); 

    for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
      if ((op->getAuxFESpaces())[k] != componentSpaces[i] ||
	  (op->getAuxFESpaces())[k] != componentSpaces[j]) {
	op->setNeedDualTraverse(true);
	break;
      }
    }    
914
915
916
917
918
919
920
921
  }

  void ProblemVec::addVectorOperator(Operator *op, int i,
				     double *factor,
				     double *estFactor)
  {
    FUNCNAME("ProblemVec::addVectorOperator()");

Thomas Witkowski's avatar
Thomas Witkowski committed
922
923
924
925
926
927
928
929
930
    rhs->getDOFVector(i)->addOperator(op, factor, estFactor);
    traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces()); 

    for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
      if ((op->getAuxFESpaces())[j] != componentSpaces[i]) {
	op->setNeedDualTraverse(true);
	break;
      }
    }    
931
932
933
934
935
936
937
938
939
  }

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

    DirichletBC *dirichlet = new DirichletBC(type, 
					     b, 
940
					     componentSpaces[system]);
941
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
942
943
      if (systemMatrix && (*systemMatrix)[system][i]) {
	(*systemMatrix)[system][i]->getBoundaryManager()->addBoundaryCondition(dirichlet);
944
945
      }
    }
946

Thomas Witkowski's avatar
Thomas Witkowski committed
947
948
    if (rhs)
      rhs->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
949

Thomas Witkowski's avatar
Thomas Witkowski committed
950
951
    if (solution)
      solution->getDOFVector(system)->getBoundaryManager()->addBoundaryCondition(dirichlet);
952
953
954
955
956
957
958
959
960
  }

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

    NeumannBC *neumann = 
      new NeumannBC(type, n, 
961
962
		    componentSpaces[row], 
		    componentSpaces[col]);
Thomas Witkowski's avatar
Thomas Witkowski committed
963
964
    if (rhs)
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
965
966
967
968
969
970
971
972
973
974
  }

  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, 
975
976
		  componentSpaces[row], 
		  componentSpaces[col]);
Thomas Witkowski's avatar
Thomas Witkowski committed
977
978
    if (rhs)
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
979

Thomas Witkowski's avatar
Thomas Witkowski committed
980
981
    if (systemMatrix && (*systemMatrix)[row][col]) {
      (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
982
983
984
985
986
987
988
    }
  }

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

989
    FiniteElemSpace *feSpace = componentSpaces[row];
990
991
992

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

Thomas Witkowski's avatar
Thomas Witkowski committed
993
994
    if (systemMatrix && (*systemMatrix)[row][col]) 
      (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);
995

Thomas Witkowski's avatar
Thomas Witkowski committed
996
997
    if (rhs) 
      rhs->getDOFVector(row)->getBoundaryManager()->
998
999
1000
	addBoundaryCondition(periodic);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1001
1002
  void ProblemVec::assembleOnOneMesh(FiniteElemSpace *feSpace, 
				     Flag assembleFlag,
Thomas Witkowski's avatar
Thomas Witkowski committed
1003
1004