diff --git a/problems/gradient-flow-curve-shortening.parset b/problems/gradient-flow-curve-shortening.parset
index 0eccf12e79723a6e58897744c7916294ecd1b9ed..2412d9f1566b0a46d2680ffbd23d2643f20c7150 100644
--- a/problems/gradient-flow-curve-shortening.parset
+++ b/problems/gradient-flow-curve-shortening.parset
@@ -5,7 +5,7 @@
structuredGrid = true
lower = 0
upper = 1
-elements = 1
+elements = 200
# Number of grid levels
numLevels = 1
@@ -15,13 +15,13 @@ numLevels = 1
#############################################
# Tolerance of the trust region solver
-tolerance = 1e-12
+tolerance = 1e-8
# Max number of steps of the trust region solver
-maxTrustRegionSteps = 0
+maxTrustRegionSteps = 100
# Initial trust-region radius
-initialTrustRegionRadius = 1
+initialTrustRegionRadius = 0.25
# Number of multigrid iterations per trust-region step
numIt = 200
@@ -53,3 +53,9 @@ energy = harmonic
# Inverse stereographic projection
initialIterate = initial-curve
+
+# Number of time steps
+numTimeSteps = 50
+
+# Time step size
+timeStepSize = 2e-4
diff --git a/problems/initial-curve.py b/problems/initial-curve.py
index 538f40c14de17b1a3e272f13b114e28f44f36b4d..05a29a529aaa793c6c7249a8819352d0698acbd1 100644
--- a/problems/initial-curve.py
+++ b/problems/initial-curve.py
@@ -1,6 +1,9 @@
import math
# Wiggly initial curve on the sphere in R^3 for a curve-shortening flow simulation
+#
+# Note: This method does not have to return unit vectors. Any non-zero vectors
+# will do; they are normalized after reading anyway.
def f(x):
- return [math.sin(0.5*math.pi*x[0]), math.cos(0.5*math.pi*x[0]), 0]
+ return [math.sin(0.5*math.pi*x[0]), math.cos(0.5*math.pi*x[0]), math.sin(10*math.pi*x[0])]
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 58536e58fd7ae6af7d53730163f562f641ae4565..596bf93ad431b87ad24332259c0976cbc457824d 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -31,14 +31,15 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
-#include <dune/gfe/rotation.hh>
#include <dune/gfe/unitvector.hh>
-#include <dune/gfe/realtuple.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
-#include <dune/gfe/harmonicenergystiffness.hh>
#include <dune/gfe/geodesicfeassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/globalgeodesicfefunction.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
+#include <dune/gfe/harmonicenergystiffness.hh>
+#include <dune/gfe/l2distancesquaredenergy.hh>
+#include <dune/gfe/weightedsumenergy.hh>
// grid dimension
const int dim = 1;
@@ -85,6 +86,8 @@ int main (int argc, char *argv[]) try
// read problem settings
const int numLevels = parameterSet.get<int>("numLevels");
+ const double timeStepSize = parameterSet.get<double>("timeStepSize");
+ const int numTimeSteps = parameterSet.get<int>("numTimeSteps");
// read solver settings
const double tolerance = parameterSet.get<double>("tolerance");
@@ -171,10 +174,18 @@ int main (int argc, char *argv[]) try
// Assembler using ADOL-C
typedef TargetSpace::rebind<adouble>::other ATargetSpace;
- std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > localEnergy;
- localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, ATargetSpace>);
- LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(localEnergy.get());
+ auto l2DistanceSquaredEnergy = std::make_shared<L2DistanceSquaredEnergy<FEBasis, ATargetSpace> >();
+
+ std::vector<std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > > addends(2);
+ addends[0] = std::make_shared<HarmonicEnergyLocalStiffness<FEBasis, ATargetSpace> >();
+ addends[1] = l2DistanceSquaredEnergy;
+
+ std::vector<double> weights = {1.0, 1.0/(2*timeStepSize)};
+
+ auto sumEnergy = std::make_shared< WeightedSumEnergy<FEBasis, ATargetSpace> >(addends, weights);
+
+ LocalGeodesicFEADOLCStiffness<FEBasis,TargetSpace> localGFEADOLCStiffness(sumEnergy.get());
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(feBasis, &localGFEADOLCStiffness);
@@ -198,30 +209,40 @@ int main (int argc, char *argv[]) try
false); // instrumentation
///////////////////////////////////////////////////////
- // Solve!
+ // Time loop
///////////////////////////////////////////////////////
- solver.setInitialIterate(x);
- solver.solve();
+ auto previousTimeStep = x;
- x = solver.getSol();
+ for (int i=0; i<numTimeSteps; i++)
+ {
+ auto previousTimeStepFct = std::make_shared<GlobalGeodesicFEFunction<FufemFEBasis,TargetSpace> >(fufemFeBasis,previousTimeStep);
+ l2DistanceSquaredEnergy->origin_ = previousTimeStepFct;
- ////////////////////////////////
- // Output result
- ////////////////////////////////
+ solver.setInitialIterate(x);
+ solver.solve();
- typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
- EmbeddedVectorType xEmbedded(x.size());
- for (size_t i=0; i<x.size(); i++)
- xEmbedded[i] = x[i].globalCoordinates();
+ x = solver.getSol();
+
+ previousTimeStep = x;
- auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,
- TypeTree::hybridTreePath(),
- xEmbedded);
+ ////////////////////////////////
+ // Output result
+ ////////////////////////////////
- SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),0);
- vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::scalar, xEmbedded[0].size()));
- vtkWriter.write("gradientflow_result");
+ typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
+ EmbeddedVectorType xEmbedded(x.size());
+ for (size_t i=0; i<x.size(); i++)
+ xEmbedded[i] = x[i].globalCoordinates();
+
+ auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,
+ TypeTree::hybridTreePath(),
+ xEmbedded);
+
+ SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),0);
+ vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::scalar, xEmbedded[0].size()));
+ vtkWriter.write("gradientflow_result_" + std::to_string(i));
+ }
return 0;
}