LocalIntegralEnergy: Make LocalInterpolationRule a template parameter

Rather than using preprocessor switches.  That's too fragile
and not extensible.

dune/gfe/assemblers/localintegralenergy.hh

@@ -11,11 +11,6 @@
 #include <dune/gfe/densities/localdensity.hh>
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
-#ifdef PROJECTED_INTERPOLATION
-#include <dune/gfe/localprojectedfefunction.hh>
-#else
-#include <dune/gfe/localgeodesicfefunction.hh>
-#endif
 
 #include <dune/elasticity/materials/localdensity.hh>
 
@@ -24,9 +19,10 @@ namespace Dune::GFE {
   /** \brief An energy given as an integral over a density
    *
    * \tparam Basis The scalar finite element basis used to construct the interpolation rule
+   * \tparam LocalInterpolationRule The rule that turns coefficients into functions
    * \tparam TargetSpace The space that the geometric finite element function maps into
    */
-  template<class Basis, class TargetSpace>
+  template<class Basis, class LocalInterpolationRule, class TargetSpace>
   class LocalIntegralEnergy
     : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
   {
@@ -78,13 +74,9 @@ namespace Dune::GFE {
       const auto& deformationLocalFiniteElement = localView.tree().child(_0,0).finiteElement();
       const auto& orientationLocalFiniteElement = localView.tree().child(_1,0).finiteElement();
 
-#ifdef PROJECTED_INTERPOLATION
-      using LocalDeformationGFEFunctionType = GFE::LocalProjectedFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<RT,gridDim> >;
-      using LocalOrientationGFEFunctionType = GFE::LocalProjectedFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<RT,gridDim> >;
-#else
-      using LocalDeformationGFEFunctionType = LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<RT,gridDim> >;
-      using LocalOrientationGFEFunctionType = LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<RT,gridDim> >;
-#endif
+      using LocalDeformationGFEFunctionType = typename std::tuple_element<0, LocalInterpolationRule>::type;
+      using LocalOrientationGFEFunctionType = typename std::tuple_element<1, LocalInterpolationRule>::type;
+
       LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
       LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localOrientationConfiguration);
 

src/film-on-substrate.cc

@@ -508,7 +508,13 @@ int main (int argc, char *argv[]) try
 
     using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
 
-    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ActiveRigidBodyMotion> >(elasticDensity);
+    // Select which type of geometric interpolation to use
+    using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
+    using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+
+    using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
+
+    auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensity);
     auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,*neumannFunctionPtr);
     auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<
         decltype(stressFreeShellFunction), CompositeBasis, RealTuple<ValueType,dim>, Rotation<ValueType,dim> > >(

test/cosseratcontinuumtest.cc

@@ -27,6 +27,7 @@
 #include <dune/gfe/assemblers/sumenergy.hh>
 #include <dune/gfe/assemblers/localintegralenergy.hh>
 #include <dune/gfe/densities/bulkcosseratdensity.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
 #include <dune/gfe/mixedriemanniantrsolver.hh>
 #include <dune/gfe/neumannenergy.hh>
 
@@ -203,8 +204,15 @@ int main (int argc, char *argv[])
 
   GFE::SumEnergy<CompositeBasis, RealTuple<adouble,dim>,Rotation<adouble,dim> > sumEnergy;
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<adouble,dim>, Rotation<adouble,dim> > >(neumannBoundary,neumannFunction);
+
+  // Select which type of geometric interpolation to use
+  using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
+  using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+
+  using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
+
   auto bulkCosseratDensity = std::make_shared<GFE::BulkCosseratDensity<FieldVector<double,dim>,adouble> >(parameters);
-  auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ARigidBodyMotion> >(bulkCosseratDensity);
+  auto bulkCosseratEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ARigidBodyMotion> >(bulkCosseratDensity);
   sumEnergy.addLocalEnergy(bulkCosseratEnergy);
   sumEnergy.addLocalEnergy(neumannEnergy);
 

test/filmonsubstratetest.cc

@@ -335,7 +335,13 @@ int main (int argc, char *argv[])
 
   using ActiveRigidBodyMotion = GFE::ProductManifold<RealTuple<adouble,dim>, Rotation<adouble,dim> >;
 
-  auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, ActiveRigidBodyMotion> >(elasticDensity);
+  // Select which type of geometric interpolation to use
+  using LocalDeformationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(deformationFEBasis.localView().tree().finiteElement()), RealTuple<adouble,dim> >;
+  using LocalOrientationInterpolationRule = LocalGeodesicFEFunction<dim, GridType::ctype, decltype(orientationFEBasis.localView().tree().finiteElement()), Rotation<adouble,dim> >;
+
+  using LocalInterpolationRule = std::tuple<LocalDeformationInterpolationRule,LocalOrientationInterpolationRule>;
+
+  auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<CompositeBasis, LocalInterpolationRule, ActiveRigidBodyMotion> >(elasticDensity);
   auto neumannEnergy = std::make_shared<GFE::NeumannEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,dim> > >(neumannBoundary,neumannFunction);
 
   auto surfaceCosseratEnergy = std::make_shared<GFE::SurfaceCosseratEnergy<