Proper boundary handling

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");
     }