Do not use shared pointers needlessly

dune/tectonic/globalruinanonlinearity.hh

@@ -18,26 +18,25 @@ template <class MatrixType, class VectorType>
 class GlobalRuinaNonlinearity
     : public GlobalNonlinearity<MatrixType, VectorType> {
 private:
-  typedef shared_ptr<BlockVector<FieldVector<double, 1>> const> dataptr;
+  typedef BlockVector<FieldVector<double, 1>> const &dataref;
 
 public:
   using GlobalNonlinearity<MatrixType, VectorType>::dim;
 
-  GlobalRuinaNonlinearity(dataptr nodalIntegrals, dataptr a, dataptr mu,
-                          dataptr eta, dataptr normalStress, dataptr b,
-                          dataptr state, dataptr L, double h)
-      // TODO: can we get the size from another place?
-      : restrictions(nodalIntegrals->size()) {
+  GlobalRuinaNonlinearity(dataref nodalIntegrals, dataref a, dataref mu,
+                          dataref eta, dataref normalStress, dataref b,
+                          dataref state, dataref L, double h)
+      : restrictions(nodalIntegrals.size()) {
     auto trivialNonlinearity = make_shared<LocalNonlinearity<dim> const>(
         make_shared<TrivialFunction const>());
     for (size_t i = 0; i < restrictions.size(); ++i) {
       restrictions[i] =
-          (*nodalIntegrals)[i] == 0
+          nodalIntegrals[i] == 0
               ? trivialNonlinearity
               : make_shared<LocalNonlinearity<dim> const>(
                     make_shared<RuinaFunction const>(
-                        (*nodalIntegrals)[i], (*a)[i], (*mu)[i], (*eta)[i],
-                        (*normalStress)[i], (*b)[i], (*state)[i], (*L)[i], h));
+                        nodalIntegrals[i], a[i], mu[i], eta[i], normalStress[i],
+                        b[i], state[i], L[i], h));
     }
   }
 

src/assemblers.cc

@@ -56,33 +56,31 @@ template <class MatrixType, class VectorType>
 Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
 assemble_nonlinearity(
     Dune::ParameterTree const &parset,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-        nodalIntegrals,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
-    double h) {
-  auto const size = (*nodalIntegrals).size();
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h) {
+  auto const size = nodalIntegrals.size();
 
   typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
   // {{{ Assemble terms for the nonlinearity
-  auto mu = Dune::make_shared<SingletonVectorType>(size);
-  *mu = parset.get<double>("mu");
+  SingletonVectorType mu(size);
+  mu = parset.get<double>("mu");
 
-  auto normalStress = Dune::make_shared<SingletonVectorType>(size);
-  *normalStress = parset.get<double>("normalstress");
+  SingletonVectorType normalStress(size);
+  normalStress = parset.get<double>("normalstress");
 
   switch (parset.get<Config::model>("model")) {
     case Config::Exponential: {
-      auto a = Dune::make_shared<SingletonVectorType>(size);
-      *a = parset.get<double>("ruina.a");
+      SingletonVectorType a(size);
+      a = parset.get<double>("ruina.a");
 
-      auto eta = Dune::make_shared<SingletonVectorType>(size);
-      *eta = parset.get<double>("eta");
+      SingletonVectorType eta(size);
+      eta = parset.get<double>("eta");
 
-      auto b = Dune::make_shared<SingletonVectorType>(size);
-      *b = parset.get<double>("ruina.b");
+      SingletonVectorType b(size);
+      b = parset.get<double>("ruina.b");
 
-      auto L = Dune::make_shared<SingletonVectorType>(size);
-      *L = parset.get<double>("ruina.L");
+      SingletonVectorType L(size);
+      L = parset.get<double>("ruina.L");
 
       return Dune::make_shared<
           Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(

src/assemblers.hh

@@ -26,9 +26,7 @@ template <class MatrixType, class VectorType>
 Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
 assemble_nonlinearity(
     Dune::ParameterTree const &parset,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-        nodalIntegrals,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
-    double h);
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h);
 
 #endif

src/assemblers_tmpl.cc

@@ -33,7 +33,5 @@ template Dune::shared_ptr<
     Dune::GlobalNonlinearity<MatrixType, VectorType> const>
 assemble_nonlinearity<MatrixType, VectorType>(
     Dune::ParameterTree const &parset,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-        nodalIntegrals,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
-    double h);
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
+    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state, double h);

src/one-body-sample.cc

@@ -198,8 +198,7 @@ int main(int argc, char *argv[]) {
     VectorType u_diff(finestSize);
     u_diff = 0.0; // Has to be zero!
 
-    auto alpha = Dune::make_shared<SingletonVectorType>(finestSize);
-    *alpha = alpha_old;
+    SingletonVectorType alpha(alpha_old);
 
     SingletonVectorType vonMisesStress;
     VectorType rhs;
@@ -268,7 +267,7 @@ int main(int argc, char *argv[]) {
              ++state_fpi) {
           auto myGlobalNonlinearity =
               assemble_nonlinearity<MatrixType, VectorType>(
-                  parset.sub("boundary.friction"), nodalIntegrals, alpha, h);
+                  parset.sub("boundary.friction"), *nodalIntegrals, alpha, h);
           MyConvexProblemType const myConvexProblem(stiffnessMatrix,
                                                     *myGlobalNonlinearity, rhs);
 
@@ -295,11 +294,10 @@ int main(int argc, char *argv[]) {
               switch (parset.get<Config::state_model>(
                   "boundary.friction.state.model")) {
                 case Config::Dieterich:
-                  (*alpha)[i] =
-                      state_update_dieterich(h, unorm / L, alpha_old[i]);
+                  alpha[i] = state_update_dieterich(h, unorm / L, alpha_old[i]);
                   break;
                 case Config::Ruina:
-                  (*alpha)[i] = state_update_ruina(h, unorm / L, alpha_old[i]);
+                  alpha[i] = state_update_ruina(h, unorm / L, alpha_old[i]);
                   break;
               }
             }
@@ -318,7 +316,7 @@ int main(int argc, char *argv[]) {
         if (parset.get<bool>("writeEvolution")) {
           double out;
           neumannFunction.evaluate(time, out);
-          octave_writer << (*alpha)[first_frictional_node][0] << " "
+          octave_writer << alpha[first_frictional_node][0] << " "
                         << u[first_frictional_node][0] * 1e6 << " " << out
                         << std::endl;
         }
@@ -329,7 +327,7 @@ int main(int argc, char *argv[]) {
         if (parset.get<bool>("printVelocitySteppingComparison")) {
           double const v = u_diff[first_frictional_node].two_norm() / h / L;
 
-          double const euler = (*alpha)[first_frictional_node];
+          double const euler = alpha[first_frictional_node];
           double direct;
           if (run < 120) {
             direct = euler;
@@ -352,7 +350,7 @@ int main(int argc, char *argv[]) {
         // Record the coefficient of friction at a fixed node
         if (parset.get<bool>("printCoefficient")) {
           double const V = u_diff[first_frictional_node].two_norm();
-          double const state = (*alpha)[first_frictional_node];
+          double const state = alpha[first_frictional_node];
 
           coefficient_writer << (mu + a * std::log(V * eta) +
                                  b * (state - std::log(eta * L))) << std::endl;
@@ -360,7 +358,7 @@ int main(int argc, char *argv[]) {
       }
 
       u += u_diff;
-      alpha_old = *alpha;
+      alpha_old = alpha;
 
       // Compute von Mises stress and write everything to a file
       if (parset.get<bool>("writeVTK")) {
@@ -372,7 +370,7 @@ int main(int argc, char *argv[]) {
             .assemble(localStressAssembler, vonMisesStress, true);
 
         writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>(
-            p1Basis, u, *alpha, p0Basis, vonMisesStress, leafView,
+            p1Basis, u, alpha, p0Basis, vonMisesStress, leafView,
             (boost::format("obs%d") % run).str());
       }
     }