ExtendedProblemStat.h 17.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_EXTENDED_PROBLEM_STAT_H
#define EXTENSIONS_EXTENDED_PROBLEM_STAT_H
21
22

#include "AMDiS.h"
23
#include "SingularDirichletBC2.h"
24
25
26
27
28
29
#if defined NONLIN_PROBLEM
  #include "nonlin/ProblemNonLin.h"
#endif

using namespace AMDiS;

Praetorius, Simon's avatar
Praetorius, Simon committed
30
31
32
const Flag INIT_EXACT_SOLUTION = 0X10000L;
const Flag UPDATE_PERIODIC_BC  = 0X20000L;
const Flag UPDATE_DIRICHLET_BC = 0X40000L;
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

#if defined NONLIN_PROBLEM
typedef ProblemNonLin ProblemStat_;
#else
typedef ProblemStat ProblemStat_;
#endif
class ExtendedProblemStat : public ProblemStat_
{
public:
  
  ExtendedProblemStat(std::string nameStr, ProblemIterationInterface *problemIteration = NULL)
  :
  #if defined NONLIN_PROBLEM
  ProblemStat_(nameStr)
  #else
  ProblemStat_(nameStr, problemIteration)
  #endif
  , oldMeshChangeIdx(0)
  {
    exactSolutions.resize(nComponents);
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = NULL;
  }
56
57
58
59
60
61
62
63
64
65
66
67
  
  template<typename SubProblemType>
  ExtendedProblemStat(std::string nameStr, SubProblemType *subProblem)
  :
  ProblemStat_(nameStr, subProblem)
  , oldMeshChangeIdx(0)
  {
    exactSolutions.resize(nComponents);
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = NULL;
  }
  
68
  
69
  //////////////////////////////////////////////////////////////////////////////
70
71
  ~ExtendedProblemStat()
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
72
    for (int i = 0; i < nComponents; ++i) {
73
74
75
76
      if (exactSolutions[i]) {
	delete exactSolutions[i];
	exactSolutions[i] = NULL;
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
77
78
79
80
    }
    
    for (size_t i = 0; i < DirichletBcDataList.size(); i++)
      delete DirichletBcDataList[i];
81
82
    
  }
83

84
  
85
  //////////////////////////////////////////////////////////////////////////////
86
87
88
  void initialize(Flag initFlag,
		  ProblemStatSeq *adoptProblem = NULL,
		  Flag adoptFlag = INIT_NOTHING)
Praetorius, Simon's avatar
Praetorius, Simon committed
89
  {    
90
    ProblemStat_::initialize(initFlag, adoptProblem, adoptFlag);
91
92
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution");
93
  }
94
  
95
96
97
98
99
100

  inline DOFVector<double>* getRhsVector(int i = 0)
  {
    return rhs->getDOFVector(i);
  }

101
  
102
  //////////////////////////////////////////////////////////////////////////////
103
104
105
106
107
108
109
  void setExactSolution(DOFVector<double> *dof, int component) 
  {
    TEST_EXIT(exactSolutions[component] != NULL)
      ("You have to initialize the exactSolutions vector! Set the initFlag INIT_EXACT_SOLUTION for the problem.\n");
    exactSolutions[component]->copy(*dof);
  }

110
  
111
  //////////////////////////////////////////////////////////////////////////////
112
113
114
115
116
117
118
  inline DOFVector<double> *getExactSolution(int i = 0) 
  {
    TEST_EXIT(exactSolutions[i] != NULL)
      ("You have to initialize the exactSolutions vector! Set the initFlag INIT_EXACT_SOLUTION for the problem.\n");
    return exactSolutions[i];
  }

119
  
120
  //////////////////////////////////////////////////////////////////////////////
121
122
  void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
			  bool asmMatrix, bool asmVector)
123
  {  
124
125
126
127
128
129
130
131
132
    ProblemStat_::buildAfterCoarsen(adaptInfo, flag, asmMatrix, asmVector);

    // update periodic data
    if (oldMeshChangeIdx != getMesh()->getChangeIndex()
      || flag.isSet(UPDATE_PERIODIC_BC)
      || (PeriodicBcDataList.size() > 0 && manualPeriodicBC.size() == 0)) {
      manualPeriodicBC.clear();
      std::vector<PeriodicBcData>::iterator it;
      for (it = PeriodicBcDataList.begin(); it != PeriodicBcDataList.end(); it++)
133
	it->addToList(getFeSpace(it->row), manualPeriodicBC);
134
    }
135
136
137
138

    // apply periodic boundary conditions
    for (size_t k = 0; k < manualPeriodicBC.size(); k++)
      applyPeriodicBC(manualPeriodicBC[k], asmMatrix, asmVector);
139
140
141
142
143
144
    
    // update dirichlet data
    if (oldMeshChangeIdx != getMesh()->getChangeIndex()
      || flag.isSet(UPDATE_DIRICHLET_BC)
      || (DirichletBcDataList.size() > 0 && singularDirichletBC.size() == 0)) {
      singularDirichletBC.clear();
145
      std::vector<DirichletBcDataBase*>::iterator it;
146
      for (it = DirichletBcDataList.begin(); it != DirichletBcDataList.end(); it++) {
147
	(*it)->addToList(getFeSpace((*it)->row), singularDirichletBC);
148
149
150
151
152
153
      }
    }
    
    // apply dirichlet boundary conditions
    for (size_t k = 0; k < singularDirichletBC.size(); k++)
      applyDirichletBC(singularDirichletBC[k], asmMatrix, asmVector);
Praetorius, Simon's avatar
Praetorius, Simon committed
154
    
155
156
157
158
159
160
    // update solverMatrix
    if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
      solverMatrix.setMatrix(*getSystemMatrix());
    }
  }

161
  
162
  //////////////////////////////////////////////////////////////////////////////
163
164
165
  void solve(AdaptInfo *adaptInfo,
	      bool createMatrixData = true,
	      bool storeMatrixData = false)
166
  {
167
168
169
170
171
    ProblemStat_::solve(adaptInfo, createMatrixData, storeMatrixData);

    oldMeshChangeIdx = getMesh()->getChangeIndex();
  }

172
  
173
  /// Add arbitrary boundary condition to system
174
  //////////////////////////////////////////////////////////////////////////////
175
176
177
178
179
180
181
182
183
184
  void addBoundaryCondition(BoundaryCondition *bc, int row, int col)
  {
    boundaryConditionSet = true;

    if (systemMatrix && (*systemMatrix)[row][col])
      (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(bc);
    if (rhs)
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(bc);
  }

185
  
186
187
188
  //============================================================================
  
  
189
  //////////////////////////////////////////////////////////////////////////////
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
  template<typename ValueContainer>
  void addDirichletBC(BoundaryType type, int row, int col,
		      ValueContainer *values)
  {
    BoundaryTypeContainer *bound = new BoundaryTypeContainer(type);
    DirichletBcDataList.push_back(
      new DirichletBcData<BoundaryTypeContainer, ValueContainer>(
	row, col, *bound, *values));
  }
  
  // general templatized version of addDirichletBC
  
  template<typename PosType, typename ValueContainer>
  void addDirichletBC(PosType *pos, int row, int col,
		      ValueContainer *values)
  {
    DirichletBcDataList.push_back(
      new DirichletBcData<PosType, ValueContainer>(row, col, *pos, *values));
  }
  
210
211
212
213
214
215
216
217
  /**
    * Dirichlet boundary condition at DOF-Index of vertex near to coords 'pos'.
    * Value defined by AbstractFunction, that is evaluated at 'pos'
    **/
  template<typename ValueContainer>
  void addSingularDirichletBC(WorldVector<double> &pos, int row, int col,
			      ValueContainer &values)
  {
218
219
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&pos, row, col, &values);
220
221
  }

222
  
223
224
225
226
  /**
    * Dirichlet boundary condition at DOF-Index 'idx'. Value defined by double.
    **/
  template<typename ValueContainer>
227
228
  void addDirichletBC(DofIndex* idx, int row, int col,
		      ValueContainer *values)
229
  {
230
231
232
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (MPI::COMM_WORLD.Get_rank() == 0)
    #endif
233
234
235
236
237
238
239
240
241
242
    DirichletBcDataList.push_back(
      new DirichletBcData<DofIndex, ValueContainer>(row, col, *idx, *values));
  }
  
  template<typename ValueContainer>
  void addSingularDirichletBC(DegreeOfFreedom idx, int row, int col,
			      ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    DofIndex* dofIndex = new DofIndex(idx);
Praetorius, Simon's avatar
Praetorius, Simon committed
243
    addDirichletBC(dofIndex, row, col, &values);
244
245
  }

246
  
247
248
249
250
251
252
253
254
255
256
257
  /**
    * set dirichlet boundary condition on implicitly defined boundary by zero levelset
    * of signed distance function.
    * An algebraic equation is forced in the region with dist<0
    * 
    **/
  template<typename ValueContainer>
  void addImplicitDirichletBC(DOFVector<double> &signedDist,
			      int row, int col,
			      ValueContainer &values)
  {
258
259
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&signedDist, row, col, &values);
260
261
  }
  
262
  
263
  ///
264
265
266
267
268
  template<typename ValueContainer>
  void addImplicitDirichletBC(AbstractFunction<double, WorldVector<double> > &signedDist,
			      int row, int col,
			      ValueContainer &values)
  {
269
270
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&signedDist, row, col, &values);
271
272
  }

273
  
274
275
276
277
278
279
280
281
282
283
284
  ///  
  template<typename ValueContainer>
  void addDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
		      int row, int col,
		      ValueContainer *values)
  {
    DirichletBcDataList.push_back(
      new DirichletBcData<AbstractFunction<bool, WorldVector<double> >, ValueContainer>(
	row, col, nr, *meshIndicator, *values));
  }
    
285
286
287
288
289
  template<typename ValueContainer>
  void addManualDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    int row, int col,
			    ValueContainer &values)
  {
290
291
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(nr, meshIndicator, row, col, &values);
292
293
294
  }
  

295
  ///
296
297
298
299
300
  template<typename ValueContainer>
  void addManualDirichletBC(AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    int row, int col,
			    ValueContainer &values)
  {
301
302
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(meshIndicator, row, col, &values);
303
304
305
  }
  
  
306
307
308
  // ===========================================================================
  
  
309
  //////////////////////////////////////////////////////////////////////////////
310
311
312
313
314
315
316
  void addManualPeriodicBC(int row,
			    BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
  {
    PeriodicBcDataList.push_back(PeriodicBcData(row, nr, meshIndicator, periodicMap));
  }

317
  
318
  //////////////////////////////////////////////////////////////////////////////
319
320
321
322
323
324
325
  void addManualPeriodicBC(int row,
			    AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
  {
    PeriodicBcDataList.push_back(PeriodicBcData(row, meshIndicator, periodicMap));
  }

326
  
327
328
329
330
331
332
333
334
335
336
  /// write Systemmatrix to file in matlab-format
  void writeMatrix(std::string filename)
  {
    mtl::io::matrix_market_ostream out(filename);
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
    solverMatrix.setMatrix(*getSystemMatrix());
    out << solverMatrix.getMatrix();
    out.close();
  }
  
337
  
338
339
340
341
342
343
344
345
346
347
  /// Returns the name of the problem
  inline string getComponentName(int comp = 0)
  {
    TEST_EXIT(comp < static_cast<int>(componentNames.size()) && comp >= 0)
      ("invalid component number\n");
    return componentNames[comp];
  }
    
protected:
  // traverse matrix rows and set unity row where dirichlet condition shall be applied.
Praetorius, Simon's avatar
Praetorius, Simon committed
348
349
350
  void applyDirichletBC(size_t row_, size_t col_, 
			DegreeOfFreedom idx, double value, 
			bool asmMatrix = true, bool asmVector = true)
351
352
353
354
355
356
357
358
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

    Matrix::size_type idx_= idx;
    bool value1set = false;

    if (asmMatrix) {
359
360
      typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
      typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
361
362
363
364
365
366
367
368
369
370
      
      // if matrix-block for the identity row is NULL, create new DOFMatrix
      if (getSystemMatrix(row_, col_) == NULL) {
	TEST_EXIT(row_ != col_)("should have been created already\n");
	(*systemMatrix)[row_][col_] = new DOFMatrix(getFeSpace(row_), getFeSpace(col_), "");
	(*systemMatrix)[row_][col_]->setCoupleMatrix(true);
	(*systemMatrix)[row_][col_]->getBoundaryManager()->
	  setBoundaryConditionMap((*systemMatrix)[row_][row_]->getBoundaryManager()->
				  getBoundaryConditionMap());
      }
371

372
373
      size_t numCols = static_cast<size_t>(getNumComponents());
      for (size_t col = 0; col < numCols; col++) {
374
375
376
377
378
379
	if (getSystemMatrix(row_, col) == NULL)
	  continue;

	// set Dirichlet-row in matrix
	Matrix &m = getSystemMatrix(row_, col)->getBaseMatrix();

380
381
382
	mtl::traits::row<Matrix>::type r(m);
	mtl::traits::col<Matrix>::type c(m);
	mtl::traits::value<Matrix>::type v(m);
383
384
385

	c_type cursor(begin<tag::row>(m)+idx_);
	for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
386
	  value1set = value1set || (r(*icursor) == c(*icursor) && col == col_);
387
388
389
390
391
392
	  v(*icursor, (r(*icursor) == c(*icursor) && col == col_ ? 1.0 : 0.0));
	}
      }

      if (!value1set) {
	matrix::inserter<Matrix, update_plus<double> > ins(getSystemMatrix(row_, col_)->getBaseMatrix());
393
	ins[idx_][idx_] << 1.0;
394
395
396
397
398
      }
    }
    if (asmVector)
      (*getRhsVector(row_))[idx] = value; // set Dirichlet-Value at rhs-vector
      
399
    (*solution->getDOFVector(col_))[idx] = value; // set Dirichlet-value for solution component
400
401
  }
  
402
  
403
404
405
406
  void applyDirichletBC(SingularDirichletBC &sbc, bool asmMatrix = true, bool asmVector = true)
  {
    applyDirichletBC(sbc.row, sbc.col, sbc.idx, sbc.value, asmMatrix, asmVector);
  }
407
  
408
409
410
411
412
413

  void applyPeriodicBC(size_t row, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &indices, bool asmMatrix = true, bool asmVector = true)
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

414
415
    typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
416

417
418
    size_t numCols = static_cast<size_t>(getNumComponents());
    for (size_t col = 0; col < numCols; col++) {
419
420
421
422
423
      TEST_EXIT(getSystemMatrix(row, col) != NULL)
	("SystemMatrix block (%d,%d) must not be NULL! Insert a Simple_ZOT(0.0) as Workaround.\n",row,col);

      Matrix &m = getSystemMatrix(row, col)->getBaseMatrix();

424
425
426
      mtl::traits::row<Matrix>::type r(m);
      mtl::traits::col<Matrix>::type c(m);
      mtl::traits::value<Matrix>::type v(m);
427
428
429
430
431
432
      
      std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values;
      row_values.resize(indices.size());

      // erase the rows for all first indices
      for (size_t i = 0; i < indices.size(); i++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
433
434
435
	if (indices[i].first == indices[i].second)
	  continue;
	
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
	if (asmMatrix) {
	  c_type cursor(begin<tag::row>(m)+indices[i].first);
	  for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
	    row_values[i].push_back(std::make_pair(c(*icursor), v(*icursor)));
	    v(*icursor, 0.0);
	  }
	}
	if (asmVector) {
	  (*(getRhsVector(row)))[indices[i].second] += (*(getRhsVector(row)))[indices[i].first];
	  (*(getRhsVector(row)))[indices[i].first] = 0.0;
	}
      }

      // add periodic associations of first and second indices, but only in the diagonal blocks
      if (asmMatrix) {
	matrix::inserter<Matrix, update_plus<double> > ins(m);
	if (row == col) {
	  for (size_t i = 0; i < indices.size(); i++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
454
455
	    if (indices[i].first == indices[i].second)
	      continue;
456
457
458
459
460
	    ins[indices[i].first][indices[i].first] << 1.0;
	    ins[indices[i].first][indices[i].second] << -1.0;
	  }
	}
	for (size_t i = 0; i < indices.size(); i++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
461
462
	  if (indices[i].first == indices[i].second)
	    continue;
463
464
465
466
467
468
469
	  for (size_t j = 0; j < row_values[i].size(); j++) {
	    ins[indices[i].second][row_values[i][j].first] << row_values[i][j].second;
	  }
	}
      }
    }
  }
470
  
471
472
473
474
475
476

  void applyPeriodicBC(ManualPeriodicBC &pbc, bool asmMatrix = true, bool asmVector = true)
  {
    applyPeriodicBC(pbc.row, pbc.indices, asmMatrix, asmVector);
  }

477
  
478
479
480
481
482
483
484
485
  void applyAlgebraicEquation(size_t row,
			      std::vector<DegreeOfFreedom> &row_idx,
			      std::vector<std::vector<double> > &coefficients,
			      std::vector<double> &rhs)
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

486
487
    typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
488

489
    size_t numCols = static_cast<size_t>(getNumComponents());
490
491
    TEST_EXIT(row_idx.size() == coefficients.size() && row_idx.size() == rhs.size() && rhs.size()>0)
      ("rhs_idx, coefficients and rhs must have the same size and size >! 0\n");
492
    TEST_EXIT(coefficients[0].size() == numCols)
493
494
      ("You have to give coefficients for all variables\n");

495
    for (size_t col = 0; col < numCols; col++) {
496
497
498
499
500
      TEST_EXIT(getSystemMatrix(row, col) != NULL)
	("SystemMatrix block (%d,%d) must not be NULL! Insert a Simple_ZOT(0.0) as Workaround.\n",row,col);

      Matrix &m = getSystemMatrix(row, col)->getBaseMatrix();

501
502
503
      mtl::traits::row<Matrix>::type r(m);
      mtl::traits::col<Matrix>::type c(m);
      mtl::traits::value<Matrix>::type v(m);
504
505
506
507
508
509
510
511
512
513
514

      // erase the rows for all row-indices and set rhs values
      for (size_t i = 0; i < coefficients.size(); i++) {
	c_type cursor(begin<tag::row>(m)+row_idx[i]);
	for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
	  v(*icursor, (r(*icursor) == c(*icursor) ? coefficients[i][col] : 0.0));
	}
	(*(getRhsVector(row)))[row_idx[i]] = rhs[i];
      }
    }
  }
515
  
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539

  void applyAlgebraicEquation(size_t row,
			      DegreeOfFreedom &row_idx0,
			      std::vector<double> &coefficients0,
			      double rhs0)
  {
    std::vector<DegreeOfFreedom> row_idx; row_idx.push_back(row_idx0);
    std::vector<std::vector<double> > coefficients; coefficients.push_back(coefficients0);
    std::vector<double> rhs; rhs.push_back(rhs0);

    applyAlgebraicEquation(row, row_idx, coefficients, rhs);
  }

private:
  
  int oldMeshChangeIdx;
  std::vector<DOFVector<double>*> exactSolutions;

  // data for periodic boundary conditions
  std::vector<ManualPeriodicBC> manualPeriodicBC;
  std::vector<PeriodicBcData> PeriodicBcDataList;

  // data for dirichlet boundary conditions
  std::vector<SingularDirichletBC> singularDirichletBC;
540
  std::vector<DirichletBcDataBase*> DirichletBcDataList;
541
542
543
544

  std::map<const FiniteElemSpace*, bool> feSpaceVisited;
  
};
545
#endif // EXTENSIONS_EXTENDED_PROBLEM_STAT_H