[Cleanup] Do everything on the leaf grid

src/sand-wedge.cc

@@ -111,10 +111,10 @@ int main(int argc, char *argv[]) {
 
     auto const refinements = parset.get<size_t>("grid.refinements");
     grid->globalRefine(refinements);
-    size_t const fineVertexCount = grid->size(grid->maxLevel(), dims);
 
     using GridView = Grid::LeafGridView;
     GridView const leafView = grid->leafView();
+    size_t const leafVertexCount = leafView.size(dims);
 
     auto myFaces = gridConstructor.constructFaces(leafView);
 
@@ -123,18 +123,18 @@ int main(int argc, char *argv[]) {
 
     // Frictional Boundary
     BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
+    Dune::BitSetVector<1> frictionalNodes(leafVertexCount);
     frictionalBoundary.getVertices(frictionalNodes);
 
     // Surface
     BoundaryPatch<GridView> const &surface = myFaces.upper;
-    Dune::BitSetVector<1> surfaceNodes(fineVertexCount);
+    Dune::BitSetVector<1> surfaceNodes(leafVertexCount);
     surface.getVertices(surfaceNodes);
 
     // Dirichlet Boundary
-    Dune::BitSetVector<dims> noNodes(fineVertexCount);
-    Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
-    for (size_t i = 0; i < fineVertexCount; ++i) {
+    Dune::BitSetVector<dims> noNodes(leafVertexCount);
+    Dune::BitSetVector<dims> dirichletNodes(leafVertexCount);
+    for (size_t i = 0; i < leafVertexCount; ++i) {
       if (myFaces.right.containsVertex(i))
         dirichletNodes[i][0] = true;
 
@@ -200,7 +200,7 @@ int main(int argc, char *argv[]) {
                                   _relativeTime);
       _ell += gravityFunctional;
     };
-    Vector ell(fineVertexCount);
+    Vector ell(leafVertexCount);
     computeExternalForces(0.0, ell);
 
     // {{{ Initial conditions
@@ -235,16 +235,16 @@ int main(int argc, char *argv[]) {
     };
 
     // Solve the stationary problem
-    Vector u_initial(fineVertexCount);
+    Vector u_initial(leafVertexCount);
     u_initial = 0.0;
     solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
                        parset.sub("u0.solver"));
-    Vector ur_initial(fineVertexCount);
+    Vector ur_initial(leafVertexCount);
     ur_initial = 0.0;
 
-    ScalarVector alpha_initial(fineVertexCount);
+    ScalarVector alpha_initial(leafVertexCount);
     alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
-    ScalarVector normalStress(fineVertexCount);
+    ScalarVector normalStress(leafVertexCount);
     myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
                                      body.getYoungModulus(),
                                      body.getPoissonRatio(), u_initial);
@@ -254,21 +254,21 @@ int main(int argc, char *argv[]) {
         frictionalBoundary, frictionInfo, normalStress);
     myGlobalFriction->updateAlpha(alpha_initial);
 
-    Vector v_initial(fineVertexCount);
+    Vector v_initial(leafVertexCount);
     v_initial = 0.0;
     {
       double v_initial_const;
       velocityDirichletFunction.evaluate(0.0, v_initial_const);
       assert(std::abs(v_initial_const) < 1e-14);
     }
-    Vector vr_initial(fineVertexCount);
+    Vector vr_initial(leafVertexCount);
     vr_initial = 0.0;
 
-    Vector a_initial(fineVertexCount);
+    Vector a_initial(leafVertexCount);
     a_initial = 0.0;
     {
       // We solve Ma = ell - [Au + Cv + Psi(v)]
-      Vector accelerationRHS(fineVertexCount);
+      Vector accelerationRHS(leafVertexCount);
       {
         accelerationRHS = 0.0;
         Arithmetic::addProduct(accelerationRHS, A, u_initial);
@@ -283,7 +283,7 @@ int main(int argc, char *argv[]) {
     }
     // }}}
 
-    Vector vertexCoordinates(fineVertexCount);
+    Vector vertexCoordinates(leafVertexCount);
     {
       Dune::MultipleCodimMultipleGeomTypeMapper<
           GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
@@ -354,8 +354,8 @@ int main(int argc, char *argv[]) {
         parset.get<double>("boundary.friction.V0"));
 
     Vector v = v_initial;
-    Vector v_m(fineVertexCount);
-    ScalarVector alpha(fineVertexCount);
+    Vector v_m(leafVertexCount);
+    ScalarVector alpha(leafVertexCount);
 
     auto const timeSteps = parset.get<size_t>("timeSteps.number"),
                maximumStateFPI = parset.get<size_t>("v.fpi.maximumIterations"),
@@ -379,8 +379,8 @@ int main(int argc, char *argv[]) {
       computeExternalForces(relativeTime, ell);
 
       Matrix velocityMatrix;
-      Vector velocityRHS(fineVertexCount);
-      Vector velocityIterate(fineVertexCount);
+      Vector velocityRHS(leafVertexCount);
+      Vector velocityIterate(leafVertexCount);
 
       stateUpdater->setup(tau);
       timeSteppingScheme->setup(ell, tau, relativeTime, velocityRHS,