Prototype: compute harmonic maps for functions to a nonflat manifold

[[Imported from SVN: r4025]]

Makefile.am

@@ -4,11 +4,15 @@
 LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
 AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
 
-noinst_PROGRAMS = staticrod staticrod2 rod3d dirneucoupling
+noinst_PROGRAMS = staticrod staticrod2 rod3d harmonicmaps dirneucoupling
 
 staticrod_SOURCES = staticrod.cc
 staticrod2_SOURCES = staticrod2.cc
 rod3d_SOURCES = rod3d.cc
+harmonicmaps_SOURCES = harmonicmaps.cc
+harmonicmaps_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
+harmonicmaps_LDADD    = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
+                          $(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
 
 dirneucoupling_SOURCES  = dirneucoupling.cc
 dirneucoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)

harmonicmaps.cc

+#include <config.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/amirameshreader.hh>
+
+#include <dune-solvers/solvers/iterativesolver.hh>
+#include <dune-solvers/norms/energynorm.hh>
+
+#include "src/roddifference.hh"
+#include "src/rodwriter.hh"
+#include "src/rotation.hh"
+#include "src/geodesicfeassembler.hh"
+#include "src/rodsolver.hh"
+
+// grid dimension
+const int dim = 3;
+
+// Image space of the geodesic fe functions
+typedef Rotation<3,double> TargetSpace;
+
+// Tangent vector of the image space
+const int blocksize = TargetSpace::TangentVector::size;
+
+using namespace Dune;
+
+int main (int argc, char *argv[]) try
+{
+    typedef std::vector<TargetSpace> SolutionType;
+
+    // parse data file
+    ConfigParser parameterSet;
+    if (argc==2)
+        parameterSet.parseFile(argv[1]);
+    else
+        parameterSet.parseFile("harmonicmaps.parset");
+
+    // read solver settings
+    const int numLevels                   = parameterSet.get<int>("numLevels");
+    const double tolerance                = parameterSet.get<double>("tolerance");
+    const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
+    const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+    const int multigridIterations         = parameterSet.get<int>("numIt");
+    const int nu1                         = parameterSet.get<int>("nu1");
+    const int nu2                         = parameterSet.get<int>("nu2");
+    const int mu                          = parameterSet.get<int>("mu");
+    const int baseIterations              = parameterSet.get<int>("baseIt");
+    const double mgTolerance              = parameterSet.get<double>("mgTolerance");
+    const double baseTolerance            = parameterSet.get<double>("baseTolerance");
+    const bool instrumented               = parameterSet.get<bool>("instrumented");
+    std::string resultPath                = parameterSet.get("resultPath", "");
+
+    // read problem settings
+    std::string path                = parameterSet.get<std::string>("path");
+    std::string gridFile            = parameterSet.get<std::string>("gridFile");
+
+    // ///////////////////////////////////////
+    //    Create the grid
+    // ///////////////////////////////////////
+    typedef UGGrid<dim> GridType;
+    GridType grid;
+    grid.setRefinementType(GridType::COPY);
+    AmiraMeshReader<GridType>::read(grid, path + gridFile);
+
+    grid.globalRefine(numLevels-1);
+
+    SolutionType x(grid.size(dim));
+
+    // //////////////////////////
+    //   Initial solution
+    // //////////////////////////
+
+    for (int i=0; i<x.size(); i++) {
+        x[i] = Rotation<3,double>::identity();
+    }
+
+    // backup for error measurement later
+    SolutionType initialIterate = x;
+
+    // /////////////////////////////////////////
+    //   Read Dirichlet values
+    // /////////////////////////////////////////
+
+    BitSetVector<1> allNodes(grid.size(dim));
+    allNodes.setAll();
+    LeafBoundaryPatch<GridType> dirichletBoundary(grid, allNodes);
+
+    BitSetVector<blocksize> dirichletNodes(grid.size(1));
+    for (int i=0; i<dirichletNodes.size(); i++)
+        dirichletNodes[i] = dirichletBoundary.containsVertex(i);
+    
+    // ////////////////////////////////////////////////////////////
+    //   Create an assembler for the Harmonic Energy Functional
+    // ////////////////////////////////////////////////////////////
+    GeodesicFEAssembler<GridType::LeafGridView,TargetSpace> assembler(grid.leafView());
+
+    // /////////////////////////////////////////////////
+    //   Create a Riemannian trust-region solver
+    // /////////////////////////////////////////////////
+#if 0
+    RodSolver<GridType> rodSolver;
+    rodSolver.setup(grid, 
+                    &assembler,
+                    x,
+                    dirichletNodes,
+                    tolerance,
+                    maxTrustRegionSteps,
+                    initialTrustRegionRadius,
+                    multigridIterations,
+                    mgTolerance,
+                    mu, nu1, nu2,
+                    baseIterations,
+                    baseTolerance,
+                    instrumented);
+
+    // /////////////////////////////////////////////////////
+    //   Solve!
+    // /////////////////////////////////////////////////////
+
+    std::cout << "Energy: " << assembler.computeEnergy(x) << std::endl;
+    
+    rodSolver.setInitialSolution(x);
+    rodSolver.solve();
+
+    x = rodSolver.getSol();
+#endif
+    // //////////////////////////////
+    //   Output result
+    // //////////////////////////////
+#if 0
+    writeRod(x, resultPath + "rod3d.result");
+#endif
+
+#if 0
+    // //////////////////////////////////////////////////////////
+    //   Recompute and compare against exact solution
+    // //////////////////////////////////////////////////////////
+    
+    SolutionType exactSolution = x;
+
+    // //////////////////////////////////////////////////////////
+    //   Compute hessian of the rod functional at the exact solution
+    //   for use of the energy norm it creates.
+    // //////////////////////////////////////////////////////////
+
+    BCRSMatrix<FieldMatrix<double, 6, 6> > hessian;
+    MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
+    rodAssembler.getNeighborsPerVertex(indices);
+    indices.exportIdx(hessian);
+    rodAssembler.assembleMatrixFD(exactSolution, hessian);
+
+
+    double error = std::numeric_limits<double>::max();
+
+    SolutionType intermediateSolution(x.size());
+
+    // Create statistics file
+    std::ofstream statisticsFile((resultPath + "trStatistics").c_str());
+
+    // Compute error of the initial iterate
+    typedef BlockVector<FieldVector<double,6> > RodDifferenceType;
+    RodDifferenceType rodDifference = computeRodDifference(exactSolution, initialIterate);
+    double oldError = std::sqrt(computeEnergyNormSquared(rodDifference, hessian));
+
+    int i;
+    for (i=0; i<maxTrustRegionSteps; i++) {
+        
+        // /////////////////////////////////////////////////////
+        //   Read iteration from file
+        // /////////////////////////////////////////////////////
+        char iSolFilename[100];
+        sprintf(iSolFilename, "tmp/intermediateSolution_%04d", i);
+            
+        FILE* fp = fopen(iSolFilename, "rb");
+        if (!fp)
+            DUNE_THROW(IOError, "Couldn't open intermediate solution '" << iSolFilename << "'");
+        for (int j=0; j<intermediateSolution.size(); j++) {
+            fread(&intermediateSolution[j].r, sizeof(double), 3, fp);
+            fread(&intermediateSolution[j].q, sizeof(double), 4, fp);
+        }
+        
+        fclose(fp);
+
+        // /////////////////////////////////////////////////////
+        //   Compute error
+        // /////////////////////////////////////////////////////
+
+        rodDifference = computeRodDifference(exactSolution, intermediateSolution);
+        
+        error = std::sqrt(computeEnergyNormSquared(rodDifference, hessian));
+        
+
+        double convRate = error / oldError;
+
+        // Output
+        std::cout << "Trust-region iteration: " << i << "  error : " << error << ",      "
+                  << "convrate " << convRate << std::endl;
+        statisticsFile << i << "  " << error << "  " << convRate << std::endl;
+
+        if (error < 1e-12)
+          break;
+
+        oldError = error;
+        
+    }            
+#endif
+
+    // //////////////////////////////
+ } catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }