Merge branch 'feature/compute-gradient' into 'master'

Add assembleGradient() method to FE assemblers

See merge request !72

dune/elasticity/assemblers/feassembler.hh

@@ -54,6 +54,11 @@ public:
           localStiffness_(localStiffness)
     {}
 
+    /** \brief Assemble the functional gradient */
+    void assembleGradient(const VectorType& configuration,
+                          VectorType& gradient) const;
+
+
     /** \brief Assemble the tangent stiffness matrix and the functional gradient together
      *
      * This may be more efficient than computing them separately
@@ -71,6 +76,47 @@ public:
 }; // end class
 
 
+
+template <class Basis, class VectorType>
+void FEAssembler<Basis,VectorType>::
+assembleGradient(const VectorType& configuration,
+                 VectorType& gradient) const
+{
+  // create backends for multi index access
+  auto configurationBackend      = Dune::Functions::istlVectorBackend(configuration);
+  auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
+
+  gradientBackend.resize(basis_);
+
+  gradient = 0;
+
+  auto localView = basis_.localView();
+
+  for (const auto& element : elements(basis_.gridView(), partitionType))
+  {
+    localView.bind(element);
+    const auto size = localView.size();
+
+    // Extract local configuration
+    std::vector<field_type> localConfiguration(size);
+    std::vector<field_type> localGradient(size);
+
+    for (size_t i=0; i<size; i++)
+      localConfiguration[i] = configurationBackend[ localView.index(i) ];
+
+    // setup local matrix and gradient
+    localStiffness_->assembleGradient(localView, localConfiguration, localGradient);
+
+    for (size_t i=0; i<size; i++)
+    {
+      auto row = localView.index(i);
+      gradientBackend[row] += localGradient[i];
+    }
+  }
+}
+
+
+
 template <class Basis, class VectorType>
 template <class MatrixType>
 void FEAssembler<Basis,VectorType>::

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -41,6 +41,14 @@ public:
     virtual double energy (const LocalView& localView,
                            const std::vector<double>& localConfiguration) const;
 
+    /** \brief Assemble the local gradient at the current position
+    *
+    *     This uses the automatic differentiation toolbox ADOL_C.
+    */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<double>& localConfiguration,
+                                  std::vector<double>& localGradient);
+
     /** \brief Assemble the local stiffness matrix at the current position
      *
      *     This uses the automatic differentiation toolbox ADOL_C.
@@ -88,6 +96,31 @@ energy(const LocalView& localView,
 }
 
 
+
+template<class LocalView>
+void LocalADOLCStiffness<LocalView>::
+assembleGradient(const LocalView& localView,
+                 const std::vector<double>& localConfiguration,
+                 std::vector<double>& localGradient
+                )
+{
+  int rank = Dune::MPIHelper::getCommunication().rank();
+  // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
+  energy(localView, localConfiguration);
+
+  /////////////////////////////////////////////////////////////////
+  // Compute the energy gradient
+  /////////////////////////////////////////////////////////////////
+
+  // Compute the actual gradient
+  size_t nDoubles = localConfiguration.size();
+
+  // Compute gradient
+  localGradient.resize(nDoubles);
+  gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
+}
+
+
 template<class LocalView>
 void LocalADOLCStiffness<LocalView>::
 assembleGradientAndHessian(const LocalView& localView,
@@ -97,24 +130,19 @@ assembleGradientAndHessian(const LocalView& localView,
                           )
 {
     int rank = Dune::MPIHelper::getCommunication().rank();
-    // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
-    energy(localView, localConfiguration);
 
     /////////////////////////////////////////////////////////////////
     // Compute the energy gradient
     /////////////////////////////////////////////////////////////////
 
-    // Compute the actual gradient
-    size_t nDoubles = localConfiguration.size();
-
-    // Compute gradient
-    localGradient.resize(nDoubles);
-    gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
+    // this tapes the energy computation
+    assembleGradient(localView, localConfiguration, localGradient);
 
     /////////////////////////////////////////////////////////////////
     // Compute Hessian
     /////////////////////////////////////////////////////////////////
 
+    size_t nDoubles = localConfiguration.size();
     localHessian.setSize(nDoubles,nDoubles);
 
     // allocate raw doubles: ADOL-C requires "**double" arguments

dune/elasticity/assemblers/localfestiffness.hh

@@ -16,6 +16,13 @@ class LocalFEStiffness
 
 public:
 
+    /** \brief Assemble the local gradient at the current position
+    */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<RT>& localConfiguration,
+                                  std::vector<RT>& localGradient
+                                 ) = 0;
+
     /** \brief Assemble the local gradient and tangent matrix at the current position
     */
     virtual void assembleGradientAndHessian(const LocalView& localView,

dune/elasticity/assemblers/localhyperdualstiffness.hh

@@ -8,7 +8,7 @@
 
 namespace Dune::Elasticity {
 
-/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers 
+/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers
  */
 template<class LocalView>
 class LocalHyperDualStiffness
@@ -27,6 +27,15 @@ public:
     virtual double energy (const LocalView& localView,
                            const std::vector<double>& localConfiguration) const;
 
+
+    /** \brief Assemble the local gradient at the current position
+     *
+     *     This uses the automatic differentiation using hyperdual numbers
+     */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<double>& localConfiguration,
+                                  std::vector<double>& localGradient);
+
     /** \brief Assemble the local stiffness matrix at the current position
      *
      *     This uses the automatic differentiation using hyperdual numbers
@@ -36,7 +45,7 @@ public:
                                             std::vector<double>& localGradient,
                                             Dune::Matrix<double>& localHessian);
 
-    std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_; 
+    std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_;
 
 };
 
@@ -50,10 +59,10 @@ energy(const LocalView& localView,
 
     std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size());
     hyperdual energy;
-    
+
     for (size_t i=0; i<localConfiguration.size(); i++)
-        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);        
-        
+        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
+
     energy = localEnergy_->energy(localView,localHyperDualConfiguration);
 
     pureEnergy = energy.real();
@@ -62,6 +71,37 @@ energy(const LocalView& localView,
 }
 
 
+
+template<class LocalView>
+void LocalHyperDualStiffness<LocalView>::
+assembleGradient(const LocalView& localView,
+                 const std::vector<double>& localConfiguration,
+                 std::vector<double>& localGradient
+)
+{
+    size_t nDoubles = localConfiguration.size();
+
+    localGradient.resize(nDoubles);
+
+    // Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
+    std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
+    for (size_t i=0; i<nDoubles; i++)
+        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
+
+    // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
+    // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
+    for (size_t i=0; i<nDoubles; i++)
+    {
+        localHyperDualConfiguration[i].seteps1(1.0);
+
+        hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
+        localGradient[i] = temp.eps1();
+
+        localHyperDualConfiguration[i].seteps1(0.0);
+    }
+}
+
+
 template<class LocalView>
 void LocalHyperDualStiffness<LocalView>::
 assembleGradientAndHessian(const LocalView& localView,
@@ -70,7 +110,7 @@ assembleGradientAndHessian(const LocalView& localView,
                            Dune::Matrix<double>& localHessian
                           )
 {
-    size_t nDoubles = localConfiguration.size(); 
+    size_t nDoubles = localConfiguration.size();
 
     localGradient.resize(nDoubles);
 
@@ -79,28 +119,28 @@ assembleGradientAndHessian(const LocalView& localView,
     std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
     for (size_t i=0; i<nDoubles; i++)
        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
-    
+
     // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
     // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
-    for (size_t i=0; i<nDoubles; i++) 
+    for (size_t i=0; i<nDoubles; i++)
     {
       localHyperDualConfiguration[i].seteps1(1.0);
       localHyperDualConfiguration[i].seteps2(1.0);
-    
+
       hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
-      localGradient[i] = temp.eps1();                   
+      localGradient[i] = temp.eps1();
       localHessian[i][i] = temp.eps1eps2();
-      
+
       localHyperDualConfiguration[i].seteps2(0.0);
 
-      for (size_t j=i+1; j<nDoubles; j++) 
-      { 
+      for (size_t j=i+1; j<nDoubles; j++)
+      {
         localHyperDualConfiguration[j].seteps2(1.0);
         localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2();
-        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian 
+        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian
         localHyperDualConfiguration[j].seteps2(0.0);
       }
-        
+
       localHyperDualConfiguration[i].seteps1(0.0);
     }
 

test/localhyperdualstiffnesstest.cc

@@ -85,10 +85,10 @@ TestSuite hyperdualEqualsADOL_C()
 
   Timer timer;
   hyperdualAssembler.assembleGradientAndHessian(x,hyperdualGradient,hyperdualHessian);
-  std::cout <<  dim <<"D: hyperdual assembler took " << timer.elapsed() << " sec." << std::endl;
+  std::cout <<  dim <<"D: hyperdual hessian and gradient assembler took " << timer.elapsed() << " sec." << std::endl;
   timer.reset();
   adolcAssembler.assembleGradientAndHessian(x,adolcGradient,adolcHessian);
-  std::cout << dim << "D: ADOL-C assembler took " << timer.elapsed() << " sec." << std::endl;
+  std::cout << dim << "D: ADOL-C hessian and gradient assembler took " << timer.elapsed() << " sec." << std::endl;
 
 
   hyperdualHessian -= adolcHessian;
@@ -97,6 +97,18 @@ TestSuite hyperdualEqualsADOL_C()
   // check the relative error
   t.check(hyperdualHessian.frobenius_norm()/adolcHessian.frobenius_norm() < 1e-12 && hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
 
+
+  timer.reset();
+  hyperdualAssembler.assembleGradient(x,hyperdualGradient);
+  std::cout <<  dim <<"D: hyperdual gradient assembler took " << timer.elapsed() << " sec." << std::endl;
+  timer.reset();
+  adolcAssembler.assembleGradient(x,adolcGradient);
+  std::cout << dim << "D: ADOL-C gradient assembler took " << timer.elapsed() << " sec." << std::endl;
+  hyperdualGradient -= adolcGradient;
+
+  // check the relative error
+  t.check(hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
+
   return t;
 
 }