ProblemStat.h 22.2 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * 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.
18
 *
19
 ******************************************************************************/
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
#include <functional>
31
#include "AMDiS_fwd.h"
32
#include "ProblemStatBase.h"
33
#include "Initfile.h"
34
35
36
#include "Boundary.h"
#include "MatrixVector.h"
#include "StandardProblemIteration.h"
37
#include "io/ElementFileWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
38
#include "ComponentTraverseInfo.h"
39
#include "AbstractFunction.h"
40
#include "solver/SolverMatrix.h"
41
#include "SystemVector.h"
42
43
44

namespace AMDiS {

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

51

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
73
    /// Destructor
Thomas Witkowski's avatar
Thomas Witkowski committed
74
    virtual ~ProblemStatSeq();
75

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

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

84
85
86
    /// Used in \ref initialize().
    virtual void createRefCoarseManager();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
146
147
148
149
150
    /// assemble all operators of matrix and vector side
    void assemble(AdaptInfo* adaptInfo)
    {
      buildAfterCoarsen(adaptInfo, 0, true, true);
    }
151

152
153
    bool dualMeshTraverseRequired();

154
    void dualAssemble(AdaptInfo *adaptInfo, Flag flag,
155
		      bool asmMatrix = true, bool asmVector = true);
156

157
    /// Returns nr of components \ref nComponents
Praetorius, Simon's avatar
Praetorius, Simon committed
158
    virtual int getNumComponents() const
159
160
    {
      return nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
161
    }
162

163
    /// Returns nr of additional components \ref nAddComponents
Praetorius, Simon's avatar
Praetorius, Simon committed
164
    virtual int getNumAddComponents() const
165
166
167
    {
      return nAddComponents;
    }
168

169
    /// Writes output files. TODO: Make obsolete.
170
171
    void writeFiles(AdaptInfo *adaptInfo, bool force);

172
173
174
    /// Writes output files.
    void writeFiles(AdaptInfo &adaptInfo, bool force);

Thomas Witkowski's avatar
Thomas Witkowski committed
175
    /// Interpolates fct to \ref solution.
176
    void interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct);
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 A.
183
    void addMatrixOperator(Operator &op, int i, int j,
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
    /// Adds an operator to \ref rhs.
191
    void addVectorOperator(Operator &op, int i,
192
			   double *factor = NULL, double *estFactor = NULL);
193

194
    /// Adds a Dirichlet boundary condition, where the rhs is given by an
195
    /// abstract function.
196
    virtual void addDirichletBC(BoundaryType type, int row, int col,
197
				AbstractFunction<double, WorldVector<double> > *b);
198
199


200
#if AMDIS_HAS_CXX11
201
    /// Adds a Dirichlet boundary condition, where the rhs is given by an
202
203
204
205
    /// lambda function or a std::function object
    virtual void addDirichletBC(BoundaryType type, int row, int col,
				std::function<double(WorldVector<double>)> b);
#endif
206

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

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

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

    /// Adds Robin boundary condition.
223
    virtual void addRobinBC(BoundaryType type, int row, int col,
224
225
226
			    AbstractFunction<double, WorldVector<double> > *n,
			    AbstractFunction<double, WorldVector<double> > *r);

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

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

237
    /// Adds a periodic boundary condition.
238
239
    virtual void addPeriodicBC(BoundaryType type, int row, int col);

240
    /// add boundary operator to matrix side
241
    virtual void addBoundaryMatrixOperator(BoundaryType type,
242
243
          Operator *op, int row, int col);

244
    virtual void addBoundaryMatrixOperator(BoundaryType type,
245
246
247
248
249
250
          Operator &op, int row, int col)
    {
      addBoundaryMatrixOperator(type, &op, row, col);
    }

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

254
    virtual void addBoundaryVectorOperator(BoundaryType type,
255
256
257
258
259
          Operator &op, int row)
    {
      addBoundaryVectorOperator(type, &op, row);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
266
267

    ///
268
269
270
271
    void assembleBoundaryConditions(DOFVector<double> *rhs,
				    DOFVector<double> *solution,
				    Mesh *mesh,
				    Flag assembleFlag);
272

273

274
    /** \name getting methods
275
     * \{
276
277
     */

Thomas Witkowski's avatar
Thomas Witkowski committed
278
    /// Returns \ref solution.
279
280
281
    inline SystemVector* getSolution()
    {
      return solution;
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    }
283
284
285
286
287
    
    inline SystemVector const* getSolution() const
    {
      return solution;
    }
288

289
    inline DOFVector<double>* getSolution(int i)
290
291
292
293
    {
      return solution->getDOFVector(i);
    }

294
295
296
297
298
    inline DOFVector<double> const* getSolution(int i) const
    {
      return solution->getDOFVector(i);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
299
    /// Returns \ref rhs.
300
301
302
    inline SystemVector* getRhs()
    {
      return rhs;
Thomas Witkowski's avatar
Thomas Witkowski committed
303
    }
304

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

Thomas Witkowski's avatar
Thomas Witkowski committed
310
    /// Returns \ref systemMatrix.
311
312
313
    inline Matrix<DOFMatrix*> *getSystemMatrix()
    {
      return systemMatrix;
Thomas Witkowski's avatar
Thomas Witkowski committed
314
    }
315

Thomas Witkowski's avatar
Thomas Witkowski committed
316
    /// Returns a pointer to the corresponding DOFMatrix.
317
    inline DOFMatrix* getSystemMatrix(int row, int col)
318
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
321
322
      return (*systemMatrix)[row][col];
    }

    /// Returns mesh of given component
323
    inline Mesh* getMesh(int comp = 0)
324
    {
325
      FUNCNAME("ProblemStatSeq::getMesh()");
Thomas Witkowski's avatar
Thomas Witkowski committed
326
      TEST_EXIT(comp < static_cast<int>(componentMeshes.size()) && comp >= 0)
327
	("invalid component number\n");
328
      return componentMeshes[comp];
Thomas Witkowski's avatar
Thomas Witkowski committed
329
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
330
331
332
333
334
335
336
337
    
    Mesh const* getMesh(int comp = 0) const
    {
      FUNCNAME("ProblemStatSeq::getMesh()");
      TEST_EXIT(comp < static_cast<int>(componentMeshes.size()) && comp >= 0)
	("invalid component number\n");
      return componentMeshes[comp];
    }
338

Thomas Witkowski's avatar
Thomas Witkowski committed
339
    /// Returns \ref meshes
340
    inline std::vector<Mesh*> getMeshes()
341
    {
342
      return meshes;
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    }
344

Thomas Witkowski's avatar
Thomas Witkowski committed
345
    /// Returns \ref feSpace_.
Praetorius, Simon's avatar
Praetorius, Simon committed
346
    inline const FiniteElemSpace* getFeSpace(int comp = 0) const
347
    {
348
      FUNCNAME("ProblemStatSeq::getFeSpace()");
349
      TEST_EXIT(comp < static_cast<int>(componentSpaces.size()) && comp >= 0)
350
	("invalid component number: %d\n", comp);
351
      return componentSpaces[comp];
Thomas Witkowski's avatar
Thomas Witkowski committed
352
    }
353

Thomas Witkowski's avatar
Thomas Witkowski committed
354
    /// Returns \ref feSpaces.
355
356
357
    inline std::vector<const FiniteElemSpace*>& getFeSpaces()
    {
      return feSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    }
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    /// Returns \ref componentSpaces;
361
    inline std::vector<const FiniteElemSpace*>& getComponentSpaces()
362
    {
363
      return componentSpaces;
364
365
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
366
    /// Returns \ref estimator.
367
368
369
    inline std::vector<Estimator*> getEstimators()
    {
      return estimator;
Thomas Witkowski's avatar
Thomas Witkowski committed
370
    }
371

Thomas Witkowski's avatar
Thomas Witkowski committed
372
    /// Returns \ref estimator.
373
374
375
    inline Estimator* getEstimator(int comp = 0)
    {
      return estimator[comp];
Thomas Witkowski's avatar
Thomas Witkowski committed
376
    }
377

Thomas Witkowski's avatar
Thomas Witkowski committed
378
    /// Returns \ref refinementManager.
379
380
381
    inline RefinementManager* getRefinementManager(int comp = 0)
    {
      return refinementManager;
Thomas Witkowski's avatar
Thomas Witkowski committed
382
    }
383

Thomas Witkowski's avatar
Thomas Witkowski committed
384
    /// Returns \ref refinementManager.
385
386
387
    inline CoarseningManager* getCoarseningManager(int comp = 0)
    {
      return coarseningManager;
Thomas Witkowski's avatar
Thomas Witkowski committed
388
    }
389

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    /// Returns \ref solver.
391
392
393
    inline LinearSolverInterface* getSolver()
    {
      return solver;
Thomas Witkowski's avatar
Thomas Witkowski committed
394
    }
395

Thomas Witkowski's avatar
Thomas Witkowski committed
396
    /// Returns \ref marker.
397
398
399
    inline Marker *getMarker(int comp = 0)
    {
      return marker[comp];
Thomas Witkowski's avatar
Thomas Witkowski committed
400
    }
401

Thomas Witkowski's avatar
Thomas Witkowski committed
402
    /// Returns \ref marker.
403
404
405
    inline std::vector<Marker*> getMarkers()
    {
      return marker;
Thomas Witkowski's avatar
Thomas Witkowski committed
406
    }
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    /// Returns the name of the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
409
    inline std::string getName() const override
410
411
    {
      return name;
Thomas Witkowski's avatar
Thomas Witkowski committed
412
    }
413

Praetorius, Simon's avatar
Praetorius, Simon committed
414
    /// Returns the name of the problem
Praetorius, Simon's avatar
Praetorius, Simon committed
415
    inline std::string getComponentName(int comp = 0) const
Praetorius, Simon's avatar
Praetorius, Simon committed
416
417
418
419
420
421
    {
      TEST_EXIT(comp < static_cast<int>(componentNames.size()) && comp >= 0)
	("invalid component number\n");
      return componentNames[comp];
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
422
    /// Returns \ref useGetBound.
Praetorius, Simon's avatar
Praetorius, Simon committed
423
    inline bool getBoundUsed() const
424
425
    {
      return useGetBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
426
    }
427

428
    /// Returns \ref info.
Praetorius, Simon's avatar
Praetorius, Simon committed
429
    int getInfo() const
430
431
432
433
    {
      return info;
    }

434
    /// Returns \ref deserialized;
Praetorius, Simon's avatar
Praetorius, Simon committed
435
    bool isDeserialized() const
436
437
438
439
    {
      return deserialized;
    }

440
441
442
    /** \} */

    /** \name setting methods
443
     * \{
444
445
     */

Thomas Witkowski's avatar
Thomas Witkowski committed
446
    /// Sets \ref estimator.
447
448
449
    inline void setEstimator(std::vector<Estimator*> est)
    {
      estimator = est;
Thomas Witkowski's avatar
Thomas Witkowski committed
450
    }
451

Thomas Witkowski's avatar
Thomas Witkowski committed
452
    /// Sets the FE space for the given component.
453
    inline void setFeSpace(const FiniteElemSpace *feSpace, int comp = 0)
454
    {
455
      feSpaces[comp] = feSpace;
Thomas Witkowski's avatar
Thomas Witkowski committed
456
    }
457

458
459
460
461
462
    void setFeSpaces(std::vector<FiniteElemSpace const*> feSpaces_)
    {
      feSpaces = feSpaces_;
    }

463
    void setComponentSpace(int comp, const FiniteElemSpace *feSpace)
464
465
466
    {
      if (static_cast<int>(componentSpaces.size()) < nComponents)
        componentSpaces.resize(nComponents);
467
      TEST_EXIT(comp >= 0 && comp < nComponents + nAddComponents)
468
469
470
471
        ("Component number not in feasable range!");

      componentSpaces[comp] = feSpace;
    }
472

Thomas Witkowski's avatar
Thomas Witkowski committed
473
    /// Sets \ref estimator.
474
475
476
    inline void setEstimator(Estimator* est, int comp = 0)
    {
      estimator[comp] = est;
Thomas Witkowski's avatar
Thomas Witkowski committed
477
    }
478

Thomas Witkowski's avatar
Thomas Witkowski committed
479
    /// Sets \ref marker.
480
481
482
    inline void setMarker(Marker* mark, int comp = 0)
    {
      marker[comp] = mark;
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    }
484

Thomas Witkowski's avatar
Thomas Witkowski committed
485
    /// Sets \ref solver.
486
487
488
    inline void setSolver(LinearSolverInterface* sol)
    {
      solver = sol;
Thomas Witkowski's avatar
Thomas Witkowski committed
489
    }
490

491
492
493
494
495
    void setSolverMatrix(SolverMatrix<Matrix<DOFMatrix*> >& solverMatrix_)
    {
      solverMatrix.setMatrix(*solverMatrix_.getOriginalMat());
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
496
    ///
497
    inline void setAssembleMatrixOnlyOnce(int i = 0, int j = 0, bool value = true)
498
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
499
      assembleMatrixOnlyOnce[i][j] = value;
500
501
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
502
    ///
503
    void setExactSolutionFct(AbstractFunction<double, WorldVector<double> > *fct,
504
			     int component)
505
    {
506
507
508
      exactSolutionFcts[component] = fct;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
509
    ///
510
    AbstractFunction<double, WorldVector<double> > *getExactSolutionFct(int i = 0)
511
    {
512
      return exactSolutionFcts[i];
513
514
    }

515
    ///
516
    std::vector< AbstractFunction<double, WorldVector<double> >* > getExactSolutionFcts()
517
518
519
520
    {
      return exactSolutionFcts;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
521
    ///
522
    void setComputeExactError(bool v)
523
    {
524
525
526
      computeExactError = v;
    }

527
528
529
530
531
532
    /// Sets \ref writeAsmInfo;
    void setWriteAsmInfo(bool b)
    {
      writeAsmInfo = b;
    }

533
    void setMeshes(std::vector<Mesh*> meshes_)
534
535
536
537
538
539
540
541
542
    {
      meshes = meshes_;
      nMeshes = static_cast<int>(meshes.size());
    }

    void setComponentMesh(int comp, Mesh* mesh)
    {
      if (static_cast<int>(componentMeshes.size()) < nComponents)
        componentMeshes.resize(nComponents);
543
      TEST_EXIT(comp >= 0 && comp < nComponents + nAddComponents)
544
545
546
547
548
549
550
551
552
553
554
555
556
557
        ("Component number not in feasable range!");

      componentMeshes[comp] = mesh;
    }

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

    void setCoarseningManager(CoarseningManager *coarse)
    {
      coarseningManager = coarse;
    }
558
559
    /** \} */

560
    /// Outputs the mesh of the given component, but the values are taken from
561
    /// the residual error estimator.
562
    void writeResidualMesh(int comp, AdaptInfo *adaptInfo, std::string name);
563

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

Thomas Witkowski's avatar
Thomas Witkowski committed
567
    /// Function that implements the deserialization procedure.
568
    void deserialize(std::istream &in) override;
Thomas Witkowski's avatar
Thomas Witkowski committed
569
570

    /// Returns \ref fileWriters.
571
    std::vector<FileWriterInterface*>& getFileWriterList()
572
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
      return fileWriters;
    }
575

576
    /// Returns \ref solutionTime.
Praetorius, Simon's avatar
Praetorius, Simon committed
577
    double getSolutionTime() const
578
579
580
    {
      return solutionTime;
    }
581

582
    /// Returns \ref buildTime.
Praetorius, Simon's avatar
Praetorius, Simon committed
583
    double getBuildTime() const
584
585
586
    {
      return buildTime;
    }
587

Thomas Witkowski's avatar
Thomas Witkowski committed
588
  protected:
589
590
    /// 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
591
592
    void computeError(AdaptInfo *adaptInfo);

593
  protected:
594

Thomas Witkowski's avatar
Thomas Witkowski committed
595
    /// Name of this problem.
596
    std::string name;
597

Thomas Witkowski's avatar
Thomas Witkowski committed
598
    /// Number of problem components
599
    int nComponents;
600

601
602
    /// Number of additional components
    int nAddComponents;
603

604
605
    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
606
    std::vector<std::string> componentNames;
607

608
609
610
    /// 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
611
612
613
    int nMeshes;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
616
    /// Meshes of this problem.
617
    std::vector<Mesh*> meshes;
618

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

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

625
626
627
    /// 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
628
    ComponentTraverseInfo traverseInfo;
629

Thomas Witkowski's avatar
Thomas Witkowski committed
630
    /// Responsible for element marking.
631
    std::vector<Marker*> marker;
632

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

Thomas Witkowski's avatar
Thomas Witkowski committed
636
    /// Linear solver of this problem. Used in \ref solve().
637
    LinearSolverInterface *solver;
638

Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
    /// System vector  storing the calculated solution of the problem.
    SystemVector *solution;
641

642
    /// System vector for the right hand side
Thomas Witkowski's avatar
Thomas Witkowski committed
643
    SystemVector *rhs;
644

Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
    /// System matrix
    Matrix<DOFMatrix*> *systemMatrix;
647

648
649
650
    /// Composed system matrix
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;

651
652
653
654
655
    /// 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.
656
    std::vector<std::vector<bool> > assembleMatrixOnlyOnce;
657

658
659
660
661
    /// 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.
662
    std::vector<std::vector<bool> > assembledMatrix;
663

Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
    /// Determines whether domain boundaries should be considered at assembling.
    bool useGetBound;
666

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

670
671
672
    /// 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
673
    RefinementManager *refinementManager;
674

675
676
677
    /// 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
678
    CoarseningManager *coarseningManager;
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
    /// Info level.
    int info;
682

683
684
    /// If true, the stationary problem was deserialized from some serialization
    /// file.
685
686
    bool deserialized;

687
688
689
    /// 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.
690
    std::vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
691

692
693
    /// If true, the error is not estimated but computed from the exact solution
    /// defined by \ref exactSolutionFcts.
694
    bool computeExactError;
695
696

    /// If at least on boundary condition is set, this variable is true. It is
697
    /// used to ensure that no operators are added after boundary condition were
698
    /// set. If this would happen, boundary conditions could set wrong on off
699
    /// diagonal matrices.
700
    bool boundaryConditionSet;
701

702
    /// If true, AMDiS prints information about the assembling procedure to
703
    /// the screen.
704
    bool writeAsmInfo;
705

706
    std::map<Operator*, std::vector<OperatorPos> > operators;
707

708
    /// time needed to solve the linear system
709
    double solutionTime;
710

711
    /// time needed to assemble the linear system
712
    double buildTime;
713

714
    template <class> friend class detail::CouplingProblemStat;
715
  };
716

717
  namespace detail
718
  {
719
    template <class ProblemStatType>
720
    struct ProblemStat : public ProblemStatType,
721
			 public StandardProblemIteration
722
    {
723
      using ProblemStatType::getName;
724
      using ProblemStatType::info;
725

726
727
728
729
730
731
732
733
734
      /// Constructor
      ProblemStat(std::string nameStr,
		  ProblemIterationInterface *problemIteration = NULL)
	: ProblemStatType(nameStr, problemIteration),
	  StandardProblemIteration(this)
      { }

      /// Returns number of managed problems
      // implements StandardProblemIteration::getNumProblems()
Praetorius, Simon's avatar
Praetorius, Simon committed
735
      virtual int getNumProblems() const override
736
737
      {
	return 1;
738
739
740
741
742
743
      }

      /// 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
744
745
      {
	return this;