diff --git a/examples/old/BlockPreconditioner.hpp b/examples/old/BlockPreconditioner.hpp
deleted file mode 100644
index 423e3f7ef9743a104d902a0307a36e4544104fbd..0000000000000000000000000000000000000000
--- a/examples/old/BlockPreconditioner.hpp
+++ /dev/null
@@ -1,64 +0,0 @@
-#pragma once
-
-#include <amdis/AMDiS.hpp>
-#include <amdis/LinearAlgebra.hpp>
-
-namespace AMDiS
-{
- /// Basis preconditioner structure for block-preconditioners
- template <class Matrix, class Vector>
- class BlockPreconditioner
- : public PreconditionerInterface<Matrix, Vector>
- {
- public:
- using Self = BlockPreconditioner;
- using Super = PreconditionerInterface<Matrix, Vector>;
- using PreconBase = Super;
-
- public:
- BlockPreconditioner() = default;
-
- /// extract iranges from BlockMatrix to be used to extract sub-vectors and sub-matrices.
- virtual void init(Matrix const& A_) override
- {
- A = &A_;
-
- A_.getRanges(rows, cols);
- initialized = true;
- }
-
- virtual void solve(Vector const& b, Vector& x) const override
- {
- error_exit("Must be implemented in derived classes!\n");
- }
-
- virtual void adjoint_solve(Vector const& x, Vector& y) const override
- {
- error_exit("Must be implemented in derived classes!\n");
- }
-
- template <size_t I>
- mtl::irange const& getRowRange(const index_t<I> i = {}) const
- {
- test_exit_dbg(initialized,
- "BlockPreconditioner not initialized. Call init() before.");
- return std::get<I>(rows);
- }
-
- template <size_t I>
- mtl::irange const& getColRange(const index_t<I> i = {}) const
- {
- test_exit_dbg(initialized,
- "BlockPreconditioner not initialized. Call init() before.");
- return std::get<I>(cols);
- }
-
- protected:
- Matrix const* A = NULL;
- bool initialized = false;
-
- std::array<mtl::irange, Matrix::N()> rows;
- std::array<mtl::irange, Matrix::M()> cols;
- };
-
-} // end namespace AMDiS
diff --git a/examples/old/MTLPfcPrecon.hpp b/examples/old/MTLPfcPrecon.hpp
deleted file mode 100644
index 4f9ca3021b7c23f5359e8fee701cc5fe68cdc35b..0000000000000000000000000000000000000000
--- a/examples/old/MTLPfcPrecon.hpp
+++ /dev/null
@@ -1,125 +0,0 @@
-#pragma once
-
-#include "BlockPreconditioner.hpp"
-#include <amdis/common/Literals.hpp>
-
-namespace AMDiS
-{
- template <class Matrix, class Vector>
- class MTLPfcPrecon
- : public BlockPreconditioner<Matrix, Vector>
- {
- public:
- using Self = MTLPfcPrecon;
- using Super = BlockPreconditioner<Matrix, Vector>;
- using PreconBase = typename Super::PreconBase;
-
- using MTLMatrix = BaseMatrix_t<Matrix>;
- using MTLVector = Vector;
- using LinearSolverType = LinearSolverInterface<MTLMatrix, MTLVector>;
-
- public:
- struct Creator : CreatorInterface<PreconBase>
- {
- virtual std::shared_ptr<PreconBase> create() override
- {
- precon = std::make_shared<Self>();
- return precon;
- }
-
- std::shared_ptr<Self> precon = NULL;
- };
-
- public:
- MTLPfcPrecon()
- : Super()
- {
- std::string solverNameM = "cg",
- solverNameMpL = "cg",
- solverNameMpL2 = "cg";
- Parameters::get("precon_pfc->M->solver", solverNameM);
- Parameters::get("precon_pfc->MpL->solver", solverNameMpL);
- Parameters::get("precon_pfc->MpL2->solver", solverNameMpL2);
-
- CreatorInterfaceName<LinearSolverType>* creator;
-
- creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameM, "precon_pfc->M->solver") );
- solverM = creator->create("precon_pfc->M");
-
- creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL, "precon_pfc->MpL->solver") );
- solverMpL = creator->create("precon_pfc->MpL");
-
- creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL2, "precon_pfc->MpL2->solver") );
- solverMpL2 = creator->create("precon_pfc->MpL2");
- }
-
- virtual ~MTLPfcPrecon()
- {
- exit();
- }
-
- void setData(double* tau_, double M0_ = 1.0)
- {
- tau = tau_;
- M0 = M0_;
- }
-
- // implementation in MTLPfcPrecon.inc.hpp
- virtual void init(Matrix const& A) override;
-
- virtual void exit() override {}
-
- // implementation in MTLPfcPrecon.inc.hpp
- virtual void solve(MTLVector const& b, MTLVector& x) const override;
-
-
- template <typename VectorX, typename VectorB>
- void solveM(VectorX& x, const VectorB& b) const
- {
- SolverInfo solverInfo("precon_pfc->M");
- solverM->solve(getM(), x, b, solverInfo);
- }
-
- template <typename VectorX, typename VectorB>
- void solveMpL(VectorX& x, const VectorB& b) const
- {
- SolverInfo solverInfo("precon_pfc->MpL");
- solverMpL->solve(MpL, x, b, solverInfo);
- }
-
- template <typename VectorX, typename VectorB>
- void solveMpL2(VectorX& x, const VectorB& b) const
- {
- SolverInfo solverInfo("precon_pfc->MpL2");
- solverMpL2->solve(MpL2, x, b, solverInfo);
- }
-
- MTLMatrix const& getM() const { return matM ? *matM : (*Super::A)(_2, _2); } // < u, v >
- MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(_1, _0); } // < M0*tau*grad(u), grad(v) >
- MTLMatrix const& getL() const { return matL ? *matL : (*Super::A)(_2, _1); } // < grad(u), grad(v) >
-
- double getTau() const { return *tau; }
-
- protected:
- double* tau = NULL;
- double M0 = 1.0;
-
- MTLMatrix* matM = NULL;
- MTLMatrix* matL = NULL;
- MTLMatrix* matL0 = NULL;
-
- MTLMatrix MpL;
- MTLMatrix MpL2;
-
- mutable MTLVector y0;
- mutable MTLVector y1;
- mutable MTLVector tmp;
-
- std::shared_ptr<LinearSolverType> solverM;
- std::shared_ptr<LinearSolverType> solverMpL;
- std::shared_ptr<LinearSolverType> solverMpL2;
- };
-
-} // end namespace AMDiS
-
-#include "MTLPfcPrecon.inc.hpp"
diff --git a/examples/old/MTLPfcPrecon.inc.hpp b/examples/old/MTLPfcPrecon.inc.hpp
deleted file mode 100644
index 3b15726247fbce0401aba5988c685980fdec55e1..0000000000000000000000000000000000000000
--- a/examples/old/MTLPfcPrecon.inc.hpp
+++ /dev/null
@@ -1,64 +0,0 @@
-namespace AMDiS
-{
-
- template <class Matrix, class Vector>
- void MTLPfcPrecon<Matrix, Vector>::init(Matrix const& A)
- {
- test_exit(tau, "Preconditioner not initialized!");
- Super::init(A);
-
- using size_type = typename MTLMatrix::size_type;
- size_type n = num_rows(getM());
-
- double delta = std::sqrt(M0 * getTau());
- msg("delta = {}", delta);
-
- // helper-matrix MpL = M + delta*L
- MpL.change_dim(n, n);
- MpL = getM() + (1.0/delta) * getL0();
-
- // helper-matrix MpL = M + sqrt(delta)*L
- MpL2.change_dim(n, n);
- MpL2 = getM() + std::sqrt(delta) * getL();
-
- // temporary variables
- y0.change_dim(num_rows(getM()));
- y1.change_dim(num_rows(getM()));
- tmp.change_dim(num_rows(getM()));
- }
-
-
- template <class Matrix, class Vector>
- void MTLPfcPrecon<Matrix, Vector>::solve(Vector const& b, Vector& x) const
- {
- x.change_dim(num_rows(b));
-
- Vector x0(x[Super::getColRange(0_c)]);
- Vector x1(x[Super::getColRange(1_c)]);
- Vector x2(x[Super::getColRange(2_c)]);
-
- const Vector b0(b[Super::getRowRange(0_c)]);
- const Vector b1(b[Super::getRowRange(1_c)]);
- const Vector b2(b[Super::getRowRange(2_c)]);
-
- double delta = std::sqrt(M0 * getTau());
-
- solveM(y0, b0); // M*y0 = b0
- y1 = getL0() * y0; // y1 = K*y0
- tmp = b1 - y1; // tmp := b1 - tau*y1
-
- solveMpL(y1, tmp); // (M + delta*K) * y1 = tmp
- x0 = y0 + (1.0/delta)*y1; // x0 = y0 + (1/delta)*y1
-
- tmp = getM() * y1; // tmp := M*y1
- solveMpL2(y1, tmp); // (M+eps*sqrt(delta)K) * y1 = tmp
- tmp = getM() * y1; // tmp := M*y1
- solveMpL2(x1, tmp); // (M+eps*sqrt(delta)K) * x1 = tmp
-
- x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
-
- tmp = b2 - getL() * x1; // tmp := b2 - K*x1
- solveM(x2, tmp);
- }
-
-} // end namespace AMDiS
diff --git a/examples/old/PfcPrecon.hpp b/examples/old/PfcPrecon.hpp
deleted file mode 100644
index 1b8b9fe17f5a279a60a12915ebade4006d5a5ec9..0000000000000000000000000000000000000000
--- a/examples/old/PfcPrecon.hpp
+++ /dev/null
@@ -1,110 +0,0 @@
-#pragma once
-
-#include <tuple>
-#include <cmath>
-
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-#include <dune/istl/multitypeblockmatrix.hh>
-#include <dune/istl/multitypeblockvector.hh>
-#include <dune/istl/preconditioner.hh>
-
-#include <amdis/common/Mpl.hpp>
-
-namespace AMDiS {
-
-template <class Matrix, class Vector>
-class PfcPrecon
- : public Dune::Preconditioner<Vector, Vector>
-{
- using SubMatrix = std::decay_t<decltype( std::declval<Matrix>()[_0][_0] )>;
- using SubVector = std::decay_t<decltype( std::declval<Vector>()[_0] )>;
-
-public:
- PfcPrecon(Matrix const& matrix, double* tauPtr, double M0)
- : matrix(matrix)
- , tauPtr(tauPtr)
- , M0(M0)
- , matL0(matrix[_1][_0])
- , matL (matrix[_2][_1])
- , matM (matrix[_2][_2])
- {}
-
- virtual void pre(Vector& x, Vector& b) override
- {
- double delta = std::sqrt(M0 * (*tauPtr));
-
- matMpL = matM;
- matMpL.axpy(1.0/delta, matL0); // => MpL = M + 1/delta * L0
-
- matMpL2 = matM;
- matMpL2.axpy(std::sqrt(delta), matL);
-
- y0.resize(matM.N());
- y1.resize(matM.N());
- tmp.resize(matM.N());
- }
-
- virtual void apply(Vector& x, const Vector& b) override
- {
- double delta = std::sqrt(M0 * (*tauPtr));
-
- solve(matM, y0, b[_0]); // M*y0 = b0
- matL0.mv(y0, y1); // y1 = K*y0
- tmp = b[_1];
- tmp-= y1; // tmp := b1 - tau*y1
-
- solve(matMpL, y1, tmp); // (M + delta*K) * y1 = tmp
- x[_0] = y0;
- x[_0].axpy(1.0/delta, y1); // x0 = y0 + (1/delta)*y1
-
- matM.mv(y1, tmp); // tmp := M*y1
- solve(matMpL2, y1, tmp); // (M+eps*sqrt(delta)K) * y1 = tmp
- matM.mv(y1, tmp); // tmp := M*y1
- solve(matMpL2, x[_1], tmp); // (M+eps*sqrt(delta)K) * x1 = tmp
-
- x[_0].axpy(-1.0/delta, x[_1]); // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
-
- matL.mv(x[_1], y1);
- tmp = b[_2];
- tmp-= y1; // tmp := b2 - K*x1
- solve(matM, x[_2], tmp);
- }
-
- virtual void post(Vector& x) override
- {
-
- }
-
-
-private:
-
- template <class Mat, class Vec>
- void solve(Mat const& A, Vec& x, Vec b)
- {
- Dune::MatrixAdapter<Mat, Vec, Vec> op(A);
- Dune::SeqJac<Mat, Vec, Vec> precon(A, 1, 1.0);
- Dune::CGSolver<Vec> solver(op, precon, 1.e-3, 20, 0);
-
- Dune::InverseOperatorResult statistics;
- solver.apply(x, b, statistics);
- }
-
-
-private:
-
- Matrix const& matrix;
- double* tauPtr;
- double M0;
-
- SubMatrix const& matL0;
- SubMatrix const& matL;
- SubMatrix const& matM;
-
- SubMatrix matMpL;
- SubMatrix matMpL2;
-
- SubVector y0, y1, tmp;
-};
-
-} // end namespace AMDiS
diff --git a/examples/old/pfc.cc b/examples/old/pfc.cc
deleted file mode 100644
index 33a32e1f37a2daee4c10dc8f9f603cc2300351a8..0000000000000000000000000000000000000000
--- a/examples/old/pfc.cc
+++ /dev/null
@@ -1,109 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
-#include <iostream>
-#include <ctime>
-#include <cmath>
-
-#include <amdis/AMDiS.hpp>
-#include <amdis/AdaptInstationary.hpp>
-#include <amdis/ProblemStat.hpp>
-#include <amdis/ProblemInstat.hpp>
-#include <amdis/Terms.hpp>
-
-#include "MTLPfcPrecon.hpp"
-
-using namespace AMDiS;
-
-// using Mesh = Dune::YaspGrid<AMDIS_DIM>;
-using Mesh = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
-
-using PfcParam = LagrangeBasis<Mesh::LeafGridView, 2, 2, 2>;
-
-
-// 3 component with polynomial degree 1
-using PfcProblem = ProblemStat<PfcParam>;
-using PfcProblemInstat = ProblemInstat<PfcParam>;
-
-using MatrixType = PfcProblem::SystemMatrixType::MultiMatrix;
-using VectorType = PfcProblem::SystemVectorType::BaseVector;
-
-using Precon = MTLPfcPrecon<MatrixType,VectorType>;
-
-int main(int argc, char** argv)
-{
- AMDiS::init(argc, argv);
-
- auto creator = new Precon::Creator;
- CreatorMap<typename Precon::PreconBase>::addCreator("pfc", creator);
-
- PfcProblem prob("pfc");
- prob.initialize(INIT_ALL);
-
- PfcProblemInstat probInstat("pfc", prob);
- probInstat.initialize(INIT_UH_OLD);
-
- AdaptInfo adaptInfo("adapt", prob.getNumComponents());
-
- double r = -0.4;
- double psi_mean = -0.3;
- double M0 = 1.0;
-
- auto& Mu = prob.getSolution<0>();
- auto& Psi = prob.getSolution<1>();
- auto& Nu = prob.getSolution<2>();
-
- using Op = PfcProblem::ElementOperator;
-
- // < [-(1+r) - 3*psi^2]*u, v > + < 2*grad(u), grad(v) >
- auto opLhs01_ = Op::zot( -(1.0 + r) );
- opLhs01_.addZOT( valueOfFunc(Psi, [](auto psi) { return -3.0 * Math::pow<2>(psi); }, 2) );
- opLhs01_.addSOT( 2.0 );
- auto opLhs01 = std::make_pair(opLhs01_, true);
-
- // < -2*psi^3, v >
- auto opRhs0 = Op::zot(valueOfFunc(Psi, [](auto psi) { return -2.0 * Math::pow<3>(psi); }, 3));
-
- // < psi, v >
- auto opRhs1 = Op::zot(valueOf(Psi));
-
- // < u, v >
- auto opM = std::make_pair(Op::zot(1.0), false);
-
- // < grad(u), grad(v) >
- auto opL = std::make_pair(Op::sot(1.0), false);
-
- // < M0*tau * grad(u), grad(v) >
- double* tauPtr = adaptInfo.getTimestepPtr();
- auto opL0 = std::make_pair(Op::sot([M0, tauPtr](auto /*x*/) { return M0 * (*tauPtr); }), false);
-
-
- // === Define the linear system ===
-
- prob.addMatrixOperator({{ {{0,0}, opM}, {{0,1}, opLhs01}, {{0,2}, opL},
- {{1,0}, opL0}, {{1,1}, opM},
- {{2,1}, opL}, {{2,2}, opM}}});
-
- prob.addVectorOperator({{ {0, opRhs0},
- {1, opRhs1} }});
-
- // === set the initial solution ===
- Mu = 0.0;
- Nu = 0.0;
- std::srand( 12345 );
- Psi.interpol([psi_mean](auto const& x) { return 0.5*(std::rand()/double(RAND_MAX) - 0.5) + psi_mean; });
-
- // === configure preconditioner, if necessary ===
- if (creator->precon) {
- creator->precon->setData(tauPtr, M0);
- }
-
- AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
- adapt.adapt();
-
- msg("|mu| = {}, |nu| = {}, |psi| = {}",
- two_norm(Mu.getVector()), two_norm(Nu.getVector()), two_norm(Psi.getVector()));
-
- AMDiS::finalize();
- return 0;
-}