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
  */
-template<class F, class BV, class Addend1, class Addend2>
+template<class F, class BV, class... Addends>
 class SumFunctionalConstrainedLinearization
 {
 public:
@@ -18,28 +24,34 @@ 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))),
+  : addends_(std::apply([&](auto&&... fi) {
+      return std::make_tuple(Addends(fi, ignore)...);
+    }, f.functions())),
+    ignore_(ignore),
     truncationTolerance_(1e-10)
   {}
 
   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
@@ -49,10 +61,9 @@ public:
     truncationFlags_ = ignore_;
 
     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
@@ -67,13 +78,13 @@ public:
         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])
@@ -105,9 +116,11 @@ public:
 
 private:
 
-  const BitVector& ignore_;
+  // The ConstrainedLinearization objects that are being summed
+  std::tuple<Addends...> addends_;
 
-  std::tuple<Addend1,Addend2> addends_;
+  // Mark degrees of freedom that should be truncated
+  const BitVector& ignore_;
 
   double truncationTolerance_;