From c8774a97b99907448e22a8bb467588ffac21c77a Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sat, 22 Feb 2014 21:42:56 +0100
Subject: [PATCH] [Cleanup] Inherit from BNGSP
---
dune/tectonic/myblockproblem.hh | 69 ++++++++++++++++++---------------
src/one-body-sample.cc | 2 +-
2 files changed, 39 insertions(+), 32 deletions(-)
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index f2baf72a..119ea0e7 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -11,6 +11,7 @@
#include <dune/solvers/common/interval.hh>
#include <dune/solvers/computeenergy.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
#include "ellipticenergy.hh"
#include "globalnonlinearity.hh"
@@ -19,13 +20,18 @@
/** \brief Base class for problems where each block can be solved with a
* modified gradient method */
-template <class ConvexProblem> class MyBlockProblem {
+template <class ConvexProblem>
+class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
+private:
+ typedef BlockNonlinearGSProblem<ConvexProblem> BNGSP;
+
public:
- using ConvexProblemType = ConvexProblem;
- using VectorType = typename ConvexProblem::VectorType;
- using MatrixType = typename ConvexProblem::MatrixType;
- using LocalVector = typename ConvexProblem::LocalVectorType;
- using LocalMatrix = typename ConvexProblem::LocalMatrixType;
+ using typename BNGSP::ConvexProblemType;
+ using typename BNGSP::LocalMatrixType;
+ using typename BNGSP::LocalVectorType;
+ using typename BNGSP::MatrixType;
+ using typename BNGSP::VectorType;
+
size_t static const block_size = ConvexProblem::block_size;
size_t static const coarse_block_size = block_size;
@@ -35,7 +41,7 @@ template <class ConvexProblem> class MyBlockProblem {
struct Linearization {
size_t static const block_size = coarse_block_size;
- using LocalMatrix = typename MyBlockProblem<ConvexProblem>::LocalMatrix;
+ using LocalMatrix = typename MyBlockProblem<ConvexProblem>::LocalMatrixType;
using MatrixType = Dune::BCRSMatrix<typename Linearization::LocalMatrix>;
using VectorType =
Dune::BlockVector<Dune::FieldVector<double, Linearization::block_size>>;
@@ -48,9 +54,10 @@ template <class ConvexProblem> class MyBlockProblem {
Dune::BitSetVector<Linearization::block_size> truncation;
};
- MyBlockProblem(Dune::ParameterTree const &parset,
- ConvexProblem const &problem)
- : parset(parset), problem(problem), localBisection() // NOTE: defaults
+ MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
+ : BNGSP(parset, problem),
+ parset_(parset),
+ localBisection() // TODO
{}
std::string getOutput(bool header = false) const {
@@ -92,8 +99,8 @@ template <class ConvexProblem> class MyBlockProblem {
MyDirectionalConvexFunction<
GlobalNonlinearity<MatrixType, VectorType>> const
- psi(computeDirectionalA(problem.A, v),
- computeDirectionalb(problem.A, problem.f, u, v), problem.phi, u, v);
+ psi(computeDirectionalA(problem_.A, v),
+ computeDirectionalb(problem_.A, problem_.f, u, v), problem_.phi, u, v);
Dune::Solvers::Interval<double> D;
psi.subDiff(0, D);
@@ -115,7 +122,7 @@ template <class ConvexProblem> class MyBlockProblem {
linearization.truncation.resize(u.size());
linearization.truncation.unsetAll();
for (size_t i = 0; i < u.size(); ++i) {
- if (problem.phi.regularity(i, u[i]) > 1e8) { // TODO: Make customisable
+ if (problem_.phi.regularity(i, u[i]) > 1e8) { // TODO: Make customisable
linearization.truncation[i] = true;
continue;
}
@@ -126,31 +133,31 @@ template <class ConvexProblem> class MyBlockProblem {
}
// construct sparsity pattern for linearization
- Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
- indices.import(problem.A);
- problem.phi.addHessianIndices(indices);
+ Dune::MatrixIndexSet indices(problem_.A.N(), problem_.A.M());
+ indices.import(problem_.A);
+ problem_.phi.addHessianIndices(indices);
// construct matrix from pattern and initialize it
indices.exportIdx(linearization.A);
linearization.A = 0.0;
- // compute quadratic part of hessian (linearization.A += problem.A)
- for (size_t i = 0; i < problem.A.N(); ++i) {
- auto const end = problem.A[i].end();
- for (auto it = problem.A[i].begin(); it != end; ++it)
+ // compute quadratic part of hessian (linearization.A += problem_.A)
+ for (size_t i = 0; i < problem_.A.N(); ++i) {
+ auto const end = problem_.A[i].end();
+ for (auto it = problem_.A[i].begin(); it != end; ++it)
linearization.A[i][it.index()] += *it;
}
// compute nonlinearity part of hessian
- problem.phi.addHessian(u, linearization.A);
+ problem_.phi.addHessian(u, linearization.A);
// compute quadratic part of gradient
linearization.b.resize(u.size());
- problem.A.mv(u, linearization.b);
- linearization.b -= problem.f;
+ problem_.A.mv(u, linearization.b);
+ linearization.b -= problem_.f;
// compute nonlinearity part of gradient
- problem.phi.addGradient(u, linearization.b);
+ problem_.phi.addGradient(u, linearization.b);
// -grad is needed for Newton step
linearization.b *= -1.0;
@@ -182,14 +189,14 @@ template <class ConvexProblem> class MyBlockProblem {
/** \brief Constructs and returns an iterate object */
IterateObject getIterateObject() {
- return IterateObject(parset, localBisection, problem);
+ return IterateObject(parset_, localBisection, problem_);
}
private:
- Dune::ParameterTree const &parset;
+ Dune::ParameterTree const &parset_;
// problem data
- ConvexProblem const &problem;
+ using BNGSP::problem_;
Bisection const localBisection;
@@ -220,7 +227,7 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
}
/** \brief Update the i-th block of the current iterate */
- void updateIterate(LocalVector const &ui, size_t i) {
+ void updateIterate(LocalVectorType const &ui, size_t i) {
u[i] = ui;
return;
}
@@ -231,11 +238,11 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
* \param ignore Set of degrees of freedom to leave untouched
*/
void solveLocalProblem(
- LocalVector &ui, size_t m,
+ LocalVectorType &ui, size_t m,
typename Dune::BitSetVector<block_size>::const_reference ignore) {
{
- LocalMatrix const *localA = nullptr;
- LocalVector localb(problem.f[m]);
+ LocalMatrixType const *localA = nullptr;
+ LocalVectorType localb(problem.f[m]);
auto const end = problem.A[m].end();
for (auto it = problem.A[m].begin(); it != end; ++it) {
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 96f64c53..0e8cb785 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -361,7 +361,7 @@ int main(int argc, char *argv[]) {
myGlobalNonlinearity->updateLogState(_alpha);
// NIT: Do we really need to pass u here?
- typename NonlinearFactory::ConvexProblem const convexProblem(
+ typename NonlinearFactory::ConvexProblem convexProblem(
1.0, velocityMatrix, *myGlobalNonlinearity, velocityRHS,
_velocityIterate);
typename NonlinearFactory::BlockProblem velocityProblem(parset,
--
GitLab