new application to test the convergence order of geodesic finite elements on a...

new application to test the convergence order of geodesic finite elements on a harmonic maps problem

[[Imported from SVN: r5764]]

Makefile.am

@@ -9,7 +9,7 @@ SUBDIRS = m4 src test
 LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
 AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
 
-noinst_PROGRAMS = rodobstacle rod3d harmonicmaps dirneucoupling rod-eoc 
+noinst_PROGRAMS = rodobstacle rod3d harmonicmaps harmonicmaps-eoc dirneucoupling rod-eoc 
 
 rodobstacle_SOURCES = rodobstacle.cc
 rod3d_SOURCES = rod3d.cc
@@ -21,6 +21,11 @@ harmonicmaps_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $
 harmonicmaps_LDADD    = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
                           $(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
 
+harmonicmaps_eoc_SOURCES = harmonicmaps-eoc.cc
+harmonicmaps_eoc_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
+harmonicmaps_eoc_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)
 dirneucoupling_LDADD    = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \

harmonicmaps-eoc.cc

+#include <config.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/onedgrid.hh>
+#include <dune/grid/../../doc/grids/gridfactory/structuredgridfactory.hh>
+
+#include <dune/ag-common/functionspacebases/p1nodalbasis.hh>
+#include <dune/ag-common/assemblers/operatorassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/massassembler.hh>
+
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+#include "src/unitvector.hh"
+#include "src/harmonicenergystiffness.hh"
+#include "src/geodesicfeassembler.hh"
+#include "src/riemanniantrsolver.hh"
+#include "src/rodrefine.hh"
+#include "src/rodwriter.hh"
+
+// grid dimension
+const int dim = 2;
+
+typedef UnitVector<3> TargetSpace;
+typedef std::vector<TargetSpace> SolutionType;
+
+const int blocksize = TargetSpace::TangentVector::size;
+
+using namespace Dune;
+using std::string;
+
+template <class GridType>
+void solve (const shared_ptr<GridType>& grid,
+            SolutionType& x, 
+            int numLevels,
+            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 Dirichlet values
+    // /////////////////////////////////////////
+
+    BitSetVector<1> allNodes(grid->size(dim));
+    allNodes.setAll();
+    LeafBoundaryPatch<GridType> dirichletBoundary(*grid, allNodes);
+
+    BitSetVector<blocksize> dirichletNodes(grid->size(dim));
+    for (int i=0; i<dirichletNodes.size(); i++)
+        dirichletNodes[i] = dirichletBoundary.containsVertex(i);
+
+    // //////////////////////////
+    //   Initial solution
+    // //////////////////////////
+
+    x.resize(grid->size(dim));
+
+        FieldVector<double,3> yAxis(0);
+    yAxis[1] = 1;
+
+    typename GridType::LeafGridView::template Codim<dim>::Iterator vIt    = grid->template leafbegin<dim>();
+    typename GridType::LeafGridView::template Codim<dim>::Iterator vEndIt = grid->template leafend<dim>();
+
+    for (; vIt!=vEndIt; ++vIt) {
+        int idx = grid->leafIndexSet().index(*vIt);
+
+        FieldVector<double,3> v;
+        FieldVector<double,2> pos = vIt->geometry().corner(0);
+        FieldVector<double,3> axis;
+        axis[0] = pos[0];  axis[1] = pos[1]; axis[2] = 1;
+        Rotation<3,double> rotation(axis, pos.two_norm()*M_PI*1.5);
+
+        if (dirichletNodes[idx][0]) {
+//             FieldMatrix<double,3,3> rMat;
+//             rotation.matrix(rMat);
+//             v = rMat[2];
+            v[0] = std::sin(pos[0]*M_PI);
+            v[1] = 0;
+            v[2] = std::cos(pos[0]*M_PI);
+        } else {
+            v[0] = 1;
+            v[1] = 0;
+            v[2] = 0;
+        }            
+
+        x[idx] = v;
+
+    }
+
+
+    // ////////////////////////////////////////////////////////////
+    //   Create an assembler for the Harmonic Energy Functional
+    // ////////////////////////////////////////////////////////////
+
+    HarmonicEnergyLocalStiffness<typename GridType::LeafGridView,TargetSpace> harmonicEnergyLocalStiffness;
+
+    GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace> assembler(grid->leafView(),
+                                                                      &harmonicEnergyLocalStiffness);
+    
+    // ///////////////////////////////////////////
+    //   Create a solver for the rod problem
+    // ///////////////////////////////////////////
+
+    RiemannianTrustRegionSolver<GridType,TargetSpace> solver;
+
+    solver.setup(*grid, 
+                 &assembler,
+                 x,
+                 dirichletNodes,
+                 tolerance,
+                 maxTrustRegionSteps,
+                 initialTrustRegionRadius,
+                 multigridIterations,
+                 innerTolerance,
+                 1, 3, 3,
+                 100,     // iterations of the base solver
+                 1e-8,    // base tolerance
+                 false);  // instrumentation
+
+    // /////////////////////////////////////////////////////
+    //   Solve!
+    // /////////////////////////////////////////////////////
+
+    solver.setInitialSolution(x);
+    solver.solve();
+
+    x = solver.getSol();
+}
+
+int main (int argc, char *argv[]) try
+{
+    // parse data file
+    ConfigParser parameterSet;
+    if (argc==2)
+        parameterSet.parseFile(argv[1]);
+    else
+        parameterSet.parseFile("harmonicmaps-eoc.parset");
+
+    // read solver settings
+    const int numLevels        = parameterSet.get<int>("numLevels");
+    const int baseIterations      = parameterSet.get<int>("baseIt");
+    const double baseTolerance    = parameterSet.get<double>("baseTolerance");
+
+    const int numBaseElements = parameterSet.get<int>("numBaseElements");
+    
+    // /////////////////////////////////////////
+    //   Read Dirichlet values
+    // /////////////////////////////////////////
+
+    
+
+
+    // ///////////////////////////////////////////////////////////
+    //   First compute the 'exact' solution on a very fine grid
+    // ///////////////////////////////////////////////////////////
+
+    typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
+
+    //    Create the reference grid
+
+    array<unsigned int,dim> elements;
+    elements.fill(4);
+    shared_ptr<GridType> referenceGrid = StructuredGridFactory<GridType>::createSimplexGrid(FieldVector<double,dim>(0),
+                                                                                            FieldVector<double,dim>(1),
+                                                                                            elements);
+    referenceGrid->globalRefine(numLevels-1);
+
+    //  Solve the rod Dirichlet problem
+    SolutionType referenceSolution;
+    solve(referenceGrid, referenceSolution, numLevels, parameterSet);
+
+
+    // //////////////////////////////////////////////////////////////////////
+    //   Compute mass matrix and laplace matrix to emulate L2 and H1 norms
+    // //////////////////////////////////////////////////////////////////////
+#if 0
+    typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
+    FEBasis basis(referenceGrid->leafView());
+    OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
+
+    LaplaceAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> laplaceLocalAssembler;
+    MassAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> massMatrixLocalAssembler;
+
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+    ScalarMatrixType laplace, massMatrix;
+
+    operatorAssembler.assemble(laplaceLocalAssembler, laplace);
+    operatorAssembler.assemble(massMatrixLocalAssembler, massMatrix);
+#endif
+    // ///////////////////////////////////////////////////////////
+    //   Compute on all coarser levels, and compare
+    // ///////////////////////////////////////////////////////////
+    
+    for (int i=1; i<=numLevels; i++) {
+
+        array<unsigned int,dim> elements;
+        elements.fill(numBaseElements);
+        shared_ptr<GridType> grid = StructuredGridFactory<GridType>::createSimplexGrid(FieldVector<double,dim>(0),
+                                                                                       FieldVector<double,dim>(1),
+                                                                                       elements);
+
+        grid->globalRefine(i-1);
+
+        // compute again
+        SolutionType solution;
+        solve(grid, solution, i, parameterSet);
+
+#if 0
+        // Prolong solution to the very finest grid
+        for (int j=i; j<numLevels; j++)
+            globalRodRefine(grid, solution);
+
+        std::stringstream numberAsAscii;
+        numberAsAscii << i;
+        writeRod(solution, "rodGrid_" + numberAsAscii.str());
+
+        assert(referenceSolution.size() == solution.size());
+#endif
+#if 0
+        BlockVector<TargetSpace::TangentVector> difference = computeGeodesicDifference(solution,referenceSolution);
+
+        H1SemiNorm< BlockVector<TargetSpace::TangentVector> > h1Norm(laplace);
+        H1SemiNorm< BlockVector<TargetSpace::TangentVector> > l2Norm(massMatrix);
+
+        // Compute max-norm difference
+        std::cout << "Level: " << i-1 
+                  << ",   max-norm error: " << difference.infinity_norm()
+                  << std::endl;
+
+        std::cout << "Level: " << i-1 
+                  << ",   L2 error: " << l2Norm(difference)
+                  << std::endl;
+
+        std::cout << "Level: " << i-1 
+                  << ",   H1 error: " << h1Norm(difference)
+                  << std::endl;
+#endif
+    }    
+        
+ } catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+ }

harmonicmaps-eoc.parset

+# Number of grid levels
+numLevels = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-12
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+# Initial trust-region radius
+initialTrustRegionRadius = 1
+
+# Number of multigrid iterations per trust-region step
+numIt = 200
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the inner solver
+innerTolerance = 1e-5
+
+# Tolerance of the base grid solver
+baseTolerance = -1
+
+############################
+#   Problem specifications
+############################
+
+# Number of elements of the rod grid
+numBaseElements = 4
+