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

@@ -11,6 +11,11 @@ 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... Addends>
 class SumFunctionalConstrainedLinearization
@@ -29,10 +34,11 @@ public:
   : addends_(std::apply([&](auto&&... fi) {
       return std::make_tuple(Addends(fi, ignore)...);
     }, f.functions())),
-    ignore_(ignore),
-    truncationTolerance_(1e-10)
+    ignore_(ignore)
   {}
 
+  /** \brief Pre-compute derivative information at the configuration x
+   */
   void bind(const Vector& x)
   {
     Hybrid::forEach(addends_, [&](auto&& addend) {
@@ -57,9 +63,10 @@ public:
     // 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++)
       Hybrid::forEach(addends_, [&](auto&& addend) {
         truncationFlags_[i] |= addend.truncated()[i];
@@ -69,9 +76,14 @@ public:
     // 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)
@@ -86,29 +98,41 @@ public:
 #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_;
@@ -122,8 +146,6 @@ private:
   // Mark degrees of freedom that should be truncated
   const BitVector& ignore_;
 
-  double truncationTolerance_;
-
   Vector negativeGradient_;
   Matrix hessian_;
   BitVector truncationFlags_;