Play around with gridviews and the dirichlet nodes

src/one-body-sample.cc

@@ -25,6 +25,8 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/common/numproc.hh> // Solver::FULL
 
+#include <dune/grid/common/mcmgmapper.hh>
+
 #include <exception>
 
 int const dim = 2;
@@ -83,7 +85,41 @@ int main() {
     BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f);
     EnergyNorm<OperatorType, VectorType> energyNorm(blockGSStep);
 
-    Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(false);
+    Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(
+        grid.size(grid.maxLevel(), dim), false);
+
+    { // Play around with the boundary
+      typedef GridType::LeafGridView GridView;
+      GridView leafView = grid.leafView();
+      typedef GridView::Codim<dim>::Iterator VertexLeafIterator;
+
+      typedef Dune::MultipleCodimMultipleGeomTypeMapper<
+          GridView, Dune::MCMGVertexLayout> VertexMapper;
+      VertexMapper myVertexMapper(leafView);
+
+      size_t bounding_nodes = 0;
+      for (VertexLeafIterator it = leafView.begin<dim>();
+           it != leafView.end<dim>(); ++it) {
+        Dune::GeometryType gt = it->type();
+        assert(it->geometry().corners() == 1);
+        SmallVector coordinates = it->geometry().corner(0);
+        bool bounding(false);
+        for (int i = 0; i < dim; ++i) {
+          if (coordinates[i] == end_points[i] || coordinates[i] == 0) {
+            bounding = true;
+            break;
+          }
+        }
+        if (bounding) {
+          ++bounding_nodes;
+          size_t id = myVertexMapper.map(*it);
+          std::cout << "Ignoring id #" << id << std::endl;
+          ignoreNodes[id] = true;
+        }
+      }
+      std::cout << "Number of bounding nodes: " << bounding_nodes << std::endl;
+    }
+
     blockGSStep.ignoreNodes_ = &ignoreNodes;
 
     LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations,