Cosserat-Continuum-Nonplanar

@osander: Can you have a look here?

This MR cleans up the usage of MIXED_SPACE in cosserat-continuum.

Now all combinations of dim = dimworld or dim != dimworld and MIXED_SPACE = 0 or MIXED_SPACE = 1 at least compile.

Also, this MR points out the one combination that still can't be solved: Problems with dim != dimworld do not work for MIXED_SPACE = 1!

As soon as I need this for locking tests I will add this; I will have to change Nonplanarcosseratshellenergy so it can be assembled with the GFEMixedAssembler.

dune/gfe/nonplanarcosseratshellenergy.hh

@@ -9,6 +9,8 @@
 
 #include <dune/fufem/boundarypatch.hh>
 
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
 #include <dune/gfe/localenergy.hh>
 #include <dune/gfe/localgeodesicfefunction.hh>
 #include <dune/gfe/rigidbodymotion.hh>
@@ -21,7 +23,8 @@
 #include <dune/localfunctions/lagrange/lfecache.hh>
 #endif
 
-template<class Basis, int dim, class field_type=double>
+template<class Basis, int dim, class field_type=double, class GridFunction = 
+  Dune::Functions::DiscreteGlobalBasisFunction<Basis,std::vector<Dune::FieldVector<double, Basis::GridView::dimensionworld>> > >
 class NonplanarCosseratShellEnergy
   : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >
 {
@@ -43,10 +46,12 @@ public:
    * \param parameters The material parameters
    */
   NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
+                               const GridFunction* gridFunctionPtr,
                                const BoundaryPatch<GridView>* neumannBoundary,
                                const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
                                const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
-  : neumannBoundary_(neumannBoundary),
+  : gridFunctionPtr_(gridFunctionPtr),
+    neumannBoundary_(neumannBoundary),
     neumannFunction_(neumannFunction),
     volumeLoad_(volumeLoad)
   {
@@ -111,6 +116,9 @@ public:
   /** \brief Curvature parameters */
   double b1_, b2_, b3_;
 
+  /** \brief The geometry used for assembling */
+  const GridFunction* gridFunctionPtr_;
+
   /** \brief The Neumann boundary */
   const BoundaryPatch<GridView>* neumannBoundary_;
 
@@ -121,9 +129,9 @@ public:
   const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
 };
 
-template <class Basis, int dim, class field_type>
-typename NonplanarCosseratShellEnergy<Basis,dim,field_type>::RT
-NonplanarCosseratShellEnergy<Basis,dim,field_type>::
+template <class Basis, int dim, class field_type, class GridFunction>
+typename NonplanarCosseratShellEnergy<Basis, dim, field_type, GridFunction>::RT
+NonplanarCosseratShellEnergy<Basis,dim,field_type, GridFunction>::
 energy(const typename Basis::LocalView& localView,
        const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
 {
@@ -220,9 +228,14 @@ energy(const typename Basis::LocalView& localView,
     // and localGeometry.global maps them to the world coordinates
     // we want to take the derivative of the normal field on the element in world coordinates
     Dune::CurvedGeometry<DT, gridDim, dimworld, Dune::CurvedGeometryTraits<DT, Dune::LagrangeLFECache<DT,DT,gridDim>>> curvedGeometry(referenceElement(element),
-      [localGeometry=element.geometry()](const auto& local) {
-        return localGeometry.global(local);
-      }, 1); //order = 1
+      [this,element](const auto& local) {
+        if (not gridFunctionPtr_) {
+          return element.geometry().global(local);
+        }
+        auto localGridFunction = localFunction(*gridFunctionPtr_);
+        localGridFunction.bind(element);
+        return localGridFunction(local);
+      }, 2); //order = 2
 
     // Second fundamental form: The derivative of the normal field
     auto normalDerivative = curvedGeometry.normalGradient(quad[pt].position());