ProblemStat.hpp 6.91 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
#pragma once

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

#include <dune/istl/matrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>

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

#include <dune/functions/functionspacebases/pqknodalbasis.hh>

#include "Basic.hpp"
#include "Flag.hpp"
#include "Timer.hpp"

namespace AMDiS
{  
  // forward declarations
  class Operator;
  
  
  template <int _dim, int _dimworld = _dim, int _degree0 = 1, int... _degrees>
  struct ProblemParametersBase
  {
    static constexpr int dim = _dim;
    static constexpr int dimworld = _dimworld;
    static constexpr int nComponents = sizeof...(_degrees) + 1;
    
    template <class... Bs>
    using FeSpaceTuple = std::tuple<Bs...>;
    
    template <class Mesh, int deg>
    using Lagrange = Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
    
    // default values
    using Mesh = AlbertaGrid<dim, dimworld>;
    using FeSpaces = FeSpaceTuple<Lagrange<Mesh, _degree0>, 
				  Lagrange<Mesh, _degrees>...>;
  };
  

  template <class Param>
  class ProblemStat : public ProblemStatBase
  {
    using FeSpaces = typename Param::FeSpaces;
    using Mesh     = typename Param::Mesh;
    using MeshView = typename Mesh::LeafGridView;
    
    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

    using ElementVector = BlockVector<FieldVector<double,1> >;
    using ElementMatrix = Matrix<FieldMatrix<double,1,1> >;
    
    
    /// Dimension of the mesh
    static constexpr int dim = Param::dim;
    
    /// Dimension of the world
    static constexpr int dow = Param::dimworld;
    
    /// Number of problem components
    static constexpr int nComponents = Param::nComponents;
    
  public:
    template <size_t I>
    using FeSpace = std::tuple_element_t<I, FeSpaces>;
    
    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
    
    using DOFVector = BlockVector<FieldVector<double,1> >;
    using DOFMatrix = BCRSMatrix<FieldMatrix<double,1,1> >;
    
    using SystemVector = FieldVector<DOFVector*, nComponents>;
    using SystemMatrix = FieldMatrix<DOFMatrix*, nComponents, nComponents>;
    
    using OperatorType = Operator<MeshView>;
    
  public:
    /// Constructor
    /**
     * Parameters read by ProblemStatSeq, with name 'PROB'
     *   PROB->names:                  \ref componentNames
     **/
    ProblemStat(std::string name)
      : name(name)
    {
      Parameters::get(name + "->names", componentNames);
    }
    
    
    /// Initialisation of the problem.
    /**
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
    void initialize(Flag initFlag, ProblemStat* adoptProblem, Flag adoptFlag);
    

    /// 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>
    // requires( Dune::Functions::Concept::isCallable<Predicate, WorldVector>() &&
    // 	   	 Dune::Functions::Concept::isCallable<Values, WorldVector>() )
    void addDirichletBC(Predicate const& predicate, 
			int row, int col, 
			Values const& values);

    
    /// Implementation of \ref ProblemStatBase::solve
    virtual void solve(AdaptInfo& adaptInfo, 
		       bool createMatrixData, 
		       bool storeMatrixData) override;
    
		       
    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo, 
				   Flag flag, 
				   bool asmMatrix, 
				   bool asmVector) override;
    

    /// Returns nr of components \ref nComponents
    constexpr int getNumComponents() const
    {
      return nComponents;
    }
    
    
    /// Return the I'th finite element space
    template <size_t I = 0>
    FeSpace<I> const& getFeSpace() const
    {
      return std::get<I>(feSpaces);
    }

    
    /// Return the mesh
    Mesh const& getMesh() const
    {
      return *mesh;
    }
    
    
    /// Return the i'th solution component
    DOFVector& getSolution(int i)
    {
      return *solution[i];
    }
    
    
    /// Return the name of the problem. Implementation of \ref ProblemStatBase::getName
    virtual std::string getName() const override
    {
      return name;
    }
    
  protected:
    // Reads a filename from the init file and creates the mesh object with 
    // respect to this file.
    void createMesh()
    {
      // TODO: Creation from meshname must be generalized to other meshes than AlbertaGrid.
      meshName = "";
      Parameters::get(name + "->mesh", meshName);
      TEST_EXIT(!meshName.empty(),
	"No mesh name specified for '" << name << "->mesh'!\n");

      mesh = std::make_shared<Mesh>(meshName);
      meshView = std::make_shared<MeshView>(mesh->leafGridView());
    }
    
    //
    void createFeSpaces()
    {
      feSpaces = std::make_shared<FeSpaces>(construct_tuple<FeSpaces>(meshView));
    }
    
    //
    template <class RowBasis, class ColBasis>
    void getOccupationPattern(const RowBasis& rowBasis, 
			      const ColBasis& colBasis, 
			      MatrixIndexSet& nb);
    
    //
    template <class RowView, class ColView>
    void getElementMatrix(RowView const& rowView,
			  ColView const& colView,
			  ElementMatrix& elementMatrix,
			  std::list<std::shared_ptr<OperatorType>>& operators);
    
    //
    template <class RowView>
    void getElementVector(RowView const& rowView,
			  ElementVector& elementVector,
			  std::list<std::shared_ptr<OperatorType>>& operators);
    
    //
    template <class RowBasis, class ColBasis>
    void assemble(std::pair<int, int> row_col,
		  const RowBasis& rowBasis, 
		  const ColBasis& colBasis, 
		  DOFMatrix* matrix,
		  DOFVector* rhs);
    
  private:
    /// 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
    
    SystemMatrix systemMatrix;
    SystemVector solution;
    SystemVector rhs;
    
    IdxPairList<DirichletBC<WorldVector>>	dirichletBc;
    IdxPairList<OperatorType>			matrixOperators;
    IdxList<OperatorType>			vectorOperators;
  };
  
}

#include "ProblemStat.inc.hpp"