ProblemStat.hpp 8.56 KB
Newer Older
1
2
3
4
5
6
7
8
9
#pragma once

#include <tuple>
#include <string>
#include <memory>
#include <list>
#include <map>

#include <dune/istl/bcrsmatrix.hh>
10
11
12
#include <dune/istl/bvector.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
13
14
15
16

#include <dune/grid/common/grid.hh>

#include "Basic.hpp"
17
#include "DirichletBC.hpp"
18
#include "FileWriter.hpp"
19
#include "Flag.hpp"
20
21
#include "Initfile.hpp"
#include "Operator.hpp"
22
#include "Preconditioner.hpp"
23
#include "ProblemStatBase.hpp"
24
#include "SystemMatrix.hpp"
25
#include "SystemVector.hpp"
26
27
#include "Timer.hpp"

Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
28
29
namespace AMDiS {
    
30
31
    

32
template <class Traits>
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
33
class ProblemStat : public ProblemStatBase
34
{
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
35
public:
36
      
37
38
    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
39
40
41
42
43
    using MeshView = typename Mesh::LeafGridView;
    
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

44
45
    using ElementVector = Dune::BlockVector<Dune::FieldVector<double,1> >;
    using ElementMatrix = Dune::Matrix<Dune::FieldMatrix<double,1,1> >;
46
47
48
    
    
    /// Dimension of the mesh
49
    static constexpr int dim = Traits::dim;
50
51
    
    /// Dimension of the world
52
    static constexpr int dow = Traits::dimworld;
53
54
    
    /// Number of problem components
55
    static constexpr int nComponents = Traits::nComponents;
56
    
57
    
58
59
60
61
62
    template <size_t I>
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
    
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
    
63
64
    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;
65
66
67
    
    using OperatorType = Operator<MeshView>;
    
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
68
69
70
public:
    /// \brief Constructor. Takes the name of the problem that is used to 
    /// access values correpsonding to this püroblem in the parameter file.
71
72
73
74
75
76
77
    /**
     * Parameters read by ProblemStatSeq, with name 'PROB'
     *   PROB->names:                  \ref componentNames
     **/
    ProblemStat(std::string name)
      : name(name)
    {
78
79
80
	Parameters::get(name + "->names", componentNames);
	for (int i = componentNames.size(); i < nComponents; ++i)
	    componentNames.push_back("solution_" + std::to_string(i));
81
82
83
    }
    
    
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
84
    /// \brief Initialisation of the problem.
85
86
87
88
    /**
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
89
90
91
    void initialize(Flag initFlag, 
		    ProblemStat* adoptProblem = NULL, 
		    Flag adoptFlag = INIT_NOTHING);
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
    

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

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

    /// Adds a Dirichlet boundary condition
    template <class Predicate, class Values>
    void addDirichletBC(Predicate const& predicate, 
			int row, int col, 
			Values const& values);

    
    /// Implementation of \ref ProblemStatBase::solve
    virtual void solve(AdaptInfo& adaptInfo, 
114
115
		       bool createMatrixData = true, 
		       bool storeMatrixData = false) override;
116
117
118
119
120
		       
    using BasePreconditioner = Dune::Preconditioner<typename SystemVectorType::MultiVector,
						    typename SystemVectorType::MultiVector>;
    virtual void solveAdvanced(AdaptInfo& adaptInfo, 
			       BasePreconditioner& precon);
121
122
123
124
125
    
		       
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo, 
				   Flag flag, 
126
127
128
129
130
131
				   bool asmMatrix = true, 
				   bool asmVector = true) override;
    

    /// Writes output files.
    void writeFiles(AdaptInfo& adaptInfo, double time = -1.0);
132
133
134
    

    /// Returns nr of components \ref nComponents
135
    static constexpr size_t getNumComponents()
136
137
138
139
140
141
142
    {
      return nComponents;
    }
    
    
    /// Return the I'th finite element space
    template <size_t I = 0>
143
    FeSpace<I> const& getFeSpace(const index_<I> = index_<I>()) const
144
145
146
147
148
149
    {
	return std::get<I>(*feSpaces);
    }
    
    
    std::shared_ptr<FeSpaces> getFeSpaces()
150
    {
151
	return feSpaces;
152
153
154
155
    }

    
    /// Return the mesh
156
    std::shared_ptr<Mesh> getMesh()
157
    {
158
	return mesh;
159
160
161
162
    }
    
    
    /// Return the i'th solution component
163
    template <size_t I = 0>
164
    decltype(auto) getSolution(const index_<I> _i = index_<I>())
165
    {
166
	return (*solution)[_i];
167
168
169
    }
    
    
170
171
172
173
174
175
176
    decltype(auto) getSolution()	{ return *solution; }
    decltype(auto) getSolution() const	{ return *solution; }
    
    decltype(auto) getSystemMatrix()	   { return *systemMatrix; }
    decltype(auto) getSystemMatrix() const { return *systemMatrix; }
    
    
177
178
179
180
181
182
    /// Return the name of the problem. Implementation of \ref ProblemStatBase::getName
    virtual std::string getName() const override
    {
      return name;
    }
    
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
183
public:
184
185
186
187
    // Reads a filename from the init file and creates the mesh object with 
    // respect to this file.
    void createMesh()
    {
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
188
189
190
191
192
193
194
195
	// TODO: Creation from meshname must be generalized to other meshes than AlbertaGrid.
	meshName = "";
	Parameters::get(name + "->mesh", meshName);
	AMDIS_TEST_EXIT(!meshName.empty(),
	    "No mesh name specified for '" << name << "->mesh'!");
	
	std::string macro_filename = "";
	Parameters::get(meshName + "->macro file name", macro_filename);
196

Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
197
198
	mesh = std::make_shared<Mesh>(macro_filename);
	meshView = std::make_shared<MeshView>(mesh->leafGridView());
199
200
201
202
203
    }
    
    //
    void createFeSpaces()
    {
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
204
205
	feSpaces = std::make_shared<FeSpaces>(construct_tuple<FeSpaces>(*meshView));
	
206
	systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
207
	solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
208
	rhs = std::make_shared<SystemVectorType>(*feSpaces, std::vector<std::string>(nComponents, "rhs"));
209
210
211
212
213
214
    }
    
    //
    template <class RowBasis, class ColBasis>
    void getOccupationPattern(const RowBasis& rowBasis, 
			      const ColBasis& colBasis, 
215
			      Dune::MatrixIndexSet& nb);
216
217
218
    
    //
    template <class RowView, class ColView>
219
    bool getElementMatrix(RowView const& rowView,
220
221
222
223
224
225
			  ColView const& colView,
			  ElementMatrix& elementMatrix,
			  std::list<std::shared_ptr<OperatorType>>& operators);
    
    //
    template <class RowView>
226
    bool getElementVector(RowView const& rowView,
227
228
229
			  ElementVector& elementVector,
			  std::list<std::shared_ptr<OperatorType>>& operators);
    
230
231
232
233
234
235
236
237
238
239
240
241
    
    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
			  RowIndexSet const& rowIndexSet,
			  ColIndexSet const& colIndexSet,
			  ElementMatrix const& elementMatrix);
    
    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
			  IndexSet const& indexSet,
			  ElementVector const& elementvector);
    
242
    //
243
    template <class RowBasis, class ColBasis, class Matrix, class Vector>
244
245
246
    void assemble(std::pair<int, int> row_col,
		  const RowBasis& rowBasis, 
		  const ColBasis& colBasis, 
247
		  Matrix* matrix,
248
		  Vector* rhs);
249
    
250
251
252
253
254
255
256
257
private:
    template <size_t I = 0>
    typename FeSpace<I>::NodeFactory& 
    getNodeFactory(const index_<I> = index_<I>())
    {
	return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
    }
    
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
258
private:
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
    /// Name of this problem.
    std::string name;

    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
    std::vector<std::string> componentNames;
    
    /// Mesh of this problem.
    // TODO: generalize to multi-mesh problems
    std::shared_ptr<Mesh> mesh;
    std::string meshName;
    
    ///
    std::shared_ptr<MeshView> meshView;
    
    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;
    
    /// FE spaces of this problem.
    std::shared_ptr<FeSpaces> feSpaces; // eventuell const
    
280
281
    std::shared_ptr<FileWriter<Traits>> filewriter;
    
282
    std::shared_ptr<SystemMatrixType> systemMatrix;
283
284
    std::shared_ptr<SystemVectorType> solution;
    std::shared_ptr<SystemVectorType> rhs;
285
286
287
288
    
    IdxPairList<DirichletBC<WorldVector>>	dirichletBc;
    IdxPairList<OperatorType>			matrixOperators;
    IdxList<OperatorType>			vectorOperators;
Praetorius, Simon's avatar
cleanup    
Praetorius, Simon committed
289
};
290
  
291
292
  
#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
293
294
295
296
//   extern template class ProblemStat<ProblemStatTraits<1>>;
  extern template class ProblemStat<ProblemStatTraits<2,2,1>>;
  extern template class ProblemStat<ProblemStatTraits<2,2,2>>;
//   extern template class ProblemStat<ProblemStatTraits<3>>;
297
298
#endif
  
299
300
301
302
}

#include "ProblemStat.inc.hpp"