Call leafView() only once; pass gridview only

src/one-body-sample.cc

@@ -40,25 +40,23 @@
 
 int const dim = 2;
 
-template <class GridType>
-void setup_boundary(GridType const &grid,
+template <class GridView>
+void setup_boundary(GridView const &gridView,
                     Dune::FieldVector<double, dim> const &end_points,
                     Dune::BitSetVector<dim> &ignoreNodes,
                     Dune::BitSetVector<1> &neumannNodes,
                     Dune::BitSetVector<1> &frictionalNodes) {
-  typedef typename GridType::LeafGridView GridView;
-  GridView const leafView = grid.leafView();
   typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;
 
   typedef Dune::MultipleCodimMultipleGeomTypeMapper<
       GridView, Dune::MCMGVertexLayout> VertexMapper;
-  VertexMapper const myVertexMapper(leafView);
+  VertexMapper const myVertexMapper(gridView);
 
   size_t dirichlet_nodes = 0;
   size_t neumann_nodes = 0;
   size_t frictional_nodes = 0;
-  for (VertexLeafIterator it = leafView.template begin<dim>();
-       it != leafView.template end<dim>(); ++it) {
+  for (VertexLeafIterator it = gridView.template begin<dim>();
+       it != gridView.template end<dim>(); ++it) {
     Dune::GeometryType const gt = it->type();
     assert(it->geometry().corners() == 1);
     Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
@@ -136,7 +134,7 @@ int main() {
     size_t const solver_maxIterations = 100;
     double const solver_tolerance = 1e-4;
 
-    // Set up grid
+    // {{{ Set up grid
     typedef Dune::YaspGrid<dim> GridType;
     Dune::FieldVector<double, dim> const end_points(
         1); // nth dimension (zero-indexed) goes from 0 to end_points[n]
@@ -147,9 +145,13 @@ int main() {
         0);                                  // zero overlap (whatever that is)
     grid.globalRefine(refinements);
 
+    typedef GridType::LeafGridView GridView;
+    GridView const leafView = grid.leafView();
+    // }}}
+
     // Set up nodal basis
     typedef P1NodalBasis<GridType::LeafGridView, double> P1Basis;
-    P1Basis const p1Basis(grid.leafView());
+    P1Basis const p1Basis(leafView);
 
     // Assemble elastic force on the body
     StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
@@ -163,8 +165,8 @@ int main() {
     Dune::BitSetVector<1> neumannNodes(grid.size(grid.maxLevel(), dim), false);
     Dune::BitSetVector<1> frictionalNodes(grid.size(grid.maxLevel(), dim),
                                           false);
-    setup_boundary<GridType>(grid, end_points, ignoreNodes, neumannNodes,
-                             frictionalNodes);
+    setup_boundary<GridType::LeafGridView>(leafView, end_points, ignoreNodes,
+                                           neumannNodes, frictionalNodes);
 
     VectorType u1(grid.size(grid.maxLevel(), dim));
     u1 = 0;
@@ -172,12 +174,12 @@ int main() {
 
     VectorType f(grid.size(grid.maxLevel(), dim));
     assemble_neumann<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
-        grid.leafView(), p1Basis, neumannNodes, f);
+        leafView, p1Basis, neumannNodes, f);
 
     // {{{ Assemble terms for the nonlinearity
     std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
     assemble_frictional<GridType, GridType::LeafGridView, SmallVector, P1Basis>(
-        grid.leafView(), p1Basis, frictionalNodes, nodalIntegrals);
+        leafView, p1Basis, frictionalNodes, nodalIntegrals);
 
     // TODO: populate on S_F
     std::vector<double> normalStress;