Introduce SumFunctionalConstrainedLinearization

Previously, the brittle-fracture program used a taylor-made linearization called ElasticEnergyConstrainedLinearization. However, there wasn't actually anything elasticity-specific in the code, and therefore this MR replaces said class by a generic SumFunctionalConstrainedLinearization.

dune/fracture-phasefields/sumfunctionalconstrainedlinearization.hh

 #ifndef DUNE_FRACTUREPHASEFIELDS_SUMFUNCTIONALCONSTRAINEDLINEARIZATION_HH
 #define DUNE_FRACTUREPHASEFIELDS_SUMFUNCTIONALCONSTRAINEDLINEARIZATION_HH
 
+#include <tuple>
+
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/solvers/common/resize.hh>
+
 namespace Dune::FracturePhasefields {
 
 /**
  * \brief A constrained linearization for a sum of functionals
+ *
+ * \tparam F The SumFunctional to be linearized
+ * \tparam BV Bit vector types for the degrees of freedom to be ignored
+ * \tparam Addends A list of ConstrainedLinearization implementations,
+ *     corresponding to the sum functional F
  */
-template<class F, class BV, class Addend1, class Addend2>
+template<class F, class BV, class... Addends>
 class SumFunctionalConstrainedLinearization
 {
 public:
@@ -18,86 +29,110 @@ public:
   using ConstrainedBitVector = BitVector;
 
   /** \brief Constructor
-   *
-   * \todo Having precisely two addends is hard-wired here!
    */
   SumFunctionalConstrainedLinearization(const F& f, const BitVector& ignore)
-  : ignore_(ignore),
-    addends_(std::make_tuple(Addend1(std::get<0>(f.functions()), ignore),
-                             Addend2(std::get<1>(f.functions()), ignore))),
-    truncationTolerance_(1e-10)
+  : addends_(std::apply([&](auto&&... fi) {
+      return std::make_tuple(Addends(fi, ignore)...);
+    }, f.functions())),
+    ignore_(ignore)
   {}
 
+  /** \brief Pre-compute derivative information at the configuration x
+   */
   void bind(const Vector& x)
   {
-    std::get<0>(addends_).bind(x);
-    std::get<1>(addends_).bind(x);
+    Hybrid::forEach(addends_, [&](auto&& addend) {
+      addend.bind(x);
+    });
 
     // Compute first derivatives of the functionals
-    negativeGradient_ = std::get<0>(addends_).negativeGradient();
-    negativeGradient_ += std::get<1>(addends_).negativeGradient();
+    Solvers::resizeInitializeZero(negativeGradient_, x);
+
+    Hybrid::forEach(addends_, [&](auto&& addend) {
+      negativeGradient_ += addend.negativeGradient();
+    });
 
     // Compute second derivatives of the functionals
     hessian_ = std::get<0>(addends_).hessian();
-    hessian_ += std::get<1>(addends_).hessian();
+
+    Hybrid::forEach(Hybrid::integralRange(std::integral_constant<int, 1>(),Hybrid::size(addends_)), [&](auto i) {
+      hessian_ += std::get<i>(addends_).hessian();
+    });
 
     ////////////////////////////////////////////////////
     // Determine which dofs to truncate
     ////////////////////////////////////////////////////
 
-    // The ignore_ field contains the Dirichlet dofs
+    // Truncate all degrees of freedom explicitly requested in the constructor
     truncationFlags_ = ignore_;
 
+    // Also truncate a degree of freedom if one of the addends wants it truncated
     for (std::size_t i=0; i<truncationFlags_.size(); i++)
-    {
-      truncationFlags_[i] |= std::get<0>(addends_).truncated()[i];
-      truncationFlags_[i] |= std::get<1>(addends_).truncated()[i];
-    }
+      Hybrid::forEach(addends_, [&](auto&& addend) {
+        truncationFlags_[i] |= addend.truncated()[i];
+      });
 
     /////////////////////////////////////////////////////
     // Do the actual truncation
     /////////////////////////////////////////////////////
 
-    // Upper left diagonal block
     for (std::size_t i=0; i<x.size(); ++i)
     {
+      // Truncate the gradient
+      for (std::size_t j=0; j<x[i].size(); ++j)
+        if (truncationFlags_[i][j])
+          negativeGradient_[i][j] = 0;
+
+      // Truncate the Hesse matrix
       auto it = hessian_[i].begin();
       auto end = hessian_[i].end();
       for (; it!=end; ++it)
         for (std::size_t j=0; j<x[i].size(); ++j)
           for (std::size_t k=0; k<x[i].size(); ++k)
             if (truncationFlags_[i][j] || truncationFlags_[it.index()][k])
-	      {
+            {
               (*it)[j][k] = 0;
 #ifndef GAUSSSEIDEL
-	      if (j==k && i == it.index())
-		(*it)[j][k] = 1;
+              if (j==k && i == it.index())
+                (*it)[j][k] = 1;
 #endif
-	      }
+            }
 
-      for (std::size_t j=0; j<x[i].size(); ++j)
-        if (truncationFlags_[i][j])
-          negativeGradient_[i][j] = 0;
     }
   }
 
-  void extendCorrection(ConstrainedVector& cv, Vector& v) const
+  /** \brief Fill the truncated degrees of freedom in a vector with zeros
+   *
+   * \param cv Input vector
+   * \param[out] v Output vector: will be a copy of cv, with the truncated
+   *   degrees of freedom overwritten by zero
+   */
+  void extendCorrection(const ConstrainedVector& cv, Vector& v) const
   {
     for (std::size_t i=0; i<v.size(); ++i)
       for (std::size_t j=0; j<v[i].size(); ++j)
         v[i][j] = (truncationFlags_[i][j]) ? 0 : cv[i][j];
   }
 
+  /** \brief Access to which degrees of freedom have been truncated */
   const BitVector& truncated() const
   {
     return truncationFlags_;
   }
 
+  /** \brief The negative gradient of the sum functional
+   *
+   * Only call this after a call to bind()!
+   */
   const auto& negativeGradient() const
   {
     return negativeGradient_;
   }
 
+  /** \brief The Hesse matrix of the sum functional
+   *
+   * Only call this after a call to bind()!
+   */
   const auto& hessian() const
   {
     return hessian_;
@@ -105,11 +140,11 @@ public:
 
 private:
 
-  const BitVector& ignore_;
-
-  std::tuple<Addend1,Addend2> addends_;
+  // The ConstrainedLinearization objects that are being summed
+  std::tuple<Addends...> addends_;
 
-  double truncationTolerance_;
+  // Mark degrees of freedom that should be truncated
+  const BitVector& ignore_;
 
   Vector negativeGradient_;
   Matrix hessian_;