diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 248f3ae74a35d1c718a123c74f24338055919914..2ab6ddc9eee2ec47aa310a1893859a0470371deb 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -6,8 +6,11 @@ #include <dune/common/exceptions.hh> +#include <dune/common/bitsetvector.hh> #include <dune/common/fmatrix.hh> +#include <dune/common/fvector.hh> #include <dune/istl/bcrsmatrix.hh> +#include <dune/istl/bvector.hh> #include <dune/grid/yaspgrid.hh> @@ -15,13 +18,23 @@ #include <dune/fufem/assemblers/operatorassembler.hh> #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh> +//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh> + +#include <dune/solvers/iterationsteps/blockgsstep.hh> +#include <dune/solvers/solvers/loopsolver.hh> +#include <dune/solvers/norms/energynorm.hh> +#include <dune/solvers/common/numproc.hh> // Solver::FULL + #include <exception> int const dim = 2; int main() { try { - typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, dim, dim>> OperatorType; + typedef Dune::FieldVector<double, dim> SmallVector; + typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix; + typedef Dune::BCRSMatrix<SmallMatrix> OperatorType; + typedef Dune::BlockVector<SmallVector> VectorType; // FIXME: Random values double const E = 1e8; @@ -52,6 +65,32 @@ int main() { OperatorType stiffnessMatrix; globalAssembler.assemble(localStiffness, stiffnessMatrix); + VectorType f(grid.size(grid.maxLevel(), dim)); + // Fill right-hand side with semi-random numbers + for (size_t i = 0; i < f.size(); ++i) + for (size_t j = 0; j < dim; ++j) + f[i][j] = i + j; + + VectorType u(grid.size(grid.maxLevel(), dim)); + // Fill initial guess with other semi-random numbers + for (size_t i = 0; i < f.size(); ++i) + for (size_t j = 0; j < dim; ++j) + f[i][j] = i * j; + + // TODO: Why does blockGSStep even provide a default constructor? + BlockGSStep<OperatorType, VectorType> blockGSStep(stiffnessMatrix, u, f); + EnergyNorm<OperatorType, VectorType> energyNorm(blockGSStep); + + // TODO: Why does blockGSStep not warn about this? + // Why is this not documented? + Dune::BitSetVector<VectorType::block_type::dimension> ignoreNodes(false); + blockGSStep.ignoreNodes_ = &ignoreNodes; + + LoopSolver<VectorType> solver(&blockGSStep, 1000, // maxIterations + 1e-5, // tolerance + &energyNorm, Solver::FULL); + + solver.solve(); return 0; } catch (Dune::Exception& e) {