[Cleanup] Construct boundary from faces

src/one-body-sample.cc

@@ -137,42 +137,52 @@ int main(int argc, char *argv[]) {
     GridView const leafView = grid->leafView();
     // }}}
 
-    // Set up the boundary
-    Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
-    Dune::BitSetVector<dims> noNodes(fineVertexCount);
-    BoundaryPatch<GridView> const neumannBoundary(leafView);
+    // Set up faces
+    BoundaryPatch<GridView> lowerFace(leafView);
+    BoundaryPatch<GridView> upperFace(leafView);
+    {
+      auto const isClose = [](double a,
+                              double b) { return std::abs(a - b) < 1e-14; };
+      lowerFace.insertFacesByProperty([&](
+          typename GridView::Intersection const &in) {
+        return isClose(lowerLeft[1], in.geometry().center()[1]);
+      });
+
+      upperFace.insertFacesByProperty([&](
+          typename GridView::Intersection const &in) {
+        return isClose(upperRight[1], in.geometry().center()[1]);
+      });
+    }
 
     Vector vertexCoordinates(fineVertexCount);
-    Dune::BitSetVector<1> frictionalNodes(fineVertexCount, false);
     {
       Dune::MultipleCodimMultipleGeomTypeMapper<
           GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
       for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
-        assert(it->geometry().corners() == 1);
-        size_t const id = vertexMapper.map(*it);
-        vertexCoordinates[id] = it->geometry().corner(0);
-        auto const &localCoordinates = vertexCoordinates[id];
-
-        // lower face
-        if (localCoordinates[1] == lowerLeft[1]) {
-          frictionalNodes[id] = true;
-          velocityDirichletNodes[id][1] = true;
-        }
+        auto const geometry = it->geometry();
+        assert(geometry.corners() == 1);
+        vertexCoordinates[vertexMapper.map(*it)] = geometry.corner(0);
+      }
+    }
 
-        // upper face
-        else if (localCoordinates[1] == upperRight[1])
-          velocityDirichletNodes[id] = true;
+    // Neumann boundary
+    BoundaryPatch<GridView> const neumannBoundary(leafView);
 
-        // right face (except for both corners)
-        else if (localCoordinates[0] == upperRight[0])
-          ;
+    // Frictional Boundary
+    BoundaryPatch<GridView> const &frictionalBoundary = lowerFace;
+    Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
+    frictionalBoundary.getVertices(frictionalNodes);
 
-        // left face (except for both corners)
-        else if (localCoordinates[0] == lowerLeft[0])
-          ;
-      }
+    // Dirichlet Boundary
+    Dune::BitSetVector<dims> noNodes(fineVertexCount);
+    Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
+    for (size_t i = 0; i < fineVertexCount; ++i) {
+      if (lowerFace.containsVertex(i))
+        dirichletNodes[i][1] = true;
+
+      if (upperFace.containsVertex(i))
+        dirichletNodes[i] = true;
     }
-    BoundaryPatch<GridView> const frictionalBoundary(leafView, frictionalNodes);
 
     // Set up functions for time-dependent boundary conditions
     using Function = Dune::VirtualFunction<double, double>;