Only write out data for a central node

src/one-body-sample.cc

@@ -86,7 +86,8 @@ void setup_boundary(GridView const &gridView,
                     Dune::BitSetVector<dim> &ignoreNodes,
                     Dune::BitSetVector<1> &neumannNodes,
                     Dune::BitSetVector<1> &frictionalNodes,
-                    GridCorner const &lowerLeft, GridCorner const &upperRight) {
+                    GridCorner const &lowerLeft, GridCorner const &upperRight,
+                    size_t &specialNode) {
   typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
 
   Dune::MultipleCodimMultipleGeomTypeMapper<
@@ -97,6 +98,22 @@ void setup_boundary(GridView const &gridView,
     assert(it->geometry().corners() == 1);
     Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
     size_t const id = myVertexMapper.map(*it);
+
+    // Find the center of the lower face
+    switch (dim) {
+      case 3:
+        if (coordinates[2] != lowerLeft[2])
+          break;
+      // fallthrough
+      case 2:
+        if (coordinates[1] == lowerLeft[1] &&
+            std::abs(coordinates[0] - lowerLeft[0]) < 1e-8)
+          specialNode = id;
+        break;
+      default:
+        assert(false);
+    }
+
     if (coordinates[1] == upperRight[1]) // upper face
       ignoreNodes[id] = true;
     else if (coordinates[1] == lowerLeft[1]) { // lower face
@@ -234,11 +251,14 @@ int main(int argc, char *argv[]) {
     EnergyNorm<MatrixType, VectorType> energyNorm(stiffnessMatrix);
 
     // Set up the boundary
+    size_t specialNode = finestSize;
     Dune::BitSetVector<dim> ignoreNodes(finestSize, false);
     Dune::BitSetVector<1> neumannNodes(finestSize, false);
     Dune::BitSetVector<1> frictionalNodes(finestSize, false);
     setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes,
-                   lowerLeft, upperRight);
+                   lowerLeft, upperRight, specialNode);
+    assert(specialNode != finestSize);
+    assert(frictionalNodes[specialNode][0]);
 
     auto const nodalIntegrals =
         assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
@@ -425,26 +445,16 @@ int main(int argc, char *argv[]) {
       alpha_old = alpha;
 
       { // Write all kinds of data
-        for (size_t i = 0; i < frictionalNodes.size(); ++i)
-          if (frictionalNodes[i][0])
-            state_writer << alpha[i][0] << " ";
-        state_writer << std::endl;
-
-        for (size_t i = 0; i < frictionalNodes.size(); ++i)
-          if (frictionalNodes[i][0])
-            displacement_writer << u[i][0] << " ";
-        displacement_writer << std::endl;
-
-        for (size_t i = 0; i < frictionalNodes.size(); ++i)
-          if (frictionalNodes[i][0])
-            velocity_writer << ud[i][0] << " ";
-        velocity_writer << std::endl;
-
-        for (size_t i = 0; i < frictionalNodes.size(); ++i)
-          if (frictionalNodes[i][0])
-            coefficient_writer << mu + a *std::log(ud[i].two_norm() / V0) +
-                                      b * (alpha[i] + std::log(V0 / L)) << " ";
-        coefficient_writer << std::endl;
+        state_writer << alpha[specialNode][0] << " " << std::endl;
+
+        displacement_writer << u[specialNode][0] << " " << std::endl;
+
+        velocity_writer << ud[specialNode][0] << " " << std::endl;
+
+        coefficient_writer << mu +
+                                  a *std::log(ud[specialNode].two_norm() / V0) +
+                                  b * (alpha[specialNode] + std::log(V0 / L))
+                           << " " << std::endl;
       }
 
       // Compute von Mises stress and write everything to a file