Do not hardwire the local interpolation rule in the energy implementations

Make them a parameter instead.  This allows to switch more easily between
geodesic finite elements and projected finite elements.

dune/gfe/chiralskyrmionenergy.hh

@@ -5,7 +5,6 @@
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/gfe/localgeodesicfestiffness.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
 
 namespace Dune {
 
@@ -16,7 +15,7 @@ namespace GFE {
  * The energy is discussed in:
  * - Christof Melcher, "Chiral skyrmions in the plane", Proc. of the Royal Society, online DOI DOI: 10.1098/rspa.2014.0394
  */
-template<class Basis, class field_type>
+template<class Basis, class LocalInterpolationRule, class field_type>
 class ChiralSkyrmionEnergy
 : public LocalGeodesicFEStiffness<Basis,UnitVector<field_type,3> >
 {
@@ -49,9 +48,9 @@ public:
   field_type kappa_;
 };
 
-template <class Basis, class field_type>
-typename ChiralSkyrmionEnergy<Basis, field_type>::RT
-ChiralSkyrmionEnergy<Basis, field_type>::
+template <class Basis, class LocalInterpolationRule, class field_type>
+typename ChiralSkyrmionEnergy<Basis, LocalInterpolationRule, field_type>::RT
+ChiralSkyrmionEnergy<Basis, LocalInterpolationRule, field_type>::
 energy(const typename Basis::LocalView& localView,
        const std::vector<TargetSpace>& localConfiguration) const
 {
@@ -61,8 +60,7 @@ energy(const typename Basis::LocalView& localView,
 
   const auto element = localView.element();
   const auto& localFiniteElement = localView.tree().finiteElement();
-  typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
-  LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localConfiguration);
+  LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
 
   int quadOrder = (element.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
                                                : localFiniteElement.localBasis().order() * 2 * gridDim;
@@ -84,13 +82,13 @@ energy(const typename Basis::LocalView& localView,
     double weight = quad[pt].weight() * integrationElement;
 
     // The value of the function
-    auto value = localGeodesicFEFunction.evaluate(quadPos);
+    auto value = localInterpolationRule.evaluate(quadPos);
 
     // The derivative of the local function defined on the reference element
-    typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
+    typename LocalInterpolationRule::DerivativeType referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos,value);
 
     // The derivative of the function defined on the actual element
-    typename LocalGFEFunctionType::DerivativeType derivative(0);
+    typename LocalInterpolationRule::DerivativeType derivative(0);
 
     for (size_t comp=0; comp<referenceDerivative.N(); comp++)
       jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);

dune/gfe/harmonicenergystiffness.hh

@@ -5,9 +5,8 @@
 #include <dune/geometry/quadraturerules.hh>
 
 #include "localgeodesicfestiffness.hh"
-#include "localgeodesicfefunction.hh"
 
-template<class Basis, class TargetSpace>
+template<class Basis, class LocalInterpolationRule, class TargetSpace>
 class HarmonicEnergyLocalStiffness
     : public LocalGeodesicFEStiffness<Basis,TargetSpace>
 {
@@ -27,17 +26,16 @@ public:
 
 };
 
-template <class Basis, class TargetSpace>
-typename HarmonicEnergyLocalStiffness<Basis, TargetSpace>::RT
-HarmonicEnergyLocalStiffness<Basis, TargetSpace>::
+template <class Basis, class LocalInterpolationRule, class TargetSpace>
+typename HarmonicEnergyLocalStiffness<Basis, LocalInterpolationRule, TargetSpace>::RT
+HarmonicEnergyLocalStiffness<Basis, LocalInterpolationRule, TargetSpace>::
 energy(const typename Basis::LocalView& localView,
        const std::vector<TargetSpace>& localSolution) const
 {
     RT energy = 0;
 
     const auto& localFiniteElement = localView.tree().finiteElement();
-    typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
-    LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+    LocalInterpolationRule localInterpolationRule(localFiniteElement,localSolution);
 
     int quadOrder = (localFiniteElement.type().isSimplex()) ? (localFiniteElement.localBasis().order()-1) * 2
                                                  : localFiniteElement.localBasis().order() * 2 * gridDim;
@@ -58,10 +56,10 @@ energy(const typename Basis::LocalView& localView,
         double weight = quad[pt].weight() * integrationElement;
 
         // The derivative of the local function defined on the reference element
-        auto referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
+        auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos);
 
         // The derivative of the function defined on the actual element
-        typename LocalGFEFunctionType::DerivativeType derivative(0);
+        typename LocalInterpolationRule::DerivativeType derivative(0);
 
         for (size_t comp=0; comp<referenceDerivative.N(); comp++)
             jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);

src/harmonicmaps.cc

@@ -35,6 +35,8 @@
 #include <dune/gfe/rotation.hh>
 #include <dune/gfe/unitvector.hh>
 #include <dune/gfe/realtuple.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
 #include <dune/gfe/localgeodesicfeadolcstiffness.hh>
 #include <dune/gfe/harmonicenergystiffness.hh>
 #include <dune/gfe/chiralskyrmionenergy.hh>
@@ -189,20 +191,23 @@ int main (int argc, char *argv[]) try
     //   Create an assembler for the Harmonic Energy Functional
     // ////////////////////////////////////////////////////////////
 
-    // Assembler using ADOL-C
     typedef TargetSpace::rebind<adouble>::other ATargetSpace;
+    typedef LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace> LocalInterpolationRule;
+    //typedef GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace> LocalInterpolationRule;
+
+    // Assembler using ADOL-C
     std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > localEnergy;
 
     std::string energy = parameterSet.get<std::string>("energy");
     if (energy == "harmonic")
     {
 
-      localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, ATargetSpace>);
+      localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, LocalInterpolationRule, ATargetSpace>);
 
     } else if (energy == "chiral_skyrmion")
     {
 
-      localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, adouble>(parameterSet.sub("energyParameters")));
+      localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, LocalInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
 
     } else
       DUNE_THROW(Exception, "Unknown energy type '" << energy << "'");