diff --git a/src/07-poisson-problem-functions.cc b/src/07-poisson-problem-functions.cc index f75882fd12900eb0d71e7e581072629a89f1a676..428a36a84186ac58efc180fcab76bda3fd0023ef 100644 --- a/src/07-poisson-problem-functions.cc +++ b/src/07-poisson-problem-functions.cc @@ -190,36 +190,50 @@ int main(int argc, char** argv) using Basis = typename Dune::Functions::LagrangeBasis<GridView,4>; Basis basis(gridView); + std::vector<bool> isBoundary; + isBoundary.resize(basis.size(), false); + computeBoundaryVertices(gridView, isBoundary, basis); + assemblePoissonProblem(gridView, A, rhs, rhsFunction, basis); - x.resize(basis.size(), 0); + x.resize(basis.size()); + x = 0; auto dirichletFunction = [](auto x) { return sin(x[0] * 2.0*M_PI)*.1; }; - Dune::Functions::interpolate(basis, x, dirichletFunction); - - std::vector<bool> isBoundary; - isBoundary.resize(basis.size(), false); - computeBoundaryVertices(gridView, isBoundary, basis); - - for(std::size_t i=0; i<rhs.size(); ++i) - { - if (isBoundary[i]) - { - for(auto& entry: A[i]) - entry = 0; - A[i][i] = 1; - rhs[i] = x[i]; + Dune::Functions::interpolate(basis, x, dirichletFunction, isBoundary); + + A.mmv(x, rhs); + for (std::size_t row=0; row<A.N(); row++) { + if (isBoundary[row]) { + rhs[row] = x[row]; + // loop over nonzero matrix entries in current row + auto columnIt = A[row].begin(); + auto columnEndIt = A[row].end(); + for (; columnIt!=columnEndIt; ++columnIt) + if (row==columnIt.index()) + *columnIt = 1.0; + else + *columnIt = 0.0; + } + else { + // loop over nonzero matrix entries in current row + auto columnIt = A[row].begin(); + auto columnEndIt = A[row].end(); + for (; columnIt!=columnEndIt; ++columnIt) + if (isBoundary[columnIt.index()]) + *columnIt = 0; } } + Dune::MatrixAdapter<Matrix,Vector,Vector> op(A); - Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,1.0); Dune::SeqSSOR<Matrix,Vector,Vector> preconditioner(A,1,1.0); +// Dune::SeqILDL<Matrix,Vector,Vector> preconditioner(A,1.0); Dune::CGSolver<Vector> cg(op, preconditioner, // preconditione - 1e-4, // desired residual reduction factor + 1e-10, // desired residual reduction factor 200, // maximum number of iterations 2); // verbosity of the solver @@ -236,7 +250,7 @@ int main(int argc, char** argv) vtkWriter.addVertexData(xFunction, Dune::VTK::FieldInfo("solution", Dune::VTK::FieldInfo::Type::scalar, 1)); vtkWriter.write("07-solution"); } { - Dune::SubsamplingVTKWriter<GridView> vtkWriter(gridView,4); + Dune::SubsamplingVTKWriter<GridView> vtkWriter(gridView,Dune::refinementLevels(4)); vtkWriter.addVertexData(xFunction, Dune::VTK::FieldInfo("solution", Dune::VTK::FieldInfo::Type::scalar, 1)); vtkWriter.write("07-solution-subsampled"); }