diff --git a/src/myblockproblem.hh b/src/myblockproblem.hh index 284a71f59f5a120fc6271c793e75c7a5e9486857..27d006b4935bdffc419e9806853458add4ed0c07 100644 --- a/src/myblockproblem.hh +++ b/src/myblockproblem.hh @@ -4,7 +4,7 @@ #include <dune/common/bitsetvector.hh> #include <dune/tnnmg/problem-classes/bisection.hh> -#include <dune/tnnmg/problem-classes/convexproblem.hh> +//#include <dune/tnnmg/problem-classes/convexproblem.hh> #include <dune/tnnmg/problem-classes/nonlinearity.hh> #include <dune/tnnmg/problem-classes/onedconvexfunction.hh> @@ -14,21 +14,21 @@ /** \brief Base class for problems where each block can be solved with a scalar * Gauss-Seidel method */ -template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem { +template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem { public: - typedef ConvexProblemTypeTEMPLATE ConvexProblemType; - typedef typename ConvexProblemType::NonlinearityType NonlinearityType; - typedef typename ConvexProblemType::VectorType VectorType; - typedef typename ConvexProblemType::MatrixType MatrixType; - typedef typename ConvexProblemType::LocalVectorType LocalVectorType; - typedef typename ConvexProblemType::LocalMatrixType LocalMatrixType; + typedef MyConvexProblemTypeTEMPLATE MyConvexProblemType; + typedef typename MyConvexProblemType::NonlinearityType NonlinearityType; + typedef typename MyConvexProblemType::VectorType VectorType; + typedef typename MyConvexProblemType::MatrixType MatrixType; + typedef typename MyConvexProblemType::LocalVectorType LocalVectorType; + typedef typename MyConvexProblemType::LocalMatrixType LocalMatrixType; - static const int block_size = ConvexProblemType::block_size; + static const int block_size = MyConvexProblemType::block_size; /** \brief Solves one local system using a scalar Gauss-Seidel method */ class IterateObject; - MyBlockProblem(ConvexProblemType& problem) : problem(problem) { + MyBlockProblem(MyConvexProblemType& problem) : problem(problem) { bisection = Bisection(0.0, 1.0, 1e-15, true, 1e-14); }; @@ -37,15 +37,15 @@ template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem { private: // problem data - ConvexProblemType& problem; + MyConvexProblemType& problem; // commonly used minimization stuff Bisection bisection; }; /** \brief Solves one local system using a scalar Gauss-Seidel method */ -template <class ConvexProblemTypeTEMPLATE> -class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject { +template <class MyConvexProblemTypeTEMPLATE> +class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject { friend class MyBlockProblem; protected: @@ -54,7 +54,7 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject { * \param problem The problem including quadratic part and nonlinear/nonsmooth * part */ - IterateObject(const Bisection& bisection, ConvexProblemType& problem) + IterateObject(const Bisection& bisection, MyConvexProblemType& problem) : problem(problem), bisection(bisection), local_J(1.0, 0.0, problem.phi, 0, 0) {}; @@ -126,7 +126,7 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject { private: // problem data - ConvexProblemType& problem; + MyConvexProblemType& problem; // commonly used minimization stuff Bisection bisection; diff --git a/src/myconvexproblem.hh b/src/myconvexproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..b7c06e836ebc35d91fe51fe3fbf7b7925585e7ba --- /dev/null +++ b/src/myconvexproblem.hh @@ -0,0 +1,55 @@ +// Based on dune/tnnmg/problem-classes/convexproblem.hh + +#include <dune/istl/bcrsmatrix.hh> + +#include <dune/tnnmg/problem-classes/nonlinearity.hh> + +/** \brief General convex problem for a Truncated Nonsmooth Newton Multigrid + (TNNMG) solver + + \tparam NonlinearityTypeTEMPLATE The type used to implement the nonlinearity + \tparam MatrixTypeTEMPLATE The type used for the matrix of the quadratic + part +*/ +template <class NonlinearityTypeTEMPLATE = Nonlinearity< + Dune::FieldVector<double, 1>, Dune::FieldMatrix<double, 1, 1>>, + class MatrixTypeTEMPLATE = + Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>> +class MyConvexProblem { +public: + typedef NonlinearityTypeTEMPLATE NonlinearityType; + + typedef typename NonlinearityType::VectorType VectorType; + typedef MatrixTypeTEMPLATE MatrixType; + typedef typename NonlinearityType::LocalVectorType LocalVectorType; + typedef typename MatrixType::block_type LocalMatrixType; + + static const int block_size = NonlinearityType::block_size; + + /** \brief Constructor with the problem components + + \param a A scalar factor in front of the quadratic part (the quadratic + part includes a factor of 1/2 already) + \param A The matrix of the quadratic part + \param am A scalar factor in front of the optional rank-one matrix + \param Am A rank-one matrix given by a single vector. The matrix is + AmAm^T + \param phi The nonlinearity + \param f The linear functional + \param u The solution vector + */ + MyConvexProblem(double a, const MatrixType& A, double am, + const VectorType& Am, NonlinearityType& phi, + const VectorType& f, VectorType& u) + : a(a), A(A), am(am), Am(Am), phi(phi), f(f), u(u) {}; + + double a; + const MatrixType& A; + double am; + const VectorType& Am; + NonlinearityType& phi; + + const VectorType& f; + + VectorType& u; +}; diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc index 066d9dfe726187cfaf5c9cd59a16882955f5ea9f..fd86fae31c01d33480d8a2939d5c89f0e696e246 100644 --- a/src/one-body-sample.cc +++ b/src/one-body-sample.cc @@ -28,13 +28,14 @@ #include <dune/grid/common/mcmgmapper.hh> -#include <dune/tnnmg/problem-classes/convexproblem.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 "myconvexproblem.hh" #include "myblockproblem.hh" int const dim = 2; @@ -141,13 +142,14 @@ int main() { typedef ZeroNonlinearity<SmallVector, SmallMatrix> NonlinearityType; NonlinearityType phi; - typedef ConvexProblem<NonlinearityType, OperatorType> ConvexProblemType; + typedef MyConvexProblem<NonlinearityType, OperatorType> + MyConvexProblemType; VectorType unused_vector(grid.size(grid.maxLevel(), dim), 0); - ConvexProblemType myConvexProblem(1, stiffnessMatrix, 0, unused_vector, - phi, f, u1); + MyConvexProblemType myConvexProblem(1, stiffnessMatrix, 0, unused_vector, + phi, f, u1); - typedef MyBlockProblem<ConvexProblemType> MyBlockProblemType; + typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType; MyBlockProblemType myBlockProblem(myConvexProblem); GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep;