diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index e6a3eee246f5960d57f1a46e7e1e290a9f36861d..0dfd443e5f726fa2cb75a79b23fde91b1749101e 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -19,7 +19,7 @@ #include <dune/fufem/assemblers/operatorassembler.hh> #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh> -//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh> +#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh> #include <dune/solvers/iterationsteps/blockgsstep.hh> #include <dune/solvers/solvers/loopsolver.hh> @@ -28,9 +28,15 @@ #include <dune/grid/common/mcmgmapper.hh> +#include <dune/tnnmg/problem-classes/convexproblem.hh> +#include <dune/tnnmg/problem-classes/nonlinearity.hh> +#include <dune/tnnmg/nonlinearities/zerononlinearity.hh> + #include <exception> #include <iostream> +#include "myblockproblem.hh" + int const dim = 2; int main() { @@ -131,6 +137,36 @@ int main() { std::cout << "Number of bounding nodes: " << bounding_nodes << std::endl; std::cout << "Number of extremal nodes: " << extremal_nodes << std::endl; } + + VectorType copy_of_u = u; + + { // experiment with convex problems and the like + + typedef ZeroNonlinearity<SmallVector, SmallMatrix> NonlinearityType; + NonlinearityType phi; + + typedef ConvexProblem<NonlinearityType, OperatorType> ConvexProblemType; + + VectorType unused_vector(grid.size(grid.maxLevel(), dim), 0); + ConvexProblemType myConvexProblem(1, stiffnessMatrix, 0, unused_vector, + phi, f, copy_of_u); + + typedef MyBlockProblem<ConvexProblemType> MyBlockProblemType; + MyBlockProblemType myBlockProblem(myConvexProblem); + + GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep; + nonlinearGSStep.ignoreNodes_ = &ignoreNodes; + nonlinearGSStep.setProblem(copy_of_u, myBlockProblem); + + // FIXME: Does this make any sense? + EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix); + LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations, + solver_tolerance, &energyNorm, + // Solver::QUIET); + Solver::FULL); + solver.solve(); + } + // TODO: Why does blockGSStep even provide a default constructor? BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f); blockGSStep.ignoreNodes_ = &ignoreNodes; @@ -138,10 +174,17 @@ int main() { // FIXME: Does this make any sense? EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix); LoopSolver<VectorType> solver(&blockGSStep, solver_maxIterations, - solver_tolerance, &energyNorm, Solver::FULL); + solver_tolerance, &energyNorm, + // Solver::QUIET); + Solver::FULL); solver.solve(); - return 0; + + VectorType diff = u; + diff -= copy_of_u; + std::cout << diff << std::endl; + std::cout << std::endl; + std::cout << diff.infinity_norm() << std::endl; } catch (Dune::Exception& e) { Dune::derr << "Dune reported error: " << e << std::endl;