Start working on a semilinear functional

dune/tnnmg/functionals/semilinearfunctional.hh

+#pragma once
+#include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/solvers/common/interval.hh>
+
+#include <memory>
+
+namespace Dune::TNNMG {
+
+namespace Impl {
+
+/** Compute Euclidian scalar product */
+template<typename V>
+auto
+dot(const V& v, const V& w)
+{
+  // TODO: Make this more general for vector types that do not implement
+  // operator*.
+  // Maybe, dune-matrix-vector also has this.
+  return v * w;
+}
+
+template<typename V>
+auto
+zeroCopy(const V& v)
+{
+  auto ret = autoCopy(v);
+  ret = 0.;
+  return ret;
+}
+
+/** Apply a function f componentwise to a vector v and store it in u, i.e.
+ *  u_i = f(v_i)
+ */
+template<typename F, typename V>
+void
+applyComponentwise(F&& f, V& u, const V& v)
+{
+  namespace H = Dune::Hybrid;
+  H::ifElse(
+    Dune::IsNumber<std::decay_t<V>>(),
+    [&](auto&& id) { u = f(id(v)); },
+    [&](auto id) {
+      H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+        applyComponentwise(f, H::elementAt(u, i), H::elementAt(v, i));
+      });
+    });
+}
+
+template<typename F, typename V, typename S>
+void
+directionalDomain(const V& origin,
+                  const V& direction,
+                  const F& phi,
+                  S& lower,
+                  S& upper)
+{
+  namespace H = Dune::Hybrid;
+  if constexpr (Dune::IsNumber<std::decay_t<V>>()) {
+    if (direction == 0) {
+      return;
+    }
+    // TODO check this
+    else if (direction > 0) {
+      lower = std::max(lower, (phi.domain()[0] - origin) / direction);
+      upper = std::min(upper, (phi.domain()[1] - origin) / direction);
+    } else {
+      lower = std::max(lower, (phi.domain()[1] - origin) / direction);
+      upper = std::min(upper, (phi.domain()[0] - origin) / direction);
+    }
+  } else {
+    H::forEach(H::integralRange(H::size(origin)), [&](auto i) {
+      directionalDomain(
+        H::elementAt(origin, i), H::elementAt(direction, i), phi, lower, upper);
+    });
+  }
+}
+
+template<typename F, typename V, typename S>
+void
+addSubdiffs(const V& origin,
+            const V& direction,
+            const V& weights,
+            const F& phi,
+            double t,
+            S& interval)
+{
+  // TODO: Check this
+  namespace H = Dune::Hybrid;
+  if constexpr (Dune::IsNumber<std::decay_t<V>>()) {
+    auto local_subdiff = phi.subDifferential(origin + t * direction);
+    local_subdiff *= weights * direction; // TODO stimmt das?
+    interval += local_subdiff;
+  } else {
+    H::forEach(H::integralRange(H::size(origin)), [&](auto i) {
+      addSubdiffs(H::elementAt(origin, i),
+                  H::elementAt(direction, i),
+                  H::elementAt(weights, i),
+                  phi,
+                  t,
+                  interval);
+    });
+  }
+}
+
+auto
+reals()
+{
+  return Dune::Solvers::Interval<double>(
+    { std::numeric_limits<double>::lowest(),
+      std::numeric_limits<double>::max() });
+}
+
+template<typename V, typename F>
+auto
+computeFunctionalValue(const V& weights, V pos, const F& phi)
+{
+  applyComponentwise(phi, pos, pos);
+  return dot(weights, pos);
+}
+}
+
+template<typename V, typename ScalarFunction, typename R = double>
+class SemilinearfunctionalDirectionalRestriction
+{
+public:
+  using Vector = V;
+  using Range = R;
+  template<typename FF>
+  SemilinearfunctionalDirectionalRestriction(const Vector& weights,
+                                             const Vector& origin,
+                                             const Vector& direction,
+                                             FF&& phi)
+    : weights_(weights)
+    , origin_(origin)
+    , direction_(direction)
+    , phi_(std::forward<FF>(phi))
+  {}
+
+  Range operator()(const Range& t) const
+  {
+    // This is of course a PITA having to copy vectors
+    // every single time. Well, what can you do...
+
+    // TODO: origin_ part is missing!
+    DUNE_THROW(Dune::Exception, "Dont use this yet");
+
+    auto val = direction_;
+    val *= t;
+    return Impl::computeFunctionalValue(weights_, std::move(val), phi_);
+  }
+
+  auto domain() const
+  {
+    auto dom = Impl::reals();
+    Impl::directionalDomain(origin_, direction_, phi_, dom[0], dom[1]);
+
+    return dom;
+  }
+
+  auto subDifferential(double t) const
+  {
+    // TODO:
+
+    /*
+    Okay, we need to compute the subdifferential of
+    \sum_i w_i * phi(x_i + t*c_i)
+    in t.
+    This should be
+    \sum_i w_i * 6phi(x_i + t*c_i)*c_i, right?
+    */
+    auto sd = Dune::Solvers::Interval<double>({ 0., 0. });
+    Impl::addSubdiffs(origin_, direction_, weights_, phi_, t, sd);
+
+    return sd;
+  }
+
+  // TODO: Provide some read-only access to data such that a clever
+  // person could write a sophisticated optimization routine for her/his
+  // particular problem.
+
+private:
+  Vector weights_;
+  Vector origin_;
+  Vector direction_;
+  ScalarFunction phi_;
+};
+
+template<typename V, typename ScalarFunction, typename R = double>
+class Semilinearfunctional
+{
+public:
+  using Vector = V;
+  using Range = R;
+  template<typename FF>
+  Semilinearfunctional(const Vector& weights, FF&& phi)
+    : weights_(weights)
+    , origin_(Impl::zeroCopy(weights_))
+    , phi_(std::forward<FF>(phi))
+  {}
+
+  template<typename FF>
+  Semilinearfunctional(const Vector& weights, const Vector& origin, FF&& phi)
+    : weights_(weights)
+    , origin_(origin)
+    , phi_(std::forward<FF>(phi))
+  {}
+
+  void updateOrigin() {}
+
+  // TODO what should this actually do?
+  void updateOrigin(std::size_t) {}
+
+  auto domain() const { return phi_.domain(); }
+
+  auto subDifferential(double pos) const
+  {
+    auto shiftedPos = origin_ + pos;
+    return phi_.subDifferential(shiftedPos);
+  }
+
+  Range operator()(Vector x) const
+  {
+    // TODO copying x is expensive. We should supply both vectors to the compute
+    // function.
+    x += origin_;
+    return Impl::computeFunctionalValue(weights_, std::move(x), phi_);
+  }
+
+  friend auto shift(const Semilinearfunctional& f, const Vector& origin)
+  {
+    auto new_origin = f.origin_;
+    new_origin += origin;
+    return Semilinearfunctional(f.weights_, new_origin, f.phi_);
+  }
+
+  friend auto derivative(const Semilinearfunctional& f)
+  {
+    // TODO: is this too much to ask for?
+    // auto phi_prime = derivative(f.phi_);
+    DUNE_THROW(Dune::NotImplemented, "Derivative not yet implemented");
+    return 42;
+  }
+
+  friend auto directionalRestriction(const Semilinearfunctional& f,
+                                     const Vector& origin,
+                                     const Vector& direction)
+  {
+    // TODO What happens if f already has a nontrivial origin? And if so,
+    // what SHOULD happen?
+    return SemilinearfunctionalDirectionalRestriction<decltype(f.weights_),
+                                                      decltype(f.phi_)>(
+      f.weights_, origin, direction, f.phi_);
+  }
+
+  template<typename Index>
+  friend auto coordinateRestriction(const Semilinearfunctional& f,
+                                    const Index& i)
+  {
+    using BT = std::decay_t<decltype(f.weights_[0])>;
+    // TODO: This will of course leave out the constant contributions in the
+    // other components. But since the minimization is independent from these,
+    // why bother? The coupling between components should be handled in the
+    // quadratic part.
+    return Semilinearfunctional<BT, decltype(f.phi_)>(
+      f.weights_[i], f.origin_[i], f.phi_);
+  }
+
+  const Vector& weights() const { return weights_; }
+
+  /** TODO: Doc me! */
+  const ScalarFunction& phi() const { return phi_; }
+
+private:
+  Vector weights_;
+  Vector origin_;
+  ScalarFunction phi_;
+};
+}
\ No newline at end of file

dune/tnnmg/functionals/semilinearfunctionallinearization.hh

+#pragma once
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/matrix-vector/algorithm.hh>
+#include <type_traits>
+
+namespace Dune::TNNMG {
+
+namespace Helper {
+/* In this namespace, we have a bunch of helper functions that are
+actually could be applied in other contexts, too. Some are even part of open MR
+that gathered dust by now.
+*/
+
+/** \brief For two vectors of equal size (and blocking), apply a scalar function
+ * componentswise such that u_i = f(v_i).
+ */
+template<typename F, typename V>
+void
+applyComponentwise(F&& f, V& u, const V& v)
+{
+  namespace H = Dune::Hybrid;
+  H::ifElse(
+    Dune::IsNumber<std::decay_t<V>>(),
+    [&](auto&& id) { u = f(id(v)); },
+    [&](auto id) {
+      H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+        applyComponentwise(f, H::elementAt(u, i), H::elementAt(v, i));
+      });
+    });
+}
+
+/** Multiply the i-th component of u by the i-th component of v, i.e.
+ *  u_i <- u_i*v_i
+ */
+template<typename V>
+void
+multiplyComponentwise(V& u, const V& v)
+{
+  namespace H = Dune::Hybrid;
+  H::ifElse(
+    Dune::IsNumber<std::decay_t<V>>(),
+    [&](auto&& id) { u *= id(v); },
+    [&](auto id) {
+      H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+        multiplyComponentwise(H::elementAt(u, i), H::elementAt(v, i));
+      });
+    });
+}
+
+template<class M, class V>
+void
+writeDiagonal(M& m, V&& v)
+{
+  namespace H = Dune::Hybrid;
+  H::ifElse(
+    Dune::IsNumber<std::decay_t<V>>(),
+    [&](auto&& id) { m = id(v); },
+    [&](auto&& id) {
+      H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+        writeDiagonal(H::elementAt(H::elementAt(m, i), i), H::elementAt(v, i));
+      });
+    });
+}
+template<class NV, class NBV, class F>
+void
+determineTruncation(const NV& x, NBV&& truncationFlags, const F& function)
+{
+  namespace H = Dune::Hybrid;
+  if constexpr (IsNumber<NV>()) {
+    truncationFlags = function.should_truncate(x);
+  } else {
+    assert(x.size() == truncationFlags.size());
+    H::forEach(H::integralRange(H::size(x)), [&](auto&& i) {
+      determineTruncation(x[i], truncationFlags[i], function); // todo element_at
+    });
+  }
+}
+
+template<class NV, class NBV>
+void
+truncateVector(NV& x, const NBV& truncationFlags)
+{
+  namespace H = Dune::Hybrid;
+  if constexpr (IsNumber<NV>()) {
+    if (truncationFlags)
+      x = 0;
+  } else {
+    H::forEach(H::integralRange(H::size(x)),
+               [&](auto&& i) { truncateVector(x[i], truncationFlags[i]); });
+  }
+}
+
+template<class NM, class RBV, class CBV>
+void
+truncateMatrix(NM& A,
+               const RBV& rowTruncationFlags,
+               const CBV& colTruncationFlags)
+{
+  namespace H = Dune::Hybrid;
+  using namespace Dune::MatrixVector;
+  if constexpr (IsNumber<NM>()) {
+    if (rowTruncationFlags or colTruncationFlags)
+      A = 0;
+  } else {
+    H::forEach(H::integralRange(H::size(rowTruncationFlags)), [&](auto&& i) {
+      auto&& Ai = A[i];
+      sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+        truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+      });
+    });
+  }
+}
+
+}
+
+/** A constrained linearization for the SemilinearFunctional */
+template<typename F,
+         typename ScalarFunction,
+         typename M,
+         typename V,
+         typename BV>
+class SemilinearFunctionalConstrainedLinearization
+{
+
+public:
+  using BitVector = BV;
+  using Matrix = M;
+  using Vector = V;
+  using ConstrainedMatrix = Matrix;
+  using ConstrainedVector = Vector;
+  using ConstrainedBitVector = BitVector;
+
+  SemilinearFunctionalConstrainedLinearization(const F& functional,
+                                               const BitVector& ignore)
+    : f_(functional)
+    , ignore_(ignore)
+  {}
+
+  /** \brief Set a matrix that serves as a template for the Hessian.
+   *
+   * This needs to be a matrix of correct size and sparsity pattern (including
+   * at least the diagonal entries).
+   * TODO: I don't like neither the name nor the fact of having a setter
+   * function at all. This should rather be done in the constructor.
+   */
+  template<typename MM>
+  void setMatrixDummy(const MM& dummy)
+  {
+    // mat_ = std::forward<MM>(dummy);
+    mat_ = dummy;
+  }
+
+  void bind(const Vector& x)
+  {
+    // Prepare the negative gradient
+    ngrad_ = x;
+
+    // Assuming that the supplied function does not explode if it is evaluated
+    // in a point where it's actually not defined.
+    Helper::applyComponentwise(derivative(f_.phi()), ngrad_, x);
+    Helper::multiplyComponentwise(ngrad_, f_.weights());
+    ngrad_ *= -1;
+
+    mat_ = 0;
+    // compute diagonal entries
+    auto diag = x;
+    Helper::applyComponentwise(derivative(derivative(f_.phi())), diag, x);
+    Helper::multiplyComponentwise(diag, f_.weights());
+    // now, write them into the diagonal of the Hessian
+    Helper::writeDiagonal(mat_, diag);
+
+    truncationFlags_ = ignore_;
+    Helper::determineTruncation(x, truncationFlags_, f_.phi());
+    Helper::truncateMatrix(mat_, truncationFlags_, truncationFlags_);
+    Helper::truncateVector(ngrad_, truncationFlags_);
+  }
+
+  void extendCorrection(ConstrainedVector& cv, Vector& v) const
+  {
+    v = cv;
+    Helper::truncateVector(v, truncationFlags_);
+  }
+
+  const BitVector& truncated() const { return truncationFlags_; }
+
+  const auto& negativeGradient() const { return ngrad_; }
+
+  const auto& hessian() const { return mat_; }
+
+private:
+  const F& f_;
+  Matrix mat_;
+  Vector ngrad_;
+  const BitVector& ignore_;
+  BitVector truncationFlags_;
+};
+}
\ No newline at end of file

dune/tnnmg/localsolvers/functionproperties.hh

+#pragma once
+#include <algorithm>
+#include <dune/common/concept.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/solvers/common/interval.hh>
+#include <tuple>
+#include <type_traits>
+
+namespace Dune {
+
+namespace Concept {
+
+struct HasDomain
+{
+  template<typename C>
+  auto require(C&& c) -> decltype(c.domain());
+};
+
+struct HasSubdifferential
+{
+  template<typename C>
+  auto require(C&& c) -> decltype(c.subDifferential(0.));
+};
+
+}
+
+namespace TNNMG {
+
+using namespace Dune::Concept;
+
+template<typename Functional,
+         typename = std::enable_if_t<models<HasDomain, Functional>()>>
+auto
+domain(const Functional& f)
+{
+  return f.domain();
+}
+
+/** For a tuple of intervals, return their intersection.
+ *
+ * If the intersection is empty, the function will throw!
+ */
+template<typename... I>
+auto
+intervalIntersection(const std::tuple<I...>& intervals)
+{
+  // We want the greatest interval such that all functionals are well-defined.
+  auto left_max = std::numeric_limits<double>::lowest();
+  auto right_min = std::numeric_limits<double>::max();
+
+  auto shrink_intervals = [&](const auto&... interval) {
+    ((left_max = std::max(left_max, interval[0])), ...);
+    ((right_min = std::min(right_min, interval[1])), ...);
+  };
+  std::apply(shrink_intervals, intervals);
+
+  // In case there is no such interval (not even a point), we throw an
+  // exception.
+  // TODO: Is this the most reasonable thing to do here?
+  if (right_min < left_max)
+    DUNE_THROW(Dune::Exception,
+               "The intersection over all intervals is empty.");
+  return Dune::Solvers::Interval<double>({ left_max, right_min });
+}
+
+template<typename... F>
+auto
+domain(const std::tuple<F...>& f)
+{
+  auto domains = std::apply(
+    [](const auto&... x) { return std::make_tuple(TNNMG::domain(x)...); }, f);
+
+  return intervalIntersection(domains);
+}
+
+template<typename Functional,
+         typename = std::enable_if_t<models<HasSubdifferential, Functional>()>>
+auto
+subDifferential(Functional&& f, double pos)
+{
+  return f.subDifferential(pos);
+}
+
+template<typename... F>
+auto
+subDifferential(const std::tuple<F...>& f, double pos)
+{
+  // Compute each subdifferential and sum them up.
+  return std::apply(
+    [pos](const auto&... x) { return (TNNMG::subDifferential(x, pos) += ...); },
+    f);
+}
+
+}
+}
\ No newline at end of file

dune/tnnmg/projections/semilineardefectprojection.hh

+#pragma once
+
+#include <dune/tnnmg/functionals/semilinearfunctional.hh> // constains the helper functions
+
+namespace Dune {
+namespace TNNMG {
+
+namespace Helper {
+template<typename F, typename V>
+void
+projectInDomain(F&& f, const V& x, V& c)
+{
+
+  // TODO:
+  // This operators only on the semilinear functional.
+  // For the sumfunctional, we need to ignore the other functionals.
+  //
+  // Plus, what should happen if more than one function restrict the domain?
+  // Likely, one should implement the machinery for the sumfunctional! (Maybe
+  // check the fracture phasefield module, Carsten might have done that before).
+  namespace H = Dune::Hybrid;
+  if constexpr (Dune::IsNumber<std::decay_t<V>>()) {
+    // auto point = f.phi().domain().projectIn(x + c);
+    auto point = std::get<1>(f.functions()).phi().domain().projectIn(x + c);
+    c = point - x;
+  } else {
+    H::forEach(H::integralRange(H::size(x)), [&](auto i) {
+      projectInDomain(f, H::elementAt(x, i), H::elementAt(c, i));
+    });
+  }
+  // namespace H = Dune::Hybrid;
+  // H::ifElse(
+  //   Dune::IsNumber<std::decay_t<V>>(),
+  //   [&](auto&& id) { u = f(id(v)); },
+  //   [&](auto id) {
+  //     H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+  //       applyComponentwise(f, H::elementAt(u, i), H::elementAt(v, i));
+  //     });
+  //   });
+}
+}
+
+struct SemilinearFunctionalDefectProjection
+{
+  template<typename F, typename V>
+  void operator()(const F& functional, const V& x, V& correction)
+  {
+    Helper::projectInDomain(functional, x, correction);
+  }
+};
+
+} // end namespace TNNMG
+} // end namespace Dune
\ No newline at end of file