ProblemVec.cc 49.3 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <sstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
2
#include <boost/lexical_cast.hpp>
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include "ProblemVec.h"
#include "RecoveryEstimator.h"
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "FileWriter.h"
#include "CoarseningManager.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
17
#include "DualTraverse.h"
18
19
20
21
22
23
#include "Mesh.h"
#include "OEMSolver.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
#include "Flag.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "TraverseParallel.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include "VtkWriter.h"
27
#include "ValueReader.h"
28
#include "ProblemVecDbg.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
29
#include "Debug.h"
30
31
32

namespace AMDiS {

33
34
  using boost::lexical_cast;

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

41
    // === create meshes ===
Thomas Witkowski's avatar
Thomas Witkowski committed
42
    if (meshes.size() != 0) { 
43
44
45
      WARNING("meshes already created\n");
    } else {
      if (initFlag.isSet(CREATE_MESH) || 
Thomas Witkowski's avatar
Thomas Witkowski committed
46
	  (!adoptFlag.isSet(INIT_MESH) &&
Thomas Witkowski's avatar
Thomas Witkowski committed
47
	   (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
48
	createMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
49

50
      if (adoptProblem && 
Thomas Witkowski's avatar
Thomas Witkowski committed
51
	  (adoptFlag.isSet(INIT_MESH) ||
52
53
	   adoptFlag.isSet(INIT_SYSTEM) ||
	   adoptFlag.isSet(INIT_FE_SPACE))) {
Thomas Witkowski's avatar
Thomas Witkowski committed
54
	meshes = adoptProblem->getMeshes();
Thomas Witkowski's avatar
Thomas Witkowski committed
55
	componentMeshes = adoptProblem->componentMeshes;
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
	refinementManager = adoptProblem->refinementManager;
	coarseningManager = adoptProblem->coarseningManager;
58
59
60

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

69
70
71
      }
    }

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

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

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

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

99
100
101
      }
    }

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

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

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

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

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

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

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

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

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

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

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

	// If AMDiS is compied for parallel computations, the global refinements are
	// ignored here. Later, each rank will add the global refinements to its 
	// private mesh.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
194
	GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d", 
195
		      &globalRefinements);
196
#endif
197

198
	// Initialize the meshes if there is no serialization file.
199
200
201
202
203
	for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
	  bool initMesh = initFlag.isSet(INIT_MESH) ||
	    (adoptProblem && adoptFlag.isSet(INIT_MESH));

	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
Thomas Witkowski's avatar
Thomas Witkowski committed
204
	    meshes[i]->initialize();	    
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); 
Thomas Witkowski's avatar
nix.    
Thomas Witkowski committed
210
	if (valueFilename.length()) {
211
212
213
214
215
216
217
218
	  ValueReader::readValue(valueFilename,
				 meshes[0],
				 solution->getDOFVector(0),
				 meshes[0]->getMacroFileInfo());
	  meshes[0]->clearMacroFileInfo();
	}

	// === do global refinements ===
219
220

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

    doOtherStuff();
  }

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

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

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

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

276
  void ProblemVec::createFESpace(DOFAdmin *admin)
277
278
279
  {
    FUNCNAME("ProblemVec::createFESpace()");

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

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

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

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

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

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

    for (int i = 0; i < nComponents; i++) {
309
      for (int j = 0; j < nComponents; j++)
310
	traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
311
      
312
      traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
313
314
315
    }

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

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

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

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

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

341
      std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]";
Thomas Witkowski's avatar
Thomas Witkowski committed
342
      rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName));
343
344

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

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

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

367

368
369
370
371
372
  void ProblemVec::createEstimator()
  {
    FUNCNAME("ProblemVec::createEstimator()");

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

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

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


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

405

406
407
408
409
  void ProblemVec::createMarker()
  {
    FUNCNAME("ProblemVec::createMarker()");

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    int nMarkersCreated = 0;
411

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

416
      if (marker[i]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
420
	nMarkersCreated++;

	// If there is more than one marker, and all components are defined
	// on the same mesh, the maximum marking has to be enabled.
421
 	if (nMarkersCreated > 1 && nMeshes == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
422
 	  marker[i]->setMaximumMarking(true);
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
					    solutionList));
    }

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

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

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

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

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

  void ProblemVec::doOtherStuff()
  {
  }

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

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

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

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

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
#endif

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

517

518
519
520
521
522
523
  void ProblemVec::estimate(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::estimate()");

    clock_t first = clock();

524
525
526
527
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

528
    if (computeExactError) {
Thomas Witkowski's avatar
Thomas Witkowski committed
529
      computeError(adaptInfo);
530
531
    } else {
      for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
532
	Estimator *scalEstimator = estimator[i];
533
534
	
	if (scalEstimator) {
535
536
	  traverseInfo.updateStatus();
	  scalEstimator->setTraverseInfo(traverseInfo);
537
	  scalEstimator->estimate(adaptInfo->getTimestep());
538

539
540
541
542
543
	  adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
	  adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
	  adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
	  adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
	} else {
544
	  WARNING("No estimator for component %d\n" , i);
545
	}
546
547
548
      }
    }

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

558

559
560
561
562
563
564
565
566
  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();

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
581

582
583
584
585
  Flag ProblemVec::refineMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::refineMesh()");

Thomas Witkowski's avatar
Thomas Witkowski committed
586
    int nMeshes = static_cast<int>(meshes.size());
587
    Flag refineFlag = 0;
588
589
590
591
    for (int i = 0; i < nMeshes; i++)
      if (adaptInfo->isRefinementAllowed(i))
	refineFlag |= refinementManager->refineMesh(meshes[i]);

592
593
594
    return refineFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
595

596
597
598
599
  Flag ProblemVec::coarsenMesh(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("ProblemVec::coarsenMesh()");

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

    return coarsenFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
609

610
611
612
613
  Flag ProblemVec::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ProblemVec::oneIteration()");

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
630

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

636
637
638
639
    if (dualMeshTraverseRequired()) {
      dualAssemble(adaptInfo, flag);
      return;
    }
640

641
    //    printOpenmpTraverseInfo(this, true);
642
643
    //    std::cout << "ElInfo = " << ElInfo::subElemMatrices.size() << std::endl;

Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
646
647
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


648
    clock_t first = clock();
649
650
651
652
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

653
654
655

    Flag assembleFlag = 
      flag | 
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
658
659
660
661
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
662
663
      Mesh::FILL_NEIGH;

Thomas Witkowski's avatar
Thomas Witkowski committed
664
    if (useGetBound)
665
      assembleFlag |= Mesh::FILL_BOUND;
Thomas Witkowski's avatar
Thomas Witkowski committed
666
667

  
Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
    traverseInfo.updateStatus();

670
671
672
    // Used to calculate the overall number of non zero entries.
    int nnz = 0;

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

675
676
677
678
679
680
      MSG("%d DOFs for %s\n", 
	  componentSpaces[i]->getAdmin()->getUsedSize(), 
	  componentSpaces[i]->getName().c_str());

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

681
      for (int j = 0; j < nComponents; j++) {
682

683
684
	if (writeAsmInfo)
	  std::cout << "-------" << i << " " << j << "----------------" << std::endl;
685

686
687
688
689
	// 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).
690
	DOFMatrix *matrix = (asmMatrix ? (*systemMatrix)[i][j] : NULL);
691

692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
	if (writeAsmInfo && matrix) {
	  for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
	       it != matrix->getOperatorsEnd(); ++it) {
	    Assembler *assembler = (*it)->getAssembler();
	    if (assembler->getZeroOrderAssembler()) 
	      std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
	    if (assembler->getFirstOrderAssembler(GRD_PSI)) 
	      std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
	    if (assembler->getFirstOrderAssembler(GRD_PHI)) 
	      std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
	    if (assembler->getSecondOrderAssembler()) 
	      std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
	  }
	}

707
708
709
	if (matrix) 
	  matrix->calculateNnz();
	
710
711
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
712
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
713
	  assembleMatrix = false;
714
715
716
717
718
719
720
	} else if (matrix) {
	  matrix->getBaseMatrix().
	    change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
		       componentSpaces[j]->getAdmin()->getUsedSize());

	  set_to_zero(matrix->getBaseMatrix());
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722

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

726
727
728
	// 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.
729
730
731
732
	if (!assembleMatrix && i != j) {
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();

733
	  continue;
734
	}
735

736
737
738
	if (assembleMatrix && matrix->getBoundaryManager())
	  matrix->getBoundaryManager()->initMatrix(matrix);

739
740
741
742
743
744
745
746
	// The simplest case: either the right hand side has no operaters, no aux
	// fe spaces, or all aux fe spaces are equal to the row and col fe space.
	assembleOnOneMesh(componentSpaces[i],
			  assembleFlag,
			  assembleMatrix ? matrix : NULL,
			  ((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
      
	assembledMatrix[i][j] = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
747

Thomas Witkowski's avatar
Thomas Witkowski committed
748
	
749
750
	if (assembleMatrix)
	  matrix->finishInsertion();
Thomas Witkowski's avatar
Thomas Witkowski committed
751

752
753
754
755
 	if (assembleMatrix && matrix->getBoundaryManager())
 	  matrix->getBoundaryManager()->exitMatrix(matrix);	
	
	if (matrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
756
	  nnz += matrix->getBaseMatrix().nnz();
757
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
758

759
760
      // And now assemble boundary conditions on the vectors
      assembleBoundaryConditions(rhs->getDOFVector(i),
Thomas Witkowski's avatar
Thomas Witkowski committed
761
762
763
 				 solution->getDOFVector(i),
 				 componentMeshes[i],
 				 assembleFlag);     
764
    }
765

766
767
768
    if (asmMatrix) {
      solverMatrix.setMatrix(*systemMatrix);      
      createPrecon();
Thomas Witkowski's avatar
Thomas Witkowski committed
769

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

773
774
775
776
777
778
779
780
#ifdef _OPENMP
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
		  TIME_USED(first, clock()), omp_get_wtime() - wtime);
#else
    INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 
		  TIME_USED(first, clock()));    
#endif    
  }
781

782
783
784
785
786
787
788
789
790
791
792
793
794
795
796

  bool ProblemVec::dualMeshTraverseRequired()
  {
    FUNCNAME("ProblemVec::dualMeshTraverseRequired()");

    TEST_EXIT(meshes.size() <= 2)("More than two mneshes are not supported yet!\n");    

    return (meshes.size() == 2);
  }
  

  void ProblemVec::dualAssemble(AdaptInfo *adaptInfo, Flag flag)
  {
    FUNCNAME("ProblemVec::dualAssemble()");

Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
799
800
    for (unsigned int i = 0; i < meshes.size(); i++)
      meshes[i]->dofCompress();


801
802
803
804
805
    clock_t first = clock();
#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
806
    
807
808
809
810
    Flag assembleFlag = 
      flag | 
      (*systemMatrix)[0][0]->getAssembleFlag() | 
      rhs->getDOFVector(0)->getAssembleFlag()  |
811
812
813
814
      Mesh::CALL_LEAF_EL                        | 
      Mesh::FILL_COORDS                         |
      Mesh::FILL_DET                            |
      Mesh::FILL_GRD_LAMBDA |
815
816
817
818
      Mesh::FILL_NEIGH;

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

820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
    traverseInfo.updateStatus();

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

    for (int i = 0; i < nComponents; i++) {

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

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

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

	if (writeAsmInfo)
	  std::cout << "-------" << i << " " << j << "----------------" << std::endl;

	// Only if this variable is true, the current matrix will be assembled.	
	bool assembleMatrix = true;
	// The DOFMatrix which should be assembled (or not, if assembleMatrix
	// will be set to false).
	DOFMatrix *matrix = (*systemMatrix)[i][j];

	if (writeAsmInfo && matrix) {
	  for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
	       it != matrix->getOperatorsEnd(); ++it) {
	    Assembler *assembler = (*it)->getAssembler();
	    if (assembler->getZeroOrderAssembler()) 
	      std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
	    if (assembler->getFirstOrderAssembler(GRD_PSI)) 
	      std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
	    if (assembler->getFirstOrderAssembler(GRD_PHI)) 
	      std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
	    if (assembler->getSecondOrderAssembler()) 
	      std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
Thomas Witkowski's avatar
Thomas Witkowski committed
856
	  }
857
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
858

859
860
861
862
863
864
865
866
867
868
869
	if (matrix) 
	  matrix->calculateNnz();
	
	// If the matrix was assembled before and it is marked to be assembled
	// only once, it will not be assembled.
	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
	  assembleMatrix = false;
	} else if (matrix) {
	  matrix->getBaseMatrix().
	    change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
		       componentSpaces[j]->getAdmin()->getUsedSize());
870

871
872
	  set_to_zero(matrix->getBaseMatrix());
	  matrix->startInsertion(matrix->getNnz());
873

874
875
876
	  if (matrix->getBoundaryManager())
	    matrix->getBoundaryManager()->initMatrix(matrix);
	}
877

878
879
880
	// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
	if (!matrix || !assembleMatrix)
	  assembleMatrix = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
881

882
883
884
885
886
887
888
889
	// If the matrix should not be assembled, the rhs vector has to be considered.
	// This will be only done, if i == j. So, if both is not true, we can jump
	// to the next matrix.
	if (!assembleMatrix && i != j)
	  if (matrix)
	    nnz += matrix->getBaseMatrix().nnz();
      }
    }
890

891
    TEST_EXIT(meshes.size() == 2)("There is something wrong!\n");
892

893
894
895
896
897
898
899
900
901
902
903
904
    const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
    BoundaryType *bound = 
      useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;

    DualTraverse dualTraverse;
    DualElInfo dualElInfo;
    int oldElIndex0 = -1;
    int oldElIndex1 = -1;
    dualTraverse.setFillSubElemMat(true, basisFcts);
    bool cont = dualTraverse.traverseFirst(meshes[0], meshes[1], -1, -1, 
					   assembleFlag, assembleFlag, dualElInfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
905
    while (cont) {     
906
907
908
      bool newEl0 = (dualElInfo.rowElInfo->getElement()->getIndex() != oldElIndex0);
      bool newEl1 = (dualElInfo.colElInfo->getElement()->getIndex() != oldElIndex1);
      oldElIndex0 = dualElInfo.rowElInfo->getElement()->getIndex();
Thomas Witkowski's avatar
Thomas Witkowski committed
909
      oldElIndex1 = dualElInfo.colElInfo->getElement()->getIndex();      
910
911
912
913
914
915
916
917
918

      for (int i = 0; i < nComponents; i++) {
	for (int j = 0; j < nComponents; j++) {
	  DOFMatrix *matrix = (*systemMatrix)[i][j];

	  if (!matrix)
	    continue;

	  if (traverseInfo.eqSpaces(i, j)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
919

920
921
922
923
924
	    ElInfo *elInfo = NULL;	    
	    if (componentMeshes[i] == meshes[0] && newEl0)
	      elInfo = dualElInfo.rowElInfo;	    
	    if (componentMeshes[i] == meshes[1] && newEl1)
	      elInfo = dualElInfo.colElInfo;
Thomas Witkowski's avatar
Thomas Witkowski committed
925
926
	  

927
928
929
930
	    if (elInfo != NULL) {
	      if (useGetBound)
		basisFcts->getBound(elInfo, bound);
	      
Thomas Witkowski's avatar
Thomas Witkowski committed
931
	      if (matrix)
932
		matrix->assemble(1.0, elInfo, bound);
Thomas Witkowski's avatar
Thomas Witkowski committed
933
	      
934
935
936
937
	      if (i == j)
		rhs->getDOFVector(i)->assemble(1.0, elInfo, bound);
	    }

Thomas Witkowski's avatar
Thomas Witkowski committed
938
939
940
941
942
943
944
	    ElInfo *mainElInfo, *auxElInfo;
	    if (traverseInfo.getRowFeSpace(i)->getMesh() == meshes[0]) {
	      mainElInfo = dualElInfo.rowElInfo;
	      auxElInfo = dualElInfo.colElInfo;
	    } else {
	      mainElInfo = dualElInfo.colElInfo;
	      auxElInfo = dualElInfo.rowElInfo;
945
946
	    }

Thomas Witkowski's avatar
Thomas Witkowski committed
947
948
949
950
951
952
	    if (useGetBound && mainElInfo != elInfo)
	      basisFcts->getBound(mainElInfo, bound);
	    
 	    if (traverseInfo.difAuxSpace(i) && i == j)
 	      rhs->getDOFVector(i)->assemble2(1.0, mainElInfo, auxElInfo,
 					      dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);	    
953

Thomas Witkowski's avatar
Thomas Witkowski committed
954
955
956
957
958
959
 	    if (traverseInfo.difAuxSpace(i, j) && matrix)
 	      matrix->assemble2(1.0, mainElInfo, auxElInfo,
 				dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);

	    if (matrix && matrix->getBoundaryManager())
	      matrix->getBoundaryManager()->fillBoundaryConditions(mainElInfo, matrix);	
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977

	  } else {
	    TEST_EXIT_DBG(traverseInfo.getStatus(i, j) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX)
	      ("Not yet supported!\n");

	    ElInfo *rowElInfo, *colElInfo;
	    if (componentMeshes[i] == meshes[0]) {
	      rowElInfo = dualElInfo.rowElInfo;
	      colElInfo = dualElInfo.colElInfo;
	    } else {
	      rowElInfo = dualElInfo.colElInfo;
	      colElInfo = dualElInfo.rowElInfo;
	    }
	    
	    if (useGetBound)
	      basisFcts->getBound(rowElInfo, bound);	      
	    
	    if (matrix) {
Thomas Witkowski's avatar
Thomas Witkowski committed
978
979
 	      matrix->assemble(1.0, rowElInfo, colElInfo, 
 			       dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
980
981
982
983
984
	      
	      if (matrix->getBoundaryManager())
		matrix->getBoundaryManager()->fillBoundaryConditions(rowElInfo, matrix);
	    }
	    
Thomas Witkowski's avatar
Thomas Witkowski committed
985
986
 	    if (i == j) 
 	      rhs->getDOFVector(i)->assemble(1.0, rowElInfo, bound);    
987
	  }