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
#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
46
47
48
49
50
  struct OperatorPos 
  {
    int row, col;
    double *factor, *estFactor;
  };

51

52
53
54
55
56
57
58
  /**
   * \ingroup Problem 
   *
   * \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
106
    /// 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.
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

146
147
    bool dualMeshTraverseRequired();

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

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

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

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

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

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

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

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

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

188
189
    /// Adds a Dirichlet boundary condition, where the rhs is given by an 
    /// abstract function.
190
    virtual void addDirichletBC(BoundaryType type, int row, int col,
191
				AbstractFunction<double, WorldVector<double> > *b);
192
193
				
				
Praetorius, Simon's avatar
Praetorius, Simon committed
194
#if HAS_CXX11
195
196
197
198
199
    /// 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
200

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

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

211
212
213
214
215
216
    /// 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.
217
218
219
220
    virtual void addRobinBC(BoundaryType type, int row, int col, 
			    AbstractFunction<double, WorldVector<double> > *n,
			    AbstractFunction<double, WorldVector<double> > *r);

221
222
223
224
225
226
227
228
229
230
    /// 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);

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

234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
    /// 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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
260
261

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

 
268
269
270
271
    /** \name getting methods
     * \{ 
     */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

416
417
418
419
420
421
    /** \} */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

      componentMeshes[comp] = mesh;
    }

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

    void setCoarseningManager(CoarseningManager *coarse)
    {
      coarseningManager = coarse;
    }
534
535
    /** \} */

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
564
  protected:
565
566
    /// 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
567
568
    void computeError(AdaptInfo *adaptInfo);

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

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

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

584
585
586
    /// 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
587
588
589
    int nMeshes;

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

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

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

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

601
602
603
    /// 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
604
605
606
    ComponentTraverseInfo traverseInfo;
    
    /// Responsible for element marking.
607
    std::vector<Marker*> marker;
608

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

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

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

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

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

624
625
626
    /// Composed system matrix
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;

627
628
629
630
631
    /// 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.
632
    std::vector<std::vector<bool> > assembleMatrixOnlyOnce;
633

634
635
636
637
    /// 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.
638
    std::vector<std::vector<bool> > assembledMatrix;
639

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

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

646
647
648
    /// 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
649
    RefinementManager *refinementManager;
650

651
652
653
    /// 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
654
    CoarseningManager *coarseningManager;
655
  
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
    /// Info level.
    int info;
658

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

663
664
665
    /// 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.
666
    std::vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
667

668
669
    /// If true, the error is not estimated but computed from the exact solution
    /// defined by \ref exactSolutionFcts.
670
    bool computeExactError;
671
    
672
673
674
675
    /// 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.
676
    bool boundaryConditionSet;
677

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

682
    std::map<Operator*, std::vector<OperatorPos> > operators;
683
    
684
    /// time needed to solve the linear system
685
    double solutionTime;
686
687
    
    /// time needed to assemble the linear system
688
    double buildTime;
689
690
    
    template <class> friend class detail::CouplingProblemStat;
691
  };
692
693
  
  namespace detail
694
  {
695
696
697
    template <class ProblemStatType>
    struct ProblemStat : public ProblemStatType, 
			 public StandardProblemIteration
698
    {
699
      using ProblemStatType::getName;
700
      using ProblemStatSeq::info;
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
      
      /// 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;
743
#endif
744
745
746
}

#endif