diff --git a/dune/elasticity/materials/neumannenergy.hh b/dune/elasticity/materials/neumannenergy.hh index ea62bdcbd147be512da772bf4013ad6df7431b29..6a11c498a836068fa5d59ce40d35d668374be695 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 adf18afa768097fd76f4fb08e8b0ba015943e430..8116391303bd35d05d6699f301da03af2166fc00 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);