a method to compute the gradient by fd approximation

[[Imported from SVN: r1524]]

src/rodassembler.cc

@@ -1183,6 +1183,106 @@ assembleGradient(const std::vector<Configuration>& sol,
 
 }
 
+template <class GridType>
+void RodAssembler<GridType>::
+assembleGradientFD(const std::vector<Configuration>& sol,
+                 Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
+{
+    double eps = 1e-2;
+    if (grad.size() != sol.size())
+        grad.resize(sol.size());
+
+    // ///////////////////////////////////////////////////////////
+    //   Compute gradient by finite-difference approximation
+    // ///////////////////////////////////////////////////////////
+    std::vector<Configuration> forwardSolution = sol;
+    std::vector<Configuration> backwardSolution = sol;
+
+    // ////////////////////////////////////////////////////
+    //   Create local assembler
+    // ////////////////////////////////////////////////////
+
+    Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
+    Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
+    RodLocalStiffness<GridType,double> localStiffness(K, A);
+
+    // //////////////////////////////////////////////////////////////////
+    //   Store pointers to all elements so we can access them by index
+    // //////////////////////////////////////////////////////////////////
+
+    const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet();
+
+    std::vector<EntityPointer> elements(grid_->size(0),grid_->template leafbegin<0>());
+    
+    typename GridType::template Codim<0>::LeafIterator eIt    = grid_->template leafbegin<0>();
+    typename GridType::template Codim<0>::LeafIterator eEndIt = grid_->template leafend<0>();
+
+    for (; eIt!=eEndIt; ++eIt)
+        elements[indexSet.index(*eIt)] = eIt;
+
+    std::vector<Configuration> localReferenceConfiguration(2);
+    std::vector<Configuration> localSolution(2);
+
+    // ///////////////////////////////////////////////////////////////
+    //   Loop over all blocks of the outer matrix
+    // ///////////////////////////////////////////////////////////////
+    for (size_t i=0; i<sol.size(); i++) {
+
+        for (int j=0; j<6; j++) {
+
+            double forwardEnergy = 0;
+            double backwardEnergy = 0;
+            
+            infinitesimalVariation(forwardSolution[i], eps, j);
+            infinitesimalVariation(backwardSolution[i], -eps, j);
+            
+            if (i != grid_->size(1)-1) {
+                
+                localReferenceConfiguration[0] = referenceConfiguration_[i];
+                localReferenceConfiguration[1] = referenceConfiguration_[i+1];
+                
+                localSolution[0] = forwardSolution[i];
+                localSolution[1] = forwardSolution[i+1];
+                
+                forwardEnergy += localStiffness.energy(elements[i], localSolution, localReferenceConfiguration);
+                
+                localSolution[0] = backwardSolution[i];
+                localSolution[1] = backwardSolution[i+1];
+                
+                backwardEnergy += localStiffness.energy(elements[i], localSolution, localReferenceConfiguration);
+                
+            } 
+            
+            if (i != 0) {
+                
+                localReferenceConfiguration[0] = referenceConfiguration_[i-1];
+                localReferenceConfiguration[1] = referenceConfiguration_[i];
+                
+                localSolution[0] = forwardSolution[i-1];
+                localSolution[1] = forwardSolution[i];
+                
+                forwardEnergy += localStiffness.energy(elements[i-1], localSolution, localReferenceConfiguration);
+                
+                localSolution[0] = backwardSolution[i-1];
+                localSolution[1] = backwardSolution[i];
+                
+                backwardEnergy += localStiffness.energy(elements[i-1], localSolution, localReferenceConfiguration);
+                
+            } 
+            
+            // Second derivative
+            grad[i][j] = (forwardEnergy - backwardEnergy) / (2*eps);
+            
+            forwardSolution[i]  = sol[i];
+            backwardSolution[i] = sol[i];
+            
+        }
+
+    }
+
+
+}
+
 
 template <class GridType>
 double RodAssembler<GridType>::

src/rodassembler.hh

@@ -108,13 +108,6 @@ public:
         return u;
     }
         
-
-    //! should allow to assmble boundary conditions only
-//     template<typename Tag>
-//     void assembleBoundaryCondition (const Entity& e, int k=1)
-//     {
-//     }
-    
 };
 
 
@@ -253,6 +246,9 @@ public:
         void assembleGradient(const std::vector<Configuration>& sol,
                               Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
 
+        void assembleGradientFD(const std::vector<Configuration>& sol,
+                                Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
+
         /** \brief Compute the energy of a deformation state */
         double computeEnergy(const std::vector<Configuration>& sol) const;