diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc index 50326584720c881b0116bc713e44bd0eaf4d62e1..52cd4665f543f40f4382da76ee63638a995b2016 100644 --- a/src/sand-wedge.cc +++ b/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,