diff --git a/dune/solvers/test/cgsteptest.cc b/dune/solvers/test/cgsteptest.cc
index 6386cfe7fe19fe4d265a36e8c15fe844e3bf9f00..792fe56ceffaeaa446b30b2a3fb1caed2c357e68 100644
--- a/dune/solvers/test/cgsteptest.cc
+++ b/dune/solvers/test/cgsteptest.cc
@@ -61,52 +61,24 @@ struct CGTestSuite
bool passed = true;
- typedef typename Dune::FieldMatrix<double, 1, 1> MatrixBlock;
- typedef typename Dune::FieldVector<double, 1> VectorBlock;
- typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
- typedef typename Dune::BlockVector<VectorBlock> Vector;
- typedef typename Dune::BitSetVector<1> BitVector;
-
- typedef typename GridType::LevelGridView GridView;
-
- const GridView gridView = grid.levelGridView(grid.maxLevel());
-
- Matrix A;
- constructPQ1Pattern(gridView, A);
- A=0.0;
- assemblePQ1Stiffness(gridView, A);
- EnergyNorm<Matrix, Vector> energyNorm(A);
-
- Vector u(A.N()), u_sol(A.N());
- u = 0;
- for (size_t i=0; i<u.size(); ++i)
- u_sol[i] = (1.0*rand())/RAND_MAX;
-
- Vector rhs(A.N());
- A.mv(u_sol, rhs);
-
- BitVector ignore(A.N());
- ignore.unsetAll();
- markBoundaryDOFs(gridView, ignore);
-
- // set Dirichlet values
- for (size_t i=0; i<u.size(); ++i)
- if (ignore[i][0])
- u[i] = u_sol[i];
+ using Problem =
+ SymmetricSampleProblem<1, typename GridType::LevelGridView>;
+ Problem p(grid.levelGridView(grid.maxLevel()));
const auto verbosity = Solver::REDUCED;
const bool relativeErrors = false;
- Analyser<Vector> analyser(u_sol, energyNorm, solveTol);
+ Analyser<typename Problem::Vector> analyser(p.u_ex, p.energyNorm,
+ solveTol);
{
std::string const testCase = "CGSolver ";
- Vector u_copy = u;
- Vector rhs_copy = rhs;
+ typename Problem::Vector u_copy = p.u;
+ typename Problem::Vector rhs_copy = p.rhs;
- CGSolver<Matrix,Vector> solver(&A, &u_copy, &rhs_copy, nullptr,
- maxIterations, stepTol, &energyNorm,
- verbosity, relativeErrors);
+ CGSolver<typename Problem::Matrix, typename Problem::Vector> solver(
+ &p.A, &u_copy, &rhs_copy, nullptr, maxIterations, stepTol,
+ &p.energyNorm, verbosity, relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
@@ -116,16 +88,19 @@ struct CGTestSuite
{
std::string const testCase = "CGSolver, preconditioned ";
- Vector u_copy = u;
- Vector rhs_copy = rhs;
+ typename Problem::Vector u_copy = p.u;
+ typename Problem::Vector rhs_copy = p.rhs;
- using BGS = BlockGSStep<Matrix,Vector,BitVector>;
+ using BGS =
+ BlockGSStep<typename Problem::Matrix, typename Problem::Vector,
+ typename Problem::BitVector>;
BGS blockgs(BGS::Direction::SYMMETRIC);
- blockgs.ignoreNodes_ = new BitVector(u.size(), false);
+ blockgs.ignoreNodes_ =
+ new typename Problem::BitVector(p.u.size(), false);
- CGSolver<Matrix,Vector> solver(&A, &u_copy, &rhs_copy, &blockgs,
- maxIterations, stepTol, &energyNorm,
- verbosity, relativeErrors);
+ CGSolver<typename Problem::Matrix, typename Problem::Vector> solver(
+ &p.A, &u_copy, &rhs_copy, &blockgs, maxIterations, stepTol,
+ &p.energyNorm, verbosity, relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
@@ -135,12 +110,15 @@ struct CGTestSuite
{
std::string const testCase = "Dune::Solvers::CGStep ";
- Vector u_copy = u;
- Vector rhs_copy = rhs;
+ typename Problem::Vector u_copy = p.u;
+ typename Problem::Vector rhs_copy = p.rhs;
- Dune::Solvers::CGStep<Matrix,Vector> cgStep(A, u_copy, rhs_copy);
- ::LoopSolver<Vector> solver(&cgStep, maxIterations, stepTol, &energyNorm,
- verbosity, relativeErrors);
+ Dune::Solvers::CGStep<typename Problem::Matrix,
+ typename Problem::Vector> cgStep(p.A, u_copy,
+ rhs_copy);
+ ::LoopSolver<typename Problem::Vector> solver(
+ &cgStep, maxIterations, stepTol, &p.energyNorm, verbosity,
+ relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
@@ -150,16 +128,23 @@ struct CGTestSuite
{
std::string const testCase = "Dune::Solvers::CGStep, preconditioned";
- Vector u_copy = u;
- Vector rhs_copy = rhs;
+ typename Problem::Vector u_copy = p.u;
+ typename Problem::Vector rhs_copy = p.rhs;
- using BGS = BlockGSStep<Matrix,Vector,BitVector>;
+ using BGS =
+ BlockGSStep<typename Problem::Matrix, typename Problem::Vector,
+ typename Problem::BitVector>;
BGS blockgs(BGS::Direction::SYMMETRIC);
- blockgs.ignoreNodes_ = new BitVector(u.size(), false);
-
- Dune::Solvers::CGStep<Matrix,Vector> cgStep(A, u_copy, rhs_copy, blockgs);
- ::LoopSolver<Vector> solver(&cgStep, maxIterations, stepTol, &energyNorm,
- verbosity, relativeErrors);
+ blockgs.ignoreNodes_ =
+ new typename Problem::BitVector(p.u.size(), false);
+
+ Dune::Solvers::CGStep<typename Problem::Matrix,
+ typename Problem::Vector> cgStep(p.A, u_copy,
+ rhs_copy,
+ blockgs);
+ ::LoopSolver<typename Problem::Vector> solver(
+ &cgStep, maxIterations, stepTol, &p.energyNorm, verbosity,
+ relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
diff --git a/dune/solvers/test/common.hh b/dune/solvers/test/common.hh
index 7c7b6badeb13dc50046288cec3a4b7dbf1d9da80..1c585060b19fb0117c8a65b643d92cceb8b85fbe 100644
--- a/dune/solvers/test/common.hh
+++ b/dune/solvers/test/common.hh
@@ -348,6 +348,61 @@ void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
}
}
+template <size_t blocksize, class GridView>
+class SymmetricSampleProblem {
+public:
+ typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+ typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
+ typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
+ typedef typename Dune::BlockVector<VectorBlock> Vector;
+ typedef typename Dune::BitSetVector<blocksize> BitVector;
+ typedef EnergyNorm<Matrix, Vector> Norm;
+
+ SymmetricSampleProblem(GridView const& gridView) {
+ constructPQ1Pattern(gridView, A);
+ A = 0.0;
+ assemblePQ1Stiffness(gridView, A);
+
+ energyNorm.setMatrix(&A);
+
+ ignore.resize(A.N());
+ ignore.unsetAll();
+ markBoundaryDOFs(gridView, ignore);
+
+ u.resize(A.N());
+ u_ex.resize(A.N());
+ rhs.resize(A.N());
+ randomize();
+ }
+
+ void randomize() {
+ // Set up a random 'solution'
+ u = 0;
+ for (size_t i = 0; i < u.size(); ++i)
+ for (size_t j = 0; j < blocksize; ++j) {
+ // in case of no dirichlet values you can make the
+ // matrix pd hereby
+ // A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
+ u_ex[i][j] = (1.0 * rand()) / RAND_MAX;
+ }
+
+ // Construct right hand side corresponding to the 'solution'
+ A.mv(u_ex, rhs);
+
+ // set dirichlet values
+ for (size_t i = 0; i < u.size(); ++i)
+ for (size_t j = 0; j < blocksize; ++j)
+ if (ignore[i][j])
+ u[i][j] = u_ex[i][j];
+ }
+
+ Matrix A;
+ Vector u;
+ Vector u_ex;
+ Vector rhs;
+ BitVector ignore;
+ Norm energyNorm;
+};
template<class GridType, class TestSuite>
bool checkWithGrid(TestSuite& suite, int uniformRefinements)
diff --git a/dune/solvers/test/multigridtest.cc b/dune/solvers/test/multigridtest.cc
index 86c8e71fb096af17d597f857f1ad36a2c5d0ce59..4529b250b7284941dc1df10c894918cb045b2fc9 100644
--- a/dune/solvers/test/multigridtest.cc
+++ b/dune/solvers/test/multigridtest.cc
@@ -43,45 +43,13 @@ struct MultigridTestSuite
bool passed = true;
- typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
- typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
- typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
- typedef typename Dune::BlockVector<VectorBlock> Vector;
- typedef typename Dune::BitSetVector<blocksize> BitVector;
+ using Problem =
+ SymmetricSampleProblem<1, typename GridType::LevelGridView>;
+ Problem p(grid.levelGridView(grid.maxLevel()));
- typedef typename GridType::LevelGridView GridView;
-
-
- const GridView gridView = grid.levelGridView(grid.maxLevel());
-
- Matrix A;
- constructPQ1Pattern(gridView, A);
- A=0.0;
- assemblePQ1Stiffness(gridView, A);
-
- Vector u(A.N()), u_ex(A.N());
- u = 0;
- for (size_t i=0; i<u.size(); ++i)
- for (size_t j=0; j<blocksize; ++j)
- {
- u_ex[i][j] = (1.0*rand())/RAND_MAX;
-// A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]); // in case of no dirichlet values you can make the matrix pd hereby
- }
-
- Vector rhs(A.N());
- A.mv(u_ex, rhs);
-
- BitVector ignore(A.N());
- ignore.unsetAll();
- markBoundaryDOFs(gridView, ignore);
-
- // set dirichlet values
- for (size_t i=0; i<u.size(); ++i)
- for (size_t j=0; j<blocksize; ++j)
- if (ignore[i][j])
- u[i][j] = u_ex[i][j];
-
- typedef EnergyNorm<Matrix, Vector> Norm;
+ typedef typename Problem::Vector Vector;
+ typedef typename Problem::Matrix Matrix;
+ typedef typename Problem::BitVector BitVector;
typedef ::LoopSolver<Vector> Solver;
typedef MultigridStep<Matrix, Vector, BitVector> MGStep;
typedef MultigridTransfer<Vector, BitVector, Matrix> Transfer;
@@ -98,71 +66,51 @@ struct MultigridTestSuite
transfer[i] = t;
}
- // we want to control errors with respect to the energy norm induced by the matrix
- Norm energynorm(A);
-
// set up smoothers and basesolver
Smoother smoother, basesolver_step;
- Norm basenorm(basesolver_step);
+ typename Problem::Norm basenorm(basesolver_step);
Solver basesolver(&basesolver_step, 1, 1e-10, &basenorm, Solver::QUIET);
// create iteration step and set all kinds of data
MGStep mgStep;
mgStep.setTransferOperators(transfer);
- mgStep.setProblem(A,u,rhs);
+ mgStep.setProblem(p.A,p.u,p.rhs);
mgStep.setSmoother(&smoother);
mgStep.setMGType(1,3,3);
- mgStep.ignoreNodes_ = &ignore;
+ mgStep.ignoreNodes_ = &p.ignore;
mgStep.basesolver_ = &basesolver;
// create loop solver
- Solver solver(&mgStep, maxIterations, tol, &energynorm, Solver::FULL);
+ Solver solver(&mgStep, maxIterations, tol, &p.energyNorm, Solver::FULL);
// solve problem
solver.preprocess();
solver.solve();
- u = mgStep.getSol();
+ p.u = mgStep.getSol();
- if (energynorm.diff(u,u_ex)>tol*10)
+ if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
{
- std::cout << "The MultigridStep did not produce a satisfactory result. ||u-u_ex||=" << energynorm.diff(u,u_ex) << std::endl;
- std::cout << "||u_ex||=" << energynorm(u_ex) << std::endl;
+ std::cout << "The MultigridStep did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
+ std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
passed = false;
}
-
- // now solve with the same operator but different rhs (and dirichlet values)
-
- // create new exact solution
- for (size_t i=0; i<u.size(); ++i)
- for (size_t j=0; j<blocksize; ++j)
- u_ex[i][j] = (1.0*rand())/RAND_MAX;
-
- // ... and its corresponding rhs
- A.mv(u_ex, rhs);
-
- // set initial iterate (remember that the iterationStep member points to u)
- u = 0.0;
- // set dirichlet values
- for (size_t i=0; i<u.size(); ++i)
- for (size_t j=0; j<blocksize; ++j)
- if (ignore[i][j])
- u[i][j] = u_ex[i][j];
+ p.randomize();
// reset rhs in mgStep without calling preprocess subsequently
- mgStep.setRhs(rhs);
+ mgStep.setRhs(p.rhs);
solver.solve();
- u = mgStep.getSol();
+ p.u = mgStep.getSol();
- if (energynorm.diff(u,u_ex)>tol*10)
+ if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
{
- std::cout << "The MultigridStep did not produce a satisfactory result after calling setRhs(). ||u-u_ex||=" << energynorm.diff(u,u_ex) << std::endl;
- std::cout << "||u_ex||=" << energynorm(u_ex) << std::endl;
+ std::cout << "The MultigridStep did not produce a satisfactory result after calling setRhs(). ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
+ std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
passed = false;
}
diff --git a/dune/solvers/test/umfpacksolvertest.cc b/dune/solvers/test/umfpacksolvertest.cc
index 6b4722308b2ee7993f883ed21f1de9daa7d1dc99..e70a37bbcf7dd52e1b12a739c533e1174f941ef5 100644
--- a/dune/solvers/test/umfpacksolvertest.cc
+++ b/dune/solvers/test/umfpacksolvertest.cc
@@ -34,52 +34,14 @@ struct UMFPackSolverTestSuite
bool passed = true;
- typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
- typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
- typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
- typedef typename Dune::BlockVector<VectorBlock> Vector;
- typedef typename Dune::BitSetVector<blocksize> BitVector;
- typedef EnergyNorm<Matrix, Vector> Norm;
-
- const auto gridView = grid.leafGridView();
-
- Matrix A;
- constructPQ1Pattern(gridView, A);
- A=0.0;
- assemblePQ1Stiffness(gridView, A);
-
- // Set up a random 'solution'
- Vector u(A.N()), u_ex(A.N());
- u = 0;
- for (size_t i=0; i<u.size(); ++i)
- for (size_t j=0; j<blocksize; ++j)
- {
- u_ex[i][j] = (1.0*rand())/RAND_MAX;
-// A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]); // in case of no dirichlet values you can make the matrix pd hereby
- }
-
- // Construct right hand side corresponding to the 'solution'
- Vector rhs(A.N());
- A.mv(u_ex, rhs);
-
- BitVector ignore(A.N());
- ignore.unsetAll();
- markBoundaryDOFs(gridView, ignore);
-
- // set dirichlet values
- for (size_t i=0; i<u.size(); ++i)
- for (size_t j=0; j<blocksize; ++j)
- if (ignore[i][j])
- u[i][j] = u_ex[i][j];
-
- // we want to control errors with respect to the energy norm induced by the matrix
- Norm energynorm(A);
-
- typedef Solvers::UMFPackSolver<Matrix,Vector> Solver;
-
- // create loop solver
- Solver solver(A,u,rhs);
- solver.ignoreNodes_ = &ignore;
+ using Problem =
+ SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>;
+ Problem p(grid.leafGridView());
+
+ typedef Solvers::UMFPackSolver<typename Problem::Matrix,
+ typename Problem::Vector> Solver;
+ Solver solver(p.A,p.u,p.rhs);
+ solver.ignoreNodes_ = &p.ignore;
// solve problem
solver.preprocess();
@@ -87,10 +49,10 @@ struct UMFPackSolverTestSuite
//std::cout << u_ex << std::endl;
- if (energynorm.diff(u,u_ex)>tol*10)
+ if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
{
- std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << energynorm.diff(u,u_ex) << std::endl;
- std::cout << "||u_ex||=" << energynorm(u_ex) << std::endl;
+ std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
+ std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
passed = false;
}