diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index 0b8059cd4f76bacb82cb79ad620d850636ca24a8..f0f3495f52760ba457e04954462a5f2bc8ee7d21 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -17,7 +17,7 @@ #include <dune/mmdg/nonconformingp1vtkfunction.hh> -template<class Grid, class GridView, class Mapper, class Problem> +template<class GridView, class Mapper, class Problem> class DG { public: @@ -37,10 +37,8 @@ public: using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>; //constructor - DG (const Grid& grid, const GridView& gridView, const Mapper& mapper, const Problem& problem) : - //NOTE: why does this not work? - //grid_(grid), gridView_(grid.leafGridView()), mapper_(*(new Mapper(gridView_, Dune::mcmgElementLayout()))), problem_(problem), dof_((1 + dim) * gridView_.size(0)) - grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem), + DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) : + gridView_(gridView), mapper_(mapper), problem_(problem), dof_((1 + dim) * gridView.size(0)) { //initialize stiffness matrix A, load vector b and solution vector d @@ -106,7 +104,6 @@ public: return std::sqrt(error); } - const Grid& grid_; const GridView& gridView_; const Mapper& mapper_; //element mapper for gridView_ const Problem& problem_; //the DG problem to be solved @@ -114,10 +111,9 @@ public: const int dof_; //degrees of freedom protected: - DG (const Grid& grid, const GridView& gridView, const Mapper& mapper, - const Problem& problem, const int dof) : - grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem), - dof_(dof) + DG (const GridView& gridView, const Mapper& mapper, const Problem& problem, + const int dof) : + gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof) { //initialize stiffness matrix A, load vector b and solution vector d A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO @@ -205,7 +201,7 @@ protected: { if constexpr (dim != 1) { - if (grid_.isInterface(intersct)) + if (gridView_.grid().isInterface(intersct)) { continue; } diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index 62d2d0b0c325eeec680f291656c60b170ff1c9e5..ccaef4d5945090dba159bf9a51f19f741a2d4813 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -3,15 +3,15 @@ #include <dune/mmdg/dg.hh> -template<class Grid, class GridView, class Mapper, class InterfaceGridView, +template<class GridView, class Mapper, class InterfaceGridView, class InterfaceMapper, class Problem> -class MMDG : public DG<Grid, GridView, Mapper, Problem> +class MMDG : public DG<GridView, Mapper, Problem> { public: static constexpr int dim = GridView::dimension; static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1); - using Base = DG<Grid, GridView, Mapper, Problem>; + using Base = DG<GridView, Mapper, Problem>; using Scalar = typename Base::Scalar; //NOTE using Base::Scalar using Matrix = typename Base::Matrix; using Vector = typename Base::Vector; @@ -24,10 +24,10 @@ public: Dune::DynamicVector<Scalar>>; //constructor - MMDG (const Grid& grid, const GridView& gridView, const Mapper& mapper, + MMDG (const GridView& gridView, const Mapper& mapper, const InterfaceGridView& iGridView, const InterfaceMapper& iMapper, const Problem& problem) : - Base(grid, gridView, mapper, problem, + Base(gridView, mapper, problem, (1 + dim) * gridView.size(0) + dim * iGridView.size(0)), iGridView_(iGridView), iMapper_(iMapper) {} @@ -59,7 +59,7 @@ public: const QRuleInterface& qrule = QRulesInterface::rule(iSimplex, quadratureOrder); const InterfaceBasis iFrame = interfaceFrame( - Base::grid_.asIntersection(iElem).centerUnitOuterNormal()); + Base::gridView_.grid().asIntersection(iElem).centerUnitOuterNormal()); //determine L2-error on the interface element iElem using a //quadrature rule @@ -147,7 +147,7 @@ private: { const int iElemIdx = iMapper_.index(iElem); const auto& iGeo = iElem.geometry(); - const auto& intersct = Base::grid_.asIntersection(iElem); + const auto& intersct = Base::gridView_.grid().asIntersection(iElem); //get basis of the interface coordinate system for evaluation of //basis function phi_iElem,i with i=1,...,dim-1 @@ -667,7 +667,7 @@ private: //get basis of the interface coordinate system for evaluation of //basis function phi_iElem,i with i=1,...,dim-1 const InterfaceBasis iFrame = interfaceFrame( - Base::grid_.asIntersection(iElem).centerUnitOuterNormal()); + Base::gridView_.grid().asIntersection(iElem).centerUnitOuterNormal()); for (int k = 0; k < iGeo.corners(); k++) { diff --git a/src/convergenceTest.py b/src/convergenceTest.py index 20488584f11ec2e91b2a03284aaf3d7a02c26c4b..92ddfa4ad6d118cb0599ee057cfe85acbd450e79 100644 --- a/src/convergenceTest.py +++ b/src/convergenceTest.py @@ -31,7 +31,7 @@ def writeDGF(N, dim): file.close() -# reads error, run time and number of elements, i.e. the output from the +# reads error, run time and maximum edge length, i.e. the output from the # excuted discontinuous Galerkin scheme def readData(): file = open("convergenceData", "r") @@ -54,12 +54,12 @@ if not exists("plots"): mkdir("plots") # number of runs from which the run times are taken as mean value -numberOfMeans = 3 +numberOfMeans = 1 # numbers of grid elements per direction numElemsPerDir = [1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130] -dimensions = [1,2] +dimensions = [1, 2] for dim in dimensions: # storage for errors, run times and largestEdgeLength diff --git a/src/dg.cc b/src/dg.cc index 60d8760f522eb8dfc2ad2f556b48cb37b913dd03..744b48cede991d649b5568e08eff1cceb0ff58ab 100644 --- a/src/dg.cc +++ b/src/dg.cc @@ -44,7 +44,7 @@ int main(int argc, char** argv) using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>; using Coordinate = Dune::FieldVector<double, dim>; using Problem = DGProblem<Coordinate>; - using DG = DG<Grid, GridView, Mapper, Problem>; + using DG = DG<GridView, Mapper, Problem>; //get parameters from file "parameterDG.ini" Dune::ParameterTree pt; @@ -83,8 +83,7 @@ int main(int argc, char** argv) "parameterDG.ini (use 'poisson' or TODO)"); } - //DG<GridView, Mapper, Problem> dg(gridView, mapper, *problem); - DG dg(grid, gridView, mapper, *problem); + DG dg(gridView, mapper, *problem); dg(mu); @@ -92,9 +91,6 @@ int main(int argc, char** argv) const double error = dg.computeL2error(); const double timeTotal = ClockHelper::elapsedRuntime(timeInitial); - //write error, runtime and number of grid elements to file - std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc); - //determine maximum edge length double largestEdgeLength = 0.; for (const auto& edge : edges(gridView)) @@ -102,9 +98,12 @@ int main(int argc, char** argv) largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume()); } + //write error, runtime and maximum edge length + std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc); + if (errorFile.is_open()) { - errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength;//gridView.size(0); + errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength; errorFile.close(); } else diff --git a/src/mmdg.cc b/src/mmdg.cc index 47c2e830efe66b890e75adcb269603fdfa4bb92f..2ed9097c7d5ef1f2c80fc36f169bd2220663a6ae 100644 --- a/src/mmdg.cc +++ b/src/mmdg.cc @@ -44,7 +44,7 @@ int main(int argc, char** argv) typename Dune::MultipleCodimMultipleGeomTypeMapper<IGridView>; using Coordinate = Dune::FieldVector<double, dim>; using Problem = CoupledDGProblem<Coordinate>; - using MMDG = MMDG<Grid, GridView, Mapper, IGridView, IMapper, Problem>; + using MMDG = MMDG<GridView, Mapper, IGridView, IMapper, Problem>; //get parameters from file "parameterDG.ini" Dune::ParameterTree pt; @@ -92,7 +92,7 @@ int main(int argc, char** argv) const Mapper mapper(gridView, Dune::mcmgElementLayout()); const IMapper iMapper(iGridView, Dune::mcmgElementLayout()); - MMDG mmdg(grid, gridView, mapper, iGridView, iMapper, *problem); + MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem); mmdg(mu, xi); //print total run time