ProblemStat.h 21.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * 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 Vey, Thomas Witkowski, Andreas Naumann, 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.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22

23
/** \file ProblemStat.h */
24

25
26
#ifndef AMDIS_PROBLEM_STAT_H
#define AMDIS_PROBLEM_STAT_H
27

Thomas Witkowski's avatar
Thomas Witkowski committed
28
29
#include <vector>
#include <list>
30
31
32
#if __cplusplus > 199711L
#include <functional>
#endif
33
#include "AMDiS_fwd.h"
34
#include "ProblemStatBase.h"
35
#include "Initfile.h"
36
37
38
#include "Boundary.h"
#include "MatrixVector.h"
#include "StandardProblemIteration.h"
39
#include "io/ElementFileWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
40
#include "ComponentTraverseInfo.h"
41
#include "AbstractFunction.h"
42
#include "solver/SolverMatrix.h"
43
#include "SystemVector.h"
44
45
46

namespace AMDiS {

47
48
49
50
51
52
  struct OperatorPos 
  {
    int row, col;
    double *factor, *estFactor;
  };

53

54
55
56
57
58
59
60
  /**
   * \ingroup Problem 
   *
   * \brief
   * This class defines the stationary problem definition in sequential
   * computations. For parallel computations, see \ref ParallelProblemStatBase.
   **/
61
  class ProblemStatSeq : public ProblemStatBase
62
  {
63
  protected:
64
    // Defines a mapping type from dof indices to world coordinates.
65
    typedef std::map<DegreeOfFreedom, WorldVector<double> > DofToCoord;
66
67

    // Defines a mapping type from dof indices to world coordinates.
68
    typedef std::map<WorldVector<double>, DegreeOfFreedom> CoordToDof;
69

70
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
71
    /// Constructor
72
    ProblemStatSeq(std::string nameStr,
73
		   ProblemIterationInterface *problemIteration = NULL);
74

Thomas Witkowski's avatar
Thomas Witkowski committed
75
    /// Destructor
Thomas Witkowski's avatar
Thomas Witkowski committed
76
    virtual ~ProblemStatSeq();
77

Thomas Witkowski's avatar
Thomas Witkowski committed
78
    /// Initialisation of the problem.
79
    virtual void initialize(Flag initFlag,
80
			    ProblemStatSeq *adoptProblem = NULL,
81
82
			    Flag adoptFlag = INIT_NOTHING);

Thomas Witkowski's avatar
Thomas Witkowski committed
83
    /// Used in \ref initialize().
84
85
    virtual void createMesh();

86
87
88
    /// Used in \ref initialize().
    virtual void createRefCoarseManager();

Thomas Witkowski's avatar
Thomas Witkowski committed
89
    /// Used in \ref initialize().
90
    virtual void createFeSpace(DOFAdmin *admin);
91

Thomas Witkowski's avatar
Thomas Witkowski committed
92
    /// Used in \ref initialize().
93
94
    virtual void createMatricesAndVectors();

Thomas Witkowski's avatar
Thomas Witkowski committed
95
    /// Used in \ref initialize().
96
97
    virtual void createSolver();

Thomas Witkowski's avatar
Thomas Witkowski committed
98
    /// Used in \ref initialize().
99
100
    virtual void createEstimator();

Thomas Witkowski's avatar
Thomas Witkowski committed
101
    /// Used in \ref initialize().
102
103
    virtual void createMarker();

Thomas Witkowski's avatar
Thomas Witkowski committed
104
    /// Used in \ref initialize().
105
106
    virtual void createFileWriter();

107
108
    /// Used in \ref initialize(). This function is deprecated and should not 
    /// be used anymore. There is no guarantee that it will work in the future.
109
110
    virtual void doOtherStuff();

111
112
    /// Implementation of ProblemStatBase::solve(). Deligates the solving
    /// of problems system to \ref solver.
113
114
    void solve(AdaptInfo *adaptInfo,
	       bool createMatrixData = true,
115
	       bool storeMatrixData = false) override;
116

117
118
    /// Implementation of ProblemStatBase::estimate(). Deligates the estimation
    /// to \ref estimator.
119
    void estimate(AdaptInfo *adaptInfo) override;
120

121
122
    /// Implementation of ProblemStatBase::markElements().
    /// Deligated to \ref adapt.
123
    Flag markElements(AdaptInfo *adaptInfo) override;
124

125
126
    /// Implementation of ProblemStatBase::refineMesh(). Deligated to the
    /// RefinementManager of \ref adapt.
127
    Flag refineMesh(AdaptInfo *adaptInfo) override;
128

129
130
    /// Implementation of ProblemStatBase::coarsenMesh(). Deligated to the
    /// CoarseningManager of \ref adapt.
131
    Flag coarsenMesh(AdaptInfo *adaptInfo) override;
132

133
134
    /// Implementation of ProblemStatBase::buildBeforeRefine().
    /// Does nothing here.
135
    void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) override {}
136

137
138
    /// Implementation of ProblemStatBase::buildBeforeCoarsen().
    /// Does nothing here.
139
    void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) override {}
140

141
142
143
    /// Implementation of ProblemStatBase::buildAfterCoarsen().
    /// Assembles \ref A and \ref rhs. With the last two parameters, assembling
    /// can be restricted to matrices or vectors only.
144
    void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
145
				   bool assembleMatrix = true,
146
				   bool assembleVector = true) override;
147

148
149
    bool dualMeshTraverseRequired();

150
151
    void dualAssemble(AdaptInfo *adaptInfo, Flag flag, 
		      bool asmMatrix = true, bool asmVector = true);
152

153
154
    /// Returns nr of components \ref nComponents
    virtual int getNumComponents()
155
    { 
156
      return nComponents; 
Thomas Witkowski's avatar
Thomas Witkowski committed
157
    }
158
    
159
    /// Returns nr of additional components \ref nAddComponents
160
161
162
163
    virtual int getNumAddComponents()
    {
      return nAddComponents;
    }
164

165
    /// Writes output files. TODO: Make obsolete.
166
167
    void writeFiles(AdaptInfo *adaptInfo, bool force);

168
169
170
    /// Writes output files.
    void writeFiles(AdaptInfo &adaptInfo, bool force);

Thomas Witkowski's avatar
Thomas Witkowski committed
171
    /// Interpolates fct to \ref solution.
172
    void interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct);
173

174
    /// Adds an operator to \ref A.
175
    void addMatrixOperator(Operator *op, int i, int j,
176
			   double *factor = NULL, double *estFactor = NULL);
177

178
    /// Adds an operator to \ref A.
179
    void addMatrixOperator(Operator &op, int i, int j,
180
			   double *factor = NULL, double *estFactor = NULL);
181

182
    /// Adds an operator to \ref rhs.
183
    void addVectorOperator(Operator *op, int i,
184
			   double *factor = NULL, double *estFactor = NULL);
185

186
    /// Adds an operator to \ref rhs.
187
    void addVectorOperator(Operator &op, int i,
188
			   double *factor = NULL, double *estFactor = NULL);
189

190
191
    /// Adds a Dirichlet boundary condition, where the rhs is given by an 
    /// abstract function.
192
    virtual void addDirichletBC(BoundaryType type, int row, int col,
193
				AbstractFunction<double, WorldVector<double> > *b);
194
195
196
197
198
199
200
201
				
				
#if __cplusplus > 199711L
    /// Adds a Dirichlet boundary condition, where the rhs is given by an 
    /// lambda function or a std::function object
    virtual void addDirichletBC(BoundaryType type, int row, int col,
				std::function<double(WorldVector<double>)> b);
#endif
202

203
204
    /// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
    /// vector.
205
206
207
    virtual void addDirichletBC(BoundaryType type, int row, int col,
				DOFVector<double> *vec);

208
209
    /// Adds a Neumann boundary condition, where the rhs is given by an
    /// abstract function.
210
211
212
    virtual void addNeumannBC(BoundaryType type, int row, int col, 
			      AbstractFunction<double, WorldVector<double> > *n);

213
214
215
216
217
218
    /// Adds a Neumann boundary condition, where the rhs is given by an DOF
    /// vector.
    virtual void addNeumannBC(BoundaryType type, int row, int col, 
			      DOFVector<double> *n);

    /// Adds Robin boundary condition.
219
220
221
222
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    AbstractFunction<double, WorldVector<double> > *n,
			    AbstractFunction<double, WorldVector<double> > *r);

223
224
225
226
227
228
229
230
231
232
    /// Adds Robin boundary condition.
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    DOFVector<double> *n,
			    DOFVector<double> *r);

    /// Adds Robin boundary condition.
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    Operator *n,
			    Operator *r);

233
    /// Adds a periodic boundary condition.
234
235
    virtual void addPeriodicBC(BoundaryType type, int row, int col);

236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
    /// add boundary operator to matrix side
    virtual void addBoundaryMatrixOperator(BoundaryType type, 
          Operator *op, int row, int col);

    virtual void addBoundaryMatrixOperator(BoundaryType type, 
          Operator &op, int row, int col)
    {
      addBoundaryMatrixOperator(type, &op, row, col);
    }

    /// add boundary operator to rhs (vector) side
    virtual void addBoundaryVectorOperator(BoundaryType type, 
          Operator *op, int row);

    virtual void addBoundaryVectorOperator(BoundaryType type, 
          Operator &op, int row)
    {
      addBoundaryVectorOperator(type, &op, row);
    }

256
257
    /// This function assembles a DOFMatrix and a DOFVector for the case,
    /// the meshes from row and col FE-space are equal.
258
    void assembleOnOneMesh(const FiniteElemSpace *feSpace, 
Thomas Witkowski's avatar
Thomas Witkowski committed
259
			   Flag assembleFlag,
Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
			   DOFMatrix *matrix, DOFVector<double> *vector);

Thomas Witkowski's avatar
Thomas Witkowski committed
262
263

    ///
264
265
266
267
    void assembleBoundaryConditions(DOFVector<double> *rhs,
				    DOFVector<double> *solution,
				    Mesh *mesh,
				    Flag assembleFlag);
268
269

 
270
271
272
273
    /** \name getting methods
     * \{ 
     */

Thomas Witkowski's avatar
Thomas Witkowski committed
274
    /// Returns \ref solution.
275
276
    inline SystemVector* getSolution() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
      return solution; 
    }
279

280
    inline DOFVector<double>* getSolution(int i)
281
282
283
284
    {
      return solution->getDOFVector(i);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
285
    /// Returns \ref rhs.
286
    inline SystemVector* getRhs() 
287
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
      return rhs; 
    }
290

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

Thomas Witkowski's avatar
Thomas Witkowski committed
296
    /// Returns \ref systemMatrix.
297
298
    inline Matrix<DOFMatrix*> *getSystemMatrix() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
      return systemMatrix; 
    }
301

Thomas Witkowski's avatar
Thomas Witkowski committed
302
    /// Returns a pointer to the corresponding DOFMatrix.
303
304
    inline DOFMatrix* getSystemMatrix(int row, int col) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
      return (*systemMatrix)[row][col];
    }

    /// Returns mesh of given component
309
    inline Mesh* getMesh(int comp = 0) 
310
    {
311
      FUNCNAME("ProblemStatSeq::getMesh()");
Thomas Witkowski's avatar
Thomas Witkowski committed
312
      TEST_EXIT(comp < static_cast<int>(componentMeshes.size()) && comp >= 0)
313
	("invalid component number\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
314
      return componentMeshes[comp]; 
Thomas Witkowski's avatar
Thomas Witkowski committed
315
    }
316

Thomas Witkowski's avatar
Thomas Witkowski committed
317
    /// Returns \ref meshes
318
    inline std::vector<Mesh*> getMeshes() 
319
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
320
321
      return meshes; 
    }
322

Thomas Witkowski's avatar
Thomas Witkowski committed
323
    /// Returns \ref feSpace_.
324
    inline const FiniteElemSpace* getFeSpace(int comp = 0) 
325
    { 
326
      FUNCNAME("ProblemStatSeq::getFeSpace()");
327
      TEST_EXIT(comp < static_cast<int>(componentSpaces.size()) && comp >= 0)
328
	("invalid component number: %d\n", comp);
329
      return componentSpaces[comp]; 
Thomas Witkowski's avatar
Thomas Witkowski committed
330
    }
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
    /// Returns \ref feSpaces.
333
    inline std::vector<const FiniteElemSpace*>& getFeSpaces() 
334
    { 
335
      return feSpaces; 
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    }
337

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    /// Returns \ref componentSpaces;
339
    inline std::vector<const FiniteElemSpace*>& getComponentSpaces() 
340
    {
341
      return componentSpaces;
342
343
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
344
    /// Returns \ref estimator.
345
    inline std::vector<Estimator*> getEstimators() 
346
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
      return estimator; 
    }
349

Thomas Witkowski's avatar
Thomas Witkowski committed
350
    /// Returns \ref estimator.
351
    inline Estimator* getEstimator(int comp = 0) 
352
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
353
354
      return estimator[comp]; 
    }
355

Thomas Witkowski's avatar
Thomas Witkowski committed
356
    /// Returns \ref refinementManager.
357
    inline RefinementManager* getRefinementManager(int comp = 0) 
358
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
      return refinementManager; 
    }
361

Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Returns \ref refinementManager.
363
    inline CoarseningManager* getCoarseningManager(int comp = 0) 
364
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
      return coarseningManager; 
    }
367

Thomas Witkowski's avatar
Thomas Witkowski committed
368
    /// Returns \ref solver.
369
    inline LinearSolverInterface* getSolver() 
370
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
      return solver; 
    }
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
    /// Returns \ref marker.
375
    inline Marker *getMarker(int comp = 0) 
376
    { 
377
      return marker[comp]; 
Thomas Witkowski's avatar
Thomas Witkowski committed
378
    }
379

Thomas Witkowski's avatar
Thomas Witkowski committed
380
    /// Returns \ref marker.
381
    inline std::vector<Marker*> getMarkers() 
382
    { 
383
      return marker; 
Thomas Witkowski's avatar
Thomas Witkowski committed
384
    }
385

Thomas Witkowski's avatar
Thomas Witkowski committed
386
    /// Returns the name of the problem
387
    inline std::string getName() override
388
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
      return name; 
    }
391

Praetorius, Simon's avatar
Praetorius, Simon committed
392
    /// Returns the name of the problem
393
    inline std::string getComponentName(int comp = 0)
Praetorius, Simon's avatar
Praetorius, Simon committed
394
395
396
397
398
399
    {
      TEST_EXIT(comp < static_cast<int>(componentNames.size()) && comp >= 0)
	("invalid component number\n");
      return componentNames[comp];
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
400
    /// Returns \ref useGetBound.
401
402
    inline bool getBoundUsed() 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
      return useGetBound; 
    }
405

406
407
408
409
410
411
    /// Returns \ref info.
    int getInfo()
    {
      return info;
    }

412
413
414
415
416
417
    /// Returns \ref deserialized;
    bool isDeserialized()
    {
      return deserialized;
    }

418
419
420
421
422
423
    /** \} */

    /** \name setting methods
     * \{ 
     */

Thomas Witkowski's avatar
Thomas Witkowski committed
424
    /// Sets \ref estimator.
425
    inline void setEstimator(std::vector<Estimator*> est) 
426
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
      estimator = est; 
    }
429

Thomas Witkowski's avatar
Thomas Witkowski committed
430
    /// Sets the FE space for the given component.
431
    inline void setFeSpace(const FiniteElemSpace *feSpace, int comp = 0) 
432
    {
433
      feSpaces[comp] = feSpace;
Thomas Witkowski's avatar
Thomas Witkowski committed
434
    }
435

436
437
438
439
440
    void setFeSpaces(std::vector<FiniteElemSpace const*> feSpaces_)
    {
      feSpaces = feSpaces_;
    }

441
    void setComponentSpace(int comp, const FiniteElemSpace *feSpace)
442
443
444
    {
      if (static_cast<int>(componentSpaces.size()) < nComponents)
        componentSpaces.resize(nComponents);
445
      TEST_EXIT(comp >= 0 && comp < nComponents + nAddComponents)
446
447
448
449
450
        ("Component number not in feasable range!");

      componentSpaces[comp] = feSpace;
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    /// Sets \ref estimator.
452
    inline void setEstimator(Estimator* est, int comp = 0) 
453
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
      estimator[comp] = est; 
    }
456

Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Sets \ref marker.
458
    inline void setMarker(Marker* mark, int comp = 0) 
459
    { 
460
      marker[comp] = mark; 
Thomas Witkowski's avatar
Thomas Witkowski committed
461
    }
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
    /// Sets \ref solver.
464
    inline void setSolver(LinearSolverInterface* sol) 
465
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
      solver = sol; 
    }
468

469
470
471
472
473
    void setSolverMatrix(SolverMatrix<Matrix<DOFMatrix*> >& solverMatrix_)
    {
      solverMatrix.setMatrix(*solverMatrix_.getOriginalMat());
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
474
    ///
475
    inline void setAssembleMatrixOnlyOnce(int i = 0, int j = 0, bool value = true) 
476
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
477
      assembleMatrixOnlyOnce[i][j] = value;
478
479
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
480
    ///
481
    void setExactSolutionFct(AbstractFunction<double, WorldVector<double> > *fct,
482
483
			     int component) 
    {
484
485
486
      exactSolutionFcts[component] = fct;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
487
    ///
488
    AbstractFunction<double, WorldVector<double> > *getExactSolutionFct(int i = 0) 
489
    {
490
      return exactSolutionFcts[i];
491
492
    }

493
    ///
494
    std::vector< AbstractFunction<double, WorldVector<double> >* > getExactSolutionFcts() 
495
496
497
498
    {
      return exactSolutionFcts;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
499
    ///
500
501
    void setComputeExactError(bool v) 
    {
502
503
504
      computeExactError = v;
    }

505
506
507
508
509
510
    /// Sets \ref writeAsmInfo;
    void setWriteAsmInfo(bool b)
    {
      writeAsmInfo = b;
    }

511
    void setMeshes(std::vector<Mesh*> meshes_)
512
513
514
515
516
517
518
519
520
    {
      meshes = meshes_;
      nMeshes = static_cast<int>(meshes.size());
    }

    void setComponentMesh(int comp, Mesh* mesh)
    {
      if (static_cast<int>(componentMeshes.size()) < nComponents)
        componentMeshes.resize(nComponents);
521
      TEST_EXIT(comp >= 0 && comp < nComponents + nAddComponents)
522
523
524
525
526
527
528
529
530
531
532
533
534
535
        ("Component number not in feasable range!");

      componentMeshes[comp] = mesh;
    }

    void setRefinementManager(RefinementManager *ref)
    {
      refinementManager = ref;
    }

    void setCoarseningManager(CoarseningManager *coarse)
    {
      coarseningManager = coarse;
    }
536
537
    /** \} */

538
539
    /// Outputs the mesh of the given component, but the values are taken from
    /// the residual error estimator. 
540
    void writeResidualMesh(int comp, AdaptInfo *adaptInfo, std::string name);
541

Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Function that implements the serialization procedure.
543
    void serialize(std::ostream &out) override;
544

Thomas Witkowski's avatar
Thomas Witkowski committed
545
    /// Function that implements the deserialization procedure.
546
    void deserialize(std::istream &in) override;
Thomas Witkowski's avatar
Thomas Witkowski committed
547
548

    /// Returns \ref fileWriters.
549
    std::vector<FileWriterInterface*>& getFileWriterList() 
550
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
      return fileWriters;
    }
553
    
554
    /// Returns \ref solutionTime.
555
556
557
558
559
    double getSolutionTime()
    {
      return solutionTime;
    }
    
560
    /// Returns \ref buildTime.
561
562
563
564
    double getBuildTime()
    {
      return buildTime;
    }
565

Thomas Witkowski's avatar
Thomas Witkowski committed
566
  protected:
567
568
    /// If the exact solution is known, the problem can compute the exact
    /// error instead of the error estimation. This is done in this function.
Thomas Witkowski's avatar
Thomas Witkowski committed
569
570
    void computeError(AdaptInfo *adaptInfo);

571
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
    
    /// Name of this problem.
574
    std::string name;
575

Thomas Witkowski's avatar
Thomas Witkowski committed
576
    /// Number of problem components
577
    int nComponents;
578
579
580
    
    /// Number of additional components
    int nAddComponents;
581

582
583
    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
584
    std::vector<std::string> componentNames;
585

586
587
588
    /// Number of problem meshes. If all components are defined on the same mesh,
    /// this number is 1. Otherwise, this variable is the number of different meshes
    /// within the problem.
Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
591
    int nMeshes;

    /// FE spaces of this problem.
592
    std::vector<const FiniteElemSpace*> feSpaces;
593

Thomas Witkowski's avatar
Thomas Witkowski committed
594
    /// Meshes of this problem.
595
    std::vector<Mesh*> meshes;
596

Thomas Witkowski's avatar
Thomas Witkowski committed
597
    /// Pointer to the fe spaces for the different problem components
598
    std::vector<const FiniteElemSpace*> componentSpaces;
599

Thomas Witkowski's avatar
Thomas Witkowski committed
600
    /// Pointer to the meshes for the different problem components
601
    std::vector<Mesh*> componentMeshes;
602

603
604
605
    /// Stores information about which meshes must be traversed to assemble the
    /// specific components. I.e., it was implemented to make use of different
    /// meshes for different components.
Thomas Witkowski's avatar
Thomas Witkowski committed
606
607
608
    ComponentTraverseInfo traverseInfo;
    
    /// Responsible for element marking.
609
    std::vector<Marker*> marker;
610

Thomas Witkowski's avatar
Thomas Witkowski committed
611
    /// Estimator of this problem. Used in \ref estimate().
612
    std::vector<Estimator*> estimator;
613

Thomas Witkowski's avatar
Thomas Witkowski committed
614
    /// Linear solver of this problem. Used in \ref solve().
615
    LinearSolverInterface *solver;
616

Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
    /// System vector  storing the calculated solution of the problem.
    SystemVector *solution;
619

Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
    /// System vector for the right hand side 
    SystemVector *rhs;
622

Thomas Witkowski's avatar
Thomas Witkowski committed
623
624
    /// System matrix
    Matrix<DOFMatrix*> *systemMatrix;
625

626
627
628
    /// Composed system matrix
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;

629
630
631
632
633
    /// Some DOFMatrices of the systemMatrix may be assembled only once (for
    /// example if they are independent of the time or older solutions). If
    /// [i][j] of this field is set to true, the corresponding DOFMatrix will
    /// be assembled only once. All other matrices will be assembled at every
    /// time step.
634
    std::vector<std::vector<bool> > assembleMatrixOnlyOnce;
635

636
637
638
639
    /// If [i][j] of this field is set to true, the corresponding DOFMatrix of
    /// the systemMatrix has been assembled at least once. This field is used
    /// to determine, if assembling of a matrix can be ommitted, if it is set
    /// to be assembled only once.
640
    std::vector<std::vector<bool> > assembledMatrix;
641

Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
    /// Determines whether domain boundaries should be considered at assembling.
    bool useGetBound;
644

Thomas Witkowski's avatar
Thomas Witkowski committed
645
    /// Writes the meshes and solution after the adaption loop.
646
    std::vector<FileWriterInterface*> fileWriters;
647

648
649
650
    /// All actions of mesh refinement are performed by refinementManager.
    /// If new refinement algorithms should be realized, one has to override
    /// RefinementManager and give one instance of it to AdaptStationary.
Thomas Witkowski's avatar
Thomas Witkowski committed
651
    RefinementManager *refinementManager;
652

653
654
655
    /// All actions of mesh coarsening are performed by coarseningManager.
    /// If new coarsening algorithms should be realized, one has to override
    /// CoarseningManager and give one instance of it to AdaptStationary.
Thomas Witkowski's avatar
Thomas Witkowski committed
656
    CoarseningManager *coarseningManager;
657
  
Thomas Witkowski's avatar
Thomas Witkowski committed
658
659
    /// Info level.
    int info;
660

661
662
    /// If true, the stationary problem was deserialized from some serialization
    /// file.
663
664
    bool deserialized;

665
666
667
    /// This vectors stores pointers to functions defining the exact solution of
    /// the problem. This may be used to compute the real error of the computed
    /// solution.
668
    std::vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
669

670
671
    /// If true, the error is not estimated but computed from the exact solution
    /// defined by \ref exactSolutionFcts.
672
    bool computeExactError;
673
    
674
675
676
677
    /// If at least on boundary condition is set, this variable is true. It is 
    /// used to ensure that no operators are added after boundary condition were
    /// set. If this would happen, boundary conditions could set wrong on off 
    /// diagonal matrices.
678
    bool boundaryConditionSet;
679

680
681
    /// If true, AMDiS prints information about the assembling procedure to 
    /// the screen.
682
    bool writeAsmInfo;
683

684
    std::map<Operator*, std::vector<OperatorPos> > operators;
685
    
686
    /// time needed to solve the linear system
687
    double solutionTime;
688
689
    
    /// time needed to assemble the linear system
690
    double buildTime;
691
692
    
    template <class> friend class detail::CouplingProblemStat;
693
  };
694
695
  
  namespace detail
696
  {
697
698
699
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType, 
			 public StandardProblemIteration
700
    {
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
      using ProblemStatType::getName;
      
      /// Constructor
      ProblemStat(std::string nameStr,
		  ProblemIterationInterface *problemIteration = NULL)
	: ProblemStatType(nameStr, problemIteration),
	  StandardProblemIteration(this)
      { }

      /// Returns number of managed problems
      // implements StandardProblemIteration::getNumProblems()
      virtual int getNumProblems() override
      { 
	return 1; 
      }

      /// Returns the problem with the given number. If only one problem
      /// is managed by this master problem, the number hasn't to be given.
      // implements StandardProblemIteration::getProblem(int)
      virtual ProblemStatBase *getProblem(int number = 0) override
      { 
	return this; 
      }

      /// Determines the execution order of the single adaption steps. If adapt is
      /// true, mesh adaption will be performed. This allows to avoid mesh adaption,
      /// e.g. in timestep adaption loops of timestep adaptive strategies.
      // implements StandardProblemIteration::oneIteration(AdaptInfo*, Flag)
      virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) override
      {	
	for (int i = 0; i < ProblemStatType::getNumComponents(); i++)
	  if (adaptInfo->spaceToleranceReached(i))
	    adaptInfo->allowRefinement(false, i);
	  else
	    adaptInfo->allowRefinement(true, i);	

	return StandardProblemIteration::oneIteration(adaptInfo, toDo);
      }
    };
  }
  
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
  typedef detail::ProblemStat<ProblemStatSeq> ProblemStat;
744
#endif
745
746
747
}

#endif