.git-blame-ignore-revs

 # First use of uncrustify, with the vanilla config from dune-project.org
 d13f34469065230dd86e63733e4dc8d1f775acfb
+
+# Cleanup after having moved everything into the Dune::GFE namespace
+c2eb8dacb5885363c30c2b7373c479ed9c66ac12

.gitlab-ci.yml

@@ -7,9 +7,10 @@ before_script: &before
   - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
   - duneci-install-module https://gitlab.dune-project.org/extensions/dune-vtk.git
   - duneci-install-module https://gitlab.dune-project.org/fufem/dune-fufem.git
-  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
-  - duneci-install-module https://gitlab.mn.tu-dresden.de/iwr/dune-gmsh4
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-gmsh4
 
+#  Tests with the 2.9 release.  That's the one in Debian bookworm
+#--------------------------------------------------------------------
 dune:2.9 gcc:
   variables:
     DUNECI_BRANCH: releases/2.9
@@ -18,6 +19,76 @@ dune:2.9 gcc:
   - *before
   script: duneci-standard-test
 
+dune:2.9 dune-elasticity gcc:
+  variables:
+    DUNECI_BRANCH: releases/2.9
+  image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
+  script: duneci-standard-test
+
+#  Tests with the 2.10 release.  That's the one in Debian trixie
+#--------------------------------------------------------------------
+dune:2.10 gcc-10:
+  image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
+  variables:
+    DUNECI_BRANCH: releases/2.10
+  before_script:
+  - *before
+  script: duneci-standard-test
+
+dune:2.10 dune-elasticity gcc-10:
+  image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
+  variables:
+    DUNECI_BRANCH: releases/2.10
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
+  script: duneci-standard-test
+
+dune:2.10 dune-curvedgrid dune-foamgrid gcc-10 C++20:
+  image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
+  variables:
+    DUNECI_BRANCH: releases/2.10
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
+  script: duneci-standard-test
+
+dune:2.10 dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity gcc-10 C++20:
+  image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
+  variables:
+    DUNECI_BRANCH: releases/2.10
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
+  script: duneci-standard-test
+
+# Full tests again, but this time with clang.
+dune:2.10 dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity clang-11 C++20:
+  image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-clang-11-20
+  variables:
+    DUNECI_BRANCH: releases/2.10
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
+  script: duneci-standard-test
+
+
+#  Tests with the git master branch
+#  Some, but not all optional dependencies are made available.
+#------------------------------------------------------------------
 dune:git gcc-10 C++20:
   image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
   before_script:
@@ -30,24 +101,44 @@ dune:git clang-11 C++20:
   - *before
   script: duneci-standard-test
 
-dune:git dune-parmg dune-curvedgeometry dune-curvedgrid dune-foamgrid gcc-10 C++20:
+dune:git dune-elasticity gcc-10 C++20:
+  image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
+  script: duneci-standard-test
+
+dune:git dune-curvedgrid dune-foamgrid gcc-10 C++20:
+  image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
+  before_script:
+  - *before
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
+  script: duneci-standard-test
+
+#  Tests with all optional dependencies
+#-------------------------------------------
+dune:git dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity gcc-10 C++20:
   image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
   before_script:
   - *before
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
-  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
-  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
   - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
   script: duneci-standard-test
 
-dune:git dune-parmg dune-curvedgeometry dune-curvedgrid dune-foamgrid clang-11 C++20:
+dune:git dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity clang-11 C++20:
   image: registry.dune-project.org/docker/ci/dune:git-debian-11-clang-11-20
   before_script:
   - *before
   - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
-  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
-  - duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
   - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
+  - duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
   script: duneci-standard-test
 
 # Check for spelling mistakes in text

CHANGELOG.md

 # Master
 
+- The `SimoFoxEnergy` class (formerly known as `SimoFoxEnergyLocalStiffness`)
+  does not accept a surface load density anymore. Use the
+  `NeumannEnergy` class for assembling surface loads.
+
+- The class `LocalIntegralEnergy` now does not construct the local interpolation
+  rule itself anymore.  Rather, it expects to be given one in its constructor.
+  This will allow interpolation rules with state.
+  
+- The same holds for `CosseratRodEnergy`.  Unfortunately, that move
+  implies that `CosseratRodEnergy` gains an extra template argument
+  which specifies the GFE function used to implement the reference
+  configuration.
+
+- The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
+  now take a `dune-functions` basis as their first arguments, instead
+  of a `LocalFiniteElement`. In addition, they need to be given an object
+  of this type in the constructor.  This is the basis to be used for
+  interpolation weights (for `LocalGeodesicFEFunction`) and for
+  interpolating in the embedding space (for `LocalProjectedFEFunction`).
+
+- The entire code is now in the `Dune::GFE` namespace.
+
+- Remove the files `localgfetestfunctionbasis.hh` and `localtangentfefunction.hh`.
+  They were never used for anything at all, I and currently wouldn't know
+  what to use them for, either.
+
+- All implementations of functions have been move to a `functions` subdirectory.
+
+- The `LocalDensity` class now has a template parameter `ElementOrIntersection`,
+  which replaces the parameter `Position`.  As the name says, this parameter
+  has to be a grid element, or a grid intersection.  As it turned out,
+  just depending on the position type of the integration domain was not enough.
+
+- The `LocalDensity` class now has a `bind` method, which binds the density
+  to a given element or intersection. This is necessary, for example, to in turn
+  bind coefficient functions that the density may own.
+
+# Release 2.10 (2024-12-29)
+
+- The module `dune-elasticity` has been downgraded from a required
+  dependency to an optional dependency.  It is currently only needed
+  by the `film-on-substrate.cc` program.
+
+- `cosserat-rod.cc` (formerly `rod3d.cc`) now reads the reference configuration
+  and the initial iterate from the Python file.
+
+- The file `rod3d.cc` has been renamed to `cosserat-rod.cc`,
+  to better reflect what it does.
+
+- In the `BulkCosseratDensity` class: Make the type of curvature
+  tensor controllable at run-time (previously, ugly preprocessor switches
+  where used).  Among its parameters `BulkCosseratDensity` now expects
+  a string `curvatureType`, which has to take one of the values `norm`,
+  `curl`, or `wryness`.
+
 - Building `dune-gfe` now requires Dune 2.9 or newer. Older configurations
   have not gotten CI-tested for quite a while, anyway.
 

CMakeLists.txt

@@ -17,6 +17,13 @@ include(DuneMacros)
 
 # start a dune project with information from dune.module
 dune_project()
+
+# Create a module library
+if(dune-common_VERSION VERSION_GREATER_EQUAL 2.9.0)
+  dune_add_library(dunegfe INTERFACE EXPORT_NAME GFE NAMESPACE Dune::
+    LINK_LIBRARIES ${DUNE_LIBS})
+endif()
+
 dune_enable_all_packages()
 add_subdirectory("src")
 add_subdirectory("dune")

dune.module

@@ -4,8 +4,8 @@
 
 #Name of the module
 Module: dune-gfe
-Version: svn
+Version: 2.11-git
 Maintainer: oliver.sander@tu-dresden.de
 #depending on
-Depends: dune-grid(>=2.9) dune-uggrid dune-istl dune-localfunctions dune-functions (>=2.9) dune-solvers dune-fufem dune-elasticity (>=2.9) dune-gmsh4
-Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgrid
+Depends: dune-grid(>=2.9) dune-uggrid dune-istl dune-localfunctions dune-functions (>=2.9) dune-solvers dune-fufem dune-gmsh4 dune-vtk
+Suggests: dune-foamgrid dune-parmg dune-curvedgrid dune-elasticity

dune/gfe/CMakeLists.txt

@@ -4,20 +4,12 @@ install(FILES
         cosseratstrain.hh
         cosseratvtkreader.hh
         cosseratvtkwriter.hh
-        embeddedglobalgfefunction.hh
         filereader.hh
         geodesicdifference.hh
         geodesicfefunctionadaptor.hh
         gfedifferencenormsquared.hh
-        globalgfefunction.hh
         gramschmidtsolver.hh
-        interpolationderivatives.hh
         linearalgebra.hh
-        localgeodesicfefunction.hh
-        localgfetestfunctionbasis.hh
-        localprojectedfefunction.hh
-        localquickanddirtyfefunction.hh
-        localtangentfefunction.hh
         maxnormtrustregion.hh
         mixedriemanniantrsolver.cc
         mixedriemanniantrsolver.hh
@@ -42,3 +34,4 @@ install(FILES
 add_subdirectory("spaces")
 add_subdirectory("assemblers")
 add_subdirectory("densities")
+add_subdirectory("functions")

dune/gfe/assemblers/CMakeLists.txt

 #install headers
 install(FILES
-        cosseratenergystiffness.hh
         cosseratrodenergy.hh
         geodesicfeassembler.hh
         geodesicfeassemblerwrapper.hh

dune/gfe/assemblers/cosseratrodenergy.hh

@@ -2,6 +2,7 @@
 #define DUNE_GFE_COSSERATRODENERGY_HH
 
 #include <array>
+#include <memory>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/version.hh>
@@ -18,9 +19,10 @@
 #include <dune/gfe/spaces/realtuple.hh>
 #include <dune/gfe/spaces/rotation.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
-  template<class Basis, class LocalInterpolationRule, class RT>
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
   class CosseratRodEnergy
     : public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > >
   {
@@ -42,6 +44,14 @@ namespace Dune::GFE {
 
   public:
 
+    /** \brief The current configuration whose energy is being computed
+     */
+    const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
+
+    /** \brief The GFE function that implements the stress-free reference configuration
+     */
+    const std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration_;
+
     /** \brief The stress-free configuration
 
        The number type cannot be RT, because RT become `adouble` when
@@ -67,9 +77,25 @@ namespace Dune::GFE {
     GridView gridView_;
 
     //! Constructor
-    CosseratRodEnergy(const GridView& gridView,
+    CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+                      std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
+                      const GridView& gridView,
+                      const std::array<double,3>& K, const std::array<double,3>& A)
+      : localGFEFunction_(localGFEFunction)
+      , localReferenceConfiguration_(localReferenceConfiguration)
+      , K_(K),
+      A_(A),
+      gridView_(gridView)
+    {}
+
+    //! Constructor
+    CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
+                      ReferenceInterpolationRule&& localReferenceConfiguration,
+                      const GridView& gridView,
                       const std::array<double,3>& K, const std::array<double,3>& A)
-      : K_(K),
+      : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
+      , localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
+      , K_(K),
       A_(A),
       gridView_(gridView)
     {}
@@ -80,9 +106,39 @@ namespace Dune::GFE {
         \param E Young's modulus
         \param nu Poisson number
      */
-    CosseratRodEnergy(const GridView& gridView,
+    CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+                      std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
+                      const GridView& gridView,
                       double A, double J1, double J2, double E, double nu)
-      : gridView_(gridView)
+      : localGFEFunction_(localGFEFunction)
+      , localReferenceConfiguration_(localReferenceConfiguration)
+      , gridView_(gridView)
+    {
+      // shear modulus
+      double G = E/(2+2*nu);
+
+      K_[0] = E * J1;
+      K_[1] = E * J2;
+      K_[2] = G * (J1 + J2);
+
+      A_[0] = G * A;
+      A_[1] = G * A;
+      A_[2] = E * A;
+    }
+
+    /** \brief Constructor setting shape constants and material parameters
+        \param A The rod section area
+        \param J1, J2 The geometric moments (Flächenträgheitsmomente)
+        \param E Young's modulus
+        \param nu Poisson number
+     */
+    CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
+                      ReferenceInterpolationRule&& localReferenceConfiguration,
+                      const GridView& gridView,
+                      double A, double J1, double J2, double E, double nu)
+      : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
+      , localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
+      , gridView_(gridView)
     {
       // shear modulus
       double G = E/(2+2*nu);
@@ -172,21 +228,18 @@ namespace Dune::GFE {
 
   };
 
-  template<class Basis, class LocalInterpolationRule, class RT>
-  RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
+  RT CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
   energy(const typename Basis::LocalView& localView,
          const std::vector<TargetSpace>& localCoefficients) const
   {
-    const auto& localFiniteElement = localView.tree().finiteElement();
-    LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
-
     const auto& element = localView.element();
+    localGFEFunction_->bind(element, localCoefficients);
 
     RT energy = 0;
 
     std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
-    using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::other;
-    InactiveLocalInterpolationRule localReferenceConfiguration(localFiniteElement, localReferenceCoefficients);
+    localReferenceConfiguration_->bind(element, localReferenceCoefficients);
 
     // ///////////////////////////////////////////////////////////////////////////////
     //   The following two loops are a reduced integration scheme.  We integrate
@@ -205,10 +258,10 @@ namespace Dune::GFE {
 
       double weight = shearingQuad[pt].weight() * integrationElement;
 
-      auto strain = getStrain(localConfiguration, element, quadPos);
+      auto strain = getStrain(*localGFEFunction_, element, quadPos);
 
       // The reference strain
-      auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+      auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
 
       for (int i=0; i<3; i++)
         energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
@@ -225,10 +278,10 @@ namespace Dune::GFE {
 
       double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
 
-      auto strain = getStrain(localConfiguration, element, quadPos);
+      auto strain = getStrain(*localGFEFunction_, element, quadPos);
 
       // The reference strain
-      auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+      auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
 
       // Part II: the bending and twisting energy
       for (int i=0; i<3; i++)
@@ -240,9 +293,9 @@ namespace Dune::GFE {
   }
 
 
-  template<class Basis, class LocalInterpolationRule, class RT>
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
   template <class ReboundLocalInterpolationRule>
-  auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+  auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
   getStrain(const ReboundLocalInterpolationRule& localInterpolation,
             const Entity& element,
             const FieldVector<double,1>& pos) const
@@ -285,9 +338,9 @@ namespace Dune::GFE {
     return strain;
   }
 
-  template<class Basis, class LocalInterpolationRule, class RT>
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
   template <class Number>
-  auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+  auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
   getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
             const Entity& element,
             const FieldVector<double, 1>& pos) const
@@ -308,8 +361,8 @@ namespace Dune::GFE {
     return stress;
   }
 
-  template<class Basis, class LocalInterpolationRule, class RT>
-  void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
+  void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
   getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
             BlockVector<FieldVector<double, blocksize> >& strain) const
   {
@@ -363,8 +416,8 @@ namespace Dune::GFE {
     }
   }
 
-  template<class Basis, class LocalInterpolationRule, class RT>
-  void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
+  void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
   getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
             BlockVector<FieldVector<double, blocksize> >& stress) const
   {
@@ -386,9 +439,9 @@ namespace Dune::GFE {
     }
   }
 
-  template<class Basis, class LocalInterpolationRule, class RT>
+  template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
   template <class PatchGridView>
-  auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+  auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
   getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
                     const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const
   {

dune/gfe/assemblers/discretekirchhoffbendingenergy.hh

+#ifndef DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
+#define DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/assemblers/localenergy.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+#include <dune/gfe/tensor3.hh>
+
+#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
+
+/**
+ * \brief Assemble the discrete Kirchhoff bending energy for a single element.
+ *
+ *  The energy is essentially the harmonic energy of a
+ *  discrete Jacobian operator which is represented as a
+ *  linear combination in a P2 finite element space.
+ *  The gradient of that linear combination serves as an
+ *  discrete Hessian in the energy.
+ */
+namespace Dune::GFE
+{
+  template<class Basis, class DiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
+  class DiscreteKirchhoffBendingEnergy
+    : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+  {
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef typename TargetSpace::ctype RT;
+    typedef typename DiscreteKirchhoffFunction::DerivativeType DerivativeType;
+
+    // Grid dimension
+    constexpr static int gridDim = GridView::dimension;
+
+    // P2-LocalFiniteElement used to represent the discrete Jacobian on the current element.
+    typedef typename Dune::LagrangeSimplexLocalFiniteElement<DT, double, gridDim, 2> P2LocalFiniteElement;
+
+  public:
+    DiscreteKirchhoffBendingEnergy(DiscreteKirchhoffFunction &localFunction)
+      : localFunction_(localFunction)
+    {}
+
+    /** \brief Assemble the energy for a single element */
+    virtual RT energy (const typename Basis::LocalView& localView,
+                       const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
+
+    /** \brief Assemble the energy for a single element */
+    virtual RT energy (const typename Basis::LocalView& localView,
+                       const std::vector<TargetSpace>& localConfiguration) const override
+    {
+
+      const auto element = localView.element();
+
+      int quadOrder = 3;
+      const auto &quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
+
+      /**
+          bind the underlying Hermite basis to the current element.
+          Note: The DiscreteKirchhoffFunction object stores a vector of global
+          and local basis coefficients that are accessed for evaluation
+          methods (linear combinations of hermite basis function and local coefficients).
+          Therefore the local configuration is passed to the bind method as well since
+          these coefficient vectors need to be updated with the new local configuration
+          previous to evaluation.
+       */
+      localFunction_.bind(element,localConfiguration);
+
+      RT bendingEnergy = 0;
+
+      for (auto&& quadPoint : quad)
+      {
+        auto discreteHessian= localFunction_.evaluateDiscreteHessian(quadPoint.position());
+
+        auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+        bendingEnergy += weight * discreteHessian.frobenius_norm2();
+      }
+
+      return 0.5 * bendingEnergy;
+    }
+
+  private:
+    DiscreteKirchhoffFunction& localFunction_;
+
+    P2LocalFiniteElement lagrangeLFE_;
+  };
+
+}  // namespace Dune::GFE
+#endif

dune/gfe/assemblers/forceenergy.hh

+#ifndef DUNE_GFE_FORCEENERGY_HH
+#define DUNE_GFE_FORCEENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/gfe/assemblers/localenergy.hh>
+#include <dune/gfe/bendingisometryhelper.hh>
+
+
+/**
+ * \brief Assemble the energy of a force term (force times deformation) for a single element.
+ */
+namespace Dune::GFE
+{
+  template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
+  class ForceEnergy
+    : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+  {
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef typename TargetSpace::ctype RT;
+
+    // Grid dimension
+    constexpr static int gridDim = GridView::dimension;
+
+  public:
+    ForceEnergy(LocalDiscreteKirchhoffFunction &localFunction,
+                LocalForce &localForce)
+      : localFunction_(localFunction),
+      localForce_(localForce)
+    {}
+
+    /** \brief Assemble the energy for a single element */
+    virtual RT energy (const typename Basis::LocalView& localView,
+                       const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
+
+    /** \brief Assemble the energy for a single element */
+    virtual RT energy (const typename Basis::LocalView& localView,
+                       const std::vector<TargetSpace>& localConfiguration) const override
+    {
+      RT energy = 0;
+
+      int quadOrder = 3;
+
+      const auto &quad = QuadratureRules<double, gridDim>::rule(GeometryTypes::simplex(gridDim), quadOrder);
+
+      const auto element = localView.element();
+
+      /** bind the local force to the current element */
+      localForce_.bind(element);
+
+      /**
+          bind the underlying Hermite basis to the current element.
+          Note: The DiscreteKirchhoffFunction object stores a vector of global
+          and local basis coefficients that are accessed for evaluation
+          methods (linear comibnations of hermite basis function and local coefficients).
+          Therefore the local configuration is passed to the bind method as well since
+          these coefficient vectors need to be updated with the new local configuration
+          previous to evaluation.
+       */
+      localFunction_.bind(element,localConfiguration);
+
+      /**
+       * @brief Compute contribution of the force Term
+       *
+       * Integrate the scalar product of the force
+       * with the local deformation function 'localFunction_'
+       * over the current element.
+       */
+      RT forceTermEnergy = 0;
+
+      for (auto&& quadPoint : quad)
+      {
+        auto deformationValue = localFunction_.evaluate(quadPoint.position());
+        auto forceValue = localForce_(quadPoint.position());
+        auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+        forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
+      }
+
+      return forceTermEnergy;
+    }
+
+  private:
+    LocalDiscreteKirchhoffFunction& localFunction_;
+
+    LocalForce& localForce_;
+
+  };
+
+}  // namespace Dune::GFE
+#endif

dune/gfe/assemblers/geodesicfeassembler.hh

@@ -10,269 +10,273 @@
 
 #include <dune/solvers/common/wrapownshare.hh>
 
-/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
- */
-template <class Basis, class TargetSpace>
-class GeodesicFEAssembler {
 
-  using field_type = typename TargetSpace::field_type;
-  typedef typename Basis::GridView GridView;
-  using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
+namespace Dune::GFE
+{
 
-  //! Dimension of the grid.
-  constexpr static int gridDim = GridView::dimension;
+  /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
+   */
+  template <class Basis, class TargetSpace>
+  class GeodesicFEAssembler {
 
-  //! Dimension of a tangent space
-  constexpr static int blocksize = TargetSpace::TangentVector::dimension;
+    using field_type = typename TargetSpace::field_type;
+    typedef typename Basis::GridView GridView;
+    using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
 
-  //!
-  typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
+    //! Dimension of the grid.
+    constexpr static int gridDim = GridView::dimension;
 
+    //! Dimension of a tangent space
+    constexpr static int blocksize = TargetSpace::TangentVector::dimension;
 
-protected:
+    //!
+    typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
 
-  //! The global basis
-  const Basis basis_;
 
-  //! The local stiffness operator
-  std::shared_ptr<LocalStiffness> localStiffness_;
+  protected:
 
-public:
+    //! The global basis
+    const Basis basis_;
 
-  /** \brief Constructor for a given grid */
-  GeodesicFEAssembler(const Basis& basis)
-    : basis_(basis)
-  {}
+    //! The local stiffness operator
+    std::shared_ptr<LocalStiffness> localStiffness_;
 
-  /** \brief Constructor for a given grid */
-  template <class LocalStiffnessT>
-  GeodesicFEAssembler(const Basis& basis,
-                      LocalStiffnessT&& localStiffness)
-    : basis_(basis),
-    localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
-  {}
+  public:
 
-  /** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
-  template <class LocalStiffnessT>
-  void setLocalStiffness(LocalStiffnessT&& localStiffness) {
-    localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
-  }
-
-  /** \brief Get the local stiffness operator. */
-  const LocalStiffness& getLocalStiffness() const {
-    return *localStiffness_;
-  }
+    /** \brief Constructor for a given grid */
+    GeodesicFEAssembler(const Basis& basis)
+      : basis_(basis)
+    {}
 
-  /** \brief Get the local stiffness operator. */
-  LocalStiffness& getLocalStiffness() {
-    return *localStiffness_;
-  }
+    /** \brief Constructor for a given grid */
+    template <class LocalStiffnessT>
+    GeodesicFEAssembler(const Basis& basis,
+                        LocalStiffnessT&& localStiffness)
+      : basis_(basis),
+      localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
+    {}
 
-  /** \brief Get the basis. */
-  const Basis& getBasis() const {
-    return basis_;
-  }
+    /** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
+    template <class LocalStiffnessT>
+    void setLocalStiffness(LocalStiffnessT&& localStiffness) {
+      localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
+    }
 
-  /** \brief Assemble the tangent stiffness matrix and the functional gradient together
-   *
-   * This is more efficient than computing them separately, because you need the gradient
-   * anyway to compute the Riemannian Hessian.
-   */
-  virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
-                                          Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >& gradient,
-                                          Dune::BCRSMatrix<MatrixBlock>& hessian,
-                                          bool computeOccupationPattern=true) const;
+    /** \brief Get the local stiffness operator. */
+    const LocalStiffness& getLocalStiffness() const {
+      return *localStiffness_;
+    }
 
-  /** \brief Assemble the gradient */
-  virtual void assembleGradient(const std::vector<TargetSpace>& sol,
-                                Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
+    /** \brief Get the local stiffness operator. */
+    LocalStiffness& getLocalStiffness() {
+      return *localStiffness_;
+    }
 
-  /** \brief Compute the energy of a deformation state */
-  virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
+    /** \brief Get the basis. */
+    const Basis& getBasis() const {
+      return basis_;
+    }
 
-  //protected:
-  void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
+    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+     *
+     * This is more efficient than computing them separately, because you need the gradient
+     * anyway to compute the Riemannian Hessian.
+     */
+    virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+                                            Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >& gradient,
+                                            Dune::BCRSMatrix<MatrixBlock>& hessian,
+                                            bool computeOccupationPattern=true) const;
 
-}; // end class
+    /** \brief Assemble the gradient */
+    virtual void assembleGradient(const std::vector<TargetSpace>& sol,
+                                  Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
 
+    /** \brief Compute the energy of a deformation state */
+    virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
 
+    //protected:
+    void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
 
-template <class Basis, class TargetSpace>
-void GeodesicFEAssembler<Basis,TargetSpace>::
-getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
-{
-  auto n = basis_.size();
+  }; // end class
 
-  nb.resize(n, n);
 
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
 
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+  template <class Basis, class TargetSpace>
+  void GeodesicFEAssembler<Basis,TargetSpace>::
+  getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
   {
-    // Bind the local FE basis view to the current element
-    localView.bind(element);
+    auto n = basis_.size();
 
-    const auto& lfe = localView.tree().finiteElement();
+    nb.resize(n, n);
 
-    for (size_t i=0; i<lfe.size(); i++) {
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
 
-      for (size_t j=0; j<lfe.size(); j++) {
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+    {
+      // Bind the local FE basis view to the current element
+      localView.bind(element);
 
-        auto iIdx = localView.index(i);
-        auto jIdx = localView.index(j);
+      const auto& lfe = localView.tree().finiteElement();
 
-        nb.add(iIdx, jIdx);
+      for (size_t i=0; i<lfe.size(); i++) {
 
-      }
+        for (size_t j=0; j<lfe.size(); j++) {
 
-    }
+          auto iIdx = localView.index(i);
+          auto jIdx = localView.index(j);
 
-  }
+          nb.add(iIdx, jIdx);
 
-}
+        }
 
-template <class Basis, class TargetSpace>
-void GeodesicFEAssembler<Basis,TargetSpace>::
-assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
-                           Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > &gradient,
-                           Dune::BCRSMatrix<MatrixBlock>& hessian,
-                           bool computeOccupationPattern) const
-{
-  if (computeOccupationPattern) {
+      }
 
-    Dune::MatrixIndexSet neighborsPerVertex;
-    getNeighborsPerVertex(neighborsPerVertex);
-    neighborsPerVertex.exportIdx(hessian);
+    }
 
   }
 
-  hessian = 0;
+  template <class Basis, class TargetSpace>
+  void GeodesicFEAssembler<Basis,TargetSpace>::
+  assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+                             Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > &gradient,
+                             Dune::BCRSMatrix<MatrixBlock>& hessian,
+                             bool computeOccupationPattern) const
+  {
+    if (computeOccupationPattern) {
 
-  gradient.resize(sol.size());
-  gradient = 0;
+      Dune::MatrixIndexSet neighborsPerVertex;
+      getNeighborsPerVertex(neighborsPerVertex);
+      neighborsPerVertex.exportIdx(hessian);
 
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
+    }
 
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
-  {
-    localView.bind(element);
+    hessian = 0;
 
-    const int numOfBaseFct = localView.tree().size();
+    gradient.resize(sol.size());
+    gradient = 0;
 
-    // Extract local solution
-    std::vector<TargetSpace> localSolution(numOfBaseFct);
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
 
-    for (int i=0; i<numOfBaseFct; i++)
-      localSolution[i] = sol[localView.index(i)];
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+    {
+      localView.bind(element);
 
-    std::vector<double> localGradient(numOfBaseFct*blocksize);
-    Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
+      const int numOfBaseFct = localView.tree().size();
 
-    // setup local matrix and gradient
-    localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
+      // Extract local solution
+      std::vector<TargetSpace> localSolution(numOfBaseFct);
 
-    // Add element matrix to global stiffness matrix
-    for(int i=0; i<numOfBaseFct; i++) {
+      for (int i=0; i<numOfBaseFct; i++)
+        localSolution[i] = sol[localView.index(i)];
 
-      auto row = localView.index(i);
+      std::vector<double> localGradient(numOfBaseFct*blocksize);
+      Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
 
-      for (int j=0; j<numOfBaseFct; j++ ) {
+      // setup local matrix and gradient
+      localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
 
-        auto col = localView.index(j);
+      // Add element matrix to global stiffness matrix
+      for(int i=0; i<numOfBaseFct; i++) {
 
-        for (std::size_t ii=0; ii<blocksize; ii++)
-          for (std::size_t jj=0; jj<blocksize; jj++)
-            hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
+        auto row = localView.index(i);
 
-      }
-    }
+        for (int j=0; j<numOfBaseFct; j++ ) {
 
-    // Add local gradient to global gradient
-    for (int i=0; i<numOfBaseFct; i++)
-      for (int j=0; j<blocksize; j++)
-        gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
+          auto col = localView.index(j);
 
-  }
+          for (std::size_t ii=0; ii<blocksize; ii++)
+            for (std::size_t jj=0; jj<blocksize; jj++)
+              hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
 
-}
+        }
+      }
 
-template <class Basis, class TargetSpace>
-void GeodesicFEAssembler<Basis,TargetSpace>::
-assembleGradient(const std::vector<TargetSpace>& sol,
-                 Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
-{
-  if (sol.size()!=basis_.size())
-    DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
+      // Add local gradient to global gradient
+      for (int i=0; i<numOfBaseFct; i++)
+        for (int j=0; j<blocksize; j++)
+          gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
 
-  grad.resize(sol.size());
-  grad = 0;
+    }
 
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
+  }
 
-  // Loop over all elements
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+  template <class Basis, class TargetSpace>
+  void GeodesicFEAssembler<Basis,TargetSpace>::
+  assembleGradient(const std::vector<TargetSpace>& sol,
+                   Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
   {
-    localView.bind(element);
+    if (sol.size()!=basis_.size())
+      DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
 
-    // The number of degrees of freedom of the current element
-    const auto nDofs = localView.tree().size();
+    grad.resize(sol.size());
+    grad = 0;
 
-    // Extract local solution
-    std::vector<TargetSpace> localSolution(nDofs);
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
 
-    for (size_t i=0; i<nDofs; i++)
-      localSolution[i] = sol[localView.index(i)];
+    // Loop over all elements
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+    {
+      localView.bind(element);
 
-    // Assemble local gradient
-    std::vector<double> localGradient(nDofs*blocksize);
+      // The number of degrees of freedom of the current element
+      const auto nDofs = localView.tree().size();
 
-    localStiffness_->assembleGradient(localView, localSolution, localGradient);
+      // Extract local solution
+      std::vector<TargetSpace> localSolution(nDofs);
 
-    // Add to global gradient
-    for (size_t i=0; i<nDofs; i++)
-      for (size_t j=0; j<blocksize; j++)
-        grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
-  }
+      for (size_t i=0; i<nDofs; i++)
+        localSolution[i] = sol[localView.index(i)];
 
-}
+      // Assemble local gradient
+      std::vector<double> localGradient(nDofs*blocksize);
 
+      localStiffness_->assembleGradient(localView, localSolution, localGradient);
 
-template <class Basis, class TargetSpace>
-double GeodesicFEAssembler<Basis, TargetSpace>::
-computeEnergy(const std::vector<TargetSpace>& sol) const
-{
-  double energy = 0;
+      // Add to global gradient
+      for (size_t i=0; i<nDofs; i++)
+        for (size_t j=0; j<blocksize; j++)
+          grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
+    }
 
-  if (sol.size() != basis_.size())
-    DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
+  }
 
-  // A view on the FE basis on a single element
-  auto localView = basis_.localView();
 
-  // Loop over all elements
-  for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+  template <class Basis, class TargetSpace>
+  double GeodesicFEAssembler<Basis, TargetSpace>::
+  computeEnergy(const std::vector<TargetSpace>& sol) const
   {
-    localView.bind(element);
+    double energy = 0;
 
-    // Number of degrees of freedom on this element
-    size_t nDofs = localView.tree().size();
+    if (sol.size() != basis_.size())
+      DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
 
-    std::vector<TargetSpace> localSolution(nDofs);
+    // A view on the FE basis on a single element
+    auto localView = basis_.localView();
 
-    for (size_t i=0; i<nDofs; i++)
-      localSolution[i] = sol[localView.index(i)[0]];
+    // Loop over all elements
+    for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
+    {
+      localView.bind(element);
 
-    energy += localStiffness_->energy(localView, localSolution);
+      // Number of degrees of freedom on this element
+      size_t nDofs = localView.tree().size();
 
-  }
+      std::vector<TargetSpace> localSolution(nDofs);
 
-  return energy;
+      for (size_t i=0; i<nDofs; i++)
+        localSolution[i] = sol[localView.index(i)[0]];
 
-}
+      energy += localStiffness_->energy(localView, localSolution);
 
+    }
+
+    return energy;
+
+  }
 
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/assemblers/geodesicfeassemblerwrapper.hh

@@ -5,7 +5,8 @@
 #include <dune/common/tuplevector.hh>
 #include <dune/gfe/spaces/productmanifold.hh>
 
-namespace Dune::GFE {
+namespace Dune::GFE
+{
 
   /** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two finite element spaces
 
@@ -13,9 +14,8 @@ namespace Dune::GFE {
    */
 
   template <class Basis, class ScalarBasis, class TargetSpace>
-  class
-    GeodesicFEAssemblerWrapper {
-
+  class GeodesicFEAssemblerWrapper
+  {
     using MixedSpace0 = std::tuple_element_t<0,TargetSpace>;
     using MixedSpace1 = std::tuple_element_t<1,TargetSpace>;
 
@@ -72,7 +72,8 @@ namespace Dune::GFE {
     auto splitVector(const std::vector<TargetSpace>& sol) const;
     std::unique_ptr<MatrixType> hessianMixed_;
   }; // end class
-} //end namespace
+
+}  // namespace Dune::GFE
 
 
 template <class Basis, class ScalarBasis, class TargetSpace>

dune/gfe/assemblers/l2distancesquaredenergy.hh

@@ -5,86 +5,94 @@
 
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/gfe/globalgfefunction.hh>
 #include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/functions/globalgfefunction.hh>
+#include <dune/gfe/functions/localgeodesicfefunction.hh>
 
-template<class Basis, class TargetSpace>
-class L2DistanceSquaredEnergy
-  : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+
+namespace Dune::GFE
 {
-  // grid types
-  typedef typename Basis::GridView GridView;
-  typedef typename GridView::ctype DT;
-  typedef typename TargetSpace::ctype RT;
 
-  // some other sizes
-  constexpr static int gridDim = GridView::dimension;
+  template<class Basis, class TargetSpace>
+  class L2DistanceSquaredEnergy
+    : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
+  {
+    // grid types
+    typedef typename Basis::GridView GridView;
+    typedef typename GridView::ctype DT;
+    typedef typename TargetSpace::ctype RT;
 
-  using LocalInterpolationRule = LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, typename TargetSpace::template rebind<double>::other>;
+    // some other sizes
+    constexpr static int gridDim = GridView::dimension;
 
-public:
+    using LocalInterpolationRule = LocalGeodesicFEFunction<Basis, typename TargetSpace::template rebind<double>::other>;
 
-  // This is the function that we are computing the L2-distance to
-  std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
+  public:
 
-  /** \brief Assemble the energy for a single element */
-  RT energy (const typename Basis::LocalView& localView,
-             const std::vector<TargetSpace>& localSolution) const override
+    // This is the function that we are computing the L2-distance to
+    std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
 
-  {
-    RT energy = 0;
+    /** \brief Assemble the energy for a single element */
+    RT energy (const typename Basis::LocalView& localView,
+               const std::vector<TargetSpace>& localSolution) const override
 
-    const auto& localFiniteElement = localView.tree().finiteElement();
-    typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
-    LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
+    {
+      RT energy = 0;
 
-    const auto element = localView.element();
-    auto localOrigin = localFunction(*origin_);
-    localOrigin.bind(element);
+      typedef LocalGeodesicFEFunction<Basis, TargetSpace> LocalGFEFunctionType;
+      LocalGFEFunctionType localGeodesicFEFunction(localView.globalBasis());
 
-    // Just guessing an appropriate quadrature order
-    auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
+      const auto element = localView.element();
+      localGeodesicFEFunction.bind(element,localSolution);
 
-    const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
+      auto localOrigin = localFunction(*origin_);
+      localOrigin.bind(element);
 
-    for (size_t pt=0; pt<quad.size(); pt++)
-    {
-      // Local position of the quadrature point
-      const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
+      // Just guessing an appropriate quadrature order
+      const auto& localFiniteElement = localView.tree().finiteElement();
+      auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
+
+      const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
 
-      const double integrationElement = element.geometry().integrationElement(quadPos);
+      for (size_t pt=0; pt<quad.size(); pt++)
+      {
+        // Local position of the quadrature point
+        const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
 
-      auto weight = quad[pt].weight() * integrationElement;
+        const double integrationElement = element.geometry().integrationElement(quadPos);
 
-      // The function value
-      auto value = localGeodesicFEFunction.evaluate(quadPos);
-      auto originValue = localOrigin(quadPos);
+        auto weight = quad[pt].weight() * integrationElement;
 
-      // The derivative of the 'origin' function
-      // First: as function defined on the reference element
-      auto originReferenceDerivative = derivative(localOrigin)(quadPos);
+        // The function value
+        auto value = localGeodesicFEFunction.evaluate(quadPos);
+        auto originValue = localOrigin(quadPos);
 
-      // The derivative of the function defined on the actual element
-      auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
+        // The derivative of the 'origin' function
+        // First: as function defined on the reference element
+        auto originReferenceDerivative = derivative(localOrigin)(quadPos);
 
-      double weightFactor = originDerivative.frobenius_norm();
-      // un-comment the following line to switch off the weight factor
-      //weightFactor = 1.0;
+        // The derivative of the function defined on the actual element
+        auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
 
-      // Add the local energy density
-      energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value);
+        double weightFactor = originDerivative.frobenius_norm();
+        // un-comment the following line to switch off the weight factor
+        //weightFactor = 1.0;
 
+        // Add the local energy density
+        energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value);
+
+      }
+
+      return energy;
     }
 
-    return energy;
-  }
+    virtual RT energy (const typename Basis::LocalView& localView,
+                       const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
+    {
+      DUNE_THROW(Dune::NotImplemented, "!");
+    }
+  };
 
-  virtual RT energy (const typename Basis::LocalView& localView,
-                     const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
-  {
-    DUNE_THROW(Dune::NotImplemented, "!");
-  }
-};
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/assemblers/localenergy.hh

@@ -7,9 +7,9 @@
 
 #include <dune/gfe/spaces/productmanifold.hh>
 
-namespace Dune {
-
-  namespace GFE:: Impl
+namespace Dune::GFE
+{
+  namespace Impl
   {
     /** \brief A class exporting container types for coefficient sets
      *
@@ -34,43 +34,40 @@ namespace Dune {
       using Coefficients = std::vector<ProductManifold<Factors...> >;
       using CompositeCoefficients = TupleVector<std::vector<Factors>... >;
     };
-  }
-
-  namespace GFE {
+  }  // namespace Impl
 
-    /** \brief Base class for energies defined by integrating over one grid element */
-    template<class Basis, class TargetSpace>
-    class LocalEnergy
-    {
-    public:
-      using RT = typename TargetSpace::ctype;
 
-      /** \brief Compute the energy
-       *
-       * \param localView Local view specifying the current element and the FE space there
-       * \param coefficients The coefficients of a FE function on the current element
-       */
-      virtual RT
-      energy (const typename Basis::LocalView& localView,
-              const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
+  /** \brief Base class for energies defined by integrating over one grid element */
+  template<class Basis, class TargetSpace>
+  class LocalEnergy
+  {
+  public:
+    using RT = typename TargetSpace::ctype;
 
-      /** \brief ProductManifolds: Compute the energy from coefficients in separate containers
-       * for each factor
-       */
-      virtual RT
-      energy (const typename Basis::LocalView& localView,
-              const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
+    /** \brief Compute the energy
+     *
+     * \param localView Local view specifying the current element and the FE space there
+     * \param coefficients The coefficients of a FE function on the current element
+     */
+    virtual RT
+    energy (const typename Basis::LocalView& localView,
+            const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
 
-      /** Empty virtual default destructor
-       *
-       * To allow proper destruction of derived classes through a base class pointer
-       */
-      virtual ~LocalEnergy() = default;
+    /** \brief ProductManifolds: Compute the energy from coefficients in separate containers
+     * for each factor
+     */
+    virtual RT
+    energy (const typename Basis::LocalView& localView,
+            const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
 
-    };
+    /** Empty virtual default destructor
+     *
+     * To allow proper destruction of derived classes through a base class pointer
+     */
+    virtual ~LocalEnergy() = default;
 
-  } // namespace GFE
+  };
 
-}  // namespace Dune
+}  // namespace Dune::GFE
 
 #endif  // DUNE_GFE_LOCALENERGY_HH

dune/gfe/assemblers/localfirstordermodel.hh

@@ -20,6 +20,6 @@ namespace Dune::GFE
 
   };
 
-}  // namespace Dune
+}  // namespace Dune::GFE
 
 #endif   // DUNE_GFE_LOCALFIRSTORDERMODEL_HH

dune/gfe/assemblers/localgeodesicfestiffness.hh

@@ -64,29 +64,30 @@ namespace Dune::GFE
       // Type of the local Hessian
       using CompositeHessian = FieldMatrix<Matrix<double>,2,2>;
     };
-  }
-}
+  }  // namespace Impl
 
-template<class Basis, class TargetSpace>
-class LocalGeodesicFEStiffness
-  : public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
-{
-public:
-
-  /** \brief Assemble the local gradient and stiffness matrix at the current position
-
-   */
-  virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                                          const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
-                                          std::vector<double>& localGradient,
-                                          typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
-
-  /** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
-   */
-  virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                                          const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
-                                          std::vector<double>& localGradient,
-                                          typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
-};
+  template<class Basis, class TargetSpace>
+  class LocalGeodesicFEStiffness
+    : public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
+  {
+  public:
+
+    /** \brief Assemble the local gradient and stiffness matrix at the current position
+
+     */
+    virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
+                                            const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
+                                            std::vector<double>& localGradient,
+                                            typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
+
+    /** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
+     */
+    virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
+                                            const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
+                                            std::vector<double>& localGradient,
+                                            typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
+  };
+
+}  // namespace Dune::GFE
 
 #endif

dune/gfe/assemblers/localintegralenergy.hh

@@ -2,6 +2,7 @@
 #define DUNE_GFE_ASSEMBLERS_LOCALINTEGRALENERGY_HH
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
 
 #include <dune/geometry/quadraturerules.hh>
 
@@ -9,8 +10,8 @@
 #include <dune/gfe/densities/localdensity.hh>
 
 
-namespace Dune::GFE {
-
+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
@@ -23,6 +24,7 @@ namespace Dune::GFE {
   {
     using LocalView = typename Basis::LocalView;
     using GridView = typename LocalView::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
     using DT = typename GridView::Grid::ctype;
     using RT = typename GFE::LocalEnergy<Basis,TargetSpace>::RT;
 
@@ -30,10 +32,28 @@ namespace Dune::GFE {
 
   public:
 
-    /** \brief Constructor with a Dune::GFE::LocalDensity
+    /** \brief Constructor from a GFE function as a shared pointer
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param ld The density that will be integrated
+     */
+    LocalIntegralEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+                        const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
+      : localGFEFunction_(localGFEFunction),
+      localDensity_(ld)
+    {}
+
+    /** \brief Constructor from a GFE function r-value reference
+     *
+     * \param localGFEFunction The geometric finite element function that
+     * the density will be evaluated on
+     * \param ld The density that will be integrated
      */
-    LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> >& ld)
-      : localDensity_(ld)
+    LocalIntegralEnergy(LocalInterpolationRule&& localGFEFunction,
+                        const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
+      : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction))),
+      localDensity_(ld)
     {}
 
   private:
@@ -44,18 +64,22 @@ namespace Dune::GFE {
     {
       RT energy = 0;
 
-      if constexpr (not Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
+      if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower)
       {
-        const auto& localFiniteElement = localView.tree().finiteElement();
-        LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
+        // Bind GFE function to the current element
+        const auto& element = localView.element();
+        localGFEFunction_->bind(element,localConfiguration);
 
-        int quadOrder = (localFiniteElement.type().isSimplex())
-           ? (localFiniteElement.localBasis().order()-1) * 2
-           : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
+        // Bind density to the element
+        localDensity_->bind(element);
 
-        const auto element = localView.element();
+        // Get a suitable quadrature rule
+        const auto localGFEOrder = localGFEFunction_->localFiniteElement().localBasis().order();
+        int quadOrder = (element.type().isSimplex())
+           ? (localGFEOrder-1) * 2
+           : (localGFEOrder * gridDim - 1) * 2;
 
-        const auto& quad = QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
+        const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
 
         for (auto&& qp : quad)
         {
@@ -70,13 +94,13 @@ namespace Dune::GFE {
           {
             if (localDensity_->dependsOnDerivative())
             {
-              auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos);
+              auto [value, derivative] = localGFEFunction_->evaluateValueAndDerivative(quadPos);
               derivative = derivative * geometryJacobianInverse;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
             }
             else
             {
-              auto value = localInterpolationRule.evaluate(quadPos);
+              const auto value = localGFEFunction_->evaluate(quadPos);
               typename LocalInterpolationRule::DerivativeType dummyDerivative;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
             }
@@ -86,7 +110,7 @@ namespace Dune::GFE {
             if (localDensity_->dependsOnDerivative())
             {
               typename TargetSpace::CoordinateType dummyValue;
-              auto derivative = localInterpolationRule.evaluateDerivative(quadPos);
+              auto derivative = localGFEFunction_->evaluateDerivative(quadPos);
               derivative = derivative * geometryJacobianInverse;
               energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
             }
@@ -97,6 +121,11 @@ namespace Dune::GFE {
           }
         }
       }
+      else
+      {
+        // You need a scalar basis or a power basis when calling this method.
+        std::abort();
+      }
 
       return energy;
     }
@@ -106,27 +135,26 @@ namespace Dune::GFE {
     {
       RT energy = 0;
 
-      if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
+      if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold
+                    && Basis::LocalView::Tree::isComposite
+                    && gridDim==GridView::dimensionworld) // TODO: Implement the case gridDim!=dimworld
       {
         static_assert(TargetSpace::size() == 2,
                       "LocalIntegralEnergy only implemented for product spaces with two factors!");
 
-        const auto& element = localView.element();
-
         using namespace Indices;
 
-        // composite Basis: grab the finite element of the first child
-        const auto& localFiniteElement0 = localView.tree().child(_0,0).finiteElement();
-        const auto& localFiniteElement1 = localView.tree().child(_1,0).finiteElement();
-
-        using LocalGFEFunctionType0 = typename std::tuple_element<0, LocalInterpolationRule>::type;
-        using LocalGFEFunctionType1 = typename std::tuple_element<1, LocalInterpolationRule>::type;
+        // Bind GFE functions to the current element
+        const auto& element = localView.element();
+        std::get<0>(*localGFEFunction_).bind(element, coefficients[_0]);
+        std::get<1>(*localGFEFunction_).bind(element, coefficients[_1]);
 
-        LocalGFEFunctionType0 localGFEFunction0(localFiniteElement0, coefficients[_0]);
-        LocalGFEFunctionType1 localGFEFunction1(localFiniteElement1, coefficients[_1]);
+        // Bind density to the element
+        localDensity_->bind(element);
 
-        int quadOrder = (element.type().isSimplex()) ? localFiniteElement0.localBasis().order()
-                                                 : localFiniteElement0.localBasis().order() * gridDim;
+        // Get a suitable quadrature rule
+        const auto feOrder = localView.tree().child(_0,0).finiteElement().localBasis().order();
+        int quadOrder = (element.type().isSimplex()) ? feOrder : feOrder * gridDim;
 
         const auto& quad = QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
 
@@ -147,19 +175,19 @@ namespace Dune::GFE {
           // A value is needed either directly, or for computing the derivative.
           TargetSpace value;
           if (localDensity_->dependsOnValue(0) || localDensity_->dependsOnDerivative(0))
-            value[_0] = localGFEFunction0.evaluate(quadPos);
+            value[_0] = std::get<0>(*localGFEFunction_).evaluate(quadPos);
           if (localDensity_->dependsOnValue(1) || localDensity_->dependsOnDerivative(1))
-            value[_1] = localGFEFunction1.evaluate(quadPos);
+            value[_1] = std::get<1>(*localGFEFunction_).evaluate(quadPos);
 
           // Compute the derivatives of the interpolation function factors
-          typename LocalGFEFunctionType0::DerivativeType derivative0;
-          typename LocalGFEFunctionType1::DerivativeType derivative1;
+          typename std::tuple_element_t<0, LocalInterpolationRule>::DerivativeType derivative0;
+          typename std::tuple_element_t<1, LocalInterpolationRule>::DerivativeType derivative1;
 
           if (localDensity_->dependsOnDerivative(0))
-            derivative0 = localGFEFunction0.evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
+            derivative0 = std::get<0>(*localGFEFunction_).evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
 
           if (localDensity_->dependsOnDerivative(1))
-            derivative1 = localGFEFunction1.evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
+            derivative1 = std::get<1>(*localGFEFunction_).evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
 
           // Copy the two derivatives into a joint matrix object
           // TODO: I am not sure about this.  May the densities should get the
@@ -177,12 +205,19 @@ namespace Dune::GFE {
                                                                     derivative);
         }
       }
+      else
+        DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis or gridDim!=dimworld");
 
       return energy;
     }
 
   protected:
-    const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> > localDensity_;
+
+    // The value and derivative of this function are evaluated at the quadrature points,
+    // and given to the density.
+    const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
+
+    const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> > localDensity_;
   };
 
 }  // namespace Dune::GFE