clean up

dune/tnnmg/functionals/semilinearfunctional.hh

@@ -237,14 +237,6 @@ public:
     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)

dune/tnnmg/functionals/semilinearfunctionallinearization.hh

@@ -2,6 +2,7 @@
 #include <dune/common/hybridutilities.hh>
 #include <dune/common/typetraits.hh>
 #include <dune/matrix-vector/algorithm.hh>
+#include <dune/tnnmg/functionals/semilinearfunctional.hh>
 #include <type_traits>
 
 namespace Dune::TNNMG {
@@ -12,24 +13,6 @@ 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
  */
@@ -38,14 +21,13 @@ 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));
-      });
+  if constexpr (Dune::IsNumber<std::decay_t<V>>()) {
+    u *= v;
+  } else {
+    H::forEach(H::integralRange(H::size(v)), [&](auto i) {
+      multiplyComponentwise(H::elementAt(u, i), H::elementAt(v, i));
     });
+  }
 }
 
 template<class M, class V>
@@ -53,15 +35,15 @@ 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));
-      });
+  if constexpr (Dune::IsNumber<std::decay_t<V>>()) {
+    m = v;
+  } else {
+    H::forEach(H::integralRange(H::size(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)
@@ -72,7 +54,8 @@ determineTruncation(const NV& x, NBV&& truncationFlags, const F& function)
   } else {
     assert(x.size() == truncationFlags.size());
     H::forEach(H::integralRange(H::size(x)), [&](auto&& i) {
-      determineTruncation(x[i], truncationFlags[i], function); // todo element_at
+      determineTruncation(
+        H::elementAt(x, i), H::elementAt(truncationFlags, i), function);
     });
   }
 }
@@ -86,8 +69,9 @@ truncateVector(NV& x, const NBV& truncationFlags)
     if (truncationFlags)
       x = 0;
   } else {
-    H::forEach(H::integralRange(H::size(x)),
-               [&](auto&& i) { truncateVector(x[i], truncationFlags[i]); });
+    H::forEach(H::integralRange(H::size(x)), [&](auto&& i) {
+      truncateVector(H::elementAt(x, i), H::elementAt(truncationFlags, i));
+    });
   }
 }
 
@@ -104,9 +88,11 @@ truncateMatrix(NM& A,
       A = 0;
   } else {
     H::forEach(H::integralRange(H::size(rowTruncationFlags)), [&](auto&& i) {
-      auto&& Ai = A[i];
+      auto&& Ai = H::elementAt(A, i);
       sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
-        truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+        truncateMatrix(Aij,
+                       H::elementAt(rowTruncationFlags, i),
+                       H::elementAt(colTruncationFlags, i));
       });
     });
   }
@@ -158,14 +144,14 @@ public:
 
     // 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);
+    Impl::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);
+    Impl::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);