diff --git a/dirneucoupling.cc b/dirneucoupling.cc index cbb9da1bcd0c95cac4d0ef3aba9b2a352bf8ac74..a576f9407011cb04ba680e13b7796ff91ad9d5fd 100644 --- a/dirneucoupling.cc +++ b/dirneucoupling.cc @@ -14,7 +14,7 @@ #include <dune/common/configparser.hh> #include <dune/ag-common/multigridstep.hh> -#include <dune/ag-common/iterativesolver.hh> +#include <dune/ag-common/solvers/loopsolver.hh> #include <dune/ag-common/projectedblockgsstep.hh> #ifdef HAVE_IPOPT #include <dune/ag-common/quadraticipopt.hh> @@ -331,7 +331,7 @@ int main (int argc, char *argv[]) try EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep); - IterativeSolver<VectorType> solver(&multigridStep, + LoopSolver<VectorType> solver(&multigridStep, // IPOpt doesn't like to be started in the solution (numLevels!=1) ? multigridIterations : 1, mgTolerance, @@ -347,7 +347,7 @@ int main (int argc, char *argv[]) try multigridStep.mgTransfer_.resize(toplevel); for (int i=0; i<multigridStep.mgTransfer_.size(); i++){ - TruncatedMGTransfer<VectorType>* newTransferOp = new TruncatedMGTransfer<VectorType>; + TruncatedCompressedMGTransfer<VectorType>* newTransferOp = new TruncatedCompressedMGTransfer<VectorType>; newTransferOp->setup(grid,i,i+1); multigridStep.mgTransfer_[i] = newTransferOp; } diff --git a/src/rodsolver.cc b/src/rodsolver.cc index 894daace796c3f14f6bbb4d3e013a0f6639b109a..23868b614784d429bc26e3120af393a5f9617047 100644 --- a/src/rodsolver.cc +++ b/src/rodsolver.cc @@ -17,34 +17,13 @@ #include "configuration.hh" #include "quaternion.hh" +#include "maxnormtrustregion.hh" // for debugging #include "../test/fdcheck.hh" // Number of degrees of freedom: // 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod -//const int blocksize = 6; - -template <class GridType> -void RodSolver<GridType>:: -setTrustRegionObstacles(double trustRegionRadius, - std::vector<BoxConstraint<field_type,blocksize> >& trustRegionObstacles) -{ - for (int j=0; j<trustRegionObstacles.size(); j++) { - - for (int k=0; k<blocksize; k++) { - - trustRegionObstacles[j].lower(k) = -trustRegionRadius; - trustRegionObstacles[j].upper(k) = trustRegionRadius; - - } - - } - - //std::cout << trustRegionObstacles << std::endl; -// exit(0); -} - template <class GridType> void RodSolver<GridType>::setup(const GridType& grid, @@ -91,12 +70,11 @@ void RodSolver<GridType>::setup(const GridType& grid, EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep); - IterativeSolver<CorrectionType>* baseSolver - = new IterativeSolver<CorrectionType>(baseSolverStep, - baseIt_, - baseTolerance_, - baseEnergyNorm, - Solver::QUIET); + LoopSolver<CorrectionType>* baseSolver = new LoopSolver<CorrectionType>(baseSolverStep, + baseIt_, + baseTolerance_, + baseEnergyNorm, + Solver::QUIET); // Make pre and postsmoothers TrustRegionGSStep<MatrixType, CorrectionType>* presmoother = new TrustRegionGSStep<MatrixType, CorrectionType>; @@ -111,7 +89,6 @@ void RodSolver<GridType>::setup(const GridType& grid, mmgStep->postsmoother_ = postsmoother; mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType>(); mmgStep->hasObstacle_ = &hasObstacle_; - mmgStep->obstacles_ = &trustRegionObstacles_; mmgStep->verbosity_ = Solver::QUIET; // ////////////////////////////////////////////////////////////////////////////////////// @@ -129,7 +106,7 @@ void RodSolver<GridType>::setup(const GridType& grid, h1SemiNorm_ = new H1SemiNorm<CorrectionType>(**A); - mmgSolver_ = new IterativeSolver<CorrectionType>(mmgStep, + mmgSolver_ = new LoopSolver<CorrectionType>(mmgStep, multigridIterations_, qpTolerance_, h1SemiNorm_, @@ -159,11 +136,6 @@ void RodSolver<GridType>::setup(const GridType& grid, for (int i=0; i<hasObstacle_.size(); i++) hasObstacle_[i].resize(grid_->size(i, 1),true); - trustRegionObstacles_.resize(numLevels); - - for (int i=0; i<numLevels; i++) - trustRegionObstacles_[i].resize(grid_->size(i,1)); - // //////////////////////////////////// // Create the transfer operators // //////////////////////////////////// @@ -183,10 +155,10 @@ void RodSolver<GridType>::setup(const GridType& grid, template <class GridType> void RodSolver<GridType>::solve() { - using namespace Dune; - - double trustRegionRadius = initialTrustRegionRadius_; + MaxNormTrustRegion<blocksize> trustRegion(x_.size(), initialTrustRegionRadius_); + std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles(dynamic_cast<MultigridStep<MatrixType,CorrectionType>*>(mmgSolver_->iterationStep_)->numLevels_); + // ///////////////////////////////////////////////////// // Set up the log file, if requested // ///////////////////////////////////////////////////// @@ -195,7 +167,7 @@ void RodSolver<GridType>::solve() fp = fopen("statistics", "w"); if (!fp) - DUNE_THROW(IOError, "Couldn't open statistics file for writing!"); + DUNE_THROW(Dune::IOError, "Couldn't open statistics file for writing!"); } @@ -207,7 +179,7 @@ void RodSolver<GridType>::solve() if (this->verbosity_ == FULL) { std::cout << "----------------------------------------------------" << std::endl; std::cout << " Trust-Region Step Number: " << i - << ", radius: " << trustRegionRadius + << ", radius: " << trustRegion.radius() << ", energy: " << rodAssembler_->computeEnergy(x_) << std::endl; std::cout << "----------------------------------------------------" << std::endl; } @@ -225,11 +197,11 @@ void RodSolver<GridType>::solve() rhs *= -1; - // Create trust-region obstacle on maxlevel - setTrustRegionObstacles(trustRegionRadius, - trustRegionObstacles_[grid_->maxLevel()]); - dynamic_cast<MultigridStep<MatrixType,CorrectionType>*>(mmgSolver_->iterationStep_)->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1); + + trustRegionObstacles.back() = trustRegion.obstacles(); + dynamic_cast<MonotoneMGStep<MatrixType, CorrectionType>*>(mmgSolver_->iterationStep_)->obstacles_ = &trustRegionObstacles; + mmgSolver_->preprocess(); @@ -248,7 +220,7 @@ void RodSolver<GridType>::solve() if (instrumented_) { fprintf(fp, "Trust-region step: %d, trust-region radius: %g\n", - i, trustRegionRadius); + i, trustRegion.radius()); // /////////////////////////////////////////////////////////////// // Compute and measure progress against the exact solution @@ -276,7 +248,7 @@ void RodSolver<GridType>::solve() FILE* fpInt = fopen(iSolFilename, "rb"); if (!fpInt) - DUNE_THROW(IOError, "Couldn't open intermediate solution"); + DUNE_THROW(Dune::IOError, "Couldn't open intermediate solution"); for (int k=0; k<intermediateSol.size(); k++) for (int l=0; l<blocksize; l++) fread(&intermediateSol[k][l], sizeof(double), 1, fpInt); @@ -396,7 +368,7 @@ void RodSolver<GridType>::solve() // very successful iteration x_ = newIterate; - trustRegionRadius *= 2; + trustRegion.scale(2); } else if ( (oldEnergy-energy) / modelDecrease > 0.01 || std::abs(oldEnergy-energy) < 1e-12) { @@ -405,7 +377,7 @@ void RodSolver<GridType>::solve() } else { // unsuccessful iteration - trustRegionRadius /= 2; + trustRegion.scale(0.5); if (this->verbosity_ == FULL) std::cout << "Unsuccessful iteration!" << std::endl; } diff --git a/src/rodsolver.hh b/src/rodsolver.hh index d91dd4c924cff6f1eb743dd61e02ed8740d6d2f0..f7f70601e5a9def97e9f5065c933a54c9d7881cc 100644 --- a/src/rodsolver.hh +++ b/src/rodsolver.hh @@ -10,7 +10,7 @@ #include <dune/ag-common/boxconstraint.hh> #include <dune/ag-common/norms/h1seminorm.hh> -#include <dune/ag-common/iterativesolver.hh> +#include <dune/ag-common/solvers/loopsolver.hh> #include "rodassembler.hh" @@ -18,7 +18,7 @@ /** \brief Riemannian trust-region solver for 3d Cosserat rod problems */ template <class GridType> -class RodSolver : public Solver +class RodSolver : public IterativeSolver<std::vector<Configuration>, Dune::BitSetVector<6> > { const static int blocksize = 6; @@ -30,12 +30,10 @@ class RodSolver : public Solver typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType; typedef std::vector<Configuration> SolutionType; - static void setTrustRegionObstacles(double trustRegionRadius, - std::vector<BoxConstraint<field_type,blocksize> >& trustRegionObstacles); public: RodSolver() - : Solver(0,NumProc::FULL), + : IterativeSolver<std::vector<Configuration>, Dune::BitSetVector<6> >(0,100,NumProc::FULL), hessianMatrix_(NULL), h1SemiNorm_(NULL) {} @@ -104,10 +102,7 @@ protected: const RodAssembler<GridType>* rodAssembler_; /** \brief The multigrid solver */ - IterativeSolver<CorrectionType>* mmgSolver_; - - /** \brief The hierarchy of trust-region obstacles */ - std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles_; + LoopSolver<CorrectionType>* mmgSolver_; /** \brief Dummy fields containing 'true' everywhere. The multigrid step expects them :-( */ diff --git a/staticrod.cc b/staticrod.cc index af11d77c80f53743d30daf3ee86856195fe812be..8278528ca6bf493eba00b997eb5875cac2dd51df 100644 --- a/staticrod.cc +++ b/staticrod.cc @@ -10,10 +10,11 @@ #include <dune/ag-common/boundarypatch.hh> #include <dune/ag-common/projectedblockgsstep.hh> #include <dune/ag-common/solvers/mmgstep.hh> -#include <dune/ag-common/iterativesolver.hh> +#include <dune/ag-common/solvers/loopsolver.hh> #include <dune/ag-common/geomestimator.hh> #include <dune/ag-common/norms/energynorm.hh> -#include <dune/ag-common/contactobsrestrict.hh> +#include <dune/ag-common/mandelobsrestrictor.hh> +#include <dune/ag-common/transferoperators/truncatedcompressedmgtransfer.hh> #include "src/rodwriter.hh" #include "src/planarrodassembler.hh" @@ -177,7 +178,7 @@ int main (int argc, char *argv[]) try EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep); - IterativeSolver<VectorType> baseSolver(&baseSolverStep, + LoopSolver<VectorType> baseSolver(&baseSolverStep, baseIt, baseTolerance, &baseEnergyNorm, @@ -196,23 +197,23 @@ int main (int argc, char *argv[]) try multigridStep.postsmoother_ = &postsmoother; multigridStep.hasObstacle_ = &hasObstacle; multigridStep.obstacles_ = &trustRegionObstacles; - multigridStep.obstacleRestrictor_ = new ContactObsRestriction<VectorType>; + multigridStep.obstacleRestrictor_ = new MandelObstacleRestrictor<VectorType>; // Create the transfer operators multigridStep.mgTransfer_.resize(maxlevel); for (int i=0; i<multigridStep.mgTransfer_.size(); i++){ - TruncatedMGTransfer<VectorType>* newTransferOp = new TruncatedMGTransfer<VectorType>; + TruncatedCompressedMGTransfer<VectorType>* newTransferOp = new TruncatedCompressedMGTransfer<VectorType>; newTransferOp->setup(rod,i,i+1); multigridStep.mgTransfer_[i] = newTransferOp; } EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep); - IterativeSolver<VectorType> solver(&multigridStep, - numIt, - tolerance, - &energyNorm, - Solver::QUIET); + LoopSolver<VectorType> solver(&multigridStep, + numIt, + tolerance, + &energyNorm, + Solver::QUIET); // /////////////////////////////////////////////////// // Do a homotopy of the material parameters diff --git a/staticrod2.cc b/staticrod2.cc index db1ca3d4b256198042cb35a7cad78c4a3578898a..0bc2250d04f61b280714a2ffd8faa8a80717b5da 100644 --- a/staticrod2.cc +++ b/staticrod2.cc @@ -12,10 +12,11 @@ #include <dune/ag-common/boundarypatch.hh> #include <dune/ag-common/projectedblockgsstep.hh> #include <dune/ag-common/solvers/mmgstep.hh> -#include <dune/ag-common/iterativesolver.hh> +#include <dune/ag-common/solvers/loopsolver.hh> #include <dune/ag-common/geomestimator.hh> #include <dune/ag-common/norms/energynorm.hh> -#include <dune/ag-common/contactobsrestrict.hh> +#include <dune/ag-common/mandelobsrestrictor.hh> +#include <dune/ag-common/transferoperators/truncatedcompressedmgtransfer.hh> #include "src/rodwriter.hh" #include "src/planarrodassembler.hh" @@ -111,7 +112,7 @@ int main (int argc, char *argv[]) try EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep); - IterativeSolver<VectorType> baseSolver(&baseSolverStep, + LoopSolver<VectorType> baseSolver(&baseSolverStep, baseIt, baseTolerance, &baseEnergyNorm, @@ -131,12 +132,12 @@ int main (int argc, char *argv[]) try multigridStep.hasObstacle_ = &hasObstacle; multigridStep.obstacles_ = &trustRegionObstacles; multigridStep.verbosity_ = Solver::QUIET; - multigridStep.obstacleRestrictor_ = new ContactObsRestriction<VectorType>; + multigridStep.obstacleRestrictor_ = new MandelObstacleRestrictor<VectorType>; EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep); - IterativeSolver<VectorType> solver(&multigridStep, + LoopSolver<VectorType> solver(&multigridStep, numIt, tolerance, &energyNorm, @@ -231,7 +232,7 @@ int main (int argc, char *argv[]) try multigridStep.mgTransfer_.resize(toplevel); for (int i=0; i<multigridStep.mgTransfer_.size(); i++){ - TruncatedMGTransfer<VectorType>* newTransferOp = new TruncatedMGTransfer<VectorType>; + TruncatedCompressedMGTransfer<VectorType>* newTransferOp = new TruncatedCompressedMGTransfer<VectorType>; newTransferOp->setup(grid,i,i+1); multigridStep.mgTransfer_[i] = newTransferOp; }