Handle extremal nodes

src/one-body-sample.cc

@@ -43,7 +43,7 @@ int main() {
     double const nu = 0.3;
     int const refinements = 5;
     size_t const solver_maxIterations = 10000;
-    double const solver_tolerance = 1e-5;
+    double const solver_tolerance = 1e-4;
 
     typedef Dune::YaspGrid<dim> GridType;
 
@@ -98,26 +98,41 @@ int main() {
       VertexMapper myVertexMapper(leafView);
 
       size_t bounding_nodes = 0;
+      size_t extremal_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);
+        bool extremal(true);
         for (int i = 0; i < dim; ++i) {
           if (coordinates[i] == end_points[i] || coordinates[i] == 0) {
             bounding = true;
             break;
           }
         }
+        for (int i = 0; i < dim; ++i) {
+          if (coordinates[i] != end_points[i] && coordinates[i] != 0) {
+            extremal = false;
+            break;
+          }
+        }
         if (bounding) {
           ++bounding_nodes;
           size_t id = myVertexMapper.map(*it);
           std::cout << "Ignoring id #" << id << std::endl;
           ignoreNodes[id] = true;
         }
+        if (extremal) {
+          ++extremal_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;
+      std::cout << "Number of extremal nodes: " << extremal_nodes << std::endl;
     }
 
     blockGSStep.ignoreNodes_ = &ignoreNodes;