From 36a99a362007f56784a591315440610815e0eed1 Mon Sep 17 00:00:00 2001 From: Elias Pipping <elias.pipping@fu-berlin.de> Date: Tue, 30 Oct 2012 16:07:56 +0100 Subject: [PATCH] Cleanup Move code around, add a comment --- src/one-body-sample.org | 78 ++++++++++++++++++++--------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/src/one-body-sample.org b/src/one-body-sample.org index 0d39401e..b9847335 100644 --- a/src/one-body-sample.org +++ b/src/one-body-sample.org @@ -256,18 +256,13 @@ #+name: mainBody #+begin_src c++ try { - typedef SharedPointerMap - <std::string, Dune::VirtualFunction<double, double> > FunctionMap; - FunctionMap functions; - initPython(functions); + Dune::Timer timer; Dune::ParameterTree parset; Dune::ParameterTreeParser::readINITree (srcdir "/one-body-sample.parset", parset); Dune::ParameterTreeParser::readOptions(argc, argv, parset); - Dune::Timer timer; - typedef Dune::FieldVector<double, dim> SmallVector; typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix; typedef Dune::BCRSMatrix<SmallMatrix> MatrixType; @@ -276,13 +271,11 @@ auto const E = parset.get<double>("body.E"); auto const nu = parset.get<double>("body.nu"); - auto const solver_tolerance = parset.get<double>("solver.tolerance"); - - auto const refinements = parset.get<size_t>("grid.refinements"); - - auto const verbose = parset.get<bool>("verbose"); - Solver::VerbosityMode const verbosity - = verbose ? Solver::FULL : Solver::QUIET; + auto const mu = parset.get<double>("boundary.friction.mu"); + auto const a = parset.get<double>("boundary.friction.a"); + auto const b = parset.get<double>("boundary.friction.b"); + auto const V0 = parset.get<double>("boundary.friction.V0"); + auto const L = parset.get<double>("boundary.friction.L"); // {{{ Set up grid typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::nonconforming> @@ -301,6 +294,7 @@ auto grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid (lowerLeft, upperRight, elements); + auto const refinements = parset.get<size_t>("grid.refinements"); grid->globalRefine(refinements); size_t const finestSize = grid->size(grid->maxLevel(), dim); @@ -314,6 +308,22 @@ P0Basis const p0Basis(leafView); P1Basis const p1Basis(leafView); + // 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); + <<setupBoundary>>; + + // Set up functions for time-dependent boundary conditions + typedef SharedPointerMap + <std::string, Dune::VirtualFunction<double, double> > FunctionMap; + FunctionMap functions; + initPython(functions); + auto const &dirichletFunction = functions.get("dirichletCondition"); + auto const &neumannFunction = functions.get("neumannCondition"); + + // Set up normal stress, mass matrix, and gravity functional double normalStress; MatrixType massMatrix; VectorType gravityFunctional; @@ -325,18 +335,16 @@ <<computeNormalStress>>; <<computeGravitationalBodyForce>>; } + SingletonVectorType surfaceNormalStress(finestSize); + surfaceNormalStress = 0.0; + for (size_t i = 0; i < frictionalNodes.size(); ++i) + if (frictionalNodes[i][0]) + surfaceNormalStress[i] = normalStress; MatrixType stiffnessMatrix; <<assembleStiffnessMatrix>>; 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); - <<setupBoundary>>; - auto const nodalIntegrals = assemble_frictional<GridType, GridView, SmallVector, P1Basis> (leafView, p1Basis, frictionalNodes); @@ -354,35 +362,22 @@ SingletonVectorType alpha(alpha_initial); SingletonVectorType vonMisesStress; - - SingletonVectorType surfaceNormalStress(finestSize); - surfaceNormalStress = 0.0; - for (size_t i = 0; i < frictionalNodes.size(); ++i) - if (frictionalNodes[i][0]) - surfaceNormalStress[i] = normalStress; // }}} + // Set up TNNMG solver + auto const solver_tolerance = parset.get<double>("solver.tolerance"); typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType; typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType; - - // Set up TNNMG solver MySolver<dim, MatrixType, VectorType, GridType, MyBlockProblemType> mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance, *grid, ignoreNodes); - <<createWriters>>; - - auto const L = parset.get<double>("boundary.friction.L"); - auto const a = parset.get<double>("boundary.friction.a"); - auto const b = parset.get<double>("boundary.friction.b"); - auto const V0 = parset.get<double>("boundary.friction.V0"); - auto const mu = parset.get<double>("boundary.friction.mu"); - auto const timesteps = parset.get<size_t>("timeSteps"); - auto const tau = parset.get<double>("endOfTime") / timesteps; + Solver::VerbosityMode const verbosity + = parset.get<bool>("verbose") ? Solver::FULL : Solver::QUIET; - auto const &dirichletFunction = functions.get("dirichletCondition"); - auto const &neumannFunction = functions.get("neumannCondition"); + <<createWriters>>; + // Set up time steppers that auto timeSteppingScheme = initTimeStepper (parset.get<Config::scheme>("timeSteppingScheme"), dirichletFunction, ignoreNodes, massMatrix, stiffnessMatrix, u_initial, ud_initial, @@ -391,6 +386,9 @@ (parset.get<Config::state_model>("boundary.friction.state_model"), alpha_initial, frictionalNodes, L); + auto const timesteps = parset.get<size_t>("timeSteps"); + auto const tau = parset.get<double>("endOfTime") / timesteps; + auto const state_fpi_max = parset.get<size_t>("solver.tnnmg.fixed_point_iterations"); for (size_t run = 1; run <= timesteps; ++run) { @@ -412,6 +410,8 @@ timeSteppingScheme->setup(ell, tau, time, problem_rhs, problem_iterate, problem_A); + // Since the velocity explodes in the quasistatic case, use the + // displacement as a convergence criterion VectorType u_saved; for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) { <<setupAndSolveProblem>>; -- GitLab