partially revert the last commit. It contained more files than intended

[[Imported from SVN: r7868]]

dune/gfe/cosseratenergystiffness.hh

@@ -82,8 +82,6 @@ public:  // for testing
         // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
         // However, since the directors of a given unit quaternion are the _columns_ of the
         // corresponding orthogonal matrix, we need to invert the i and j indices
-        //
-        // So, if I am not mistaken, DR[i][j][k] contains \partial M_ij / \partial k
         Tensor3<double,3 , 3, 4> dd_dq;
         value.q.getFirstDerivativesOfDirectors(dd_dq);
         
@@ -138,8 +136,7 @@ public:
 
     RT curvatureEnergy(const Tensor3<double,3,3,3>& DR) const
     {
-        //return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
-        return mu_ * std::pow(L_c_ * DR.frobenius_norm(),q_);
+        return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
     }
 
     RT bendingEnergy(const Dune::FieldMatrix<double,dim,dim>& R, const Tensor3<double,3,3,3>& DR) const
@@ -183,7 +180,7 @@ public:
     /** \brief Curvature exponent */
     double q_;
     
-    const BoundaryPatch<GridView>* neumannBoundary_;
+    
 };
 
 template <class GridView, class LocalFiniteElement, int dim>
@@ -198,7 +195,7 @@ energy(const Entity& element,
     LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> localGeodesicFEFunction(localFiniteElement,
                                                                                                       localSolution);
 
-    int quadOrder = 2*gridDim;
+    int quadOrder = 1;//gridDim;
 
     const Dune::QuadratureRule<double, gridDim>& quad 
         = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
@@ -279,37 +276,7 @@ energy(const Entity& element,
             DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
 
     }
-    
-    
-    //////////////////////////////////////////////////////////////////////////////
-    //   Assemble boundary contributions
-    //////////////////////////////////////////////////////////////////////////////
-
-    for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) {
-        
-        if (not neumannBoundary_ or not neumannBoundary_->contains(*it))
-            continue;
-        
-        const Dune::QuadratureRule<double, gridDim-1>& quad 
-            = Dune::QuadratureRules<double, gridDim-1>::rule(it->type(), quadOrder);
-    
-        for (size_t pt=0; pt<quad.size(); pt++) {
-        
-            // Local position of the quadrature point
-            const Dune::FieldVector<double,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position());
-        
-            const double integrationElement = it->geometry().integrationElement(quad[pt].position());
 
-            // The value of the local function
-            RigidBodyMotion<dim> value = localGeodesicFEFunction.evaluate(quadPos);
-
-            // Constant Neumann force in z-direction
-            energy += thickness_ * 40 * value.r[2] * quad[pt].weight() * integrationElement;
-            
-        }
-        
-    }
-    
     return energy;
 }
 

dune/gfe/riemanniantrsolver.cc

@@ -62,9 +62,7 @@ setup(const GridType& grid,
     //   Create a multigrid solver
     // ////////////////////////////////
 
-#if 0
     // First create a Gauss-seidel base solver
-    
     TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
 
     EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep);
@@ -74,11 +72,6 @@ setup(const GridType& grid,
                                                                             baseTolerance,
                                                                             baseEnergyNorm,
                                                                             Solver::QUIET);
-#else
-    QuadraticIPOptSolver<MatrixType,CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>();
-    baseSolver->tolerance_ = baseTolerance;
-    baseSolver->verbosity_ = Solver::QUIET;
-#endif
 
     // Make pre and postsmoothers
     TrustRegionGSStep<MatrixType, CorrectionType>* presmoother  = new TrustRegionGSStep<MatrixType, CorrectionType>;
@@ -112,21 +105,12 @@ setup(const GridType& grid,
 
     h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
 
-#if 0
     innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
                                                                                                    innerIterations_,
                                                                                                    innerTolerance_,
                                                                                                    h1SemiNorm_,
                                                                                                  Solver::FULL));
 
-#else
-    QuadraticIPOptSolver<MatrixType,CorrectionType>* solver = new QuadraticIPOptSolver<MatrixType,CorrectionType>();
-    solver->tolerance_ = innerTolerance_;
-    solver->verbosity_ = Solver::REDUCED;
-    innerSolver_ = std::shared_ptr<QuadraticIPOptSolver<MatrixType,CorrectionType> >(solver);
-#endif
-
-    
     // Write all intermediate solutions, if requested
     if (instrumented_
         && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
@@ -206,7 +190,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
     std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles((mgStep) 
                                                                                          ? mgStep->numLevels_
-                                                                                         : 1);
+                                                                                         : 0);
 
    // /////////////////////////////////////////////////////
     //   Set up the log file, if requested
@@ -258,6 +242,8 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
         // The right hand side is the _negative_ gradient
         rhs *= -1;
 
+/*        std::cout << "rhs:\n" << rhs << std::endl;
+        std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;*/
         
         // //////////////////////////////////////////////////////////////////////
         //   Modify matrix and right-hand side to account for Dirichlet values
@@ -291,17 +277,10 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
         }
 
-        //std::cout << "rhs:\n" << rhs << std::endl;
-        //std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;
-        //printmatrix(std::cout, (*hessianMatrix_), "hessian", "--");
-        //writeMatrixToMatlab(*hessianMatrix_, "cosserat_hessian");
-        
-        //mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
-        std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->setProblem(*hessianMatrix_, corr, rhs);
+        mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
         
         trustRegionObstacles.back() = trustRegion.obstacles();
-        //mgStep->obstacles_ = &trustRegionObstacles;
-        std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->obstacles_ = &trustRegionObstacles.back();
+        mgStep->obstacles_ = &trustRegionObstacles;
 
         innerSolver_->preprocess();
         
@@ -441,7 +420,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
         }
 
         if (energy >= oldEnergy &&
-            (std::abs((oldEnergy-energy)/energy) < 1e-9 || std::abs(modelDecrease/energy) < 1e-9)) {
+            (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
             if (this->verbosity_ == NumProc::FULL)
                 std::cout << "Suspecting rounding problems" << std::endl;