new application to measure the discretization error of geodesic finite elements

[[Imported from SVN: r4128]]

Makefile.am

@@ -4,11 +4,14 @@
 LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
 AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
 
-noinst_PROGRAMS = staticrod staticrod2 rod3d harmonicmaps dirneucoupling
+noinst_PROGRAMS = staticrod staticrod2 rod3d harmonicmaps dirneucoupling rod-eoc
 
 staticrod_SOURCES = staticrod.cc
 staticrod2_SOURCES = staticrod2.cc
 rod3d_SOURCES = rod3d.cc
+
+rod_eoc_SOURCES = rod-eoc.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) \

rod-eoc.cc

+#include <config.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/onedgrid.hh>
+
+#include <dune/istl/io.hh>
+
+
+#include <dune-solvers/solvers/iterativesolver.hh>
+#include <dune-solvers/norms/energynorm.hh>
+
+#include "src/rigidbodymotion.hh"
+#include "src/geodesicdifference.hh"
+#include "src/rotation.hh"
+#include "src/rodassembler.hh"
+#include "src/riemanniantrsolver.hh"
+#include "src/rodrefine.hh"
+
+typedef Dune::OneDGrid GridType;
+
+typedef RigidBodyMotion<3> TargetSpace;
+typedef std::vector<RigidBodyMotion<3> > SolutionType;
+
+const int blocksize = TargetSpace::TangentVector::size;
+
+using namespace Dune;
+using std::string;
+
+void solve (const GridType& grid,
+            SolutionType& x, int numLevels,
+            const TargetSpace& dirichletValue,
+            ConfigParser& parameters)
+{
+    // read solver setting
+    const double innerTolerance           = parameters.get<double>("innerTolerance");
+    const double tolerance                = parameters.get<double>("tolerance");
+    const int maxTrustRegionSteps         = parameters.get<int>("maxTrustRegionSteps");
+    const double initialTrustRegionRadius = parameters.get<double>("initialTrustRegionRadius");
+    const int multigridIterations         = parameters.get<int>("numIt");
+
+    // read rod parameter settings
+    const double A               = parameters.get<double>("A");
+    const double J1              = parameters.get<double>("J1");
+    const double J2              = parameters.get<double>("J2");
+    const double E               = parameters.get<double>("E");
+    const double nu              = parameters.get<double>("nu");
+
+    //   Create a local assembler
+    RodLocalStiffness<OneDGrid::LeafGridView,double> localStiffness(grid.leafView(),
+                                                                    A, J1, J2, E, nu);
+
+
+
+
+    x.resize(grid.size(1));
+
+    // //////////////////////////
+    //   Initial solution
+    // //////////////////////////
+
+    for (int i=0; i<x.size(); i++) {
+        x[i].r[0] = 0;
+        x[i].r[1] = 0;
+        x[i].r[2] = double(i)/(x.size()-1);
+        x[i].q    = Rotation<3,double>::identity();
+    }
+
+    //  set Dirichlet value
+    x.back() = dirichletValue;
+
+    // Both ends are Dirichlet
+    BitSetVector<blocksize> dirichletNodes(grid.size(1));
+    dirichletNodes.unsetAll();
+        
+    dirichletNodes[0] = dirichletNodes.back() = true;
+    
+    // ///////////////////////////////////////////
+    //   Create a solver for the rod problem
+    // ///////////////////////////////////////////
+
+    RodAssembler<GridType> rodAssembler(grid, &localStiffness);
+
+    RiemannianTrustRegionSolver<GridType,RigidBodyMotion<3> > rodSolver;
+#if 0
+    rodSolver.setup(grid, 
+                    &rodAssembler,
+                    x,
+                    dirichletNodes,
+                    tolerance,
+                    maxTrustRegionSteps,
+                    initialTrustRegionRadius,
+                    multigridIterations,
+                    innerTolerance,
+                    mu, nu1, nu2,
+                    baseIterations,
+                    baseTolerance,
+                    false);  // instrumentation
+#else
+    rodSolver.setupTCG(grid, 
+                       &rodAssembler,
+                       x,
+                       dirichletNodes,
+                       tolerance,
+                       maxTrustRegionSteps,
+                       initialTrustRegionRadius,
+                       multigridIterations,
+                       innerTolerance,
+                       false);  // instrumentation
+#endif
+
+    // /////////////////////////////////////////////////////
+    //   Solve!
+    // /////////////////////////////////////////////////////
+
+    rodSolver.setInitialSolution(x);
+    rodSolver.solve();
+
+    x = rodSolver.getSol();
+}
+
+int main (int argc, char *argv[]) try
+{
+    // parse data file
+    ConfigParser parameterSet;
+    if (argc==2)
+        parameterSet.parseFile(argv[1]);
+    else
+        parameterSet.parseFile("rod-eoc.parset");
+
+    // read solver settings
+    const int numLevels        = parameterSet.get<int>("numLevels");
+    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 baseTolerance    = parameterSet.get<double>("baseTolerance");
+
+    const int numRodBaseElements = parameterSet.get<int>("numRodBaseElements");
+    
+    // /////////////////////////////////////////
+    //   Read Dirichlet values
+    // /////////////////////////////////////////
+
+    RigidBodyMotion<3> dirichletValue;
+    dirichletValue.r[0] = parameterSet.get<double>("dirichletValueX");
+    dirichletValue.r[1] = parameterSet.get<double>("dirichletValueY");
+    dirichletValue.r[2] = parameterSet.get<double>("dirichletValueZ");
+
+    FieldVector<double,3> axis;
+    axis[0] = parameterSet.get<double>("dirichletAxisX");
+    axis[1] = parameterSet.get<double>("dirichletAxisY");
+    axis[2] = parameterSet.get<double>("dirichletAxisZ");
+    double angle = parameterSet.get<double>("dirichletAngle");
+
+    dirichletValue.q = Rotation<3,double>(axis, M_PI*angle/180);
+
+    // ///////////////////////////////////////////////////////////
+    //   First compute the 'exact' solution on a very fine grid
+    // ///////////////////////////////////////////////////////////
+
+    //    Create the reference grid
+    
+    GridType referenceGrid(numRodBaseElements, 0, 1);
+    
+    referenceGrid.globalRefine(numLevels-1);
+
+    //  Solve the rod Dirichlet problem
+    SolutionType referenceSolution;
+    solve(referenceGrid, referenceSolution, numLevels, dirichletValue, parameterSet);
+
+
+
+    // ///////////////////////////////////////////////////////////
+    //   Compute on all coarser levels, and compare
+    // ///////////////////////////////////////////////////////////
+
+    
+    
+    for (int i=1; i<=numLevels; i++) {
+
+        GridType grid(numRodBaseElements, 0, 1);
+        grid.globalRefine(i-1);
+
+        // compute again
+        SolutionType solution;
+        solve(grid, solution, i, dirichletValue, parameterSet);
+
+        // Prolong solution to the very finest grid
+        for (int j=i; j<numLevels; j++)
+            globalRodRefine(grid, solution);
+
+        assert(referenceSolution.size() == solution.size());
+
+        // Compute max-norm difference
+        std::cout << "Level: " << i-1 
+                  << ",   max-norm error: " << computeGeodesicDifference(solution,referenceSolution).infinity_norm()
+                  << std::endl;
+
+    }    
+        
+ } catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }