diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 778effaa91120f11ee753d52aadeb13c716011ca..8feb344f72223c699bc0ad7ca13b35801efaaf1c 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -131,11 +131,21 @@ setup(const GridType& grid,
LocalMapper,
LocalMapper> matrixComm(*globalMapper_, grid_->leafGridView(), localMapper, localMapper, 0);
- auto A = std::make_shared<ScalarMatrixType>(matrixComm.reduceAdd(localA));
+#if !DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ auto
+#endif
+ A = std::make_shared<ScalarMatrixType>(matrixComm.reduceAdd(localA));
#else
- auto A = std::make_shared<ScalarMatrixType>(localA);
+#if !DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ auto
+#endif
+ A = std::make_shared<ScalarMatrixType>(localA);
#endif
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ h1SemiNorm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*A);
+#else
h1SemiNorm_ = std::make_shared<H1SemiNorm<CorrectionType> >(A);
+#endif
innerSolver_ = std::make_shared<::LoopSolver<CorrectionType> >(mmgStep,
innerIterations_,
@@ -153,12 +163,19 @@ setup(const GridType& grid,
operatorAssembler.assemble(massStiffness, localMassMatrix);
+#if !DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ auto
+#endif
#if HAVE_MPI
- auto massMatrix = std::make_shared<ScalarMatrixType>(matrixComm.reduceAdd(localMassMatrix));
+ massMatrix = std::make_shared<ScalarMatrixType>(matrixComm.reduceAdd(localMassMatrix));
#else
- auto massMatrix = std::make_shared<ScalarMatrixType>(localMassMatrix);
+ massMatrix = std::make_shared<ScalarMatrixType>(localMassMatrix);
#endif
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
+#else
l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(massMatrix);
+#endif
// Write all intermediate solutions, if requested
if (instrumented_
@@ -351,7 +368,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
localMapper,
0);
#endif
- for (int i=0; i<maxTrustRegionSteps_; i++) {
+ auto& i = statistics_.finalIteration;
+ for (i=0; i<maxTrustRegionSteps_; i++) {
/* std::cout << "current iterate:\n";
for (size_t j=0; j<x_.size(); j++)
@@ -643,4 +661,6 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
// //////////////////////////////////////////////
if (instrumented_)
fclose(fp);
+
+ statistics_.finalEnergy = oldEnergy;
}
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index add3d56d61c31853ab262531d23a93892bd08b1a..a3f5ec3d28a80c2cfdef844ccfb3d5f19e5a7670 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -92,6 +92,16 @@ class RiemannianTrustRegionSolver
typedef typename MapperFactory<typename Basis::GridView,Basis>::LocalMapper LocalMapper;
#endif
+ /** \brief Records information about the last run of the RiemannianTrustRegionSolver
+ *
+ * This is used primarily for unit testing.
+ */
+ struct Statistics
+ {
+ std::size_t finalIteration;
+
+ field_type finalEnergy;
+ };
public:
@@ -144,6 +154,8 @@ public:
SolutionType getSol() const {return x_;}
+ const Statistics& getStatistics() const {return statistics_;}
+
protected:
#if HAVE_MPI
@@ -163,7 +175,7 @@ protected:
Dune::FieldVector<double,blocksize> scaling_;
/** \brief Maximum number of trust-region steps */
- int maxTrustRegionSteps_;
+ std::size_t maxTrustRegionSteps_;
/** \brief Maximum number of multigrid iterations */
int innerIterations_;
@@ -195,6 +207,14 @@ protected:
/** \brief If set to true we log convergence speed and other stuff */
bool instrumented_;
+ /** \brief Store information about solver runs for unit testing */
+ Statistics statistics_;
+
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ std::shared_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > > A;
+ std::shared_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > > massMatrix;
+#endif
+
};
#include "riemanniantrsolver.cc"
diff --git a/dune/gfe/sumcosseratenergy.hh b/dune/gfe/sumcosseratenergy.hh
index 2011d909604194c646996056ea814ad55253d8ca..9cf25de0b61431351e1667e7bba8d9b4658e4508 100644
--- a/dune/gfe/sumcosseratenergy.hh
+++ b/dune/gfe/sumcosseratenergy.hh
@@ -27,7 +27,11 @@ public:
* \param elasticEnergy The elastic energy
* \param cosseratEnergy The cosserat energy
*/
+#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
+ SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<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
std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace>> cosseratEnergy)
: elasticEnergy_(elasticEnergy),
@@ -49,7 +53,11 @@ public:
private:
+#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
+ std::shared_ptr<LocalFEStiffness<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
std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace> > cosseratEnergy_;
};
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 832d3955894819ccfb705c659bfe8a9cf5870082..e34df722a89618c99398b706ec56d4a6ed30295a 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -304,10 +304,18 @@ int main (int argc, char *argv[]) try
orientationDirichletDofs[i][j] = true;
#else
BitSetVector<1> dirichletNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
+#else
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+#endif
BitSetVector<1> neumannNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ constructBoundaryDofs(neumannBoundary,fufemFeBasis,neumannNodes);
+#else
constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
+#endif
BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
for (size_t i=0; i<feBasis.size(); i++)
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index fd9e1fbe4fdaa05ccc06a5eac834b2d739a3b221..dbd81056278ea589e34e932c598e18bc23183136 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -14,6 +14,7 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
+#include <dune/common/version.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
@@ -23,7 +24,9 @@
#include <dune/elasticity/materials/exphenckyenergy.hh>
#include <dune/elasticity/materials/henckyenergy.hh>
+#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
+#endif
#include <dune/elasticity/materials/neohookeenergy.hh>
#include <dune/elasticity/materials/neumannenergy.hh>
#include <dune/elasticity/materials/sumenergy.hh>
@@ -59,6 +62,19 @@
const int dim = WORLD_DIM;
const int order = 1;
+#if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
+template<>
+struct Dune::MathematicalConstants<adouble>
+{
+ static const adouble pi ()
+ {
+ using std::acos;
+ static const adouble pi = acos( adouble( -1 ) );
+ return pi;
+ }
+};
+#endif
+
//differentiation method
typedef adouble ValueType;
@@ -214,16 +230,13 @@ int main (int argc, char *argv[]) try
PythonFunction<FieldVector<double,dim>, bool> pythonSurfaceShellVertices(Python::evaluate(lambda));
for (auto&& v : vertices(gridView))
{
- bool isDirichlet;
- pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
+ bool isDirichlet = pythonDirichletVertices(v.geometry().corner(0));
dirichletVertices[indexSet.index(v)] = isDirichlet;
- bool isNeumann;
- pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
+ bool isNeumann = pythonNeumannVertices(v.geometry().corner(0));
neumannVertices[indexSet.index(v)] = isNeumann;
- bool isSurfaceShell;
- pythonSurfaceShellVertices.evaluate(v.geometry().corner(0), isSurfaceShell);
+ bool isSurfaceShell = pythonSurfaceShellVertices(v.geometry().corner(0));
surfaceShellVertices[indexSet.index(v)] = isSurfaceShell;
}
@@ -236,13 +249,27 @@ int main (int argc, char *argv[]) try
BitSetVector<1> dirichletNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
+ using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
+ FufemFEBasis fufemFeBasis(feBasis);
+ constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
+#else
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+#endif
BitSetVector<1> neumannNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
+ constructBoundaryDofs(neumannBoundary,fufemFeBasis,neumannNodes);
+#else
constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
+#endif
BitSetVector<1> surfaceShellNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_FUFEM, 2, 7)
+ constructBoundaryDofs(surfaceShellBoundary,fufemFeBasis,surfaceShellNodes);
+#else
constructBoundaryDofs(surfaceShellBoundary,feBasis,surfaceShellNodes);
+#endif
BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
for (size_t i=0; i<feBasis.size(); i++)
@@ -336,7 +363,11 @@ 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_LT(DUNE_ELASTICITY, 2, 7)
+ std::shared_ptr<LocalFEStiffness<GridView,
+#else
std::shared_ptr<Elasticity::LocalEnergy<GridView,
+#endif
FEBasis::LocalView::Tree::FiniteElement,
std::vector<Dune::FieldVector<ValueType, dim>> > > elasticEnergy;
@@ -344,10 +375,12 @@ int main (int argc, char *argv[]) try
elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
ValueType> >(materialParameters);
+#if !DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
ValueType> >(materialParameters);
+#endif
if (parameterSet.get<std::string>("energy") == "neohooke")
elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
@@ -477,4 +510,4 @@ int main (int argc, char *argv[]) try
}
} catch (Exception& e) {
std::cout << e.what() << std::endl;
-}
\ No newline at end of file
+}
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index e8c84d8be0c97181fd0d8d3fa97c4853673fb7ea..78cbc34011206c5a1443bdb5aecbf52bdf257457 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -234,8 +234,12 @@ int main (int argc, char *argv[])
BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ DuneFunctionsBasis<FEBasis> fufemBasis(feBasis);
+ constructBoundaryDofs(dirichletBoundary,fufemBasis,dirichletNodes);
+#else
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
-
+#endif
// //////////////////////////
// Initial iterate
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 8b4c20ff302f3f6ba957e044fa1a5da58899bcaf..e28375730ca3df2d1322b4ba0f9d0586ba21583d 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -19,3 +19,12 @@ foreach(_test ${TESTS})
dune_add_test(SOURCES ${_test}.cc)
target_compile_options(${_test} PRIVATE -g)
endforeach()
+
+# Run distributed tests
+dune_add_test(SOURCES harmonicmaptest.cc
+ MPI_RANKS 1 4
+ TIMEOUT 600
+ CMAKE_GUARD MPI_FOUND)
+
+# Copy the example grid used for testing into the build dir
+file(COPY grids/irregular-square.msh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/grids)
diff --git a/test/adolctest.cc b/test/adolctest.cc
index bb04d8fa380a0c094b7f6510b95e04a471a9def6..be7ac34c8dfa24d1331cedd4fafc207f054537c8 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -412,6 +412,8 @@ void compareMatrices(const Matrix<FieldMatrix<double,N,N> >& matrixA, std::strin
int main (int argc, char *argv[]) try
{
+ MPIHelper::instance(argc, argv);
+
typedef std::vector<TargetSpace> SolutionType;
enum { embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension };
enum { blocksize = TargetSpace::TangentVector::dimension };
diff --git a/test/grids/irregular-square.msh b/test/grids/irregular-square.msh
new file mode 100644
index 0000000000000000000000000000000000000000..2cf3522da82264d2dd27e793e8c217e41e553972
--- /dev/null
+++ b/test/grids/irregular-square.msh
@@ -0,0 +1,42 @@
+$MeshFormat
+2.1 0 8
+$EndMeshFormat
+
+$Nodes
+17
+ 1 -5 -5 0
+ 2 -2 -5 0
+ 3 3 -5 0
+ 4 5 -5 0
+ 5 -5 -2 0
+ 6 -2.8 -1.8 0
+ 7 1 -1.5 0
+ 8 5 -2 0
+ 9 -5 4 0
+10 -3 1.5 0
+11 1 2.5 0
+12 4 3.1 0
+13 -5 5 0
+14 -1 5 0
+15 1 5 0
+16 5 5 0
+17 5 0 0
+$EndNodes
+
+$Elements
+14
+ 1 2 3 0 1 0 1 2 6
+ 2 2 3 0 1 0 2 7 6
+ 3 3 3 0 1 0 2 3 8 7
+ 4 2 3 0 1 0 3 4 8
+ 5 3 3 0 1 0 1 6 10 5
+ 6 3 3 0 1 0 6 7 11 10
+ 7 3 3 0 1 0 7 17 12 11
+ 8 2 3 0 1 0 7 8 17
+ 9 2 3 0 1 0 5 10 9
+10 3 3 0 1 0 9 10 14 13
+11 2 3 0 1 0 10 11 14
+12 3 3 0 1 0 11 12 15 14
+13 2 3 0 1 0 12 16 15
+14 2 3 0 1 0 12 17 16
+$EndElements
diff --git a/test/harmonicmaptest.cc b/test/harmonicmaptest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2dbe77293ec5fa54031f577c2bc45bbddf876644
--- /dev/null
+++ b/test/harmonicmaptest.cc
@@ -0,0 +1,205 @@
+#include <config.h>
+
+#include <array>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
+#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/harmonicenergy.hh>
+#include <dune/gfe/geodesicfeassembler.hh>
+#include <dune/gfe/riemanniantrsolver.hh>
+
+// grid dimension
+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;
+
+
+int main (int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+
+ using SolutionType = std::vector<TargetSpace>;
+
+ // read problem settings
+ const int numLevels = 4;
+
+ // read solver settings
+ const double tolerance = 1e-6;
+ const int maxTrustRegionSteps = 1000;
+ const double initialTrustRegionRadius = 1;
+ const int multigridIterations = 200;
+ const int baseIterations = 100;
+ const double mgTolerance = 1e-10;
+ const double baseTolerance = 1e-8;
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ using GridType = UGGrid<dim>;
+
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ std::shared_ptr<GridType> grid(GmshReader<GridType>::read("grids/irregular-square.msh"));
+#else
+ std::shared_ptr<GridType> grid = GmshReader<GridType>::read("grids/irregular-square.msh");
+#endif
+
+ grid->globalRefine(numLevels-1);
+
+ grid->loadBalance();
+
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ //////////////////////////////////////////////////////////////////////////////////
+ // Construct the scalar function space basis corresponding to the GFE space
+ //////////////////////////////////////////////////////////////////////////////////
+
+ using FEBasis = Functions::LagrangeBasis<GridView, order>;
+ FEBasis feBasis(gridView);
+ SolutionType x(feBasis.size());
+
+ ///////////////////////////////////////////
+ // Determine Dirichlet values
+ ///////////////////////////////////////////
+ BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+ const GridView::IndexSet& indexSet = gridView.indexSet();
+
+ // Make predicate function that computes which vertices are on the Dirichlet boundary,
+ // based on the vertex positions.
+ auto dirichletVerticesPredicate = [](FieldVector<double,dim> x)
+ {
+ return (x[0] < -4.9999 or x[0] > 4.9999 or x[1] < -4.9999 or x[1] > 4.9999);
+ };
+
+ for (auto&& vertex : vertices(gridView))
+ {
+ bool isDirichlet = dirichletVerticesPredicate(vertex.geometry().corner(0));
+ dirichletVertices[indexSet.index(vertex)] = isDirichlet;
+ }
+
+ BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+ BitSetVector<blocksize> dirichletNodes(feBasis.size(), false);
+#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
+ DuneFunctionsBasis<FEBasis> fufemBasis(feBasis);
+ constructBoundaryDofs(dirichletBoundary,fufemBasis,dirichletNodes);
+#else
+ constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+#endif
+
+ ////////////////////////////
+ // Initial iterate
+ ////////////////////////////
+
+ // The inverse stereographic projection through the north pole
+ auto initialIterateFunction = [](FieldVector<double,dim> x) -> TargetSpace::CoordinateType
+ {
+ auto normSquared = x.two_norm2();
+ return {2*x[0] / (normSquared+1), 2*x[1] / (normSquared+1), (normSquared-1)/ (normSquared+1)};
+ };
+
+ 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++)
+ x[i] = v[i];
+
+ //////////////////////////////////////////////////////////////
+ // Create an assembler for the Harmonic Energy Functional
+ //////////////////////////////////////////////////////////////
+
+ using ATargetSpace = TargetSpace::rebind<adouble>::other;
+ //using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+ using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+
+ //std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy = std::make_shared<HarmonicEnergy<FEBasis, GeodesicInterpolationRule, ATargetSpace> >();
+ std::shared_ptr<GFE::LocalEnergy<FEBasis,ATargetSpace> > localEnergy = std::make_shared<HarmonicEnergy<FEBasis, ProjectedInterpolationRule, ATargetSpace> >();
+
+ LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy.get());
+
+ GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, &localGFEADOLCStiffness);
+
+ ///////////////////////////////////////////////////
+ // Create a Riemannian trust-region solver
+ ///////////////////////////////////////////////////
+
+ RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+ solver.setup(*grid,
+ &assembler,
+ x,
+ dirichletNodes,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ 3, 3, 1, // Multigrid V-cycle
+ baseIterations,
+ baseTolerance,
+ false); // Do not write intermediate iterates
+
+ ///////////////////////////////////////////////////////
+ // Solve!
+ ///////////////////////////////////////////////////////
+
+ solver.setInitialIterate(x);
+ solver.solve();
+
+ x = solver.getSol();
+
+ std::size_t expectedFinalIteration = 12;
+ if (solver.getStatistics().finalIteration != expectedFinalIteration)
+ {
+ std::cerr << "Trust-region solver did " << solver.getStatistics().finalIteration+1
+ << " iterations, instead of the expected '" << expectedFinalIteration+1 << "'!" << std::endl;
+ return 1;
+ }
+
+ double expectedEnergy = 12.2927849;
+ if ( std::abs(solver.getStatistics().finalEnergy - expectedEnergy) > 1e-7)
+ {
+ std::cerr << std::setprecision(9);
+ std::cerr << "Final energy is " << solver.getStatistics().finalEnergy
+ << " but '" << expectedEnergy << "' was expected!" << std::endl;
+ return 1;
+ }
+
+ return 0;
+}