From c417302e37b41bc0e5a4b49efecbb4851af91373 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sat, 18 Jan 2020 15:12:09 +0100
Subject: [PATCH] Add a system test for harmonic maps
This new test computes a harmonic map from a square to a sphere.
In other words, rather than testing some specific aspect of some
class, it does a whole simulation. It should therefore cover
much more of the dune-gfe code.
---
test/CMakeLists.txt | 3 +
test/grids/irregular-square.msh | 42 +++++++
test/harmonicmaptest.cc | 195 ++++++++++++++++++++++++++++++++
3 files changed, 240 insertions(+)
create mode 100644 test/grids/irregular-square.msh
create mode 100644 test/harmonicmaptest.cc
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 8b4c20ff..78251121 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -4,6 +4,7 @@ set(TESTS
cosseratenergytest
frameinvariancetest
harmonicenergytest
+ harmonicmaptest
localgeodesicfefunctiontest
localgfetestfunctiontest
localprojectedfefunctiontest
@@ -19,3 +20,5 @@ foreach(_test ${TESTS})
dune_add_test(SOURCES ${_test}.cc)
target_compile_options(${_test} PRIVATE -g)
endforeach()
+
+file(COPY grids/irregular-square.msh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/grids)
diff --git a/test/grids/irregular-square.msh b/test/grids/irregular-square.msh
new file mode 100644
index 00000000..2cf3522d
--- /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 00000000..996588be
--- /dev/null
+++ b/test/harmonicmaptest.cc
@@ -0,0 +1,195 @@
+#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>;
+
+ std::shared_ptr<GridType> grid = GmshReader<GridType>::read("grids/irregular-square.msh");
+
+ grid->globalRefine(numLevels-1);
+
+ 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);
+ constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+
+
+ ////////////////////////////
+ // 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;
+}
--
GitLab