diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index abe48b260d2b9bc6d8fb6839c644d0d8d4441261..ee20d468d3f64745681e89044d6a150a5875a3ba 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -31,9 +31,9 @@ class CosseratVTKWriter
v1Embedded[i] = v1[i].globalCoordinates();
// Interpolate
- auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,7> >(basis1, Dune::TypeTree::hybridTreePath(), v1Embedded);
+ auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,7> >(basis1, v1Embedded);
std::vector<Dune::FieldVector<double,7> > v2Embedded;
- Dune::Functions::interpolate(basis2, Dune::TypeTree::hybridTreePath(), v2Embedded, function);
+ Dune::Functions::interpolate(basis2, v2Embedded, function);
// Copy back from R^7 into RigidBodyMotions
v2.resize(v2Embedded.size());
@@ -51,9 +51,9 @@ class CosseratVTKWriter
v1Embedded[i] = v1[i].globalCoordinates();
// Interpolate
- auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,3> >(basis1, Dune::TypeTree::hybridTreePath(), v1Embedded);
+ auto function = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,3> >(basis1, v1Embedded);
std::vector<Dune::FieldVector<double,3> > v2Embedded;
- Dune::Functions::interpolate(basis2, Dune::TypeTree::hybridTreePath(), v2Embedded, function);
+ Dune::Functions::interpolate(basis2, v2Embedded, function);
// Copy back from FieldVector to RealTuple
v2.resize(v2Embedded.size());
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 0304ca2c1b60d52267c27e0154db59f1c6e8c669..cd6febcc4a33c3a14b7fae089f0e496e0c5ab198 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -476,7 +476,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
if (instrumented_) {
- fprintf(fp, "Trust-region step: %d, trust-region radius: %g\n",
+ fprintf(fp, "Trust-region step: %ld, trust-region radius: %g\n",
i, trustRegion.radius());
// ///////////////////////////////////////////////////////////////
@@ -644,7 +644,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
if (instrumented_) {
char iFilename[100];
- sprintf(iFilename, "tmp/intermediateSolution_%04d", i);
+ sprintf(iFilename, "tmp/intermediateSolution_%04ld", i);
FILE* fpIterate = fopen(iFilename, "wb");
if (!fpIterate)
diff --git a/dune/gfe/sumcosseratenergy.hh b/dune/gfe/sumcosseratenergy.hh
index 9cf25de0b61431351e1667e7bba8d9b4658e4508..51f6241a2090e51977d9fc79f0f426bf2e49911f 100644
--- a/dune/gfe/sumcosseratenergy.hh
+++ b/dune/gfe/sumcosseratenergy.hh
@@ -29,6 +29,8 @@ public:
*/
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
+#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
+ SumCosseratEnergy(std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#else
SumCosseratEnergy(std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#endif
@@ -55,6 +57,8 @@ private:
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
+#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
+ std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#else
std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#endif
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 4c2cfa1ae0a52eafade6e51c67ceaf78d8896e74..c63d0c092e470a0c594de4d8ce58d0c4d644008f 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -22,6 +22,16 @@
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/io/file/vtk.hh>
+#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
+#include <dune/elasticity/materials/exphenckydensity.hh>
+#include <dune/elasticity/materials/henckydensity.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+#include <dune/elasticity/materials/neohookedensity.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
+#include <dune/elasticity/materials/sumenergy.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
+#include <dune/elasticity/materials/neumannenergy.hh>
+#else
#include <dune/elasticity/materials/exphenckyenergy.hh>
#include <dune/elasticity/materials/henckyenergy.hh>
#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
@@ -31,6 +41,8 @@
#include <dune/elasticity/materials/neumannenergy.hh>
#include <dune/elasticity/materials/sumenergy.hh>
#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
+#endif
+
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
@@ -79,7 +91,7 @@ struct Dune::MathematicalConstants<adouble>
typedef adouble ValueType;
using namespace Dune;
-
+#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
/** \brief A constant vector-valued function, for simple Neumann boundary values */
struct NeumannFunction
: public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,dim> >
@@ -99,25 +111,7 @@ struct NeumannFunction
FieldVector<double,dim> values_;
double homotopyParameter_;
};
-
-struct VolumeLoad
- : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
-{
- VolumeLoad(const FieldVector<double,3> values,
- double homotopyParameter)
- : values_(values),
- homotopyParameter_(homotopyParameter)
- {}
-
- void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const {
- out = 0;
- out.axpy(homotopyParameter_, values_);
- }
-
- FieldVector<double,3> values_;
- double homotopyParameter_;
-};
-
+#endif
int main (int argc, char *argv[]) try
{
@@ -241,11 +235,11 @@ int main (int argc, char *argv[]) try
}
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);
BoundaryPatch<GridView> surfaceShellBoundary(gridView, surfaceShellVertices);
if (mpiHelper.rank()==0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+ std::cout << "Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";
BitSetVector<1> dirichletNodes(feBasis.size(), false);
@@ -259,9 +253,9 @@ int main (int argc, char *argv[]) try
BitSetVector<1> neumannNodes(feBasis.size(), false);
#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
- constructBoundaryDofs(neumannBoundary,fufemFeBasis,neumannNodes);
+ constructBoundaryDofs(*neumannBoundary,fufemFeBasis,neumannNodes);
#else
- constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
+ constructBoundaryDofs(*neumannBoundary,feBasis,neumannNodes);
#endif
BitSetVector<1> surfaceShellNodes(feBasis.size(), false);
@@ -325,7 +319,10 @@ int main (int argc, char *argv[]) try
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
auto localDisplacementFunction = localFunction(displacementFunction);
- auto neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
+ FieldVector<double,dim> neumannValues {0,0,0};
+
+ if (parameterSet.hasKey("neumannValues"))
+ neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
std::cout << "Neumann values: " << neumannValues << std::endl;
// We need to subsample, because VTK cannot natively display real second-order functions
@@ -343,16 +340,19 @@ int main (int argc, char *argv[]) try
const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
- std::shared_ptr<VolumeLoad> volumeLoad;
-
+#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 8)
std::shared_ptr<NeumannFunction> neumannFunction;
- if (parameterSet.hasKey("neumannValues"))
- {
- neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,dim> >("neumannValues"),
+ neumannFunction = std::make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,dim> >("neumannValues"),
homotopyParameter);
-
- std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,dim> >("neumannValues") << std::endl;
- }
+#else
+ // A constant vector-valued function, for simple Neumann boundary values
+ std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>> neumannFunctionPtr;
+ neumannFunctionPtr = std::make_shared<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<double,dim>)>>([&](FieldVector<double,dim> ) {
+ auto nv = neumannValues;
+ nv *= (-homotopyParameter);
+ return nv;
+ });
+#endif
if (mpiHelper.rank() == 0)
{
@@ -363,6 +363,31 @@ int main (int argc, char *argv[]) try
// Assembler using ADOL-C
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
+#if DUNE_VERSION_GTE(DUNE_ELASTICITY, 2,8)
+ std::shared_ptr<Elasticity::LocalDensity<dim,ValueType>> elasticDensity;
+ if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
+ elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "neohooke")
+ elasticDensity = std::make_shared<Elasticity::NeoHookeDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "hencky")
+ elasticDensity = std::make_shared<Elasticity::HenckyDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "exphencky")
+ elasticDensity = std::make_shared<Elasticity::ExpHenckyDensity<dim,ValueType>>(materialParameters);
+ if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
+ elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(materialParameters);
+
+ if(!elasticDensity)
+ DUNE_THROW(Exception, "Error: Selected energy not available!");
+
+ auto elasticEnergy = std::make_shared<LocalIntegralEnergy<GridView, FEBasis::LocalView::Tree::FiniteElement, ValueType>>(elasticDensity);
+ auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType> >(neumannBoundary,neumannFunctionPtr);
+
+ auto elasticAndNeumann = std::make_shared<Dune::SumEnergy<GridView,
+ FEBasis::LocalView::Tree::FiniteElement,
+ ValueType>>(elasticEnergy, neumannEnergy);
+#else
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
std::shared_ptr<LocalFEStiffness<GridView,
#else
@@ -402,12 +427,12 @@ int main (int argc, char *argv[]) try
auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
- ValueType> >(&neumannBoundary,neumannFunction.get());
+ ValueType> >(neumannBoundary.get(),neumannFunction.get());
auto elasticAndNeumann = std::make_shared<SumEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
ValueType>>(elasticEnergy, neumannEnergy);
-
+#endif
using LocalEnergyBase = GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble, dim> >;
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 0ada4c171bd8e209fb95163e7a9d887552fc3a79..3879b56e0fd23777645367382f2b99439832b2e3 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -242,7 +242,6 @@ int main (int argc, char *argv[]) try
xEmbedded[i] = x[i].globalCoordinates();
auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,
- TypeTree::hybridTreePath(),
xEmbedded);
SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),refinementLevels(0));
@@ -282,7 +281,6 @@ int main (int argc, char *argv[]) try
xEmbedded[i] = x[i].globalCoordinates();
auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,
- TypeTree::hybridTreePath(),
xEmbedded);
SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),refinementLevels(0));
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 772fc06e6103f0b407841c6a42b496596f107bd5..3cfa88935e56fe0dbc0d73e4d15b308097622056 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -108,7 +108,7 @@ void fillVTKWriter(Writer& vtkWriter, const Basis& feBasis, const SolutionType&
}
else
{
- auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,TypeTree::hybridTreePath(),xEmbedded);
+ auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,xEmbedded);
vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::vector, xEmbedded[0].size()));
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
index 2dbe77293ec5fa54031f577c2bc45bbddf876644..1b11a59d11f49b89f9b64fc3f1b92e3b745d269e 100644
--- a/test/harmonicmaptest.cc
+++ b/test/harmonicmaptest.cc
@@ -35,9 +35,6 @@ const int dim = 2;
// Image space of the geodesic fe functions
typedef UnitVector<double,3> TargetSpace;
-// Tangent vector of the image space
-const int blocksize = TargetSpace::TangentVector::dimension;
-
const int order = 1;
using namespace Dune;
@@ -87,6 +84,24 @@ int main (int argc, char *argv[])
FEBasis feBasis(gridView);
SolutionType x(feBasis.size());
+ using namespace Functions::BasisFactory;
+
+ // A basis for the embedding space
+ auto powerBasis = makeBasis(
+ gridView,
+ power<TargetSpace::CoordinateType::dimension>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
+ // A basis for the tangent space
+ auto tangentBasis = makeBasis(
+ gridView,
+ power<TargetSpace::TangentVector::dimension>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
///////////////////////////////////////////
// Determine Dirichlet values
///////////////////////////////////////////
@@ -107,12 +122,12 @@ int main (int argc, char *argv[])
}
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
+ BitSetVector<TargetSpace::TangentVector::dimension> dirichletNodes(powerBasis.size(), false);
#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
DuneFunctionsBasis<FEBasis> fufemBasis(feBasis);
constructBoundaryDofs(dirichletBoundary,fufemBasis,dirichletNodes);
#else
- constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+ constructBoundaryDofs(dirichletBoundary,tangentBasis,dirichletNodes);
#endif
////////////////////////////
@@ -127,15 +142,6 @@ int main (int argc, char *argv[])
};
std::vector<TargetSpace::CoordinateType> v;
- using namespace Functions::BasisFactory;
-
- auto powerBasis = makeBasis(
- gridView,
- power<TargetSpace::CoordinateType::dimension>(
- lagrange<order>(),
- blockedInterleaved()
- ));
-
Dune::Functions::interpolate(powerBasis, v, initialIterateFunction);
for (size_t i=0; i<x.size(); i++)