Compute the local corrections in local coordinates on the tangent bundle now,...

Compute the local corrections in local coordinates on the tangent bundle now, instead of using the embedding in the surrounding space

[[Imported from SVN: r5985]]

dune/gfe/geodesicfeassembler.hh

@@ -22,7 +22,7 @@ class GeodesicFEAssembler {
     enum { gridDim = GridView::dimension };
     
     //! Dimension of a tangent space
-    enum { blocksize = TargetSpace::EmbeddedTangentVector::size };
+    enum { blocksize = TargetSpace::TangentVector::size };
     
     //!
     typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;

dune/gfe/localgeodesicfestiffness.hh

@@ -316,7 +316,7 @@ class LocalGeodesicFEStiffness <GridView,UnitVector<dim> >
 public:
     
     //! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
-    enum { blocksize = TargetSpace::EmbeddedTangentVector::size };
+    enum { blocksize = TargetSpace::TangentVector::size };
 
     /** \brief Assemble the local stiffness matrix at the current position
 
@@ -330,10 +330,17 @@ public:
 
     /** \brief Assemble the element gradient of the energy functional 
 
+    The default implementation in this class uses a finite difference approximation */
+    virtual void assembleEmbeddedGradient(const Entity& element,
+                                          const std::vector<TargetSpace>& solution,
+                                          std::vector<typename TargetSpace::EmbeddedTangentVector>& gradient) const;
+
+    /** \brief Assemble the element gradient of the energy functional 
+
     The default implementation in this class uses a finite difference approximation */
     virtual void assembleGradient(const Entity& element,
                                   const std::vector<TargetSpace>& solution,
-                                  std::vector<Dune::FieldVector<double,blocksize> >& gradient) const;
+                                  std::vector<typename TargetSpace::TangentVector>& gradient) const;
 
     // assembled data
     Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > A_;
@@ -342,11 +349,13 @@ public:
 
 template <class GridView, int dim>
 void LocalGeodesicFEStiffness<GridView, UnitVector<dim> >::
-assembleGradient(const Entity& element,
-                 const std::vector<TargetSpace>& localSolution,
-                 std::vector<Dune::FieldVector<double,blocksize> >& localGradient) const
+assembleEmbeddedGradient(const Entity& element,
+                         const std::vector<TargetSpace>& localSolution,
+                         std::vector<typename TargetSpace::EmbeddedTangentVector>& localGradient) const
 {
 
+    const int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::size;
+
     // ///////////////////////////////////////////////////////////
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
@@ -360,7 +369,7 @@ assembleGradient(const Entity& element,
 
     for (size_t i=0; i<localSolution.size(); i++) {
         
-        for (int j=0; j<blocksize; j++) {
+        for (int j=0; j<embeddedBlocksize; j++) {
             
             // The return value does not have unit norm.  But assigning it to a UnitVector object
             // will normalize it.  This amounts to an extension of the energy functional 
@@ -382,6 +391,24 @@ assembleGradient(const Entity& element,
 
 }
 
+template <class GridView, int dim>
+void LocalGeodesicFEStiffness<GridView, UnitVector<dim> >::
+assembleGradient(const Entity& element,
+                 const std::vector<TargetSpace>& localSolution,
+                 std::vector<typename TargetSpace::TangentVector>& localGradient) const
+{
+    std::vector<typename TargetSpace::EmbeddedTangentVector> embeddedLocalGradient;
+
+    // first compute the gradient in embedded coordinates
+    assembleEmbeddedGradient(element, localSolution, embeddedLocalGradient);
+
+    // transform to coordinates on the tangent space
+    localGradient.resize(embeddedLocalGradient.size());
+
+    for (size_t i=0; i<localGradient.size(); i++)
+        localSolution[i].orthonormalFrame().mv(embeddedLocalGradient[i], localGradient[i]);
+
+}
 
 template <class GridType, int dim>
 void LocalGeodesicFEStiffness<GridType,UnitVector<dim> >::

dune/gfe/riemanniantrsolver.hh

@@ -18,9 +18,9 @@
 template <class GridType, class TargetSpace>
 class RiemannianTrustRegionSolver 
     : public IterativeSolver<std::vector<TargetSpace>,
-                             Dune::BitSetVector<TargetSpace::EmbeddedTangentVector::size> >
+                             Dune::BitSetVector<TargetSpace::TangentVector::size> >
 { 
-    const static int blocksize = TargetSpace::EmbeddedTangentVector::size;
+    const static int blocksize = TargetSpace::TangentVector::size;
 
     const static int gridDim = GridType::dimension;
 

harmonicmaps.cc

@@ -45,7 +45,7 @@ typedef RealTuple<1> TargetSpace;
 #endif
 
 // Tangent vector of the image space
-const int blocksize = TargetSpace::EmbeddedTangentVector::size;
+const int blocksize = TargetSpace::TangentVector::size;
 
 using namespace Dune;
 
@@ -254,6 +254,8 @@ int main (int argc, char *argv[]) try
     writeRod(x, resultPath + "rod3d.result");
 #endif
 
+    exit(0);
+
     // //////////////////////////////////////////////////////////
     //   Recompute and compare against exact solution
     // //////////////////////////////////////////////////////////