Skip to content
Snippets Groups Projects
Commit 13babae3 authored by Patrick Jaap's avatar Patrick Jaap
Browse files

NeumannEnergy: use std::function instead of Dune::VirtualFunction 

We also use std::shared_ptr's instead of C-pointers.
parent fd0913ca
Branches
No related tags found
1 merge request!32NeumannEnergy: use std::function instead of Dune::VirtualFunction
Pipeline #28506 passed
Stage: test
Show whitespace changes
Inline Side-by-side
dune/elasticity/materials/neumannenergy.hh
+ 10
− 18
View file @ 13babae3
#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
......
src/finite-strain-elasticity.cc
+ 15
− 25
View file @ 13babae3
#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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment