Use shared pointer for nonlinearity

dune/tectonic/globallaursennonlinearity.hh

@@ -7,6 +7,7 @@
 
 #include <dune/common/fvector.hh>
 
+#include "localnonlinearity.hh"
 #include "globalnonlinearity.hh"
 #include "nicefunction.hh"
 
@@ -33,11 +34,15 @@ class GlobalLaursenNonlinearity : public Dune::GlobalNonlinearity<dim> {
 
     sigma_n [id + mu id] = sigma_n (1 + mu) id
   */
-  virtual Dune::NiceFunction *restriction(int i) const {
+  virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
     double coefficient = nodalIntegrals[i][0];
     coefficient *= normalStress[i];
     coefficient *= 1 + mu[i];
-    return new OuterFunctionType(coefficient);
+
+    shared_ptr<NiceFunction const> const func(
+        new OuterFunctionType(coefficient));
+    return shared_ptr<LocalNonlinearity<dim> const>(
+        new LocalNonlinearity<dim>(func));
   }
 
 private:

dune/tectonic/globalnonlinearity.hh

@@ -6,14 +6,15 @@
 #include <dune/common/fvector.hh>
 
 #include "nicefunction.hh"
+#include "localnonlinearity.hh"
 
 namespace Dune {
 template <int dim> class GlobalNonlinearity {
 public:
   /*
-    Return a restriction of the outer function to the i'th node. If
+    Return a restriction of the outer function to the i'th node.
   */
-  virtual NiceFunction* restriction(int i) const = 0;
+  virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const = 0;
 };
 }
 #endif

dune/tectonic/globalruinanonlinearity.hh

@@ -22,20 +22,25 @@ class GlobalRuinaNonlinearity : public Dune::GlobalNonlinearity<dim> {
         a(a),
         mu(mu),
         eta(eta),
-        normalStress(normalStress) {}
+        normalStress(normalStress),
+        trivialNonlinearity(new LocalNonlinearity<dim>(
+            shared_ptr<NiceFunction const>(new TrivialFunction))) {}
 
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  virtual Dune::NiceFunction *restriction(int i) const {
+  virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
     if (nodalIntegrals[i][0] == 0)
-      return new Dune::TrivialFunction();
-    else
-      return new Dune::RuinaFunction(nodalIntegrals[i][0], a[i], mu[i], eta[i],
-                                     normalStress[i]);
+      return trivialNonlinearity;
+
+    shared_ptr<NiceFunction const> const func(new Dune::RuinaFunction(
+        nodalIntegrals[i][0], a[i], mu[i], eta[i], normalStress[i]));
+    return shared_ptr<LocalNonlinearity<dim>>(new LocalNonlinearity<dim>(func));
   }
 
 private:
+  shared_ptr<LocalNonlinearity<dim> const> const trivialNonlinearity;
+
   std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
   std::vector<double> a;
   std::vector<double> mu;

dune/tectonic/localnonlinearity.hh

@@ -18,7 +18,7 @@ template <int dimension> class LocalNonlinearity {
   typedef FieldVector<double, dimension> VectorType;
   typedef FieldMatrix<double, dimension, dimension> MatrixType;
 
-  LocalNonlinearity(Dune::shared_ptr<NiceFunction> &func) : func_(func) {}
+  LocalNonlinearity(Dune::shared_ptr<NiceFunction const> func) : func_(func) {}
 
   double operator()(VectorType const x) const {
     double ret;
@@ -63,7 +63,7 @@ template <int dimension> class LocalNonlinearity {
   }
 
 private:
-  Dune::shared_ptr<NiceFunction> &func_;
+  Dune::shared_ptr<NiceFunction const> const func_;
 };
 }
 #endif

dune/tectonic/myblockproblem.hh

@@ -9,7 +9,6 @@
 #include <dune/tnnmg/problem-classes/onedconvexfunction.hh>
 
 #include "localnonlinearity.hh"
-#include "nicefunction.hh"
 #include "samplefunctional.hh"
 
 /** \brief Base class for problems where each block can be solved with a
@@ -118,8 +117,7 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject {
       }
       assert(localA != NULL);
 
-      Dune::shared_ptr<Dune::NiceFunction> f(problem.phi.restriction(m));
-      Dune::LocalNonlinearity<block_size> const phi(f);
+      auto phi = problem.phi.restriction(m);
       Dune::SampleFunctional<block_size> localJ(*localA, localb, phi,
                                                 ignore_component);
 

dune/tectonic/samplefunctional.hh

@@ -24,15 +24,15 @@ template <int dim> class SampleFunctional {
   typedef LocalNonlinearity<dim> NonlinearityType;
 
   SampleFunctional(SmallMatrix const &A, SmallVector const &b,
-                   NonlinearityType const &phi, int ignore = dim)
+                   shared_ptr<NonlinearityType const> phi, int ignore = dim)
       : A(A), b(b), phi(phi), ignore(ignore) {}
 
   double operator()(SmallVector const &v) const {
     SmallVector y;
-    A.mv(v, y);            //      Av
-    y /= 2;                //  1/2 Av
-    y -= b;                //  1/2 Av - b
-    return y * v + phi(v); // <1/2 Av - b,v> + H(|v|)
+    A.mv(v, y);             //      Av
+    y /= 2;                 //  1/2 Av
+    y -= b;                 //  1/2 Av - b
+    return y * v + *phi(v); // <1/2 Av - b,v> + H(|v|)
   }
 
   bool descentDirection(SmallVector const x, SmallVector &ret) const {
@@ -44,7 +44,7 @@ template <int dim> class SampleFunctional {
       smoothGradient(x, d);
 
       Interval<double> D;
-      phi.directionalSubDiff(x, d, D);
+      phi->directionalSubDiff(x, d, D);
       double const nonlinearDecline = D[1];
       double const smoothDecline = -(d * d);
       double const combinedDecline = smoothDecline + nonlinearDecline;
@@ -94,7 +94,7 @@ template <int dim> class SampleFunctional {
 
   SmallMatrix const &A;
   SmallVector const &b;
-  NonlinearityType const &phi;
+  shared_ptr<NonlinearityType const> const phi;
   int const ignore;
 
 private:
@@ -107,7 +107,7 @@ template <int dim> class SampleFunctional {
   }
 
   void upperGradient(SmallVector const &x, SmallVector &y) const {
-    phi.upperGradient(x, y);
+    phi->upperGradient(x, y);
     SmallVector z;
     smoothGradient(x, z);
     y += z;
@@ -116,7 +116,7 @@ template <int dim> class SampleFunctional {
   }
 
   void lowerGradient(SmallVector const &x, SmallVector &y) const {
-    phi.lowerGradient(x, y);
+    phi->lowerGradient(x, y);
     SmallVector z;
     smoothGradient(x, z);
     y += z;
@@ -156,7 +156,7 @@ void minimise(Functional const J, typename Functional::SmallVector &x,
       return;
 
     typedef typename Functional::NonlinearityType LocalNonlinearityType;
-    LocalNonlinearityType phi = J.phi;
+    LocalNonlinearityType phi = *J.phi; // copy, see below
     if (linesearchp) {
       // {{{ Construct a restriction of J to the line x + t * descDir