Enable caching for the laursen nonlinearity

dune/tectonic/globallaursennonlinearity.hh

@@ -23,7 +23,17 @@ class GlobalLaursenNonlinearity
       shared_ptr<BlockVector<FieldVector<double, 1>> const> mu,
       shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress,
       shared_ptr<BlockVector<FieldVector<double, 1>> const> nodalIntegrals)
-      : mu(mu), normalStress(normalStress), nodalIntegrals(nodalIntegrals) {}
+      : mu(mu),
+        normalStress(normalStress),
+        nodalIntegrals(nodalIntegrals),
+        trivialNonlinearity(
+            new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
+        restrictions(nodalIntegrals->size()) // TODO: can we get the size from
+                                             // another place?
+  {
+    for (auto &x : restrictions)
+      x = shared_ptr<LocalNonlinearity<dim> const>();
+  }
 
   /*
     Return a restriction of the outer function to the i'th node. If
@@ -40,18 +50,29 @@ class GlobalLaursenNonlinearity
     sigma_n [id + mu id] = sigma_n (1 + mu) id
   */
   virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
+    if ((*nodalIntegrals)[i][0] == 0)
+      return trivialNonlinearity;
+
+    if (restrictions[i] != nullptr)
+      return restrictions[i];
+
     double coefficient = (*nodalIntegrals)[i][0];
     coefficient *= (*normalStress)[i];
     coefficient *= 1 + (*mu)[i];
 
     auto const func = make_shared<OuterFunctionType const>(coefficient);
-    return make_shared<LocalNonlinearity<dim> const>(func);
+    restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
+    return restrictions[i];
   }
 
 private:
+  shared_ptr<LocalNonlinearity<dim> const> const trivialNonlinearity;
+
   shared_ptr<BlockVector<FieldVector<double, 1>> const> mu;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> normalStress;
   shared_ptr<BlockVector<FieldVector<double, 1>> const> nodalIntegrals;
+
+  std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
 };
 }
 #endif