The LoopSolver class has been split off IterativeSolver

[[Imported from SVN: r3413]]

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;
     }

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;
         }

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 :-( */

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

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;
         }