From 13babae3732be97924402dd8423f4c200f101060 Mon Sep 17 00:00:00 2001
From: Patrick <patrick.jaap@tu-dresden.de>
Date: Fri, 8 May 2020 14:29:36 +0200
Subject: [PATCH] NeumannEnergy: use std::function instead of
 Dune::VirtualFunction

We also use std::shared_ptr's instead of C-pointers.
---
 dune/elasticity/materials/neumannenergy.hh | 28 ++++++---------
 src/finite-strain-elasticity.cc            | 40 ++++++++--------------
 2 files changed, 25 insertions(+), 43 deletions(-)

diff --git a/dune/elasticity/materials/neumannenergy.hh b/dune/elasticity/materials/neumannenergy.hh
index ea62bdc..6a11c49 100644
--- a/dune/elasticity/materials/neumannenergy.hh
+++ b/dune/elasticity/materials/neumannenergy.hh
@@ -1,9 +1,10 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_NEUMANNENERGY_HH
 #define DUNE_ELASTICITY_MATERIALS_NEUMANNENERGY_HH
 
+#include <functional>
+
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/fufem/functions/virtualgridfunction.hh>
 #include <dune/fufem/boundarypatch.hh>
 
 #include <dune/elasticity/assemblers/localenergy.hh>
@@ -26,7 +27,7 @@ public:
    * \param parameters The material parameters
    */
   NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView>>& neumannBoundary,
-                const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,dim> >* neumannFunction)
+                std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<ctype,dim>)> neumannFunction)
   : neumannBoundary_(neumannBoundary),
     neumannFunction_(neumannFunction)
   {}
@@ -66,12 +67,7 @@ public:
             value[j] += shapeFunctionValues[i] * localConfiguration[ localView.tree().child(j).localIndex(i) ];
 
         // Value of the Neumann data at the current position
-        Dune::FieldVector<double,dim> neumannValue;
-
-        if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,dim> >*>(neumannFunction_))
-            dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,dim> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
-        else
-            neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
+        auto neumannValue = neumannFunction_( it.geometry().global(quad[pt].position()) );
 
         // Only translational dofs are affected by the Neumann force
         for (size_t i=0; i<neumannValue.size(); i++)
@@ -89,7 +85,7 @@ private:
   const std::shared_ptr<BoundaryPatch<GridView>> neumannBoundary_;
 
   /** \brief The function implementing the Neumann data */
-  const Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,dim> >* neumannFunction_;
+  std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<ctype,dim>)> neumannFunction_;
 };
 
 }  // namespace Dune::Elasticity
@@ -113,8 +109,8 @@ public:
   /** \brief Constructor with a set of material parameters
    * \param parameters The material parameters
    */
-  NeumannEnergy(const BoundaryPatch<GridView>* neumannBoundary,
-                const Dune::VirtualFunction<Dune::FieldVector<ctype,dim>, Dune::FieldVector<double,dim> >* neumannFunction)
+  NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView>>& neumannBoundary,
+                const std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<ctype,dim>)>>& neumannFunction)
   : neumannBoundary_(neumannBoundary),
     neumannFunction_(neumannFunction)
   {}
@@ -154,12 +150,8 @@ public:
             value[j] += shapeFunctionValues[i] * localConfiguration[i][j];
 
         // Value of the Neumann data at the current position
-        Dune::FieldVector<double,dim> neumannValue;
+        auto neumannValue = (*neumannFunction_)( it.geometry().global(quad[pt].position()) );
 
-        if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,dim> >*>(neumannFunction_))
-            dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,dim> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
-        else
-            neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
 
         // Only translational dofs are affected by the Neumann force
         for (size_t i=0; i<neumannValue.size(); i++)
@@ -174,11 +166,11 @@ public:
 
 private:
   /** \brief The Neumann boundary */
-  const BoundaryPatch<GridView>*
+  const std::shared_ptr<BoundaryPatch<GridView>>
     DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::NeumannEnergy") neumannBoundary_;
 
   /** \brief The function implementing the Neumann data */
-  const Dune::VirtualFunction<Dune::FieldVector<double,dim>, Dune::FieldVector<double,dim> >* neumannFunction_;
+  const std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<ctype,dim>)>> neumannFunction_;
 };
 
 }  // namespace Dune
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index adf18af..8116391 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -1,4 +1,5 @@
 #include <config.h>
+#include <functional>
 
 // Includes for the ADOL-C automatic differentiation library
 // Need to come before (almost) all others.
@@ -49,27 +50,6 @@ const int order = 1;
 
 using namespace Dune;
 
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
-    : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,dim> >
-{
-  NeumannFunction(const FieldVector<double,dim> values,
-                  double homotopyParameter)
-  : values_(values),
-    homotopyParameter_(homotopyParameter)
-  {}
-
-  void evaluate(const FieldVector<double, dim>& x, FieldVector<double,dim>& out) const
-  {
-    out = 0;
-    out.axpy(-homotopyParameter_, values_);
-  }
-
-  FieldVector<double,dim> values_;
-  double homotopyParameter_;
-};
-
-
 int main (int argc, char *argv[]) try
 {
   // initialize MPI, finalize is done automatically on exit
@@ -255,13 +235,23 @@ int main (int argc, char *argv[]) try
     // ////////////////////////////////////////////////////////////
 
     const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
-    std::shared_ptr<NeumannFunction> neumannFunction;
+
+    FieldVector<double,dim> neumannValues(0.);
     if (parameterSet.hasKey("neumannValues"))
     {
-      neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,dim> >("neumannValues"),
-                                                     homotopyParameter);
+      neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
+    }
+
+    // create simple constant valued Neumann function
+    auto neumannFunction = [&]( FieldVector<double,dim> )
+    {
+      auto nv = neumannValues;
+      nv *= (-homotopyParameter);
+      return nv;
+    };
 
     if (mpiHelper.rank()==0)
+    {
       std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,dim> >("neumannValues") << std::endl;
     }
 
@@ -293,7 +283,7 @@ int main (int argc, char *argv[]) try
 
     auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);
 
-    auto neumannEnergy = std::make_shared<Elasticity::NeumannEnergy<LocalView,adouble>>(neumannBoundary,neumannFunction.get());
+    auto neumannEnergy = std::make_shared<Elasticity::NeumannEnergy<LocalView,adouble>>(neumannBoundary,neumannFunction);
 
     auto totalEnergy = std::make_shared<Elasticity::SumEnergy<LocalView,adouble>>(elasticEnergy, neumannEnergy);
 
-- 
GitLab