diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index a410cc0351b056caa1da8a15008dd8a8776244d5..e576d2df03fa4f735c70beb067c2523db500a4bd 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -18,21 +18,23 @@ class DG public: using Scalar = typename GridView::ctype; static constexpr int dim = GridView::dimension; + using Matrix = Dune::DynamicMatrix<Scalar>; + using Vector = Dune::DynamicVector<Scalar>; //constructor DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) : - gridView_(gridView), mapper_(mapper), problem_(problem) {} + gridView_(gridView), mapper_(mapper), problem_(problem), + dof((1 + dim) * gridView.size(0)) + { + A = Matrix(dof, dof, 0.0); //initialize stiffness matrix + b = Vector(dof, 0.0); //initialize load vector + d = Vector(dof, 0.0); //initialize solution vector + } - const void operator() (const Scalar K, const Scalar mu) const + const void operator() (const Scalar K, const Scalar mu) { //NOTE: what is an appropriate sparse matrix type? -> BCRS - const int dof = (1 + dim)*gridView_.size(0); //degrees of freedom - using Matrix = Dune::DynamicMatrix<Scalar>; - using Vector = Dune::DynamicVector<Scalar>; - Matrix A(dof, dof, 0.0); //stiffness matrix - Vector b(dof, 0.0); //load vector - Vector d(dof, 0.0); //NOTE: // const int order = 3; //order of the quadrature rule @@ -122,7 +124,7 @@ public: { std::cout << "\t" << intersctGeo.corner(k); } - std::cout << "\n\n"; + std::cout << "\tnormal: " << normal << "\n\n"; //TODO: quadrature rule cannot be used for dim = 1! // const Dune::QuadratureRule<double,dim-1>& secondOrderRule = @@ -480,10 +482,16 @@ public: } const GridView& gridView_; - const Mapper& mapper_; //element mapper for gridView_ - const Problem& problem_; //the DG problem to be solved + + const int dof; //degrees of freedom + +private: + Matrix A; //stiffness matrix + Vector b; //load vector + Vector d; //solution vector + }; #endif