Use shared pointers to avoid copying

dune/tectonic/globallaursennonlinearity.hh

@@ -17,7 +17,8 @@ template <int dim, class OuterFunctionType>
 class GlobalLaursenNonlinearity : public GlobalNonlinearity<dim> {
 public:
   GlobalLaursenNonlinearity(
-      std::vector<double> const &mu, std::vector<double> const &normalStress,
+      shared_ptr<std::vector<double> const> mu,
+      shared_ptr<std::vector<double> const> normalStress,
       std::vector<FieldVector<double, 1>> const &nodalIntegrals)
       : mu(mu), normalStress(normalStress), nodalIntegrals(nodalIntegrals) {}
 
@@ -37,8 +38,8 @@ class GlobalLaursenNonlinearity : public GlobalNonlinearity<dim> {
   */
   virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
     double coefficient = nodalIntegrals[i][0];
-    coefficient *= normalStress[i];
-    coefficient *= 1 + mu[i];
+    coefficient *= (*normalStress)[i];
+    coefficient *= 1 + (*mu)[i];
 
     shared_ptr<NiceFunction const> const func(
         new OuterFunctionType(coefficient));
@@ -49,8 +50,8 @@ class GlobalLaursenNonlinearity : public GlobalNonlinearity<dim> {
 private:
   // TODO: If we're clever, we only store one vector with the precomputed
   // results
-  std::vector<double> mu;
-  std::vector<double> normalStress;
+  shared_ptr<std::vector<double> const> mu;
+  shared_ptr<std::vector<double> const> normalStress;
   std::vector<FieldVector<double, 1>> nodalIntegrals;
 };
 }

dune/tectonic/globalruinanonlinearity.hh

@@ -17,8 +17,10 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<dim> {
 public:
   GlobalRuinaNonlinearity(
       std::vector<FieldVector<double, 1>> const &nodalIntegrals,
-      std::vector<double> const &a, std::vector<double> const &mu,
-      std::vector<double> const &eta, std::vector<double> const &normalStress)
+      shared_ptr<std::vector<double> const> a,
+      shared_ptr<std::vector<double> const> mu,
+      shared_ptr<std::vector<double> const> eta,
+      shared_ptr<std::vector<double> const> normalStress)
       : nodalIntegrals(nodalIntegrals),
         a(a),
         mu(mu),
@@ -34,8 +36,9 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<dim> {
     if (nodalIntegrals[i][0] == 0)
       return trivialNonlinearity;
 
-    shared_ptr<NiceFunction const> const func(new RuinaFunction(
-        nodalIntegrals[i][0], a[i], mu[i], eta[i], normalStress[i]));
+    shared_ptr<NiceFunction const> const func(
+        new RuinaFunction(nodalIntegrals[i][0], (*a)[i], (*mu)[i], (*eta)[i],
+                          (*normalStress)[i]));
     return shared_ptr<LocalNonlinearity<dim>>(new LocalNonlinearity<dim>(func));
   }
 
@@ -43,10 +46,10 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<dim> {
   shared_ptr<LocalNonlinearity<dim> const> const trivialNonlinearity;
 
   std::vector<FieldVector<double, 1>> nodalIntegrals;
-  std::vector<double> a;
-  std::vector<double> mu;
-  std::vector<double> eta;
-  std::vector<double> normalStress;
+  shared_ptr<std::vector<double> const> a;
+  shared_ptr<std::vector<double> const> mu;
+  shared_ptr<std::vector<double> const> eta;
+  shared_ptr<std::vector<double> const> normalStress;
 };
 }
 #endif

src/one-body-sample.cc

@@ -140,26 +140,27 @@ void assemble_nonlinearity(
     Dune::shared_ptr<Dune::GlobalNonlinearity<dim> const> &myGlobalNonlinearity,
     std::vector<Dune::FieldVector<double, 1>> &nodalIntegrals) {
   // {{{ Assemble terms for the nonlinearity
-  std::vector<double> mu;
-  mu.resize(size);
-  std::fill(mu.begin(), mu.end(), parset.get<double>("boundary.friction.mu"));
-
-  std::vector<double> normalStress;
-  normalStress.resize(size);
-  std::fill(normalStress.begin(), normalStress.end(),
+  auto mu = Dune::shared_ptr<std::vector<double>>(new std::vector<double>());
+  mu->resize(size);
+  std::fill(mu->begin(), mu->end(), parset.get<double>("boundary.friction.mu"));
+
+  auto normalStress =
+      Dune::shared_ptr<std::vector<double>>(new std::vector<double>());
+  normalStress->resize(size);
+  std::fill(normalStress->begin(), normalStress->end(),
             parset.get<double>("boundary.friction.normalstress"));
 
   std::string const friction_model =
       parset.get<std::string>("boundary.friction.model");
   if (friction_model == std::string("Ruina")) {
-    std::vector<double> a;
-    a.resize(size);
-    std::fill(a.begin(), a.end(),
+    auto a = Dune::shared_ptr<std::vector<double>>(new std::vector<double>());
+    a->resize(size);
+    std::fill(a->begin(), a->end(),
               parset.get<double>("boundary.friction.ruina.a"));
 
-    std::vector<double> eta;
-    eta.resize(size);
-    std::fill(eta.begin(), eta.end(),
+    auto eta = Dune::shared_ptr<std::vector<double>>(new std::vector<double>());
+    eta->resize(size);
+    std::fill(eta->begin(), eta->end(),
               parset.get<double>("boundary.friction.eta"));
 
     auto const tmp = new Dune::GlobalRuinaNonlinearity<dim>(