From 3c17d05cb6a2d79c1969b54f441da64af3da6dc4 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 6 Dec 2015 08:03:12 +0100
Subject: [PATCH] Use TupleVector to store the iterates, instead of two
separate std::vectors
---
dune/gfe/mixedriemanniantrsolver.cc | 42 ++++++++++++----------------
dune/gfe/mixedriemanniantrsolver.hh | 21 ++++++--------
src/mixed-cosserat-continuum.cc | 43 +++++++++++++++--------------
3 files changed, 50 insertions(+), 56 deletions(-)
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index 9e770d1b..4255da00 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -35,8 +35,7 @@ template <class GridType,
void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,TargetSpace1>::
setup(const GridType& grid,
const MixedGFEAssembler<Basis0, TargetSpace0, Basis1, TargetSpace1>* assembler,
- const SolutionType0& x0,
- const SolutionType1& x1,
+ const SolutionType& x,
const Dune::BitSetVector<blocksize0>& dirichletNodes0,
const Dune::BitSetVector<blocksize1>& dirichletNodes1,
double tolerance,
@@ -55,8 +54,7 @@ setup(const GridType& grid,
grid_ = &grid;
assembler_ = assembler;
- x0_ = x0;
- x1_ = x1;
+ x_ = x;
tolerance_ = tolerance;
maxTrustRegionSteps_ = maxTrustRegionSteps;
initialTrustRegionRadius_ = initialTrustRegionRadius;
@@ -320,7 +318,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
using namespace Dune::TypeTree::Indices;
- double oldEnergy = assembler_->computeEnergy(x0_, x1_);
+ double oldEnergy = assembler_->computeEnergy(x_[_0], x_[_1]);
oldEnergy = mpiHelper.getCollectiveCommunication().sum(oldEnergy);
bool recomputeGradientHessian = true;
@@ -345,17 +343,17 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
std::cout << "----------------------------------------------------" << std::endl;
}
- CorrectionType0 corr0(x0_.size());
+ CorrectionType0 corr0(x_[_0].size());
corr0 = 0;
- CorrectionType1 corr1(x1_.size());
+ CorrectionType1 corr1(x_[_1].size());
corr1 = 0;
Dune::Timer assemblyTimer;
if (recomputeGradientHessian) {
- assembler_->assembleGradientAndHessian(x0_,
- x1_,
+ assembler_->assembleGradientAndHessian(x_[_0],
+ x_[_1],
rhs0,
rhs1,
(*hessianMatrix_)[_0][_0],
@@ -481,15 +479,14 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
- SolutionType0 newIterate0 = x0_;
- for (size_t j=0; j<newIterate0.size(); j++)
- newIterate0[j] = TargetSpace0::exp(newIterate0[j], corr0[j]);
+ SolutionType newIterate = x_;
+ for (size_t j=0; j<newIterate[_0].size(); j++)
+ newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr0[j]);
- SolutionType1 newIterate1 = x1_;
- for (size_t j=0; j<newIterate1.size(); j++)
- newIterate1[j] = TargetSpace1::exp(newIterate1[j], corr1[j]);
+ for (size_t j=0; j<newIterate[_1].size(); j++)
+ newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr1[j]);
- double energy = assembler_->computeEnergy(newIterate0,newIterate1);
+ double energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
energy = mpiHelper.getCollectiveCommunication().sum(energy);
// compute the model decrease
@@ -532,8 +529,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
if (this->verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
- x0_ = newIterate0;
- x1_ = newIterate1;
+ x_ = newIterate;
break;
}
@@ -543,8 +539,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
if ( (oldEnergy-energy) / modelDecrease > 0.9) {
// very successful iteration
- x0_ = newIterate0;
- x1_ = newIterate1;
+ x_ = newIterate;
trustRegion0.scale(2);
trustRegion1.scale(2);
@@ -556,8 +551,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
} else if ( (oldEnergy-energy) / modelDecrease > 0.01
|| std::abs(oldEnergy-energy) < 1e-12) {
// successful iteration
- x0_ = newIterate0;
- x1_ = newIterate1;
+ x_ = newIterate;
// current energy becomes 'oldEnergy' for the next iteration
oldEnergy = energy;
@@ -581,8 +575,8 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
DuneFunctionsBasis<Basis1> fufemBasis1(assembler_->basis1_);
std::stringstream iAsAscii;
iAsAscii << i+1;
- CosseratVTKWriter<GridType>::template writeMixed<DuneFunctionsBasis<Basis0>, DuneFunctionsBasis<Basis1> >(fufemBasis0,x0_,
- fufemBasis1,x1_,
+ CosseratVTKWriter<GridType>::template writeMixed<DuneFunctionsBasis<Basis0>, DuneFunctionsBasis<Basis1> >(fufemBasis0,x_[_0],
+ fufemBasis1,x_[_1],
"mixed-cosserat_iterate_" + iAsAscii.str());
if (rank==0)
diff --git a/dune/gfe/mixedriemanniantrsolver.hh b/dune/gfe/mixedriemanniantrsolver.hh
index 4807dd90..4cfa5044 100644
--- a/dune/gfe/mixedriemanniantrsolver.hh
+++ b/dune/gfe/mixedriemanniantrsolver.hh
@@ -11,6 +11,8 @@
#include <dune/grid/utility/globalindexset.hh>
+#include <dune/functions/common/tuplevector.hh>
+
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -43,8 +45,7 @@ class MixedRiemannianTrustRegionSolver
Dune::MultiTypeBlockVector<MatrixType10,MatrixType11> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize0> > CorrectionType0;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize1> > CorrectionType1;
- typedef std::vector<TargetSpace0> SolutionType0;
- typedef std::vector<TargetSpace1> SolutionType1;
+ typedef Dune::Functions::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > SolutionType;
#if 0
#ifdef SECOND_ORDER
@@ -68,8 +69,7 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const MixedGFEAssembler<Basis0, TargetSpace0, Basis1, TargetSpace1>* assembler,
- const SolutionType0& x0,
- const SolutionType1& x1,
+ const SolutionType& x,
const Dune::BitSetVector<blocksize0>& dirichletNodes0,
const Dune::BitSetVector<blocksize1>& dirichletNodes1,
double tolerance,
@@ -104,16 +104,14 @@ public:
#endif
void solve();
- void setInitialIterate(const SolutionType0& x0,
- const SolutionType1& x1)
+ void setInitialIterate(const SolutionType& x)
{
- x0_ = x0;
- x1_ = x1;
+ x_ = x;
}
- std::tuple<SolutionType0,SolutionType1> getSol() const
+ SolutionType getSol() const
{
- return std::make_tuple(x0_,x1_);
+ return x_;
}
protected:
@@ -124,8 +122,7 @@ protected:
const GridType* grid_;
/** \brief The solution vectors */
- SolutionType0 x0_;
- SolutionType1 x1_;
+ SolutionType x_;
/** \brief The initial trust-region radius in the maximum-norm */
double initialTrustRegionRadius_;
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
index 42653536..8d208ce4 100644
--- a/src/mixed-cosserat-continuum.cc
+++ b/src/mixed-cosserat-continuum.cc
@@ -22,6 +22,7 @@
#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/functions/common/tuplevector.hh>
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/fufem/boundarypatch.hh>
@@ -81,8 +82,9 @@ int main (int argc, char *argv[]) try
<< std::endl << "sys.path.append('/home/sander/dune/dune-gfe/problems/')"
<< std::endl;
- typedef std::vector<RealTuple<double,3> > DeformationSolutionType;
- typedef std::vector<Rotation<double,3> > OrientationSolutionType;
+ using namespace Dune::TypeTree::Indices;
+ typedef Dune::Functions::TupleVector<std::vector<RealTuple<double,3> >,
+ std::vector<Rotation<double,3> > > SolutionType;
// parse data file
ParameterTree parameterSet;
@@ -225,7 +227,9 @@ int main (int argc, char *argv[]) try
// Initial iterate
// //////////////////////////
- DeformationSolutionType xDisp(deformationFEBasis.indexSet().size());
+ SolutionType x;
+
+ x[_0].resize(deformationFEBasis.indexSet().size());
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
@@ -233,10 +237,10 @@ int main (int argc, char *argv[]) try
std::vector<FieldVector<double,3> > v;
::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
- for (size_t i=0; i<xDisp.size(); i++)
- xDisp[i] = v[i];
+ for (size_t i=0; i<x[_0].size(); i++)
+ x[_0][i] = v[i];
- OrientationSolutionType xOrient(orientationFEBasis.indexSet().size());
+ x[_1].resize(orientationFEBasis.indexSet().size());
#if 0
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
@@ -252,8 +256,8 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
- CosseratVTKWriter<GridType>::writeMixed<FufemDeformationFEBasis,FufemOrientationFEBasis>(fufemDeformationFEBasis,xDisp,
- fufemOrientationFEBasis,xOrient,
+ CosseratVTKWriter<GridType>::writeMixed<FufemDeformationFEBasis,FufemOrientationFEBasis>(fufemDeformationFEBasis,x[_0],
+ fufemOrientationFEBasis,x[_1],
resultPath + "mixed-cosserat_homotopy_0");
for (int i=0; i<numHomotopySteps; i++) {
@@ -306,8 +310,7 @@ int main (int argc, char *argv[]) try
OrientationFEBasis, Rotation<double,3> > solver;
solver.setup(*grid,
&assembler,
- xDisp,
- xOrient,
+ x,
deformationDirichletDofs,
orientationDirichletDofs,
tolerance,
@@ -343,28 +346,28 @@ int main (int argc, char *argv[]) try
std::vector<FieldMatrix<double,3,3> > dOV;
::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
- for (size_t j=0; j<xDisp.size(); j++)
+ for (size_t j=0; j<x[_0].size(); j++)
if (deformationDirichletNodes[j][0])
- xDisp[j] = ddV[j];
+ x[_0][j] = ddV[j];
- for (size_t j=0; j<xOrient.size(); j++)
+ for (size_t j=0; j<x[_1].size(); j++)
if (orientationDirichletNodes[j][0])
- xOrient[j].set(dOV[j]);
+ x[_1][j].set(dOV[j]);
// /////////////////////////////////////////////////////
// Solve!
// /////////////////////////////////////////////////////
- solver.setInitialIterate(xDisp,xOrient);
+ solver.setInitialIterate(x);
solver.solve();
- std::tie(xDisp,xOrient) = solver.getSol();
+ x = solver.getSol();
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
- CosseratVTKWriter<GridType>::writeMixed<FufemDeformationFEBasis,FufemOrientationFEBasis>(fufemDeformationFEBasis,xDisp,
- fufemOrientationFEBasis,xOrient,
+ CosseratVTKWriter<GridType>::writeMixed<FufemDeformationFEBasis,FufemOrientationFEBasis>(fufemDeformationFEBasis,x[_0],
+ fufemOrientationFEBasis,x[_1],
resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
}
@@ -376,9 +379,9 @@ int main (int argc, char *argv[]) try
// finally: compute the average deformation of the Neumann boundary
// That is what we need for the locking tests
FieldVector<double,3> averageDef(0);
- for (size_t i=0; i<xDisp.size(); i++)
+ for (size_t i=0; i<x[_0].size(); i++)
if (neumannNodes[i][0])
- averageDef += xDisp[i].globalCoordinates();
+ averageDef += x[_0][i].globalCoordinates();
averageDef /= neumannNodes.count();
if (mpiHelper.rank()==0)
--
GitLab