Upper face: Neumann. Lower face: Dirichlet

src/one-body-sample.cc

@@ -90,6 +90,7 @@ int main() {
 
     Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(
         grid.size(grid.maxLevel(), dim), false);
+    Dune::BitSetVector<1> neumannNodes(grid.size(grid.maxLevel(), dim), false);
     { // Play around with the boundary
       typedef GridType::LeafGridView GridView;
       GridView leafView = grid.leafView();
@@ -99,6 +100,8 @@ int main() {
           GridView, Dune::MCMGVertexLayout> VertexMapper;
       VertexMapper myVertexMapper(leafView);
 
+      size_t dirichlet_nodes = 0;
+      size_t neumann_nodes = 0;
       size_t bounding_nodes = 0;
       size_t extremal_nodes = 0;
       for (VertexLeafIterator it = leafView.begin<dim>();
@@ -120,6 +123,16 @@ int main() {
             break;
           }
         }
+        if (coordinates[0] == 0) {
+          ++dirichlet_nodes;
+          size_t id = myVertexMapper.map(*it);
+          ignoreNodes[id] = true;
+        }
+        if (coordinates[0] == end_points[0]) {
+          ++neumann_nodes;
+          size_t id = myVertexMapper.map(*it);
+          neumannNodes[id] = true;
+        }
         if (bounding) {
           ++bounding_nodes;
           size_t id = myVertexMapper.map(*it);
@@ -133,6 +146,9 @@ int main() {
           // ignoreNodes[id] = true;
         }
       }
+      std::cout << "Number of Neumann nodes: " << neumann_nodes << std::endl;
+      std::cout << "Number of Dirichlet nodes: " << dirichlet_nodes
+                << std::endl;
       std::cout << "Number of bounding nodes: " << bounding_nodes << std::endl;
       std::cout << "Number of extremal nodes: " << extremal_nodes << std::endl;
     }