Do not hardwire the local interpolation rule in the embedded global function

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

dune/gfe/embeddedglobalgfefunction.hh

@@ -6,8 +6,6 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 
-#include <dune/gfe/localgeodesicfefunction.hh>
-
 namespace Dune {
 
   namespace GFE {
@@ -17,7 +15,7 @@ namespace Dune {
  *  \tparam B  - The global basis type.
  *  \tparam TargetSpace - The manifold that this functions takes its values in.
  */
-template<class B, class TargetSpace>
+template<class B, class LocalInterpolationRule, class TargetSpace>
 class EmbeddedGlobalGFEFunction
 : public VirtualGridViewFunction<typename B::GridView, typename TargetSpace::CoordinateType>
 {
@@ -29,7 +27,6 @@ public:
     typedef typename GridView::template Codim<0>::Entity Element;
     typedef typename GridView::Grid::ctype  ctype;
 
-    typedef LocalGeodesicFEFunction<GridView::dimension, ctype, LocalFiniteElement, TargetSpace> LocalGFEFunction;
     typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
 
     //! Dimension of the grid.
@@ -82,8 +79,8 @@ public:
             localCoeff[i] = coefficients_[localIndexSet.index(i)];
 
         // create local gfe function
-        LocalGFEFunction localGFE(localView.tree().finiteElement(),localCoeff);
-        return localGFE.evaluate(local).globalCoordinates();
+        LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
+        return localInterpolationRule.evaluate(local).globalCoordinates();
     }
 
     /** \brief Evaluate the derivative of the function at local coordinates. */
@@ -110,10 +107,10 @@ public:
             localCoeff[i] = coefficients_[localIndexSet.index(i)];
 
         // create local gfe function
-        LocalGFEFunction localGFE(localView.tree().finiteElement(),localCoeff);
+        LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
 
         // use it to evaluate the derivative
-        Dune::FieldMatrix<ctype, embeddedDim, gridDim> refJac = localGFE.evaluateDerivative(local);
+        Dune::FieldMatrix<ctype, embeddedDim, gridDim> refJac = localInterpolationRule.evaluateDerivative(local);
 
         Dune::FieldMatrix<ctype, embeddedDim, gridDim> out =0.0;
         //transform the gradient

src/compute-disc-error.cc

@@ -15,8 +15,9 @@
 #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/embeddedglobalgfefunction.hh>
-#include <dune/gfe/cosseratvtkwriter.hh>
 
 // grid dimension
 const int dim = 2;
@@ -39,6 +40,10 @@ void measureDiscreteEOC(const GridView gridView,
   FEBasis feBasis(gridView);
   FEBasis referenceFEBasis(referenceGridView);
 
+  typedef LocalGeodesicFEFunction<GridView::dimension, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace> LocalInterpolationRule;
+  //typedef GFE::LocalProjectedFEFunction<GridView::dimension, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace> LocalInterpolationRule;
+  std::cout << "Using local interpolation: " << className<LocalInterpolationRule>() << std::endl;
+
   //////////////////////////////////////////////////////////////////////////////////
   //  Read the data whose error is to be measured
   //////////////////////////////////////////////////////////////////////////////////
@@ -61,7 +66,7 @@ void measureDiscreteEOC(const GridView gridView,
     x[i] = TargetSpace(embeddedX[i]);
 
   // The numerical solution, as a grid function
-  GFE::EmbeddedGlobalGFEFunction<FEBasis, TargetSpace> numericalSolution(feBasis, x);
+  GFE::EmbeddedGlobalGFEFunction<FEBasis, LocalInterpolationRule, TargetSpace> numericalSolution(feBasis, x);
 
   ///////////////////////////////////////////////////////////////////////////
   // Read the reference configuration
@@ -80,7 +85,7 @@ void measureDiscreteEOC(const GridView gridView,
     referenceX[i] = TargetSpace(embeddedReferenceX[i]);
 
   // The reference solution, as a grid function
-  GFE::EmbeddedGlobalGFEFunction<FEBasis, TargetSpace> referenceSolution(referenceFEBasis, referenceX);
+  GFE::EmbeddedGlobalGFEFunction<FEBasis, LocalInterpolationRule, TargetSpace> referenceSolution(referenceFEBasis, referenceX);
 
   /////////////////////////////////////////////////////////////////
   //   Measure the discretization error
@@ -190,6 +195,10 @@ void measureAnalyticalEOC(const GridView gridView,
   typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
   FEBasis feBasis(gridView);
 
+  typedef LocalGeodesicFEFunction<GridView::dimension, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace> LocalInterpolationRule;
+  //typedef GFE::LocalProjectedFEFunction<GridView::dimension, double, typename FEBasis::LocalView::Tree::FiniteElement, TargetSpace> LocalInterpolationRule;
+  std::cout << "Using local interpolation: " << className<LocalInterpolationRule>() << std::endl;
+
   //////////////////////////////////////////////////////////////////////////////////
   //  Read the data whose error is to be measured
   //////////////////////////////////////////////////////////////////////////////////
@@ -222,7 +231,7 @@ void measureAnalyticalEOC(const GridView gridView,
   auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
 
   // The numerical solution, as a grid function
-  GFE::EmbeddedGlobalGFEFunction<FEBasis, TargetSpace> numericalSolution(feBasis, x);
+  GFE::EmbeddedGlobalGFEFunction<FEBasis, LocalInterpolationRule, TargetSpace> numericalSolution(feBasis, x);
 
   // QuadratureRule for the integral of the L^2 error
   QuadratureRuleKey quadKey(dim,6);