Kill Laursen's nonlinearity and useless functions

dune/tectonic/globallaursennonlinearity.hh

-#ifndef DUNE_TECTONIC_GLOBAL_LAURSEN_NONLINEARITY_HH
-#define DUNE_TECTONIC_GLOBAL_LAURSEN_NONLINEARITY_HH
-
-#include <vector>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/shared_ptr.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include "globalnonlinearity.hh"
-#include "localnonlinearity.hh"
-#include "nicefunction.hh"
-
-namespace Dune {
-template <class OuterFunctionType, class VectorType, class MatrixType>
-class GlobalLaursenNonlinearity
-    : public GlobalNonlinearity<VectorType, MatrixType> {
-private:
-  typedef shared_ptr<BlockVector<FieldVector<double, 1>> const> dataptr;
-
-public:
-  using GlobalNonlinearity<VectorType, MatrixType>::dim;
-
-  GlobalLaursenNonlinearity(dataptr mu, dataptr normalStress,
-                            dataptr 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
-    mu and sigma_n denote the coefficient of friction and the normal
-    stress, respectively, at the i'th node, this function is given
-    by
-
-    sigma_n [(1/eta \bar Gamma)* + mu id]
-
-    TODO: We currently assume eta = 1
-  */
-  virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
-    if ((*nodalIntegrals)[i] == 0)
-      return trivialNonlinearity;
-
-    if (restrictions[i] != nullptr)
-      return restrictions[i];
-
-    double coefficient = (*nodalIntegrals)[i];
-    coefficient *= (*normalStress)[i];
-    coefficient *= 1 + (*mu)[i];
-
-    auto const func = make_shared<OuterFunctionType const>(coefficient);
-    restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
-    return restrictions[i];
-  }
-
-private:
-  dataptr mu;
-  dataptr normalStress;
-  dataptr nodalIntegrals;
-
-  shared_ptr<LocalNonlinearity<dim> const> const trivialNonlinearity;
-
-  std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
-};
-}
-#endif

dune/tectonic/nicefunction.hh

@@ -16,13 +16,9 @@ class NiceFunction {
   double virtual leftDifferential(double s) const = 0;
   double virtual rightDifferential(double s) const = 0;
 
-  double virtual second_deriv(double x) const {
-    DUNE_THROW(NotImplemented, "second derivative not implemented");
-  }
+  double virtual second_deriv(double x) const = 0;
 
-  double virtual regularity(double s) const {
-    DUNE_THROW(NotImplemented, "regularity not implemented");
-  }
+  double virtual regularity(double s) const = 0;
 
   // Whether H(|.|) is smooth at zero
   bool virtual smoothesNorm() const { return false; }
@@ -85,50 +81,6 @@ class RuinaFunction : public NiceFunction {
   double const rho;
 };
 
-class LinearFunction : public NiceFunction {
-public:
-  LinearFunction(double a) : coefficient(a) {}
-
-  void virtual evaluate(double const &x, double &y) const {
-    y = coefficient * x;
-  }
-
-  double virtual leftDifferential(double s) const { return coefficient; }
-
-  double virtual rightDifferential(double s) const { return coefficient; }
-
-  double virtual second_deriv(double) const { return 0; }
-
-  double virtual regularity(double s) const { return 0; }
-
-private:
-  double const coefficient;
-};
-
-template <int slope> class SampleFunction : public NiceFunction {
-public:
-  void virtual evaluate(double const &x, double &y) const {
-    y = (x < 1) ? x : (slope * (x - 1) + 1);
-  }
-
-  double virtual leftDifferential(double s) const {
-    return (s <= 1) ? 1 : slope;
-  }
-
-  double virtual rightDifferential(double s) const {
-    return (s < 1) ? 1 : slope;
-  }
-};
-
-class SteepFunction : public NiceFunction {
-public:
-  void virtual evaluate(double const &x, double &y) const { y = 100 * x; }
-
-  double virtual leftDifferential(double s) const { return 100; }
-
-  double virtual rightDifferential(double s) const { return 100; }
-};
-
 class TrivialFunction : public NiceFunction {
 public:
   void virtual evaluate(double const &x, double &y) const { y = 0; }
@@ -143,61 +95,5 @@ class TrivialFunction : public NiceFunction {
 
   bool virtual smoothesNorm() const { return true; }
 };
-
-// slope in [n-1,n] is n
-class HorribleFunction : public NiceFunction {
-public:
-  void virtual evaluate(double const &x, double &y) const {
-    double const fl = floor(x);
-    double const sum = fl * (fl + 1) / 2;
-    y = sum + (fl + 1) * (x - fl);
-  }
-
-  double virtual leftDifferential(double x) const {
-    double const fl = floor(x);
-    if (x - fl < 1e-14)
-      return fl;
-    else
-      return fl + 1;
-  }
-
-  double virtual rightDifferential(double x) const {
-    double const c = ceil(x);
-    if (c - x < 1e-14)
-      return c + 1;
-    else
-      return c;
-  }
-};
-
-// slope in [n-1,n] is log(n+1)
-class HorribleFunctionLogarithmic : public NiceFunction {
-public:
-  void virtual evaluate(double const &x, double &y) const {
-    y = 0;
-    size_t const fl = floor(x);
-    for (size_t i = 1; i <= fl;)
-      y += std::log(
-          ++i); // factorials grow to fast so we compute this incrementally
-
-    y += std::log(fl + 2) * (x - fl);
-  }
-
-  double virtual leftDifferential(double x) const {
-    double const fl = floor(x);
-    if (x - fl < 1e-14)
-      return std::log(fl + 1);
-    else
-      return std::log(fl + 2);
-  }
-
-  double virtual rightDifferential(double x) const {
-    double const c = ceil(x);
-    if (c - x < 1e-14)
-      return std::log(c + 2);
-    else
-      return std::log(c + 1);
-  }
-};
 }
 #endif

src/assemblers.cc

@@ -7,7 +7,6 @@
 #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
 #include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
 
-#include <dune/tectonic/globallaursennonlinearity.hh>
 #include <dune/tectonic/globalruinanonlinearity.hh>
 
 #include "assemblers.hh"
@@ -90,14 +89,7 @@ assemble_nonlinearity(
           nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
     }
     case Config::Laursen:
-      assert(false); // Known to be broken
-
-      // TODO: take state and h into account
-      // FIXME: We should be using a quadratic rather than a linear function
-      // here!
-      return Dune::make_shared<Dune::GlobalLaursenNonlinearity<
-          Dune::LinearFunction, VectorType, MatrixType> const>(mu, normalStress,
-                                                               nodalIntegrals);
+      assert(false);
   }
 }