ExtendedProblemStat.h 17 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
72
73
74
75
76
77
  ~ExtendedProblemStat()
  {
    for (int i = 0; i < nComponents; ++i)
      if (exactSolutions[i]) {
	delete exactSolutions[i];
	exactSolutions[i] = NULL;
      }
    
  }
78

79
  
80
81
82
83
84
  void initialize(Flag initFlag,
		  ProblemStatSeq *adoptProblem = NULL,
		  Flag adoptFlag = INIT_NOTHING)
  {
    ProblemStat_::initialize(initFlag, adoptProblem, adoptFlag);
85
86
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution");
87
  }
88
  
89
90
91
92
93
94

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

95
  
96
97
98
99
100
101
102
103
  ///
  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);
  }

104
  
105
106
107
108
109
110
111
112
  ///
  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];
  }

113
  
114
115
  void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
			  bool asmMatrix, bool asmVector)
116
  {  
117
118
119
120
121
122
123
124
125
    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++)
126
	it->addToList(getFeSpace(it->row), manualPeriodicBC);
127
    }
128
129
130
131

    // apply periodic boundary conditions
    for (size_t k = 0; k < manualPeriodicBC.size(); k++)
      applyPeriodicBC(manualPeriodicBC[k], asmMatrix, asmVector);
132
133
134
135
136
137
    
    // update dirichlet data
    if (oldMeshChangeIdx != getMesh()->getChangeIndex()
      || flag.isSet(UPDATE_DIRICHLET_BC)
      || (DirichletBcDataList.size() > 0 && singularDirichletBC.size() == 0)) {
      singularDirichletBC.clear();
138
      std::vector<DirichletBcDataBase*>::iterator it;
139
      for (it = DirichletBcDataList.begin(); it != DirichletBcDataList.end(); it++) {
140
	(*it)->addToList(getFeSpace((*it)->row), singularDirichletBC);
141
142
143
144
145
146
147
148
149
150
151
152
153
      }
    }
    
    // apply dirichlet boundary conditions
    for (size_t k = 0; k < singularDirichletBC.size(); k++)
      applyDirichletBC(singularDirichletBC[k], asmMatrix, asmVector);

    // update solverMatrix
    if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
      solverMatrix.setMatrix(*getSystemMatrix());
    }
  }

154
  
155
156
157
  void solve(AdaptInfo *adaptInfo,
	      bool createMatrixData = true,
	      bool storeMatrixData = false)
158
  {
159
160
161
162
163
    ProblemStat_::solve(adaptInfo, createMatrixData, storeMatrixData);

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

164
  
165
166
167
168
169
170
171
172
173
174
175
  /// Add arbitrary boundary condition to system
  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);
  }

176
  
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
  //============================================================================
  
  
  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));
  }
  
200
201
202
203
204
205
206
207
  /**
    * 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)
  {
208
209
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&pos, row, col, &values);
210
211
  }

212
  
213
214
215
216
  /**
    * Dirichlet boundary condition at DOF-Index 'idx'. Value defined by double.
    **/
  template<typename ValueContainer>
217
218
  void addDirichletBC(DofIndex* idx, int row, int col,
		      ValueContainer *values)
219
  {
220
221
222
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (MPI::COMM_WORLD.Get_rank() == 0)
    #endif
223
224
225
226
227
228
229
230
231
232
233
    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);
    addDirichletBC(*dofIndex, row, col, values);
234
235
  }

236
  
237
238
239
240
241
242
243
244
245
246
247
  /**
    * 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)
  {
248
249
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&signedDist, row, col, &values);
250
251
  }
  
252
  
253
  ///
254
255
256
257
258
  template<typename ValueContainer>
  void addImplicitDirichletBC(AbstractFunction<double, WorldVector<double> > &signedDist,
			      int row, int col,
			      ValueContainer &values)
  {
259
260
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&signedDist, row, col, &values);
261
262
  }

263
  
264
265
266
267
268
269
270
271
272
273
274
  ///  
  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));
  }
    
275
276
277
278
279
  template<typename ValueContainer>
  void addManualDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    int row, int col,
			    ValueContainer &values)
  {
280
281
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(nr, meshIndicator, row, col, &values);
282
283
284
  }
  

285
  ///
286
287
288
289
290
  template<typename ValueContainer>
  void addManualDirichletBC(AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    int row, int col,
			    ValueContainer &values)
  {
291
292
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(meshIndicator, row, col, &values);
293
294
295
  }
  
  
296
297
298
  // ===========================================================================
  
  
299
300
301
302
303
304
305
  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));
  }

306
  
307
308
309
310
311
312
313
  void addManualPeriodicBC(int row,
			    AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
  {
    PeriodicBcDataList.push_back(PeriodicBcData(row, meshIndicator, periodicMap));
  }

314
  
315
316
317
318
319
320
321
322
323
324
  /// 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();
  }
  
325
  
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
  /// 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.
  void applyDirichletBC(size_t row_, size_t col_, DegreeOfFreedom idx, double value, bool asmMatrix = true, bool asmVector = true)
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

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

    if (asmMatrix) {
345
346
      typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
      typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
347
348
349
350
351
352
353
354
355
356
      
      // 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());
      }
357

358
359
      size_t numCols = static_cast<size_t>(getNumComponents());
      for (size_t col = 0; col < numCols; col++) {
360
361
362
363
364
365
	if (getSystemMatrix(row_, col) == NULL)
	  continue;

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

366
367
368
	mtl::traits::row<Matrix>::type r(m);
	mtl::traits::col<Matrix>::type c(m);
	mtl::traits::value<Matrix>::type v(m);
369
370
371

	c_type cursor(begin<tag::row>(m)+idx_);
	for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
372
	  value1set = value1set || (r(*icursor) == c(*icursor) && col == col_);
373
374
375
376
377
378
	  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());
379
	ins[idx_][idx_] << 1.0;
380
381
382
383
384
      }
    }
    if (asmVector)
      (*getRhsVector(row_))[idx] = value; // set Dirichlet-Value at rhs-vector
      
385
    (*solution->getDOFVector(col_))[idx] = value; // set Dirichlet-value for solution component
386
387
  }
  
388
  
389
390
391
392
  void applyDirichletBC(SingularDirichletBC &sbc, bool asmMatrix = true, bool asmVector = true)
  {
    applyDirichletBC(sbc.row, sbc.col, sbc.idx, sbc.value, asmMatrix, asmVector);
  }
393
  
394
395
396
397
398
399

  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;

400
401
    typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
402

403
404
    size_t numCols = static_cast<size_t>(getNumComponents());
    for (size_t col = 0; col < numCols; col++) {
405
406
407
408
409
      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();

410
411
412
      mtl::traits::row<Matrix>::type r(m);
      mtl::traits::col<Matrix>::type c(m);
      mtl::traits::value<Matrix>::type v(m);
413
414
415
416
417
418
      
      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
419
420
421
	if (indices[i].first == indices[i].second)
	  continue;
	
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
	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
440
441
	    if (indices[i].first == indices[i].second)
	      continue;
442
443
444
445
446
	    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
447
448
	  if (indices[i].first == indices[i].second)
	    continue;
449
450
451
452
453
454
455
	  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;
	  }
	}
      }
    }
  }
456
  
457
458
459
460
461
462

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

463
  
464
465
466
467
468
469
470
471
  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;

472
473
    typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
474

475
    size_t numCols = static_cast<size_t>(getNumComponents());
476
477
    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");
478
    TEST_EXIT(coefficients[0].size() == numCols)
479
480
      ("You have to give coefficients for all variables\n");

481
    for (size_t col = 0; col < numCols; col++) {
482
483
484
485
486
      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();

487
488
489
      mtl::traits::row<Matrix>::type r(m);
      mtl::traits::col<Matrix>::type c(m);
      mtl::traits::value<Matrix>::type v(m);
490
491
492
493
494
495
496
497
498
499
500

      // 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];
      }
    }
  }
501
  
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525

  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;
526
  std::vector<DirichletBcDataBase*> DirichletBcDataList;
527
528
529
530

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