Make compile with the new dune-functions FE bases

[[Imported from SVN: r10101]]

dune/gfe/mixedcosseratenergy.hh

@@ -20,14 +20,16 @@
 //#define QUADRATIC_MEMBRANE_ENERGY
 
 
-template<class GridView, class DisplacementLocalFiniteElement, class OrientationLocalFiniteElement, int dim, class field_type=double>
+template<class DisplacementBasis, class OrientationBasis, int dim, class field_type=double>
 class MixedCosseratEnergy
-    : public MixedLocalGeodesicFEStiffness<GridView,
-                                           DisplacementLocalFiniteElement,RealTuple<field_type,dim>,
-                                           OrientationLocalFiniteElement,Rotation<field_type,dim> >
+    : public MixedLocalGeodesicFEStiffness<DisplacementBasis,RealTuple<field_type,dim>,
+                                           OrientationBasis,Rotation<field_type,dim> >
 {
     // grid types
-    typedef typename GridView::Grid::ctype DT;
+    typedef typename DisplacementBasis::LocalView::Tree::FiniteElement DisplacementLocalFiniteElement;
+    typedef typename OrientationBasis::LocalView::Tree::FiniteElement OrientationLocalFiniteElement;
+    typedef typename DisplacementBasis::GridView GridView;
+    typedef typename GridView::ctype DT;
     typedef field_type RT;
     typedef typename GridView::template Codim<0>::Entity Entity;
 
@@ -263,11 +265,11 @@ public:
     const Dune::VirtualFunction<Dune::FieldVector<double,gridDim>, Dune::FieldVector<double,3> >* neumannFunction_;
 };
 
-template <class GridView, class DeformationLocalFiniteElement, class OrientationLocalFiniteElement, int dim, class field_type>
-typename MixedCosseratEnergy<GridView,DeformationLocalFiniteElement,OrientationLocalFiniteElement,dim,field_type>::RT
-MixedCosseratEnergy<GridView,DeformationLocalFiniteElement,OrientationLocalFiniteElement,dim,field_type>::
+template <class DeformationBasis, class OrientationBasis, int dim, class field_type>
+typename MixedCosseratEnergy<DeformationBasis,OrientationBasis,dim,field_type>::RT
+MixedCosseratEnergy<DeformationBasis,OrientationBasis,dim,field_type>::
 energy(const Entity& element,
-       const DeformationLocalFiniteElement& deformationLocalFiniteElement,
+       const DisplacementLocalFiniteElement& deformationLocalFiniteElement,
        const std::vector<RealTuple<field_type,dim> >& localDeformationConfiguration,
        const OrientationLocalFiniteElement& orientationLocalFiniteElement,
        const std::vector<Rotation<field_type,dim> >& localOrientationConfiguration) const
@@ -278,7 +280,7 @@ energy(const Entity& element,
 
     RT energy = 0;
 
-    typedef LocalGeodesicFEFunction<gridDim, DT, DeformationLocalFiniteElement, RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
+    typedef LocalGeodesicFEFunction<gridDim, DT, DisplacementLocalFiniteElement, RealTuple<field_type,dim> > LocalDeformationGFEFunctionType;
     LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localDeformationConfiguration);
 
     typedef LocalGeodesicFEFunction<gridDim, DT, OrientationLocalFiniteElement, Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
@@ -363,20 +365,20 @@ energy(const Entity& element,
     if (not neumannFunction_)
         return energy;
 
-    for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) {
-
-        if (not neumannBoundary_ or not neumannBoundary_->contains(*it))
+    for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
+    {
+        if (not neumannBoundary_ or not neumannBoundary_->contains(it))
             continue;
 
         const Dune::QuadratureRule<DT, gridDim-1>& quad
-            = Dune::QuadratureRules<DT, gridDim-1>::rule(it->type(), quadOrder);
+            = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
 
         for (size_t pt=0; pt<quad.size(); pt++) {
 
             // Local position of the quadrature point
-            const Dune::FieldVector<DT,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position());
+            const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
 
-            const DT integrationElement = it->geometry().integrationElement(quad[pt].position());
+            const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
 
             // The value of the local function
             RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
@@ -387,7 +389,7 @@ energy(const Entity& element,
             if (dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_))
                 dynamic_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> >*>(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
             else
-                neumannFunction_->evaluate(it->geometry().global(quad[pt].position()), neumannValue);
+                neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
 
             // Only translational dofs are affected by the Neumann force
             for (size_t i=0; i<neumannValue.size(); i++)
@@ -400,5 +402,5 @@ energy(const Entity& element,
     return energy;
 }
 
-#endif   //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
+#endif   //#ifndef DUNE_GFE_MIXEDCOSSERATENERGY_HH
 

dune/gfe/mixedgfeassembler.hh

@@ -35,20 +35,16 @@ public:
     const Basis0 basis0_;
     const Basis1 basis1_;
 
-    MixedLocalGeodesicFEStiffness<GridView,
-                                  typename Basis0::LocalFiniteElement,
-                                  TargetSpace0,
-                                  typename Basis1::LocalFiniteElement,
-                                  TargetSpace1>* localStiffness_;
+    MixedLocalGeodesicFEStiffness<Basis0, TargetSpace0,
+                                  Basis1, TargetSpace1>* localStiffness_;
 
 public:
 
     /** \brief Constructor for a given grid */
     MixedGFEAssembler(const Basis0& basis0,
                       const Basis1& basis1,
-                      MixedLocalGeodesicFEStiffness<GridView,
-                                               typename Basis0::LocalFiniteElement, TargetSpace0,
-                                               typename Basis0::LocalFiniteElement, TargetSpace1>* localStiffness)
+                      MixedLocalGeodesicFEStiffness<Basis0, TargetSpace0,
+                                                    Basis1, TargetSpace1>* localStiffness)
         : basis0_(basis0),
           basis1_(basis1),
           localStiffness_(localStiffness)
@@ -94,50 +90,57 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
                  Dune::MatrixIndexSet& nb10,
                  Dune::MatrixIndexSet& nb11) const
 {
-    nb00.resize(basis0_.size(), basis0_.size());
-    nb01.resize(basis0_.size(), basis1_.size());
-    nb10.resize(basis1_.size(), basis0_.size());
-    nb11.resize(basis1_.size(), basis1_.size());
+    nb00.resize(basis0_.indexSet().size(), basis0_.indexSet().size());
+    nb01.resize(basis0_.indexSet().size(), basis1_.indexSet().size());
+    nb10.resize(basis1_.indexSet().size(), basis0_.indexSet().size());
+    nb11.resize(basis1_.indexSet().size(), basis1_.indexSet().size());
+
+    // A view on the FE basis on a single element
+    typename Basis0::LocalView localView0(&basis0_);
+    typename Basis1::LocalView localView1(&basis1_);
+    auto localIndexSet0 = basis0_.indexSet().localIndexSet();
+    auto localIndexSet1 = basis1_.indexSet().localIndexSet();
 
     // Grid view must be the same for both bases
-    ElementIterator it    = basis0_.getGridView().template begin<0,Dune::Interior_Partition>();
-    ElementIterator endit = basis0_.getGridView().template end<0,Dune::Interior_Partition>  ();
+    ElementIterator it    = basis0_.gridView().template begin<0,Dune::Interior_Partition>();
+    ElementIterator endit = basis0_.gridView().template end<0,Dune::Interior_Partition>  ();
 
     for (; it!=endit; ++it) {
 
-        const typename Basis0::LocalFiniteElement& lfe0 = basis0_.getLocalFiniteElement(*it);
-        const typename Basis1::LocalFiniteElement& lfe1 = basis1_.getLocalFiniteElement(*it);
+        // Bind the local FE basis view to the current element
+        localView0.bind(*it);
+        localView1.bind(*it);
+        localIndexSet0.bind(localView0);
+        localIndexSet1.bind(localView1);
 
-        for (size_t i=0; i<lfe0.localBasis().size(); i++) {
+        for (size_t i=0; i<localView0.size(); i++) {
 
-            int iIdx = basis0_.index(*it,i);
+            int iIdx = localIndexSet0.index(i)[0];
 
-            for (size_t j=0; j<lfe0.localBasis().size(); j++) {
-                int jIdx = basis0_.index(*it,j);
+            for (size_t j=0; j<localView0.size(); j++) {
+                int jIdx = localIndexSet0.index(j)[0];
                 nb00.add(iIdx, jIdx);
             }
 
-            for (size_t j=0; j<lfe1.localBasis().size(); j++) {
-                int jIdx = basis1_.index(*it,j);
+            for (size_t j=0; j<localView1.size(); j++) {
+                int jIdx = localIndexSet1.index(j)[0];
                 nb01.add(iIdx, jIdx);
-
             }
 
         }
 
-        for (size_t i=0; i<lfe1.localBasis().size(); i++) {
+        for (size_t i=0; i<localView1.size(); i++) {
 
-            int iIdx = basis1_.index(*it,i);
+            int iIdx = localIndexSet1.index(i)[0];
 
-            for (size_t j=0; j<lfe0.localBasis().size(); j++) {
-                int jIdx = basis0_.index(*it,j);
+            for (size_t j=0; j<localView0.size(); j++) {
+                int jIdx = localIndexSet0.index(j)[0];
                 nb10.add(iIdx, jIdx);
             }
 
-            for (size_t j=0; j<lfe1.localBasis().size(); j++) {
-                int jIdx = basis1_.index(*it,j);
+            for (size_t j=0; j<localView1.size(); j++) {
+                int jIdx = localIndexSet1.index(j)[0];
                 nb11.add(iIdx, jIdx);
-
             }
 
         }
@@ -184,70 +187,82 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
     gradient1.resize(configuration1.size());
     gradient1 = 0;
 
-    ElementIterator it    = basis0_.getGridView().template begin<0,Dune::Interior_Partition>();
-    ElementIterator endit = basis0_.getGridView().template end<0,Dune::Interior_Partition>  ();
+    // A view on the FE basis on a single element
+    typename Basis0::LocalView localView0(&basis0_);
+    typename Basis1::LocalView localView1(&basis1_);
+    auto localIndexSet0 = basis0_.indexSet().localIndexSet();
+    auto localIndexSet1 = basis1_.indexSet().localIndexSet();
+
+    ElementIterator it    = basis0_.gridView().template begin<0,Dune::Interior_Partition>();
+    ElementIterator endit = basis0_.gridView().template end<0,Dune::Interior_Partition>  ();
 
     for( ; it != endit; ++it ) {
 
-        const int nDofs0 = basis0_.getLocalFiniteElement(*it).localBasis().size();
-        const int nDofs1 = basis1_.getLocalFiniteElement(*it).localBasis().size();
+        // Bind the local FE basis view to the current element
+        localView0.bind(*it);
+        localView1.bind(*it);
+        localIndexSet0.bind(localView0);
+        localIndexSet1.bind(localView1);
+
+        const int nDofs0 = localView0.size();
+        const int nDofs1 = localView1.size();
 
         // Extract local solution
         std::vector<TargetSpace0> localConfiguration0(nDofs0);
         std::vector<TargetSpace1> localConfiguration1(nDofs1);
 
         for (int i=0; i<nDofs0; i++)
-            localConfiguration0[i] = configuration0[basis0_.index(*it,i)];
+            localConfiguration0[i] = configuration0[localIndexSet0.index(i)[0]];
 
         for (int i=0; i<nDofs1; i++)
-            localConfiguration1[i] = configuration1[basis1_.index(*it,i)];
+            localConfiguration1[i] = configuration1[localIndexSet1.index(i)[0]];
 
         std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
         std::vector<Dune::FieldVector<double,blocksize1> > localGradient1(nDofs1);
 
         // setup local matrix and gradient
         localStiffness_->assembleGradientAndHessian(*it,
-                                                    basis0_.getLocalFiniteElement(*it), localConfiguration0,
-                                                    basis1_.getLocalFiniteElement(*it), localConfiguration1,
+                                                    localView0.tree().finiteElement(), localConfiguration0,
+                                                    localView1.tree().finiteElement(), localConfiguration1,
                                                     localGradient0, localGradient1);
 
         // Add element matrix to global stiffness matrix
         for (int i=0; i<nDofs0; i++) {
 
-            int row = basis0_.index(*it,i);
+            int row = localIndexSet0.index(i)[0];
 
             for (int j=0; j<nDofs0; j++ ) {
-                int col = basis0_.index(*it,j);
+                int col = localIndexSet0.index(j)[0];
                 hessian00[row][col] += localStiffness_->A00_[i][j];
             }
 
             for (int j=0; j<nDofs1; j++ ) {
-                int col = basis1_.index(*it,j);
+                int col = localIndexSet1.index(j)[0];
                 hessian01[row][col] += localStiffness_->A01_[i][j];
             }
         }
 
         for (int i=0; i<nDofs1; i++) {
 
-            int row = basis1_.index(*it,i);
+            int row = localIndexSet1.index(i)[0];
 
             for (int j=0; j<nDofs0; j++ ) {
-                int col = basis0_.index(*it,j);
+                int col = localIndexSet0.index(j)[0];
                 hessian10[row][col] += localStiffness_->A10_[i][j];
             }
 
             for (int j=0; j<nDofs1; j++ ) {
-                int col = basis1_.index(*it,j);
+                int col = localIndexSet1.index(j)[0];
                 hessian11[row][col] += localStiffness_->A11_[i][j];
             }
         }
 
         // Add local gradient to global gradient
         for (int i=0; i<nDofs0; i++)
-            gradient0[basis0_.index(*it,i)] += localGradient0[i];
+            gradient0[localIndexSet0.index(i)[0]] += localGradient0[i];
 
         for (int i=0; i<nDofs1; i++)
-            gradient1[basis1_.index(*it,i)] += localGradient1[i];
+            gradient1[localIndexSet1.index(i)[0]] += localGradient1[i];
     }
 
 }
@@ -300,34 +315,46 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
 {
     double energy = 0;
 
-    if (configuration0.size()!=basis0_.size())
-        DUNE_THROW(Dune::Exception, "Configuration vector doesn't match the grid!");
+    if (configuration0.size()!=basis0_.indexSet().size())
+        DUNE_THROW(Dune::Exception, "Configuration vector 0 doesn't match the basis!");
+
+    if (configuration1.size()!=basis1_.indexSet().size())
+        DUNE_THROW(Dune::Exception, "Configuration vector 1 doesn't match the basis!");
 
-    if (configuration1.size()!=basis1_.size())
-        DUNE_THROW(Dune::Exception, "Configuration vector doesn't match the grid!");
+    // A view on the FE basis on a single element
+    typename Basis0::LocalView localView0(&basis0_);
+    typename Basis1::LocalView localView1(&basis1_);
+    auto localIndexSet0 = basis0_.indexSet().localIndexSet();
+    auto localIndexSet1 = basis1_.indexSet().localIndexSet();
 
-    ElementIterator it    = basis0_.getGridView().template begin<0,Dune::Interior_Partition>();
-    ElementIterator endIt = basis0_.getGridView().template end<0,Dune::Interior_Partition>();
+    ElementIterator it    = basis0_.gridView().template begin<0,Dune::Interior_Partition>();
+    ElementIterator endIt = basis0_.gridView().template end<0,Dune::Interior_Partition>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
 
+        // Bind the local FE basis view to the current element
+        localView0.bind(*it);
+        localView1.bind(*it);
+        localIndexSet0.bind(localView0);
+        localIndexSet1.bind(localView1);
+
         // Number of degrees of freedom on this element
-        size_t nDofs0 = basis0_.getLocalFiniteElement(*it).localBasis().size();
-        size_t nDofs1 = basis1_.getLocalFiniteElement(*it).localBasis().size();
+        size_t nDofs0 = localView0.size();
+        size_t nDofs1 = localView1.size();
 
         std::vector<TargetSpace0> localConfiguration0(nDofs0);
         std::vector<TargetSpace1> localConfiguration1(nDofs1);
 
         for (size_t i=0; i<nDofs0; i++)
-            localConfiguration0[i] = configuration0[basis0_.index(*it,i)];
+            localConfiguration0[i] = configuration0[localIndexSet0.index(i)[0]];
 
         for (size_t i=0; i<nDofs1; i++)
-            localConfiguration1[i] = configuration1[basis1_.index(*it,i)];
+            localConfiguration1[i] = configuration1[localIndexSet1.index(i)[0]];
 
         energy += localStiffness_->energy(*it,
-                                          basis0_.getLocalFiniteElement(*it), localConfiguration0,
-                                          basis1_.getLocalFiniteElement(*it), localConfiguration1);
+                                          localView0.tree().finiteElement(), localConfiguration0,
+                                          localView1.tree().finiteElement(), localConfiguration1);
 
     }
 

dune/gfe/mixedlocalgeodesicfestiffness.hh

@@ -7,12 +7,17 @@
 #include <dune/istl/matrix.hh>
 
 
-template<class GridView,
-         class DeformationLocalFiniteElement, class DeformationTargetSpace,
-         class OrientationLocalFiniteElement, class OrientationTargetSpace>
+template<class DeformationBasis, class DeformationTargetSpace,
+         class OrientationBasis, class OrientationTargetSpace>
 class MixedLocalGeodesicFEStiffness
 {
+    static_assert(std::is_same<typename DeformationBasis::GridView, typename OrientationBasis::GridView>::value,
+                  "DeformationBasis and OrientationBasis must be designed on the same GridView!");
+
     // grid types
+    typedef typename DeformationBasis::LocalView::Tree::FiniteElement DeformationLocalFiniteElement;
+    typedef typename OrientationBasis::LocalView::Tree::FiniteElement OrientationLocalFiniteElement;
+    typedef typename DeformationBasis::GridView GridView;
     typedef typename GridView::Grid::ctype DT;
     typedef typename DeformationTargetSpace::ctype RT;
     typedef typename GridView::template Codim<0>::Entity Entity;

dune/gfe/mixedlocalgfeadolcstiffness.hh

@@ -18,15 +18,19 @@
 
 /** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
  */
-template<class GridView,
-         class LocalFiniteElement0, class TargetSpace0,
-         class LocalFiniteElement1, class TargetSpace1>
+template<class Basis0, class TargetSpace0,
+         class Basis1, class TargetSpace1>
 class MixedLocalGFEADOLCStiffness
-    : public MixedLocalGeodesicFEStiffness<GridView,
-                                           LocalFiniteElement0,TargetSpace0,
-                                           LocalFiniteElement1,TargetSpace1>
+    : public MixedLocalGeodesicFEStiffness<Basis0,TargetSpace0,
+                                           Basis1,TargetSpace1>
 {
+    static_assert(std::is_same<typename Basis0::GridView, typename Basis1::GridView>::value,
+                  "Basis0 and Basis1 must be designed on the same GridView!");
+
     // grid types
+    typedef typename Basis0::GridView GridView;
+    typedef typename Basis0::LocalView::Tree::FiniteElement LocalFiniteElement0;
+    typedef typename Basis1::LocalView::Tree::FiniteElement LocalFiniteElement1;
     typedef typename GridView::Grid::ctype DT;
     typedef typename TargetSpace0::ctype RT;
     typedef typename GridView::template Codim<0>::Entity Entity;
@@ -48,9 +52,8 @@ public:
     enum { embeddedBlocksize0 = TargetSpace0::EmbeddedTangentVector::dimension };
     enum { embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension };
 
-    MixedLocalGFEADOLCStiffness(const MixedLocalGeodesicFEStiffness<GridView,
-                                                               LocalFiniteElement0, ATargetSpace0,
-                                                               LocalFiniteElement1, ATargetSpace1>* energy)
+    MixedLocalGFEADOLCStiffness(const MixedLocalGeodesicFEStiffness<Basis0, ATargetSpace0,
+                                                                    Basis1, ATargetSpace1>* energy)
     : localEnergy_(energy)
     {}
 
@@ -82,14 +85,14 @@ public:
                                             std::vector<typename TargetSpace0::TangentVector>& localGradient0,
                                             std::vector<typename TargetSpace1::TangentVector>& localGradient1);
 
-    const MixedLocalGeodesicFEStiffness<GridView, LocalFiniteElement0, ATargetSpace0, LocalFiniteElement1, ATargetSpace1>* localEnergy_;
+    const MixedLocalGeodesicFEStiffness<Basis0, ATargetSpace0, Basis1, ATargetSpace1>* localEnergy_;
 
 };
 
 
-template <class GridView, class LocalFiniteElement0, class TargetSpace0, class LocalFiniteElement1, class TargetSpace1>
-typename MixedLocalGFEADOLCStiffness<GridView, LocalFiniteElement0, TargetSpace0, LocalFiniteElement1, TargetSpace1>::RT
-MixedLocalGFEADOLCStiffness<GridView, LocalFiniteElement0, TargetSpace0, LocalFiniteElement1, TargetSpace1>::
+template <class Basis0, class TargetSpace0, class Basis1, class TargetSpace1>
+typename MixedLocalGFEADOLCStiffness<Basis0, TargetSpace0, Basis1, TargetSpace1>::RT
+MixedLocalGFEADOLCStiffness<Basis0, TargetSpace0, Basis1, TargetSpace1>::
 energy(const Entity& element,
        const LocalFiniteElement0& localFiniteElement0,
        const std::vector<TargetSpace0>& localConfiguration0,
@@ -198,8 +201,8 @@ assembleGradient(const Entity& element,
 //   To compute the Hessian we need to compute the gradient anyway, so we may
 //   as well return it.  This saves assembly time.
 // ///////////////////////////////////////////////////////////
-template <class GridType, class LocalFiniteElement0, class TargetSpace0, class LocalFiniteElement1, class TargetSpace1>
-void MixedLocalGFEADOLCStiffness<GridType, LocalFiniteElement0, TargetSpace0, LocalFiniteElement1, TargetSpace1>::
+template <class Basis0, class TargetSpace0, class Basis1, class TargetSpace1>
+void MixedLocalGFEADOLCStiffness<Basis0, TargetSpace0, Basis1, TargetSpace1>::
 assembleGradientAndHessian(const Entity& element,
                            const LocalFiniteElement0& localFiniteElement0,
                            const std::vector<TargetSpace0>& localConfiguration0,

dune/gfe/mixedriemanniantrsolver.cc

@@ -7,17 +7,16 @@
 #include <dune/istl/io.hh>
 
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
 
 // Using a monotone multigrid as the inner solver
 #include <dune/solvers/iterationsteps/trustregiongsstep.hh>
 #include <dune/solvers/iterationsteps/mmgstep.hh>
 #include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
-#if defined THIRD_ORDER || defined SECOND_ORDER
-#include <dune/gfe/pktop1mgtransfer.hh>
-#endif
 #include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include "maxnormtrustregion.hh"
@@ -135,11 +134,16 @@ setup(const GridType& grid,
     mmgStep1->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType1>();
     mmgStep1->verbosity_         = Solver::FULL;
 
-#if 0
     // //////////////////////////////////////////////////////////////////////////////////////
     //   Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
     // //////////////////////////////////////////////////////////////////////////////////////
+    typedef DuneFunctionsBasis<Basis0> FufemBasis0;
+    FufemBasis0 basis0(grid.leafGridView());
+
+    typedef DuneFunctionsBasis<Basis1> FufemBasis1;
+    FufemBasis1 basis1(grid.leafGridView());
 
+#if 0
     BasisType basis(grid.leafGridView());
     OperatorAssembler<BasisType,BasisType> operatorAssembler(basis, basis);
 
@@ -189,47 +193,27 @@ setup(const GridType& grid,
     mmgStep0->mgTransfer_.resize(numLevels-1);
     mmgStep1->mgTransfer_.resize(numLevels-1);
 
-    if (assembler->basis0_.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
+    if (basis0.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
     {
       if (numLevels>1) {
+        typedef typename TruncatedCompressedMGTransfer<CorrectionType0>::TransferOperatorType TransferOperatorType;
         P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
 
-        PKtoP1MGTransfer<CorrectionType0>* topTransferOp0 = new PKtoP1MGTransfer<CorrectionType0>;
-        topTransferOp0->setup(assembler->basis0_,p1Basis);
-        #if 1
-        mmgStep0->mgTransfer_.back() = topTransferOp0;
-        #else
-        // If we are on more than 1 processors, join all local transfer matrices on rank 0,
-        // and construct a single global transfer operator there.
-        typedef GlobalUniqueIndex<typename GridType::LeafGridView, gridDim> LeafP1GUIndex;
-        LeafP1GUIndex p1Index(grid_->leafGridView());
-
-        typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
-        MatrixCommunicator<GUIndex, TransferOperatorType, LeafP1GUIndex> matrixComm(*guIndex_, p1Index, 0);
-
-        mmgStep->mgTransfer_.back() = new PKtoP1MGTransfer<CorrectionType>
-        (Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(topTransferOp->getMatrix())));
-        #endif
-        for (int i=0; i<mmgStep0->mgTransfer_.size()-1; i++){
+        TransferOperatorType pkToP1TransferMatrix;
+        assembleBasisInterpolationMatrix<TransferOperatorType,
+                                         P1NodalBasis<typename GridType::LeafGridView,double>,
+                                         FufemBasis0>(pkToP1TransferMatrix,p1Basis,assembler->basis0_);
+
+        mmgStep0->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType0>;
+        Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+        dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType0>*>(mmgStep0->mgTransfer_.back())->setMatrix(topTransferOperator);
+
+        for (size_t i=0; i<mmgStep0->mgTransfer_.size()-1; i++){
           // Construct the local multigrid transfer matrix
           TruncatedCompressedMGTransfer<CorrectionType0>* newTransferOp0 = new TruncatedCompressedMGTransfer<CorrectionType0>;
           newTransferOp0->setup(*grid_,i+1,i+2);
 
-          #if 1
           mmgStep0->mgTransfer_[i] = newTransferOp0;
-          #else
-          // If we are on more than 1 processors, join all local transfer matrices on rank 0,
-          // and construct a single global transfer operator there.
-          typedef GlobalUniqueIndex<typename GridType::LevelGridView, gridDim> LevelGUIndex;
-          LevelGUIndex fineGUIndex(grid_->levelGridView(i+2));
-          LevelGUIndex coarseGUIndex(grid_->levelGridView(i+1));
-
-          typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
-          MatrixCommunicator<LevelGUIndex, TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, 0);
-
-          mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>
-          (Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix())));
-          #endif
         }
 
       }
@@ -244,35 +228,27 @@ setup(const GridType& grid,
         TruncatedCompressedMGTransfer<CorrectionType0>* newTransferOp0 = new TruncatedCompressedMGTransfer<CorrectionType0>;
         newTransferOp0->setup(*grid_,i,i+1);
 
-        #if 1
         mmgStep0->mgTransfer_[i] = newTransferOp0;
-        #else
-        // If we are on more than 1 processors, join all local transfer matrices on rank 0,
-        // and construct a single global transfer operator there.
-        typedef GlobalUniqueIndex<typename GridType::LevelGridView, gridDim> LevelGUIndex;
-        LevelGUIndex fineGUIndex(grid_->levelGridView(i+1));
-        LevelGUIndex coarseGUIndex(grid_->levelGridView(i));
-
-        typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
-        MatrixCommunicator<LevelGUIndex, TransferOperatorType> matrixComm(fineGUIndex, coarseGUIndex, 0);
-
-        mmgStep->mgTransfer_[i] = new TruncatedCompressedMGTransfer<CorrectionType>
-        (Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(newTransferOp->getMatrix())));
-        #endif
       }
 
     }
 
-    if (assembler->basis1_.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
+    if (basis1.getLocalFiniteElement(*grid.leafGridView().template begin<0>()).localBasis().order() > 1)
     {
       if (numLevels>1) {
+        typedef typename TruncatedCompressedMGTransfer<CorrectionType1>::TransferOperatorType TransferOperatorType;
         P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
 
-        PKtoP1MGTransfer<CorrectionType1>* topTransferOp1 = new PKtoP1MGTransfer<CorrectionType1>;
-        topTransferOp1->setup(assembler->basis1_,p1Basis);
-        mmgStep1->mgTransfer_.back() = topTransferOp1;
+        TransferOperatorType pkToP1TransferMatrix;
+        assembleBasisInterpolationMatrix<TransferOperatorType,
+                                         P1NodalBasis<typename GridType::LeafGridView,double>,
+                                         FufemBasis1>(pkToP1TransferMatrix,p1Basis,assembler->basis1_);
+
+        mmgStep0->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType1>;
+        Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+        dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType1>*>(mmgStep1->mgTransfer_.back())->setMatrix(topTransferOperator);
 
-        for (int i=0; i<mmgStep1->mgTransfer_.size()-1; i++){
+        for (size_t i=0; i<mmgStep1->mgTransfer_.size()-1; i++){
           // Construct the local multigrid transfer matrix
           TruncatedCompressedMGTransfer<CorrectionType1>* newTransferOp1 = new TruncatedCompressedMGTransfer<CorrectionType1>;
           newTransferOp1->setup(*grid_,i+1,i+2);
@@ -303,8 +279,8 @@ setup(const GridType& grid,
       #if 0
       hasObstacle0_.resize(guIndex_->nGlobalEntity(), true);
       #else
-      hasObstacle0_.resize(assembler->basis0_.size(), true);
-      hasObstacle1_.resize(assembler->basis1_.size(), true);
+      hasObstacle0_.resize(assembler->basis0_.indexSet().size(), true);
+      hasObstacle1_.resize(assembler->basis1_.indexSet().size(), true);
       #endif
       mmgStep0->hasObstacle_ = &hasObstacle0_;
       mmgStep1->hasObstacle_ = &hasObstacle1_;
@@ -323,15 +299,9 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
     Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc,argv);
     int rank = grid_->comm().rank();
 
-    MonotoneMGStep<MatrixType00,CorrectionType0>* mgStep0 = nullptr;
-
-    // if the inner solver is a monotone multigrid set up a max-norm trust-region
-    if (dynamic_cast<LoopSolver<CorrectionType0>*>(innerSolver_.get()))
-        mgStep0 = dynamic_cast<MonotoneMGStep<MatrixType00,CorrectionType0>*>(dynamic_cast<LoopSolver<CorrectionType0>*>(innerSolver_.get())->iterationStep_);
-
     // \todo Use global index set instead of basis for parallel computations
-    MaxNormTrustRegion<blocksize0> trustRegion0(assembler_->basis0_.size(), initialTrustRegionRadius_);
-    MaxNormTrustRegion<blocksize1> trustRegion1(assembler_->basis1_.size(), initialTrustRegionRadius_);
+    MaxNormTrustRegion<blocksize0> trustRegion0(assembler_->basis0_.indexSet().size(), initialTrustRegionRadius_);
+    MaxNormTrustRegion<blocksize1> trustRegion1(assembler_->basis1_.indexSet().size(), initialTrustRegionRadius_);
 
     std::vector<BoxConstraint<field_type,blocksize0> > trustRegionObstacles0;
     std::vector<BoxConstraint<field_type,blocksize1> > trustRegionObstacles1;
@@ -607,10 +577,12 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis0,TargetSpace0,Basis1,Target
         }
 
         // Output each iterate, to better understand what the algorithm does
+        DuneFunctionsBasis<Basis0> fufemBasis0(assembler_->basis0_);
+        DuneFunctionsBasis<Basis1> fufemBasis1(assembler_->basis1_);
         std::stringstream iAsAscii;
         iAsAscii << i+1;
-        CosseratVTKWriter<GridType>::template writeMixed<Basis0,Basis1>(assembler_->basis0_,x0_,
-                                                                        assembler_->basis1_,x1_,
+        CosseratVTKWriter<GridType>::template writeMixed<DuneFunctionsBasis<Basis0>, DuneFunctionsBasis<Basis1> >(fufemBasis0,x0_,
+                                                                        fufemBasis1,x1_,
                                                                         "mixed-cosserat_iterate_" + iAsAscii.str());
 
         if (rank==0)

dune/gfe/rodassembler.cc

@@ -11,25 +11,32 @@
 
 
 
-template <class GridView>
-void RodAssembler<GridView,3>::
+template <class Basis>
+void RodAssembler<Basis,3>::
 assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
                  Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
 {
     using namespace Dune;
 
-    if (sol.size()!=this->basis_.size())
+    if (sol.size()!=this->basis_.indexSet().size())
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
 
     grad.resize(sol.size());
     grad = 0;
 
-    ElementIterator it    = this->basis_.getGridView().template begin<0>();
-    ElementIterator endIt = this->basis_.getGridView().template end<0>();
+    // A view on the FE basis on a single element
+    typename Basis::LocalView localView(&this->basis_);
+    auto localIndexSet = this->basis_.indexSet().localIndexSet();
+
+    ElementIterator it    = this->basis_.gridView().template begin<0>();
+    ElementIterator endIt = this->basis_.gridView().template end<0>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
 
+        localView.bind(*it);
+        localIndexSet.bind(localView);
+
         // A 1d grid has two vertices
         static const int nDofs = 2;
 
@@ -37,19 +44,19 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
         std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
 
         for (int i=0; i<nDofs; i++)
-            localSolution[i] = sol[this->basis_.index(*it,i)];
+            localSolution[i] = sol[localIndexSet.index(i)[0]];
 
         // Assemble local gradient
         std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
 
         this->localStiffness_->assembleGradient(*it,
-                                                this->basis_.getLocalFiniteElement(*it),
+                                                localView.tree().finiteElement(),
                                                 localSolution,
                                                 localGradient);
 
         // Add to global gradient
         for (int i=0; i<nDofs; i++)
-            grad[this->basis_.index(*it,i)] += localGradient[i];
+            grad[localIndexSet.index(i)[0]] += localGradient[i];
 
     }
 

dune/gfe/rodassembler.hh

@@ -14,7 +14,7 @@
 
 /** \brief The FEM operator for an extensible, shearable rod in 3d
  */
-template <class GridView, int spaceDim>
+template <class Basis, int spaceDim>
 class RodAssembler
 {
     static_assert(spaceDim==2 || spaceDim==3,
@@ -23,12 +23,11 @@ class RodAssembler
 
 /** \brief The FEM operator for an extensible, shearable rod in 3d
  */
-template <class GridView>
-class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,3> >
+template <class Basis>
+class RodAssembler<Basis,3> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >
 {
+  typedef typename Basis::GridView GridView;
 
-    //typedef typename GridType::template Codim<0>::Entity EntityType;
-    //typedef typename GridType::template Codim<0>::EntityPointer EntityPointer;
     typedef typename GridView::template Codim<0>::Iterator ElementIterator;
 
         //! Dimension of the grid.  This needs to be one!
@@ -44,18 +43,18 @@ class RodAssembler<GridView,3> : public GeodesicFEAssembler<P1NodalBasis<GridVie
 
 public:
         //! ???
-    RodAssembler(const GridView &gridView,
+    RodAssembler(const Basis& basis,
                  RodLocalStiffness<GridView,double>* localStiffness)
-        : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,3> >(gridView,localStiffness)
+        : GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,localStiffness)
         {
-            std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(gridDim));
+            std::vector<RigidBodyMotion<double,3> > referenceConfiguration(basis.indexSet().size());
 
-            typename GridView::template Codim<gridDim>::Iterator it    = gridView.template begin<gridDim>();
-            typename GridView::template Codim<gridDim>::Iterator endIt = gridView.template end<gridDim>();
+            auto it    = basis.gridView().template begin<gridDim>();
+            auto endIt = basis.gridView().template end<gridDim>();
 
             for (; it != endIt; ++it) {
 
-                int idx = gridView.indexSet().index(*it);
+                int idx = basis.gridView().indexSet().index(*it);
 
                 referenceConfiguration[idx].r[0] = 0;
                 referenceConfiguration[idx].r[1] = 0;
@@ -91,10 +90,11 @@ public:
 
 /** \brief The FEM operator for a 2D extensible, shearable rod
  */
-template <class GridView>
-class RodAssembler<GridView,2> : public GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,2> >
+template <class Basis>
+class RodAssembler<Basis,2> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >
 {
 
+    typedef typename Basis::GridView GridView;
     typedef typename GridView::template Codim<0>::Entity EntityType;
     typedef typename GridView::template Codim<0>::Iterator ElementIterator;
 
@@ -118,7 +118,7 @@ public:
 
     //! ???
     RodAssembler(const GridView &gridView)
-        : GeodesicFEAssembler<P1NodalBasis<GridView>, RigidBodyMotion<double,2> >(gridView,NULL)
+        : GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >(gridView,NULL)
     {
         B = 1;
         A1 = 1;

dune/gfe/rodlocalstiffness.hh

@@ -4,16 +4,15 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/matrix.hh>
 #include <dune/geometry/quadraturerules.hh>
-#include <dune/localfunctions/lagrange/p1.hh>
 
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
 
 #include "localgeodesicfestiffness.hh"
 #include "rigidbodymotion.hh"
 
 template<class GridView, class RT>
 class RodLocalStiffness
-    : public LocalGeodesicFEStiffness<GridView, typename P1NodalBasis<GridView>::LocalFiniteElement, RigidBodyMotion<RT,3> >
+    : public LocalGeodesicFEStiffness<Dune::Functions::PQKNodalBasis<GridView,1>, RigidBodyMotion<RT,3> >
 {
     typedef RigidBodyMotion<RT,3> TargetSpace;
 
@@ -102,7 +101,7 @@ public:
                        const Dune::array<RigidBodyMotion<RT,3>, dim+1>& localSolution) const;
 
     virtual RT energy (const Entity& e,
-                       const typename P1NodalBasis<GridView>::LocalFiniteElement& localFiniteElement,
+                       const typename Dune::Functions::PQKNodalBasis<GridView,1>::LocalView::Tree::FiniteElement& localFiniteElement,
                        const std::vector<RigidBodyMotion<RT,3> >& localSolution) const
     {
         assert(localSolution.size()==2);