From bcce36a11ade75ad742b023e8b11025e78dd5622 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 14 Jan 2015 21:27:09 +0000
Subject: [PATCH] Now that we use AD, remove the implementation of the gradient
 of the harmonic energy

... so much difficult code obsolete now... AD rocks...

[[Imported from SVN: r10010]]
---
 dune/gfe/harmonicenergystiffness.hh | 122 ----------------------------
 src/harmonicmaps.cc                 |   2 -
 2 files changed, 124 deletions(-)

diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index 3a41fb1d..bd913400 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -7,10 +7,6 @@
 #include "localgeodesicfestiffness.hh"
 #include "localgeodesicfefunction.hh"
 
-#ifdef HARMONIC_ENERGY_FD_GRADIENT
-#warning Finite-difference approximation of the energy gradient
-#endif
-
 template<class GridView, class LocalFiniteElement, class TargetSpace>
 class HarmonicEnergyLocalStiffness
     : public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,TargetSpace>
@@ -33,20 +29,6 @@ public:
                const LocalFiniteElement& localFiniteElement,
                const std::vector<TargetSpace>& localSolution) const;
 
-#ifndef HARMONIC_ENERGY_FD_GRADIENT
-    // The finite difference gradient method is in the base class.
-    // If the cpp macro is not set we overload it here.
-    /** \brief Assemble the gradient of the energy functional on one element */
-    virtual void assembleEmbeddedGradient(const Entity& element,
-                                          const LocalFiniteElement& localFiniteElement,
-                                          const std::vector<TargetSpace>& solution,
-                                  std::vector<typename TargetSpace::EmbeddedTangentVector>& gradient) const;
-
-    virtual void assembleGradient(const Entity& element,
-                                  const LocalFiniteElement& localFiniteElement,
-                                  const std::vector<TargetSpace>& localSolution,
-                                  std::vector<typename TargetSpace::TangentVector>& localGradient) const;
-#endif
 };
 
 template <class GridView, class LocalFiniteElement, class TargetSpace>
@@ -101,109 +83,5 @@ energy(const Entity& element,
     return 0.5 * energy;
 }
 
-#ifndef HARMONIC_ENERGY_FD_GRADIENT
-template <class GridView, class LocalFiniteElement, class TargetSpace>
-void HarmonicEnergyLocalStiffness<GridView, LocalFiniteElement, TargetSpace>::
-assembleEmbeddedGradient(const Entity& element,
-                         const LocalFiniteElement& localFiniteElement,
-                         const std::vector<TargetSpace>& localSolution,
-                         std::vector<typename TargetSpace::EmbeddedTangentVector>& localGradient) const
-{
-    typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
-    // initialize gradient
-    localGradient.resize(localSolution.size());
-    std::fill(localGradient.begin(), localGradient.end(), typename TargetSpace::EmbeddedTangentVector(0));
-
-    // Set up local gfe function from the  local coefficients
-    typedef LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> LocalGFEFunctionType;
-    LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
-
-    // I am not sure about the correct quadrature order
-    int quadOrder = (element.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
-                                                 : (localFiniteElement.localBasis().order()-1) * 2 * gridDim;
-
-    // numerical quadrature loop
-    const Dune::QuadratureRule<double, gridDim>& quad
-        = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
-
-    for (size_t pt=0; pt<quad.size(); pt++) {
-
-        // Local position of the quadrature point
-        const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
-
-        const double integrationElement = element.geometry().integrationElement(quadPos);
-
-        const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-        double weight = quad[pt].weight() * integrationElement;
-
-        // The derivative of the local function defined on the reference element
-        typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
-
-        // The derivative of the function defined on the actual element
-        typename LocalGFEFunctionType::DerivativeType derivative;
-
-        for (size_t comp=0; comp<referenceDerivative.N(); comp++)
-            jacobianInverseTransposed.mv(referenceDerivative[comp], derivative[comp]);
-
-        // loop over all the element's degrees of freedom and compute the gradient wrt it
-        for (size_t i=0; i<localSolution.size(); i++) {
-
-            typename LocalGFEFunctionType::DerivativeOfGradientWRTCoefficientType referenceDerivativeDerivative;
-#ifdef HARMONIC_ENERGY_FD_INNER_GRADIENT
-#warning Using finite differences to compute the inner gradients!
-            localGeodesicFEFunction.evaluateFDDerivativeOfGradientWRTCoefficient(quadPos, i, referenceDerivativeDerivative);
-#else
-            localGeodesicFEFunction.evaluateDerivativeOfGradientWRTCoefficient(quadPos, i, referenceDerivativeDerivative);
-#endif
-
-            // multiply the transformation from the reference element to the actual element
-            typename LocalGFEFunctionType::DerivativeOfGradientWRTCoefficientType derivativeDerivative;
-            for (int ii=0; ii<TargetSpace::EmbeddedTangentVector::dimension; ii++)
-                for (int jj=0; jj<TargetSpace::EmbeddedTangentVector::dimension; jj++)
-                    for (int kk=0; kk<gridDim; kk++) {
-                        derivativeDerivative[ii][jj][kk] = 0;
-                        for (int ll=0; ll<gridDim; ll++)
-                            derivativeDerivative[ii][jj][kk] += referenceDerivativeDerivative[ii][jj][ll] * jacobianInverseTransposed[kk][ll];
-                    }
-
-            for (int j=0; j<derivative.rows; j++) {
-
-                for (int k=0; k<derivative.cols; k++) {
-
-                    for (int l=0; l<TargetSpace::EmbeddedTangentVector::dimension; l++)
-                        localGradient[i][l] += weight*derivative[j][k] * derivativeDerivative[l][j][k];
-
-                }
-
-            }
-
-        }
-
-    }
-
-}
-
-template <class GridView, class LocalFiniteElement, class TargetSpace>
-void HarmonicEnergyLocalStiffness<GridView, LocalFiniteElement, TargetSpace>::
-assembleGradient(const Entity& element,
-                 const LocalFiniteElement& localFiniteElement,
-                 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, localFiniteElement, 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]);
-
-}
-#endif
 #endif
 
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 02cebb22..64000706 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -2,8 +2,6 @@
 
 #include <fenv.h>
 
-//#define HARMONIC_ENERGY_FD_GRADIENT
-
 //#define UNITVECTOR2
 #define UNITVECTOR3
 //#define UNITVECTOR4
-- 
GitLab