diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index 3ed8832755acdddc5f896845e8184eda3e8dc103..73835628eec65a9616fc96423d9c4b57f6b56b9e 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -26,13 +26,13 @@ public: const void operator() (const Scalar K, const Scalar mu) const { - //NOTE: what is an appropriate sparse matrix type? + //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; + Vector d(dof, 0.0); //NOTE: // const int order = 3; //order of the quadrature rule @@ -311,6 +311,8 @@ public: } } + std::cout << A << std::endl; + //NOTE: what would be an appropiate solver here? A.solve(d, b); @@ -348,11 +350,28 @@ public: } } + Dune::DynamicVector<Scalar> + pressureTEST(gridView_.size(0), 0.0); + + for (const auto& elem : elements(gridView_)) + { + const int elemIdx = mapper_.index(elem); + const int elemIdxSLE = (dim + 1)*elemIdx; + const auto& center = elem.geometry().center(); + pressureTEST[elemIdx] = d[elemIdxSLE]; + + for (int i = 0; i < dim; i++) + { + pressureTEST[elemIdx] += d[elemIdxSLE + i + 1] * center[i]; + } + } + + //NOTE: dim = 1 - auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);}; +// auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);}; //NOTE: dim = 2 -/* auto exactSolution = [](const auto& pos) + auto exactSolution = [](const auto& pos) { double solution = 0.0; for (int m = 1; m <= 21; m = m + 2) @@ -366,7 +385,7 @@ public: solution *= -16 / pow(M_PI, 4); return solution; }; -*/ + Dune::DynamicVector<double> exactPressure(gridView_.size(0)); for (const auto& elem : elements(gridView_)) @@ -379,10 +398,10 @@ public: //create output directory if necessary mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); - auto vtkWriter = - std::make_shared<Dune::VTKWriter<GridView>>(gridView_); - Dune::VTKSequenceWriter<GridView> - vtkSqWriter(vtkWriter, "pressure", "data", "data"); + Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming); + // std::make_shared<Dune::VTKWriter<GridView>>(gridView_); +// Dune::VTKSequenceWriter<GridView> +// vtkSqWriter(vtkWriter, "pressure", "data", "data"); using VTKFunction = Dune::VTKFunction<GridView>; using P1Function = Dune::NonconformingP1VTKFunction<GridView, @@ -391,10 +410,11 @@ public: new P1Function(gridView_, pressure, "pressure"); //NOTE: why does this only work with a sequence writer? - vtkSqWriter.addVertexData( + vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(vtkOperator) ); - vtkSqWriter.addCellData(exactPressure, "exactPressure"); - vtkSqWriter.write(0.0); + vtkWriter.addCellData(exactPressure, "exactPressure"); + vtkWriter.addCellData(pressureTEST, "pressureTEST"); + vtkWriter.write("linearpressure"); } const GridView& gridView_;