diff --git a/dune/fem-dg/algorithm/caller/sub/additionaloutput.hh b/dune/fem-dg/algorithm/caller/sub/additionaloutput.hh
index 41e02603f6562bfc1922923d922279fe55fbf63b..2cb031a635c674105809b562f209807332a29361 100644
--- a/dune/fem-dg/algorithm/caller/sub/additionaloutput.hh
+++ b/dune/fem-dg/algorithm/caller/sub/additionaloutput.hh
@@ -51,17 +51,11 @@ namespace Fem
       typedef typename OutDiscreteFunctionType::DiscreteFunctionSpaceType   OutDiscreteFunctionSpaceType;
 
       typedef typename InDiscreteFunctionSpaceType::GridPartType  GridPartType;
-      typedef typename InDiscreteFunctionSpaceType::IteratorType  Iterator;
-      typedef typename Iterator::Entity                           Entity;
-      typedef typename Entity::Geometry                           Geometry;
       typedef typename InDiscreteFunctionSpaceType::DomainType    DomainType;
       typedef typename InDiscreteFunctionSpaceType::RangeType     InRangeType;
       typedef typename OutDiscreteFunctionSpaceType::RangeType    OutRangeType;
 
-      typedef typename InDiscreteFunctionType::LocalFunctionType  InLocalFuncType;
-      typedef typename OutDiscreteFunctionType::LocalFunctionType OutLocalFuncType;
-
-      const InDiscreteFunctionSpaceType& space =  alg.solution().space();
+      const auto& space =  alg.solution().space();
       solution_.clear();
 
       InRangeType cons(0.0);
@@ -69,31 +63,27 @@ namespace Fem
       OutRangeType prim(0.0);
       OutRangeType prim_bg(0.0);
 
-      const Iterator endit = space.end();
-      for( Iterator it = space.begin(); it != endit ; ++it)
+      for( const auto& entity : elements( space.gridPart() ) )
       {
-        // get entity
-        const Entity& entity = *it ;
-        const Geometry& geo = entity.geometry();
+        const auto& geo = entity.geometry();
 
         // get quadrature rule for L2 projection
         Dune::Fem::CachingQuadrature< GridPartType, 0 > quad( entity, 2*space.order()+1 );
 
-        InLocalFuncType consLF = alg.solution().localFunction( entity );
-        OutLocalFuncType primLF = solution_.localFunction( entity );
+        const auto& consLF = alg.solution().localFunction( entity );
+        const auto& primLF = solution_.localFunction( entity );
 
-        const int quadNop = quad.nop();
-        for(int qP = 0; qP < quadNop; ++qP)
+        for( const auto qp : quad )
         {
-          const DomainType& xgl = geo.global( quad.point(qP) );
-          consLF.evaluate( quad[qP], cons );
+          const auto& xgl = geo.global( qp.position() );
+          consLF.evaluate( qp, cons );
 
           // it is useful to visualize better suited quantities
           bool forVisual = true;
           alg.model().conservativeToPrimitive( tp.time(), xgl, cons, prim, forVisual );
 
-          prim *=  quad.weight(qP);
-          primLF.axpy( quad[qP] , prim );
+          prim *=  qp.weight();
+          primLF.axpy( qp, prim );
         }
       }
     }
diff --git a/dune/fem-dg/algorithm/sub/elliptic.hh b/dune/fem-dg/algorithm/sub/elliptic.hh
index 0b82878d287803a805dbb2c5a859afc9799b1550..3e6fad5cde8ed478bc2dd70c8848002429f2faa6 100644
--- a/dune/fem-dg/algorithm/sub/elliptic.hh
+++ b/dune/fem-dg/algorithm/sub/elliptic.hh
@@ -225,10 +225,8 @@ namespace Fem
 
         localre_.init(entity);
         localre_.clear();
-        IntersectionIteratorType end = df_.space().gridPart().iend( entity );
-        for( IntersectionIteratorType it = df_.space().gridPart().ibegin( entity ); it != end; ++it )
+        for (const auto& intersection : intersections( df_.space().gridPart(), entity ) )
         {
-          const IntersectionType &intersection = *it;
           if ( intersection.neighbor() && df_.space().continuous(intersection) )
           {
             if( ! intersection.conforming() )
@@ -245,7 +243,6 @@ namespace Fem
         // CACHING
         typedef typename Operator::FaceQuadratureType                               FaceQuadratureType ;
         typedef Dune::Fem::IntersectionQuadrature< FaceQuadratureType, conforming > IntersectionQuadratureType;
-        typedef typename IntersectionQuadratureType::FaceQuadratureType             QuadratureImp;
 
         const EntityType &outside = intersection.outside();
 
@@ -257,8 +254,8 @@ namespace Fem
         const int quadOrder = 2 * std::max( enOrder, nbOrder ) + 1;
 
         IntersectionQuadratureType interQuad( df_.space().gridPart(), intersection, quadOrder );
-        const QuadratureImp &quadInside  = interQuad.inside();
-        const QuadratureImp &quadOutside = interQuad.outside();
+        const auto& quadInside  = interQuad.inside();
+        const auto& quadOutside = interQuad.outside();
         const int numQuadraturePoints = quadInside.nop();
 
         // obtain all required function values on intersection
@@ -434,12 +431,8 @@ namespace Fem
       polOrderContainer_.resize();
       if ( errorEstimator_.isPadaptive() )
       {
-        typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
-        typedef typename IteratorType::Entity EntityType ;
-        const IteratorType end = space_.end();
-        for( IteratorType it = space_.begin(); it != end; ++it )
+        for( const auto& entity : elements( space_.gridPart() ) )
         {
-          const EntityType& entity = *it;
           int order = polOrderContainer_[ entity ].value();
           while (order == -1) // is a new element
           {
@@ -470,14 +463,10 @@ namespace Fem
 
       const double error = errorEstimator_.estimate( problem );
       std::cout << "ESTIMATE: " << error << std::endl;
-      typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
-      const IteratorType end = space_.end();
-      for( IteratorType it = space_.begin(); it != end; ++it )
-      {
-        const typename IteratorType::Entity &entity = *it;
-        polOrderContainer_[entity].value() =
-          errorEstimator_.newOrder( 0.98*tolerance, entity );
-      }
+
+      for( const auto& entity : elements( space_.gridPart() ) )
+        polOrderContainer_[entity].value() = errorEstimator_.newOrder( 0.98*tolerance, entity );
+
       return (error < std::abs(tolerance) ? false : errorEstimator_.mark( 0.98 * tolerance));
 #else
       return false;
diff --git a/dune/fem-dg/assemble/primalmatrix.hh b/dune/fem-dg/assemble/primalmatrix.hh
index c66313ac664a79a272776c7f4e675b4853a6e1ab..0c3d8422991b92020cdd6ae8d62c65cffa18638f 100644
--- a/dune/fem-dg/assemble/primalmatrix.hh
+++ b/dune/fem-dg/assemble/primalmatrix.hh
@@ -53,7 +53,7 @@ namespace Fem
     typedef typename DestinationType::LocalFunctionType           LocalFunctionType;
 
     typedef typename GridPartType::IntersectionIteratorType       IntersectionIteratorType;
-    typedef typename IntersectionIteratorType::Intersection       IntersectionType;
+    typedef typename GridPartType::IntersectionType               IntersectionType;
     typedef typename IntersectionType::Geometry                   IntersectionGeometryType;
 
     typedef Fem::Parameter                                        ParameterType;
@@ -258,10 +258,8 @@ namespace Fem
 
       Dune::Fem::TemporaryLocalFunction< DiscreteFunctionSpaceType > rhsLocal( space_ );
 
-      const IteratorType end = space_.end();
-      for( IteratorType it = space_.begin(); it != end; ++it )
+      for( const auto& entity : elements( space_.gridPart() ) )
       {
-        const EntityType &entity = *it;
         const GeometryType &geometry = entity.geometry();
         const double volume = geometry.volume();
 
@@ -277,16 +275,16 @@ namespace Fem
         const unsigned int numBasisFunctionsEn = baseSet.size();
 
         VolumeQuadratureType quadrature( entity, elementQuadOrder( space_.order( entity ) ) );
-        const size_t numQuadraturePoints = quadrature.nop();
-        for( size_t pt = 0; pt < numQuadraturePoints; ++pt )
+
+        for( const auto qp : quadrature )
         {
-          LocalEvaluationType local( entity, quadrature, uZero, uJacZero, pt, time_, volume );
+          LocalEvaluationType local( entity, quadrature, uZero, uJacZero, qp.index(), time_, volume );
 
-          const double weight = quadrature.weight( pt ) * geometry.integrationElement( local.position() );
+          const double weight = qp.weight() * geometry.integrationElement( local.position() );
 
           // resize of phi and dphi is done in evaluate methods
-          baseSet.evaluateAll( quadrature[ pt ], phi );
-          baseSet.jacobianAll( quadrature[ pt ], dphi );
+          baseSet.evaluateAll( qp, phi );
+          baseSet.jacobianAll( qp, dphi );
 
           RangeType arhs(0);
           // first assemble rhs (evaluate source with u=0)
@@ -347,16 +345,12 @@ namespace Fem
             // assemble rhs
             arhs     *=  weight;
             arhsdphi *= -weight;
-            rhsLocal.axpy( quadrature[ pt ], arhs, arhsdphi );
+            rhsLocal.axpy( qp, arhs, arhsdphi );
           }
         }
 
-        const IntersectionIteratorType endiit = space_.gridPart().iend( entity );
-        for ( IntersectionIteratorType iit = space_.gridPart().ibegin( entity );
-              iit != endiit ; ++ iit )
+        for (const auto& intersection : intersections(space_.gridPart(), entity) )
         {
-          const IntersectionType& intersection = *iit ;
-
           if( intersection.neighbor() && calculateFluxes_ )
           {
             if ( space_.continuous(intersection) ) continue;
@@ -378,11 +372,11 @@ namespace Fem
             const size_t numFaceQuadPoints = faceQuadInside.nop();
             resize( numFaceQuadPoints, maxNumBasisFunctions );
 
-            // evalute base functions
-            for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
+            // evaluate base functions
+            for( const auto qp : faceQuadInside )
             {
-              baseSet.evaluateAll( faceQuadInside[ pt ], phiFaceEn[pt] );
-              baseSet.jacobianAll( faceQuadInside[ pt ], dphiFaceEn[pt] );
+              baseSet.evaluateAll( qp, phiFaceEn[qp.index()] );
+              baseSet.jacobianAll( qp, dphiFaceEn[qp.index()] );
             }
 
             // storage for all flux values
@@ -503,10 +497,8 @@ namespace Fem
       const RangeType uZero(0);
       const JacobianRangeType uJacZero(0);
 
-      const IteratorType end = space_.end();
-      for( IteratorType it = space_.begin(); it != end; ++it )
+      for( const auto& entity : elements( space_.gridPart() ) )
       {
-        const EntityType &entity = *it;
         const GeometryType &geometry = entity.geometry();
         const double volume = geometry.volume();
 
@@ -515,11 +507,8 @@ namespace Fem
         const BasisFunctionSetType &baseSet = rhsLocal.baseFunctionSet();
         const unsigned int numBasisFunctionsEn = baseSet.size();
 
-        const IntersectionIteratorType endiit = space_.gridPart().iend( entity );
-        for ( IntersectionIteratorType iit = space_.gridPart().ibegin( entity );
-              iit != endiit ; ++ iit )
+        for (const auto& intersection : intersections(space_.gridPart(), entity) )
         {
-          const IntersectionType& intersection = *iit ;
 
           if( intersection.neighbor() && calculateFluxes_ )
           {
@@ -619,7 +608,6 @@ namespace Fem
       const QuadratureImp &faceQuadInside  = interQuad.inside();
       const QuadratureImp &faceQuadOutside = interQuad.outside();
 
-
       IntersectionStorage intersectionStorage( intersection,
                                                entity, neighbor,
                                                entity.geometry().volume(),
@@ -630,13 +618,14 @@ namespace Fem
       resize( numFaceQuadPoints, maxNumBasisFunctions );
 
       // evalute base functions
-      for( size_t pt = 0; pt < numFaceQuadPoints; ++pt )
+      for( const auto qp : faceQuadInside )
       {
         // resize is done in evaluateAll and jacobianAll
-        baseSet.evaluateAll( faceQuadInside[ pt ], phiFaceEn[pt] );
-        baseSet.jacobianAll( faceQuadInside[ pt ], dphiFaceEn[pt] );
-        baseSetNb.evaluateAll( faceQuadOutside[ pt ], phiFaceNb[pt] );
-        baseSetNb.jacobianAll( faceQuadOutside[ pt ], dphiFaceNb[pt] );
+        const auto idx = qp.index();
+        baseSet.evaluateAll( qp, phiFaceEn[idx] );
+        baseSet.jacobianAll( qp, dphiFaceEn[idx] );
+        baseSetNb.evaluateAll( faceQuadOutside[idx], phiFaceNb[idx] );
+        baseSetNb.jacobianAll( faceQuadOutside[idx], dphiFaceNb[idx] );
       }
 
       flux(dfSpace.gridPart(), intersectionStorage,
@@ -725,10 +714,8 @@ namespace Fem
   //    typedef typename MatrixType::LocalMatrixType LocalMatrixType;
   //    double maxSymError_diagonal = 0;
   //    double maxSymError_offdiagonal = 0;
-  //    const IteratorType end = space().end();
-  //    for( IteratorType it = space().begin(); it != end; ++it )
+  //    for( const auto& entity : elements( space().gridPart() ) )
   //    {
-  //      const EntityType& entity = *it;
   //      LocalMatrixType localOpEn = matrix_->localMatrix( entity, entity  );
   //      const BasisFunctionSetType &baseSet = localOpEn.domainBasisFunctionSet();
   //      const unsigned int numBasisFunctions = baseSet.size();
@@ -744,11 +731,8 @@ namespace Fem
   //          maxSymError_diagonal = std::max(maxSymError_diagonal,std::abs(v1-v2));
   //        }
   //      }
-  //      const IntersectionIteratorType endiit = space().gridPart().iend( entity );
-  //      for ( IntersectionIteratorType iit = space().gridPart().ibegin( entity );
-  //            iit != endiit ; ++ iit )
+  //      for (const auto& intersection : intersections(space().gridPart(), entity) )
   //      {
-  //        const IntersectionType& intersection = *iit ;
   //        if ( !intersection.neighbor() )
   //          continue;
   //        EntityPointerType ep = intersection.outside();
diff --git a/dune/fem-dg/examples/dataio/checkpointing.hh b/dune/fem-dg/examples/dataio/checkpointing.hh
index fb91f80ad615236e0a471f04fd8f1cd9c3d4a01b..22f7960eddd5ba2eb080a201f1d5fbb67a2bdee0 100644
--- a/dune/fem-dg/examples/dataio/checkpointing.hh
+++ b/dune/fem-dg/examples/dataio/checkpointing.hh
@@ -31,9 +31,6 @@ namespace Fem
       typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
       typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
 
-      typedef typename GridPartType :: template Codim< 0 > ::
-        template Partition< Dune::All_Partition > :: IteratorType IteratorType ;
-
       typedef typename DiscreteFunctionType::LocalFunctionType LocalFuncType;
 
       if( polOrd < 0 ) polOrd = 2*space.order() + 4 ;
@@ -43,20 +40,19 @@ namespace Fem
 
       double sum = 0.0;
 
-      IteratorType it    = space.gridPart().template begin< 0, Dune::All_Partition > ();
-      IteratorType endit = space.gridPart().template end< 0, Dune::All_Partition > ();
-
-      for(; it != endit ; ++it)
+      for( const auto& en : elements( space.gridPart(), Dune::Partitions::all ) )
       {
-        Dune::Fem::CachingQuadrature<GridPartType,0> quad(*it, polOrd);
-        LocalFuncType lf = discFunc.localFunction(*it);
-        for( size_t qP = 0; qP < quad.nop(); ++qP )
+        Dune::Fem::CachingQuadrature<GridPartType,0> quad( en, polOrd );
+        LocalFuncType lf = discFunc.localFunction( en );
+
+        for( const auto qp : quad )
         {
-          double det = (*it).geometry().integrationElement(quad.point(qP));
-          f.evaluate((*it).geometry().global(quad.point(qP)), ret);
-          lf.evaluate(quad[qP],phi);
+          const auto& x = qp.position();
+          double det = en.geometry().integrationElement( x );
+          f.evaluate( en.geometry().global( x ), ret);
+          lf.evaluate( qp, phi );
           RangeType diff = ret - phi ;
-          sum += det * quad.weight(qP) * ( diff * diff );
+          sum += det * qp.weight() * ( diff * diff );
         }
       }
       return std::sqrt( Dune::Fem::MPIManager::comm().sum( sum ) );
diff --git a/dune/fem-dg/examples/stokes/stokesalgorithm.hh b/dune/fem-dg/examples/stokes/stokesalgorithm.hh
index 9fac9b507c0800bacce65c001fa18b1865543c9a..2449022e65c1109d44075dc0133b592943abde83 100644
--- a/dune/fem-dg/examples/stokes/stokesalgorithm.hh
+++ b/dune/fem-dg/examples/stokes/stokesalgorithm.hh
@@ -315,12 +315,8 @@ namespace Fem
       polOrderContainer_.resize();
       if ( errorEstimator_.isPadaptive() )
       {
-        typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
-        typedef typename IteratorType::Entity EntityType ;
-        const IteratorType end = space_.end();
-        for( IteratorType it = space_.begin(); it != end; ++it )
+        for( const auto& entity : elements( space_.gridPart() ) )
         {
-          const EntityType& entity = *it;
           int order = polOrderContainer_[ entity ].value();
           while (order == -1) // is a new element
           {
diff --git a/dune/fem-dg/examples/stokes/stokesassembler.hh b/dune/fem-dg/examples/stokes/stokesassembler.hh
index 1998cd01a0ea238835f810920971d9021077ed4e..03cc7f9a312e8830ee76dcbbf845e7c118ed53ac 100644
--- a/dune/fem-dg/examples/stokes/stokesassembler.hh
+++ b/dune/fem-dg/examples/stokes/stokesassembler.hh
@@ -273,23 +273,14 @@ namespace Fem
     typedef typename DiscreteFunctionSpaceType::DomainType                DomainType;
     typedef typename DiscreteFunctionSpaceType::RangeType                 RangeType;
     typedef typename DiscreteFunctionSpaceType::JacobianRangeType         JacobianRangeType;
-    typedef typename DiscreteFunctionSpaceType::IteratorType              IteratorType;
 
     typedef typename DiscretePressureSpaceType::RangeType                 PressureRangeType;
     typedef typename DiscretePressureSpaceType::JacobianRangeType         PressureJacobianRangeType;
 
 
-    typedef typename DiscreteFunctionType::LocalFunctionType              LocalFuncType;
-    typedef typename DiscretePressureFunctionType::LocalFunctionType      LocalPressureType;
-
     // Types extracted from the underlying grid
-    typedef typename GridType::Traits::LeafIntersectionIterator           IntersectionIterator;
-
-    typedef typename GridPartType::Traits::IndexSetType                   IndexSetType;
-    typedef typename IntersectionIterator::Intersection                   IntersectionType;
+    typedef typename GridPartType::IntersectionType                       IntersectionType;
     typedef typename GridType::template Codim<0>::Entity                  EntityType ;
-    typedef typename GridType::template Codim<0>::Geometry                GeometryType;
-    //! type of quadrature to be used
 
   public:
     typedef typename OperatorTraits::VolumeQuadratureType                 VolumeQuadratureType;
@@ -364,12 +355,8 @@ namespace Fem
       pressureStabMatrix_->clear();
       pressureRhs_->clear();
 
-      typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType;
-      IteratorType end = spc_.end();
-      for(IteratorType it = spc_.begin(); it != end; ++it)
-      {
-        assembleLocal( *it );
-      }
+      for( const auto& en : elements( spc_.gridPart() ) )
+        assembleLocal( en );
 
       //important: finish build process!
       pressureGradMatrix_->communicate();
@@ -421,7 +408,7 @@ namespace Fem
 #endif
       //pressureDivMatrix_->matrix().print(std::cout);
       //std::cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n\n";
-      //pressureGradMatrix_->.matrix().print(std::cout);
+      //pressureGradMatrix_->matrix().print(std::cout);
       //std::cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n\n";
       //pressureStabMatrix_->matrix().print(std::cout);
     }
@@ -433,14 +420,12 @@ namespace Fem
 
       GridPartType& gridPart = spc_.gridPart();
 
-      const BaseFunctionSetType& bsetEn = spc_.basisFunctionSet(en);
-      const PressureBaseFunctionSetType& pressurebsetEn = pressurespc_.basisFunctionSet(en);
+      const auto& bsetEn         = spc_.basisFunctionSet(en);
+      const auto& pressurebsetEn = pressurespc_.basisFunctionSet(en);
       const size_t numBaseFunctions=bsetEn.size();
       const size_t numPBaseFunctions=pressurebsetEn.size();
 
-      //- typedefs
-      typedef typename EntityType :: Geometry Geometry;
-      const Geometry& geo=en.geometry();
+      const auto& geo = en.geometry();
 
       VolumeQuadratureType volQuad(en, volumeQuadOrd_);
 
@@ -458,18 +443,18 @@ namespace Fem
 #endif
       JacobianInverseType inv;
 
-      const size_t quadNop = volQuad.nop();
-      for (size_t l = 0; l < quadNop ; ++l)
+      for( const auto qp : volQuad )
       {
-        const typename  VolumeQuadratureType::CoordinateType& x=volQuad.point(l);
+        const auto& x = qp.position();
         inv = geo.jacobianInverseTransposed(x);
-        bsetEn.evaluateAll(volQuad[l],u);
-        bsetEn.jacobianAll(volQuad[l],du);
 
-        pressurebsetEn.evaluateAll(volQuad[l],p);
-        pressurebsetEn.jacobianAll(volQuad[l],dp);
+        bsetEn.evaluateAll(qp,u);
+        bsetEn.jacobianAll(qp,du);
 
-        double quadWeight = volQuad.weight(l)*geo.integrationElement(x);
+        pressurebsetEn.evaluateAll(qp,p);
+        pressurebsetEn.jacobianAll(qp,dp);
+
+        double weight = qp.weight()*geo.integrationElement(x);
 
         for(unsigned int k = 0;k < numBaseFunctions ;++k)
         {
@@ -482,39 +467,31 @@ namespace Fem
           for(unsigned int n = 0; n < numPBaseFunctions ; ++n)
           {
             //eval (-p_n*div u_j)
-            double PGM=p[n]*divu*quadWeight;
+            double PGM=p[n]*divu*weight;
             PGM*=-1;
             enPGrad.add(k,n,PGM);
 
             //eval -(-u_j*grad p_n)
-            double PDM =(u[k]*dp[n][0])*quadWeight;
+            double PDM =(u[k]*dp[n][0])*weight;
             enPDiv.add(n,k,PDM);
           }
         }
       }
 
-      IntersectionIterator endnit = gridPart.iend(en);
       typedef typename FaceQuadratureType :: NonConformingQuadratureType NonConformingFaceQuadratureType;
 
-      for (IntersectionIterator nit = gridPart.ibegin(en); nit != endnit; ++nit)
+      for (const auto& intersection : intersections(gridPart, en) )
       {
-        const IntersectionType& edge=*nit;
 
-        if(edge.neighbor())
+        if(intersection.neighbor())
         {
-#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
-          const EntityType& nb= edge.outside();
-#else
-         //Access to the neighbor Element
-          EntityPointerType neighEp=edge.outside();
-          const EntityType& nb=  *neighEp;
-#endif
-
-          if(edge.conforming() )
+          const EntityType& nb = intersection.outside();
+          if(intersection.conforming() )
           {
-            FaceQuadratureType faceQuadInner(gridPart,edge, faceQuadOrd_, FaceQuadratureType::INSIDE);
-            FaceQuadratureType faceQuadOuter(gridPart,edge, faceQuadOrd_, FaceQuadratureType::OUTSIDE);
-            applyNeighbor(edge, en, nb, faceQuadInner, faceQuadOuter, enPGrad, enPDiv
+            typedef Dune::Fem::IntersectionQuadrature< FaceQuadratureType, true > IntersectionQuadratureType;
+
+            IntersectionQuadratureType faceQuad(gridPart,intersection, faceQuadOrd_);
+            applyNeighbor(intersection, en, nb, faceQuad, enPGrad, enPDiv
 #if PRESSURESTABILIZATION
                 , enPStab
 #endif
@@ -522,21 +499,21 @@ namespace Fem
           }
           else
           {
-            NonConformingFaceQuadratureType  faceQuadOuter(gridPart, edge, faceQuadOrd_, NonConformingFaceQuadratureType::OUTSIDE);
-            NonConformingFaceQuadratureType  faceQuadInner(gridPart, edge, faceQuadOrd_, NonConformingFaceQuadratureType::INSIDE);
-            applyNeighbor(edge, en, nb, faceQuadInner, faceQuadOuter, enPGrad, enPDiv
+            typedef Dune::Fem::IntersectionQuadrature< FaceQuadratureType, false > IntersectionQuadratureType;
+            IntersectionQuadratureType faceQuad(gridPart,intersection, faceQuadOrd_);
+            applyNeighbor(intersection, en, nb, faceQuad, enPGrad, enPDiv
 #if PRESSURESTABILIZATION
                 , enPStab
 #endif
                 );
           }
         }
-        else if(edge.boundary())
+        else if(intersection.boundary())
         {
-          if(edge.conforming() )
+          if(intersection.conforming() )
           {
-            FaceQuadratureType faceQuadInner(gridPart,edge, faceQuadOrd_, FaceQuadratureType::INSIDE);
-            applyBoundary(edge, en, faceQuadInner, enPGrad, enPDiv
+            FaceQuadratureType faceQuadInner(gridPart,intersection, faceQuadOrd_, FaceQuadratureType::INSIDE);
+            applyBoundary(intersection, en, faceQuadInner, enPGrad, enPDiv
 #if PRESSURESTABILIZATION
                 , enPStab
 #endif
@@ -544,8 +521,8 @@ namespace Fem
           }
           else
           {
-            NonConformingFaceQuadratureType   faceQuadInner(gridPart, edge, faceQuadOrd_, NonConformingFaceQuadratureType::INSIDE);
-            applyBoundary(edge, en, faceQuadInner, enPGrad, enPDiv
+            NonConformingFaceQuadratureType   faceQuadInner(gridPart, intersection, faceQuadOrd_, NonConformingFaceQuadratureType::INSIDE);
+            applyBoundary(intersection, en, faceQuadInner, enPGrad, enPDiv
 #if PRESSURESTABILIZATION
                 , enPStab
 #endif
@@ -555,13 +532,13 @@ namespace Fem
         }
       }
     }
+
   private:
     template<class Quad,class Entity>
     void applyNeighbor(const IntersectionType &edge,
                        const Entity &en,
                        const Entity &nb,
-                       const Quad &faceQuadInner,
-                       const Quad &faceQuadOuter,
+                       const Quad &faceQuad,
                        LocalPressureGradMatType&  enPGrad,
                        LocalPressureDivMatType&   enPDiv
 #if PRESSURESTABILIZATION
@@ -569,10 +546,13 @@ namespace Fem
 #endif
                        ) const
     {
-      const BaseFunctionSetType& bsetEn = spc_.basisFunctionSet(en);
-      const BaseFunctionSetType& bsetNb = spc_.basisFunctionSet(nb);
-      const PressureBaseFunctionSetType& pressurebsetEn=pressurespc_.basisFunctionSet(en);
-      const PressureBaseFunctionSetType& pressurebsetNb=pressurespc_.basisFunctionSet(nb);
+      const auto& faceQuadInner  = faceQuad.inside();
+      const auto& faceQuadOuter  = faceQuad.outside();
+
+      const auto& bsetEn         = spc_.basisFunctionSet(en);
+      const auto& bsetNb         = spc_.basisFunctionSet(nb);
+      const auto& pressurebsetEn = pressurespc_.basisFunctionSet(en);
+      const auto& pressurebsetNb = pressurespc_.basisFunctionSet(nb);
       const size_t numBaseFunctionsEn=bsetEn.size();
       const size_t numBaseFunctionsNb=bsetNb.size();
       const size_t numPBaseFunctionsEn=pressurebsetEn.size();
@@ -593,24 +573,28 @@ namespace Fem
       PressureRangeType uNormal(0.);
       RangeType pNormal(0.);
 
-      const int quadNop = faceQuadInner.nop();
-
       LocalPressureGradMatType  nbPressGrad=pressureGradMatrix_->localMatrix(nb,en);
       LocalPressureDivMatType   nbPressDiv=pressureDivMatrix_->localMatrix(nb,en);
 
 #if PRESSURESTABILIZATION
       LocalPressureStabMatType  nbPressStab=pressureStabMatrix_->localMatrix(nb,en);
 #endif
-      for (int l = 0; l < quadNop ; ++l)
+      //for( const auto qpp : faceQuad )
+      //{
+      //  const auto& qp = qpp.first;
+      //  const auto& qpOut = qpp.second;
+      for( const auto qp : faceQuadInner )
       {
-        bsetEn.evaluateAll(faceQuadInner[l],uEn);
-        bsetNb.evaluateAll(faceQuadOuter[l],uNb);
-        pressurebsetEn.evaluateAll(faceQuadInner[l],pEn);
-        pressurebsetNb.evaluateAll(faceQuadOuter[l],pNb);
+        bsetEn.evaluateAll(qp,uEn);
+        bsetNb.evaluateAll(faceQuadOuter[qp.index()],uNb);
+        //bsetNb.evaluateAll(qpOut,uNb);
+        pressurebsetEn.evaluateAll(qp,pEn);
+        pressurebsetNb.evaluateAll(faceQuadOuter[qp.index()],pNb);
+        //pressurebsetNb.evaluateAll(qpOut,pNb);
 
-        DomainType normal=edge.integrationOuterNormal(faceQuadInner.localPoint(l));
+        DomainType normal=edge.integrationOuterNormal(qp.localPosition());
 
-        double intWeight=faceQuadInner.weight(l);
+        double weight=qp.weight();
 
         for(unsigned int j=0; j< numBaseFunctionsEn; ++j)
         {
@@ -621,7 +605,7 @@ namespace Fem
             // ******************************************************************************
             pNormal=normal;
             pNormal*=pEn[n][0];
-            double PGM_en=(uEn[j]*pNormal*intWeight);
+            double PGM_en=(uEn[j]*pNormal*weight);
             PGM_en*=0.5;
             enPGrad.add(j,n,PGM_en);
 
@@ -629,7 +613,7 @@ namespace Fem
             // -p+*u+.n+
             // ******************************************************************************
             uNormal[0]=uEn[j]*normal;
-            double PDM_en=(pEn[n]*uNormal)*intWeight;
+            double PDM_en=(pEn[n]*uNormal)*weight;
             PDM_en*=-0.5;
             enPDiv.add(n,j,PDM_en);
 
@@ -641,12 +625,12 @@ namespace Fem
               // ******************************************************************************
               for(unsigned int m=0;m<numPBaseFunctionsEn;++m)
               {
-                double PSM_nb=pEn[m]*pEn[n]*intWeight*d11_;
+                double PSM_nb=pEn[m]*pEn[n]*weight*d11_;
                 nbPressStab.add(m,n,PSM_nb);
               }
               for(unsigned int m=0;m<numPBaseFunctionsNb;++m)
               {
-                double PSM_en=pEn[n]*pNb[m]*intWeight*d11_;
+                double PSM_en=pEn[n]*pNb[m]*weight*d11_;
                 enPStab.add(m,n,PSM_en);
               }
 
@@ -664,7 +648,7 @@ namespace Fem
             // ******************************************************************************
             pNormal=normal;
             pNormal*=pNb[n][0];
-            double PGM_nb=(uEn[j]*pNormal*intWeight);
+            double PGM_nb=(uEn[j]*pNormal*weight);
             PGM_nb*=0.5;
             nbPressGrad.add(j,n,PGM_nb);
           }
@@ -679,7 +663,7 @@ namespace Fem
             // ******************************************************************************
 
             uNormal[0]=uNb[j]*normal;
-            double PDM_nb =( pEn[n]*uNormal[0])*intWeight;
+            double PDM_nb =( pEn[n]*uNormal[0])*weight;
             PDM_nb*=-0.5;
             nbPressDiv.add(n,j,PDM_nb);
 
@@ -691,12 +675,12 @@ namespace Fem
               // ******************************************************************************
               for(unsigned int m=0;m<numPBaseFunctionsEn;++m)
               {
-                double PSM_nb=pEn[m]*pEn[n]*intWeight*d11_;
+                double PSM_nb=pEn[m]*pEn[n]*weight*d11_;
                 nbPressStab.add(m,n,PSM_nb);
               }
               for(unsigned int m=0;m<numPBaseFunctionsNb;++m)
               {
-                double PSM_en=pEn[n]*pNb[m]*intWeight*d11_;
+                double PSM_en=pEn[n]*pNb[m]*weight*d11_;
                 enPStab.add(m,n,PSM_en);
               }
             }
@@ -719,8 +703,8 @@ namespace Fem
 #endif
                        ) const
     {
-      const BaseFunctionSetType& bsetEn = spc_.basisFunctionSet(en);
-      const PressureBaseFunctionSetType& pressurebsetEn=pressurespc_.basisFunctionSet(en);
+      const auto& bsetEn             = spc_.basisFunctionSet(en);
+      const auto& pressurebsetEn     = pressurespc_.basisFunctionSet(en);
       const size_t numBaseFunctions=bsetEn.size();
       const size_t numPBaseFunctions=pressurebsetEn.size();
       std::vector<RangeType> uEn(numBaseFunctions);
@@ -730,21 +714,19 @@ namespace Fem
       std::vector<PressureJacobianRangeType> dpEn(numPBaseFunctions);
 
 
-      LocalPressureType localPressure=pressureRhs_->localFunction(en);
+      auto localPressure = pressureRhs_->localFunction(en);
       RangeType dirichletValue(0.0);
       RangeType pNormal(0.0);
 
-      int quadNop=faceQuadInner.nop();
-
-      for (int l = 0; l < quadNop ; ++l)
+      for( const auto qp : faceQuadInner )
       {
-        bsetEn.evaluateAll(faceQuadInner[l],uEn);
-        pressurebsetEn.evaluateAll(faceQuadInner[l],pEn);
+        bsetEn.evaluateAll(qp,uEn);
+        pressurebsetEn.evaluateAll(qp,pEn);
 
-        DomainType normal=edge.integrationOuterNormal(faceQuadInner.localPoint(l));
+        auto normal = edge.integrationOuterNormal(qp.localPosition());
 
-        double intWeight=faceQuadInner.weight(l);
-        DomainType quadInEn=edge.geometry().global(faceQuadInner.localPoint(l));
+        double weight=qp.weight();
+        DomainType quadInEn=edge.geometry().global(qp.localPosition());
 
         problem_.template get<0>().g(quadInEn,dirichletValue);
         double pressureDirichlet;
@@ -756,7 +738,7 @@ namespace Fem
           // ******************************************************************************
           pressureDirichlet=normal*dirichletValue;
           pressureDirichlet*=pEn[n][0];
-          pressureDirichlet*=intWeight;
+          pressureDirichlet*=weight;
           //pressureDirichlet*=-1.;
           localPressure[n]+=pressureDirichlet;
         }
@@ -770,7 +752,7 @@ namespace Fem
             // ******************************************************************************
             pNormal=normal;
             pNormal*=pEn[n][0];
-            double PGM_en=(pNormal*uEn[j]) *intWeight;
+            double PGM_en=(pNormal*uEn[j]) *weight;
             enPGrad.add(j,n,PGM_en);
           }
         }
diff --git a/dune/fem-dg/misc/dgnorm.hh b/dune/fem-dg/misc/dgnorm.hh
index a6dca621b2a41ae6020ddab3efe5256a734a4f8b..d462e27e73b0e8370892b22bb985903d40974dc9 100644
--- a/dune/fem-dg/misc/dgnorm.hh
+++ b/dune/fem-dg/misc/dgnorm.hh
@@ -25,17 +25,17 @@ namespace Fem
 
   protected:
 
-    typedef typename BaseType::EntityType EntityType;
-    typedef typename EntityType::Geometry Geometry;
+    typedef typename BaseType::EntityType                   EntityType;
+    typedef typename EntityType::Geometry                   Geometry;
 
-    typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
-    typedef typename IntersectionIteratorType :: Intersection IntersectionType ;
-    typedef typename IntersectionType::Geometry IntersectionGeometryType;
+    typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
+    typedef typename GridPartType::IntersectionType         IntersectionType;
+    typedef typename IntersectionType::Geometry             IntersectionGeometryType;
 
-    typedef CachingQuadrature< GridPartType, 0 > QuadratureType;
-    typedef ElementQuadrature< GridPartType, 1 > FaceQuadratureType;
+    typedef CachingQuadrature< GridPartType, 0 >            QuadratureType;
+    typedef ElementQuadrature< GridPartType, 1 >            FaceQuadratureType;
   public:
-    typedef Integrator< QuadratureType > IntegratorType;
+    typedef Integrator< QuadratureType >                    IntegratorType;
 
 
     using BaseType::gridPart;
@@ -225,10 +225,8 @@ namespace Fem
 
     double jumpTerm = 0;
     {
-      const IntersectionIteratorType endiit = gridPart().iend( entity );
-      for ( IntersectionIteratorType iit = gridPart().ibegin( entity ); iit != endiit ; ++ iit )
+      for (const auto& intersection : intersections(gridPart(), entity) )
       {
-        const IntersectionType& intersection = *iit ;
         if( intersection.neighbor() )
         {
           const EntityType& neighbor = intersection.outside();
@@ -237,8 +235,7 @@ namespace Fem
           unsigned int nbIdx = gridPart().indexSet().index(neighbor);
           if( (enIdx < nbIdx) || (neighbor.partitionType() != Dune::InteriorEntity) )
           {
-            typedef typename IntersectionType :: Geometry IntersectionGeometry;
-            const IntersectionGeometry intersectionGeometry = intersection.geometry();
+            const auto intersectionGeometry = intersection.geometry();
 
             const double intersectionArea = intersectionGeometry.volume();
             const double heInverse = intersectionArea / std::min( geometry.volume(), geometryNb.volume() );
@@ -289,10 +286,8 @@ namespace Fem
 
     double jumpTerm = 0;
     {
-      const IntersectionIteratorType endiit = gridPart().iend( entity );
-      for ( IntersectionIteratorType iit = gridPart().ibegin( entity ); iit != endiit ; ++ iit )
+      for (const auto& intersection : intersections(gridPart(), entity) )
       {
-        const IntersectionType& intersection = *iit ;
         if( intersection.neighbor() )
         {
           const EntityType& neighbor = intersection.outside();
diff --git a/dune/fem-dg/operator/adaptation/adaptation.cc b/dune/fem-dg/operator/adaptation/adaptation.cc
index b0c55568dd494e5e5fbb026893805e727b47f935..62d591cf8039bec91cf5416a5c2a62f6b434a205 100644
--- a/dune/fem-dg/operator/adaptation/adaptation.cc
+++ b/dune/fem-dg/operator/adaptation/adaptation.cc
@@ -81,15 +81,9 @@ namespace Fem
   {
     double volume = 0;
     // type of iterator, i.e. leaf iterator
-    typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
 
-    const IteratorType endit = gridPart_.template end< 0 >();
-    for( IteratorType it = gridPart_.template begin< 0 >();
-         it != endit; ++it )
+    for( const auto& entity : elements( gridPart_ ) )
     {
-      // entity
-      const GridEntityType &entity = *it;
-
       // sum up the volume
       volume += entity.geometry().volume();
     }
@@ -337,16 +331,8 @@ namespace Fem
     // get local coarsen tolerance
     const double coarsenTol = refineTol * coarsenTheta_;
 
-    // type of iterator, i.e. leaf iterator
-    typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
-
-    const IteratorType endit = gridPart_.template end< 0 >();
-    for( IteratorType it = gridPart_.template begin< 0 >();
-         it != endit; ++it )
+    for( const auto& entity : elements( gridPart_ ) )
     {
-      // entity
-      const GridEntityType &entity = *it;
-
       // get local error indicator
       const double localIndicator = getLocalIndicator( entity );
       // get entity level
@@ -420,12 +406,8 @@ namespace Fem
 
     // count elements
     size_t count = 0;
-    // type of iterator, i.e. leaf iterator
-    typedef typename GridPartType::template Codim< 0 >::IteratorType IteratorType;
 
-    const IteratorType endit = gridPart_.template end< 0 >();
-    for( IteratorType it = gridPart_.template begin< 0 >();
-         it != endit; ++it )
+    for( const auto& entity : elements( gridPart_ ) )
     {
       ++count;
     }
diff --git a/dune/fem-dg/operator/adaptation/estimator.hh b/dune/fem-dg/operator/adaptation/estimator.hh
index b4023362bdf860dbaa4f0a839ca5fca55d160095..dc956cb6ed80389497cdca93365ec5cd664ddbd4 100644
--- a/dune/fem-dg/operator/adaptation/estimator.hh
+++ b/dune/fem-dg/operator/adaptation/estimator.hh
@@ -135,11 +135,8 @@ namespace Fem
         return;
 
       // also mark all neighbors of the actual entity for refinement
-      const IntersectionIteratorType nbend = gridPart_.iend( entity );
-      for (IntersectionIteratorType nb = gridPart_.ibegin( entity );
-           nb != nbend; ++nb)
+      for (const auto& intersection : intersections(gridPart_, entity) )
       {
-        const IntersectionType& intersection = *nb ;
         if( intersection.neighbor() )
         {
 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
@@ -231,16 +228,15 @@ namespace Fem
       LocalFunctionType lf = uh.localFunction( entity );
       const int quadOrder = ( lf.order()==0 ? 1 : lf.order() );
       ElementQuadratureType quad( entity, quadOrder );
-      const int numQuad = quad.nop();
 
       // get max and min of the indicator quantity
       double ind1LocMax = -1E100;
       double ind1LocMin =  1E100;
       double ind2LocMax = -1E100;
       double ind2LocMin =  1E100;
-      for( int qp=0; qp<numQuad; ++qp )
+      for( const auto qp : quad )
       {
-        DomainType xEn = quad.point( qp );
+        DomainType xEn = qp.position();
         lf.evaluate( xEn, val );
         const double ind1 = indicator1( entity, xEn, val );
         ind1LocMax = std::max( ind1LocMax, ind1 );
@@ -263,11 +259,8 @@ namespace Fem
       }
 
       // iterate over neighbors
-      const IntersectionIteratorType nbend = gridPart_.iend( entity );
-      for (IntersectionIteratorType nb = gridPart_.ibegin( entity );
-           nb != nbend; ++nb)
+      for (const auto& intersection : intersections(gridPart_, entity) )
       {
-        const IntersectionType& intersection = *nb ;
         if( intersection.neighbor() )
         {
           // access neighbor
@@ -282,16 +275,16 @@ namespace Fem
             // get local function on the neighbor element
             LocalFunctionType lfnb = uh.localFunction( neighbor );
             ElementQuadratureType quadNeigh( entity, quadOrder );
-            const int numQuadNe = quad.nop();
 
             // get max and min of the indicator quantity in the neighbor
             double ind1nbLocMax = -1E100;
             double ind1nbLocMin =  1E100;
             double ind2nbLocMax = -1E100;
             double ind2nbLocMin =  1E100;
-            for( int qpNe=0; qpNe<numQuadNe; ++qpNe )
+
+            for( const auto qp : quadNeigh )
             {
-              DomainType xNe = quadNeigh.point( qpNe );
+              DomainType xNe = qp.position();
               lfnb.evaluate( xNe, valnb );
               const double ind1nb = indicator1( neighbor, xNe, valnb );
               ind1nbLocMax = std::max( ind1nbLocMax, ind1nb );
@@ -336,11 +329,10 @@ namespace Fem
       ind2MaxDiff_   = 0;
 
       numberOfElements_ = 0 ;
-      const IteratorType end = dfSpace_.end();
-      for( IteratorType it = dfSpace_.begin(); it != end; ++it )
+      for( const auto& en : elements( dfSpace_.gridPart() ) )
       {
         // do local estimation
-        estimateLocal( uh, *it, indMin[ 0 ], indMax[ 0 ], indMin[ 1 ], indMax[ 1 ] );
+        estimateLocal( uh, en, indMin[ 0 ], indMax[ 0 ], indMin[ 1 ], indMax[ 1 ] );
         // count number of elements
         ++ numberOfElements_ ;
       }
diff --git a/dune/fem-dg/operator/adaptation/estimatorbase.hh b/dune/fem-dg/operator/adaptation/estimatorbase.hh
index 95ae61a5a8c4f4f1b7aad20d4f446f5505a22051..6a849ce371baaddeee04827e9a9253ba3468cdc0 100644
--- a/dune/fem-dg/operator/adaptation/estimatorbase.hh
+++ b/dune/fem-dg/operator/adaptation/estimatorbase.hh
@@ -50,7 +50,7 @@ namespace Fem
     typedef typename GridPartType :: IndexSetType                     IndexSetType;
     typedef typename GridPartType :: IntersectionIteratorType         IntersectionIteratorType;
 
-    typedef typename IntersectionIteratorType :: Intersection         IntersectionType;
+    typedef typename GridPartType :: IntersectionType                 IntersectionType;
 
     typedef typename GridPartType :: template Codim< 0 > :: EntityType        ElementType;
     typedef typename GridType :: template Codim< 0 > :: Entity                GridElementType;
@@ -102,9 +102,8 @@ namespace Fem
     virtual void estimate()
     {
       clear( indicator_ );
-      const IteratorType end = dfSpace_.end();
-      for( IteratorType it = dfSpace_.begin(); it != end; ++it )
-        estimateLocal( *it );
+      for( const auto& en : elements( dfSpace_.gridPart() ) )
+        estimateLocal( en );
     }
 
 
@@ -122,12 +121,8 @@ namespace Fem
     void mark()
     {
       // loop over all elements
-      const IteratorType end = dfSpace_.end();
-      for( IteratorType it = dfSpace_.begin(); it != end; ++it )
-      {
-        const ElementType &entity = *it;
-        markLocal( entity );
-      }
+      for( const auto& en : elements( dfSpace_.gridPart() ) )
+        markLocal( en );
     }
 
   };
diff --git a/dune/fem-dg/operator/adaptation/poissonestimator.hh b/dune/fem-dg/operator/adaptation/poissonestimator.hh
index cf66177f91e197064a27caea3170fa786a5488af..62bbb81550f244ed4d000ffe8bbafd0595dba8d4 100644
--- a/dune/fem-dg/operator/adaptation/poissonestimator.hh
+++ b/dune/fem-dg/operator/adaptation/poissonestimator.hh
@@ -191,7 +191,7 @@ namespace Fem
     typedef typename GridPartType::IndexSetType                      IndexSetType;
     typedef typename GridPartType::IntersectionIteratorType          IntersectionIteratorType;
 
-    typedef typename IntersectionIteratorType::Intersection          IntersectionType;
+    typedef typename GridPartType::IntersectionType                  IntersectionType;
 
     typedef typename GridType::template Codim< 0 >::Entity           ElementType;
     typedef typename ElementType::Geometry                           GeometryType;
@@ -338,20 +338,16 @@ namespace Fem
       }
       clear();
       henMin = 1e10;
-      const IteratorType end = dfSpace_.end();
-      for( IteratorType it = dfSpace_.begin(); it != end; ++it )
-      {
-        const ElementType &entity = *it;
 
+      for( const auto& entity : elements( dfSpace_.gridPart() ) )
+      {
         const LocalFunctionType uLocal = uh_.localFunction( entity );
         const SigmaLocalFunctionType sigmaLocal = sigma_.localFunction( entity );
 
         estimateLocal( rhs, entity, uLocal, sigmaLocal );
 
-        IntersectionIteratorType end = gridPart_.iend( entity );
-        for( IntersectionIteratorType inter = gridPart_.ibegin( entity ); inter != end; ++inter )
+        for (const auto& intersection : intersections(gridPart_, entity) )
         {
-          const IntersectionType &intersection = *inter;
           if( intersection.neighbor() )
             estimateIntersection( intersection, entity, uLocal, sigmaLocal );
           else
@@ -464,10 +460,8 @@ namespace Fem
 
       if (tolerance < 0)
       {
-        const IteratorType end = dfSpace_.end();
-        for( IteratorType it = dfSpace_.begin(); it != end; ++it )
+        for( const auto& entity : elements( dfSpace_.gridPart() ) )
         {
-          const ElementType &entity = *it;
           grid_.mark( 1, entity );
           ++marked;
         }
@@ -480,10 +474,8 @@ namespace Fem
                 tolerance*tolerance / (double)indexSet_.size( 0 );
 
         // loop over all elements
-        const IteratorType end = dfSpace_.end();
-        for( IteratorType it = dfSpace_.begin(); it != end; ++it )
+        for( const auto& entity : elements( dfSpace_.gridPart() ) )
         {
-          const ElementType &entity = *it;
           // check local error indicator
           if( (maximumStrategy_ && indicator_[ indexSet_.index( entity ) ] > localTol2)
               ||
@@ -509,18 +501,15 @@ namespace Fem
     }
     void closure()
     {
-      const IteratorType end = dfSpace_.end();
-      for( IteratorType it = dfSpace_.begin(); it != end; ++it )
+      for( const auto& entity : elements( dfSpace_.gridPart() ) )
       {
-        const ElementType &entity = *it;
         int conf = 0, nconf = 0;
         int lev = entity.level();
-        IntersectionIteratorType end = gridPart_.iend( entity );
-        for( IntersectionIteratorType inter = gridPart_.ibegin( entity ); inter != end; ++inter )
+        for (const auto& intersection : intersections(gridPart_, entity) )
         {
-          if (inter->neighbor())
+          if (intersection.neighbor())
           {
-            if (inter->outside().level() > lev)
+            if (intersection.outside().level() > lev)
               ++nconf;
             else
               ++conf;
@@ -577,23 +566,24 @@ namespace Fem
                         ( uLocal.order()*uLocal.order() );
       const int index = indexSet_.index( entity );
 
-      const int quadOrder = 2*(dfSpace_.order( entity ) )+2.;
+      const int quadOrder = 2*(dfSpace_.order( entity ) )+2;
       VolumeQuadratureType quad( entity, quadOrder );
-      const int numQuadraturePoints = quad.nop();
-      for( int qp = 0; qp < numQuadraturePoints; ++qp )
+
+      for( const auto qp : quad )
       {
+        const auto x = qp.position();
         // R_2 = h^2 |f+laplace u|^2
 
         RangeType y(0.),tmp(0.);
 
-        const DomainType global = geometry.global(quad.point(qp));
+        const DomainType global = geometry.global( x );
         rhs.f(global,y);
 
-        divergence( entity, geometry, quad, qp, h2, uLocal, sigmaLocal, tmp);
+        divergence( entity, geometry, quad, qp.index(), h2, uLocal, sigmaLocal, tmp);
 
         y+=tmp;
 
-        const double weight = quad.weight( qp ) * geometry.integrationElement( quad.point( qp ) );
+        const double weight = qp.weight() * geometry.integrationElement( x );
         indicator_[ index ] +=h2*weight * (y * y);
 
         R2_[ index ] += h2*weight * (y * y);
@@ -623,16 +613,17 @@ namespace Fem
       LocalMassMatrixType localMassMatrix( dgSpace, 2*dgSpace.order(entity) ) ;
       DGLFType udg( dgSpace );
       VolumeQuadratureType quad( entity, 2*(dgSpace.order(entity) )+1 );
-      const int numQuadraturePoints = quad.nop();
 
       udg.init( entity );
       udg.clear();
-      for( int qp = 0; qp < numQuadraturePoints; ++qp )
+
+      for( const auto qp : quad )
       {
+        const auto x = qp.position();
         RangeType uh;
-        uLocal.evaluate( quad[ qp ], uh );
-        uh *= quad.weight( qp ) * geometry.integrationElement( quad.point(qp) );
-        udg.axpy( quad[ qp ], uh );
+        uLocal.evaluate( qp, uh );
+        uh *= qp.weight() * geometry.integrationElement( x );
+        udg.axpy( qp, uh );
       }
 
       localMassMatrix.applyInverse( udg );
@@ -747,36 +738,40 @@ namespace Fem
       typedef SigmaConverterType JacobianTuple;
       typedef IntersectionQuadraturePointContext< IntersectionType, ElementType, FaceQuadratureType, RangeTuple, JacobianTuple > IntersectionLocalEvaluationType;
 
-      for( int qp = 0; qp < numQuadraturePoints; ++qp )
+      for( const auto qp : quadInside )
       {
-        DomainType unitNormal
-               = intersection.integrationOuterNormal( quadInside.localPoint( qp ) );
+        const auto idx = qp.index();
+        const auto& x = qp.localPosition();
+
+        DomainType unitNormal = intersection.integrationOuterNormal( x );
         const double integrationElement = unitNormal.two_norm();
 
         unitNormal/=integrationElement;
 
         // R_1 = h| (d u_l * n_l) + (d u_r * n_r) |^2 = h| (d u_l - d u_r) * n_l |^2
         JacobianRangeType AJacEn;
-        fluxEn[qp] /= integrationElement;
+        fluxEn[idx] /= integrationElement;
 
-        const SigmaConverterType sigmaValueEn( sigmaValuesEn[qp] );
-        IntersectionLocalEvaluationType local( intersection, inside, quadInside, uValuesEn[qp], sigmaValueEn, qp, 0, volume );
+        const SigmaConverterType sigmaValueEn( sigmaValuesEn[idx] );
+        IntersectionLocalEvaluationType local( intersection, inside, quadInside, uValuesEn[idx], sigmaValueEn, idx, 0, volume );
 
         oper_.model().diffusion(local,
-                                uValuesEn[qp], sigmaValueEn, AJacEn);
+                                uValuesEn[idx], sigmaValueEn, AJacEn);
         // note that flux=-hatA therefore we compute -hatA+Agrad u
-        AJacEn.umv( unitNormal, fluxEn[qp]);
+        AJacEn.umv( unitNormal, fluxEn[idx]);
 
         // R_orth = h^{-1} |u_l-u_r|
         RangeType jump;
-        jump=uValuesEn[qp];
-        jump-=uValuesNb[qp];
+        jump=uValuesEn[idx];
+        jump-=uValuesNb[idx];
 
-        errorInside  += quadInside.weight( qp ) *h* (fluxEn[qp] * fluxEn[qp]) *integrationElement;
-        errorInside  += quadInside.weight( qp ) *1./h* (jump * jump) *integrationElement;
+        const auto& weight = qp.weight();
 
-        R1_[ insideIndex ] += quadInside.weight( qp ) *h* (fluxEn[qp]*fluxEn[qp]) *integrationElement;
-        Rorth_[ insideIndex ] += quadInside.weight( qp ) *1./h* (jump * jump) *integrationElement;
+        errorInside  += weight *h* (fluxEn[idx] * fluxEn[idx]) *integrationElement;
+        errorInside  += weight *1./h* (jump * jump) *integrationElement;
+
+        R1_[ insideIndex ] += weight *h* (fluxEn[idx]*fluxEn[idx]) *integrationElement;
+        Rorth_[ insideIndex ] += weight *1./h* (jump * jump) *integrationElement;
       }
       if( errorInside > 0.0 )
         indicator_[ insideIndex ] +=  errorInside;
@@ -893,47 +888,52 @@ namespace Fem
       typedef SigmaConverterType  JacobianTuple;
       typedef IntersectionQuadraturePointContext< IntersectionType, ElementType, QuadratureImp, RangeTuple, JacobianTuple > IntersectionLocalEvaluationType;
 
-      for( int qp = 0; qp < numQuadraturePoints; ++qp )
+      for( const auto qp : quadInside )
       {
-        DomainType unitNormal
-               = intersection.integrationOuterNormal( quadInside.localPoint( qp ) );
+        const auto idx = qp.index();
+        const auto& x = qp.localPosition();
+
+        DomainType unitNormal = intersection.integrationOuterNormal( x );
         const double integrationElement = unitNormal.two_norm();
 
         unitNormal/=integrationElement;
 
         // R_1 = h| (d u_l * n_l) + (d u_r * n_r) |^2 = h| (d u_l - d u_r) * n_l |^2
         JacobianRangeType AJacEn,AJacNb;
-        fluxEn[qp] /= integrationElement;
-        fluxNb[qp] /= integrationElement;
+        fluxEn[idx] /= integrationElement;
+        fluxNb[idx] /= integrationElement;
 
-        const SigmaConverterType sigmaValueEn( sigmaValuesEn[qp] );
-        const IntersectionLocalEvaluationType left( intersection, inside, quadInside, uValuesEn[qp], sigmaValueEn, qp, 0, volume );
-        const SigmaConverterType sigmaValueNb( sigmaValuesNb[qp] );
-        const IntersectionLocalEvaluationType right( intersection, outside, quadInside, uValuesNb[qp], sigmaValueNb, qp, 0, volume );
+        const SigmaConverterType sigmaValueEn( sigmaValuesEn[idx] );
+        const IntersectionLocalEvaluationType left( intersection, inside, quadInside, uValuesEn[idx], sigmaValueEn, idx, 0, volume );
+        const SigmaConverterType sigmaValueNb( sigmaValuesNb[idx] );
+        const IntersectionLocalEvaluationType right( intersection, outside, quadInside, uValuesNb[idx], sigmaValueNb, idx, 0, volume );
 
         oper_.model().diffusion(left,
-                                uValuesEn[qp], sigmaValueEn, AJacEn);
+                                uValuesEn[idx], sigmaValueEn, AJacEn);
         oper_.model().diffusion(right,
-                                uValuesNb[qp], sigmaValueNb, AJacNb);
+                                uValuesNb[idx], sigmaValueNb, AJacNb);
 
         // note that flux=-hatA therefore we compute -hatA+Agrad u
-        AJacEn.umv( unitNormal, fluxEn[qp]);
-        AJacNb.umv( unitNormal, fluxNb[qp]);
+        AJacEn.umv( unitNormal, fluxEn[idx]);
+        AJacNb.umv( unitNormal, fluxNb[idx]);
 
         // R_orth = h^{-1} |u_l-u_r|
         RangeType jump;
-        jump=uValuesEn[qp];
-        jump-=uValuesNb[qp];
-
-        errorInside  += quadInside.weight( qp ) *h* (fluxEn[qp] * fluxEn[qp]) *integrationElement;
-        errorInside  += quadInside.weight( qp ) *1./h* (jump * jump) *integrationElement;
-        errorOutside += quadOutside.weight( qp ) *h* (fluxNb[qp] * fluxNb[qp]) *integrationElement;
-        errorOutside += quadOutside.weight( qp ) *1./h* (jump * jump) *integrationElement;
-
-        R1_[ insideIndex ] += quadInside.weight( qp ) *h* (fluxEn[qp]*fluxEn[qp]) *integrationElement;
-        Rorth_[ insideIndex ] += quadInside.weight( qp ) *1./h* (jump * jump) *integrationElement;
-        R1_[ outsideIndex ] += quadOutside.weight( qp ) *h* (fluxNb[qp]*fluxNb[qp]) *integrationElement;
-        Rorth_[ outsideIndex ] += quadOutside.weight( qp ) *1./h* (jump * jump) *integrationElement;
+        jump=uValuesEn[idx];
+        jump-=uValuesNb[idx];
+
+        const auto& weightIn = qp.weight();
+        const auto& weightOut = quadOutside.weight( idx );
+
+        errorInside  += weightIn *h* (fluxEn[idx] * fluxEn[idx]) *integrationElement;
+        errorInside  += weightIn *1./h* (jump * jump) *integrationElement;
+        errorOutside += weightOut *h* (fluxNb[idx] * fluxNb[idx]) *integrationElement;
+        errorOutside += weightOut *1./h* (jump * jump) *integrationElement;
+
+        R1_[ insideIndex ] += weightIn *h* (fluxEn[idx]*fluxEn[idx]) *integrationElement;
+        Rorth_[ insideIndex ] += weightIn *1./h* (jump * jump) *integrationElement;
+        R1_[ outsideIndex ] += weightOut *h* (fluxNb[idx]*fluxNb[idx]) *integrationElement;
+        Rorth_[ outsideIndex ] += weightOut *1./h* (jump * jump) *integrationElement;
       }
     }
 
diff --git a/dune/fem-dg/operator/adaptation/stokesestimator.hh b/dune/fem-dg/operator/adaptation/stokesestimator.hh
index b0fb350ffa236c0948f60cb3ae73379d8fcb2648..ca5fa4f9cdb25913fcd04d2b1b26b294a039a402 100644
--- a/dune/fem-dg/operator/adaptation/stokesestimator.hh
+++ b/dune/fem-dg/operator/adaptation/stokesestimator.hh
@@ -54,7 +54,7 @@ namespace Fem
     typedef typename GridPartType::IndexSetType                      IndexSetType;
     typedef typename GridPartType::IntersectionIteratorType          IntersectionIteratorType;
 
-    typedef typename IntersectionIteratorType::Intersection          IntersectionType;
+    typedef typename GridPartType::IntersectionType                  IntersectionType;
 
     typedef typename GridType::template Codim< 0 >::Entity           ElementType;
 #if not DUNE_VERSION_NEWER(DUNE_GRID,2,4)
@@ -103,20 +103,16 @@ namespace Fem
       Rdiv_.resize( this->indexSet_.size( 0 ));
       std::fill( Rdiv_.begin(), Rdiv_.end(), 0);
 
-      const IteratorType end = dfSpace_.end();
-      for( IteratorType it = dfSpace_.begin(); it != end; ++it )
+      for( const auto& entity : elements( dfSpace_.gridPart() ) )
       {
-        const ElementType &entity = *it;
         const LocalFunctionType uLocal = uh_.localFunction( entity );
         const SigmaLocalFunctionType sigmaLocal = sigma_.localFunction( entity );
 
         BaseType::estimateLocal( rhs, entity, uLocal, sigmaLocal );
         estimateLocal( rhs, entity, uLocal, sigmaLocal );
 
-        IntersectionIteratorType end = gridPart_.iend( entity );
-        for( IntersectionIteratorType inter = gridPart_.ibegin( entity ); inter != end; ++inter )
+        for (const auto& intersection : intersections(gridPart_, entity) )
         {
-          const IntersectionType &intersection = *inter;
           if( intersection.neighbor() )
             BaseType::estimateIntersection( intersection, entity, uLocal, sigmaLocal );
            else
@@ -179,20 +175,19 @@ namespace Fem
       ElementQuadratureType quad( entity, 2*(dfSpace_.order() + 2) );
       JacobianInverseType inv;
 
-      const int numQuadraturePoints = quad.nop();
-      for( int qp = 0; qp < numQuadraturePoints; ++qp )
+      for( const auto qp : quad )
       {
-        const typename  ElementQuadratureType::CoordinateType &x=quad.point(qp);
+        const auto& x = qp.position();
         inv = geometry.jacobianInverseTransposed(x);
         JacobianRangeType uJac(0.);
-        uLocal.jacobian(quad[qp],uJac);
+        uLocal.jacobian( qp,uJac );
 
         double divergence = 0;
         for (int i=0;i<uJac.rows;++i)
         {
           divergence += uJac[i][i];
         }
-        const double weight = quad.weight( qp ) * geometry.integrationElement( quad.point( qp ) );
+        const double weight = qp.weight() * geometry.integrationElement( x );
         indicator_[ index ] += h2*weight * (divergence * divergence);
 
         Rdiv_[ index ] += h2*weight * (divergence * divergence);
diff --git a/dune/fem-dg/operator/adaptation/utility.hh b/dune/fem-dg/operator/adaptation/utility.hh
index e1204d2e3f65c7ba880265bdefeb3579535162fa..a4f42353147274b3ad1f079b676bd38b1a3bc926 100644
--- a/dune/fem-dg/operator/adaptation/utility.hh
+++ b/dune/fem-dg/operator/adaptation/utility.hh
@@ -163,7 +163,6 @@ namespace Fem
                                         const int finestLevel )
     {
       typedef typename GridPart :: GridType GridType ;
-      typedef typename GridPart :: template Codim< 0 > :: IteratorType IteratorType;
 
       const bool hasGridHierarchy = Fem :: GridPartCapabilities :: hasGrid< GridPart > :: v;
 
@@ -176,10 +175,10 @@ namespace Fem
       double maxVolume[ 1 ] = { std::numeric_limits< double > ::min() };
 
       // if grid is not empty compute smallest and biggest volume
-      const IteratorType end = gridPart.template end< 0 >();
-      for( IteratorType it = gridPart.template begin< 0 >(); it != end; ++it )
+
+      for( const auto& en : elements( gridPart ) )
       {
-        const double volume = findCoarsestVolume( *it, hasGridHierarchy );
+        const double volume = findCoarsestVolume( en, hasGridHierarchy );
         minVolume[ 0 ] = std::min( minVolume[ 0 ], volume );
         maxVolume[ 0 ] = std::max( maxVolume[ 0 ], volume );
       }
diff --git a/dune/fem-dg/operator/dg/discretemodelcommon.hh b/dune/fem-dg/operator/dg/discretemodelcommon.hh
index 152e3e9994b9aa5214f6efcf011e09556d302d66..732ae5cbdf7ca91fbaabea4b66e53d04fb043965 100644
--- a/dune/fem-dg/operator/dg/discretemodelcommon.hh
+++ b/dune/fem-dg/operator/dg/discretemodelcommon.hh
@@ -143,7 +143,7 @@ namespace Fem
     typedef typename Traits :: DiscreteFunctionSpaceType               DiscreteFunctionSpaceType;
 
     typedef typename GridPartType :: IntersectionIteratorType          IntersectionIteratorType;
-    typedef typename IntersectionIteratorType :: Intersection          IntersectionType;
+    typedef typename GridPartType :: IntersectionType                  IntersectionType;
     typedef typename BaseType :: EntityType                            EntityType;
     typedef typename DiscreteFunctionSpaceType :: FunctionSpaceType    FunctionSpaceType;
     typedef typename FunctionSpaceType :: RangeFieldType               RangeFieldType;
diff --git a/dune/fem-dg/operator/dg/fluxdiscretemodel.hh b/dune/fem-dg/operator/dg/fluxdiscretemodel.hh
index b0aac05b90e3a01419cf89209750817c93d7997c..a7d65bb859f14e24efa5ea0150d3da01e5404be0 100644
--- a/dune/fem-dg/operator/dg/fluxdiscretemodel.hh
+++ b/dune/fem-dg/operator/dg/fluxdiscretemodel.hh
@@ -62,9 +62,9 @@ namespace Fem
     typedef typename Traits::DiscreteFunctionSpaceType::FunctionSpaceType::RangeType                          RangeType;
     typedef typename Traits::DiscreteFunctionSpaceType::FunctionSpaceType::JacobianRangeType                  JacobianRangeType;
     typedef typename Traits::GridType                              GridType;
-    typedef typename Traits::GridPartType
-              ::IntersectionIteratorType                           IntersectionIterator;
-    typedef typename IntersectionIterator::Intersection            Intersection;
+    typedef typename Traits::GridPartType                          GridPartType;
+    typedef typename GridPartType::IntersectionIteratorType        IntersectionIterator;
+    typedef typename GridPartType::IntersectionType                Intersection;
     typedef typename BaseType::EntityType                          EntityType;
 
     enum { evaluateJacobian = DiffusionFluxType::evaluateJacobian };
@@ -278,7 +278,7 @@ namespace Fem
     typedef typename Traits::GridPartType                     GridPartType;
     typedef typename Traits::GridType                         GridType;
     typedef typename GridPartType::IntersectionIteratorType   IntersectionIterator;
-    typedef typename IntersectionIterator::Intersection       Intersection;
+    typedef typename GridPartType::IntersectionType           Intersection;
     typedef typename BaseType::EntityType                     EntityType;
     typedef typename Traits::DiscreteFunctionSpaceType        DiscreteFunctionSpaceType;
 
diff --git a/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh b/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh
index 13eb7bd9224cefd3e4cfccdcf4f40b909504fab1..aeed9cea39c682ca37a0a9bea6574cc053ee83e0 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh
@@ -38,36 +38,35 @@ namespace Fem
   public:
     typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
 
-    enum { dimDomain = DiscreteFunctionSpaceType :: dimDomain };
-    enum { dimRange  = DiscreteFunctionSpaceType :: dimRange };
+    enum { dimDomain = DiscreteFunctionSpaceType::dimDomain };
+    enum { dimRange  = DiscreteFunctionSpaceType::dimRange };
     enum { dimGradRange = dimDomain * dimRange };
-    enum { polOrd = DiscreteFunctionSpaceType :: polynomialOrder };
+    enum { polOrd = DiscreteFunctionSpaceType::polynomialOrder };
 
-    typedef typename DiscreteFunctionSpaceType :: RangeFieldType        RangeFieldType;
-    typedef typename DiscreteFunctionSpaceType :: DomainFieldType       DomainFieldType;
-    typedef FieldVector< DomainFieldType, dimDomain-1 >                 FaceDomainType;
-    typedef typename DiscreteFunctionSpaceType :: DomainType            DomainType;
-    typedef typename DiscreteFunctionSpaceType :: RangeType             RangeType;
-    typedef typename DiscreteFunctionSpaceType :: JacobianRangeType     JacobianRangeType;
+    typedef typename DiscreteFunctionSpaceType::RangeFieldType        RangeFieldType;
+    typedef typename DiscreteFunctionSpaceType::DomainFieldType       DomainFieldType;
+    typedef FieldVector< DomainFieldType, dimDomain-1 >               FaceDomainType;
+    typedef typename DiscreteFunctionSpaceType::DomainType            DomainType;
+    typedef typename DiscreteFunctionSpaceType::RangeType             RangeType;
+    typedef typename DiscreteFunctionSpaceType::JacobianRangeType     JacobianRangeType;
 
 
-    typedef typename DiscreteFunctionSpaceType :: GridPartType          GridPartType;
-    typedef typename GridPartType :: IntersectionIteratorType           IntersectionIterator;
-    typedef typename IntersectionIterator :: Intersection               Intersection;
-    typedef typename GridPartType :: GridType                           GridType;
-    typedef typename DiscreteFunctionSpaceType :: EntityType            EntityType;
-    typedef typename GridPartType::template Codim< 0 >::IteratorType    IteratorType;
-    typedef typename GridPartType::IntersectionIteratorType
-                                                                        IntersectionIteratorType;
+    typedef typename DiscreteFunctionSpaceType::GridPartType          GridPartType;
+    typedef typename GridPartType::IntersectionIteratorType           IntersectionIterator;
+    typedef typename GridPartType::IntersectionType                   IntersectionType;
+    typedef typename GridPartType::GridType                           GridType;
+    typedef typename DiscreteFunctionSpaceType::EntityType            EntityType;
+    typedef typename GridPartType::template Codim< 0 >::IteratorType  IteratorType;
+    typedef typename GridPartType::IntersectionIteratorType           IntersectionIteratorType;
 
-    typedef typename BaseType :: DiscreteGradientSpaceType              DiscreteGradientSpaceType;
-    typedef typename DiscreteGradientSpaceType :: RangeType             GradientType;
-    typedef Fem::TemporaryLocalFunction< DiscreteGradientSpaceType >    LiftingFunctionType;
+    typedef typename BaseType::DiscreteGradientSpaceType              DiscreteGradientSpaceType;
+    typedef typename DiscreteGradientSpaceType::RangeType             GradientType;
+    typedef Fem::TemporaryLocalFunction< DiscreteGradientSpaceType >  LiftingFunctionType;
 
-    typedef Fem::CachingQuadrature< GridPartType, 0>                    VolumeQuadratureType ;
+    typedef Fem::CachingQuadrature< GridPartType, 0>                  VolumeQuadratureType ;
 
     typedef Fem::LocalMassMatrix
-      < DiscreteGradientSpaceType, VolumeQuadratureType >               LocalMassMatrixType;
+      < DiscreteGradientSpaceType, VolumeQuadratureType >             LocalMassMatrixType;
 
     class Lifting
     {
@@ -203,19 +202,14 @@ namespace Fem
       int maxNumOutflowFaces = 0;
       if ( useTheoryParams_ )
       {
-        const IteratorType itend = gridPart.template end<0>();
-        for( IteratorType it = gridPart.template begin<0>(); it != itend; ++it )
+        for( const auto& entity : elements( gridPart ) )
         {
-          const EntityType& entity = * it ;
           const double insideVol = entity.geometry().volume();
           int numFaces = 0;
           int numOutflowFaces = 0;
-          const IntersectionIteratorType intitend = gridPart.iend( entity );
-          for(IntersectionIteratorType intit = gridPart.ibegin( entity );
-              intit != intitend; ++intit )
-          {
-            const Intersection& intersection = * intit ;
 
+          for (const auto& intersection : intersections(gridPart, entity) )
+          {
             ++numFaces ;
             if ( intersection.neighbor() )
             {
@@ -378,7 +372,7 @@ namespace Fem
     }
 
     template <class QuadratureImp, class ArgumentTupleVector >
-    void initializeIntersection(const Intersection& intersection,
+    void initializeIntersection(const IntersectionType& intersection,
                                 const EntityType& inside,
                                 const EntityType& outside,
                                 const double time,
@@ -397,7 +391,7 @@ namespace Fem
     }
 
     template <class QuadratureImp, class ArgumentTupleVector >
-    void computeLiftings(const Intersection& intersection,
+    void computeLiftings(const IntersectionType& intersection,
                          const EntityType& inside,
                          const EntityType& outside,
                          const double time,
@@ -606,7 +600,7 @@ namespace Fem
 
 
     template <class QuadratureImp, class ArgumentTupleVector>
-    void initializeBoundary(const Intersection& intersection,
+    void initializeBoundary(const IntersectionType& intersection,
                             const EntityType& entity,
                             const double time,
                             const QuadratureImp& quadInner,
@@ -637,7 +631,7 @@ namespace Fem
 
   protected:
     template <class QuadratureImp, class ArgumentTuple, class LiftingFunction >
-    void addLifting(const Intersection& intersection,
+    void addLifting(const IntersectionType& intersection,
                     const EntityType &entity,
                     const ArgumentTuple &uTuple,
                     const RangeType& u,
@@ -648,7 +642,7 @@ namespace Fem
                     const RangeType& uRight,
                     LiftingFunction& func ) const
     {
-      IntersectionQuadraturePointContext< Intersection, EntityType, QuadratureImp, ArgumentTuple, ArgumentTuple >
+      IntersectionQuadraturePointContext< IntersectionType, EntityType, QuadratureImp, ArgumentTuple, ArgumentTuple >
         local( intersection, entity, faceQuad, uTuple, uTuple, quadPoint, time, entity.geometry().volume() );
 
       const FaceDomainType& x = faceQuad.localPoint( quadPoint );
@@ -708,11 +702,8 @@ namespace Fem
     {
       double sumFaceVolSqr  = 0.0;
 
-      const IntersectionIteratorType intitend = gridPart_.iend( entity );
-      for(IntersectionIteratorType intit = gridPart_.ibegin( entity );
-          intit != intitend; ++intit )
+      for (const auto& intersection : intersections(gridPart_, entity) )
       {
-        const Intersection& intersection = * intit ;
         const double faceVol = intersection.geometry().volume();
 
         // !!!!! forget about Neumann for now
@@ -1193,7 +1184,7 @@ namespace Fem
 
   public:
     typedef typename BaseType::GridPartType        GridPartType;
-    typedef typename BaseType::Intersection        Intersection;
+    typedef typename BaseType::IntersectionType    IntersectionType;
     typedef typename BaseType::EntityType          EntityType;
     typedef typename BaseType::RangeType           RangeType;
     typedef typename BaseType::JacobianRangeType   JacobianRangeType;
@@ -1229,7 +1220,7 @@ namespace Fem
     }
 
     template <class QuadratureImp, class ArgumentTupleVector >
-    void initializeIntersection(const Intersection& intersection,
+    void initializeIntersection(const IntersectionType& intersection,
                                 const EntityType& inside,
                                 const EntityType& outside,
                                 const double time,
diff --git a/dune/fem-dg/operator/fluxes/diffusion/fluxbase.hh b/dune/fem-dg/operator/fluxes/diffusion/fluxbase.hh
index c17f44b2c3c8cc9fe1d26d26b45fd0c94a084dc9..6f1cf14e1428ccd3586d145325d77b6e8601e979 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/fluxbase.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/fluxbase.hh
@@ -77,7 +77,7 @@ namespace Fem
 
     typedef typename DiscreteFunctionSpaceType::GridPartType         GridPartType;
     typedef typename GridPartType::IntersectionIteratorType          IntersectionIterator;
-    typedef typename IntersectionIterator::Intersection              Intersection;
+    typedef typename GridPartType::IntersectionType                  IntersectionType;
     typedef typename GridPartType::GridType                          GridType;
     typedef typename DiscreteFunctionSpaceType::EntityType           EntityType;
     enum { dimGradRange = dimDomain * dimRange };
@@ -140,7 +140,7 @@ namespace Fem
 
   protected:
     bool determineDirection( const bool areaSwitch, const double enVolume, const double nbVolume,
-                             const Intersection& intersection ) const
+                             const IntersectionType& intersection ) const
     {
       if (areaSwitch && std::abs( enVolume - nbVolume ) > 1e-8)
       {
@@ -200,7 +200,7 @@ namespace Fem
     }
 
     template <class QuadratureImp, class ArgumentTupleVector >
-    void initializeIntersection(const Intersection& intersection,
+    void initializeIntersection(const IntersectionType& intersection,
                                 const EntityType& inside,
                                 const EntityType& outside,
                                 const double time,
@@ -212,7 +212,7 @@ namespace Fem
     }
 
     template <class QuadratureImp, class ArgumentTupleVector>
-    void initializeBoundary(const Intersection& intersection,
+    void initializeBoundary(const IntersectionType& intersection,
                             const EntityType& entity,
                             const double time,
                             const QuadratureImp& quadInner,
@@ -313,7 +313,7 @@ namespace Fem
 
     typedef typename BaseType::GridPartType             GridPartType;
     typedef typename BaseType::IntersectionIterator     IntersectionIterator;
-    typedef typename BaseType::Intersection             Intersection;
+    typedef typename BaseType::IntersectionType         IntersectionType;
     typedef typename BaseType::GridType                 GridType;
     typedef typename BaseType::EntityType               EntityType;
     enum { dimGradRange = dimDomain * dimRange };
diff --git a/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh b/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
index 7ce14f62d5f9f276e8a944a3987e7b8cb20cbd22..6153916d7c62e918eaa3269b43c81d17da35750a 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
@@ -37,7 +37,7 @@ namespace Fem
 
     typedef typename DiscreteFunctionSpaceType :: GridPartType         GridPartType;
     typedef typename GridPartType :: IntersectionIteratorType          IntersectionIterator;
-    typedef typename IntersectionIterator :: Intersection              Intersection;
+    typedef typename GridPartType :: IntersectionType                  Intersection;
     typedef typename GridPartType :: GridType                          GridType;
     typedef typename DiscreteFunctionSpaceType :: EntityType           EntityType;
     enum { dimGradRange = dimDomain * dimRange };
diff --git a/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh b/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
index 98aa4f637ebd9a65bb98ad75cdced2c681e6807d..54b8304f442a9adf75db4b24e13cd04feeead3a0 100644
--- a/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
+++ b/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
@@ -37,28 +37,28 @@ namespace Fem
     // These type definitions allow a convenient access to arguments of pass.
     std::integral_constant< int, passId > uVar;
   public:
-    typedef LimiterTraits<GlobalPassTraitsImp,Model, passId > Traits;
-
-    typedef typename Traits::RangeType RangeType;
-    typedef typename Traits:: DomainType DomainType;
-    typedef typename Traits:: LocalDomainType LocalDomainType;
-    typedef typename Traits:: DomainFieldType DomainFieldType;
-    typedef typename Traits::GridType GridType;
-    typedef typename Traits::GridPartType GridPartType;
-    typedef typename Traits::JacobianRangeType JacobianRangeType;
-    typedef typename GridPartType::IntersectionIteratorType :: Intersection IntersectionType;
+    typedef LimiterTraits<GlobalPassTraitsImp,Model, passId >           Traits;
+
+    typedef typename Traits::RangeType                                  RangeType;
+    typedef typename Traits::DomainType                                 DomainType;
+    typedef typename Traits::LocalDomainType                            LocalDomainType;
+    typedef typename Traits::DomainFieldType                            DomainFieldType;
+    typedef typename Traits::GridType                                   GridType;
+    typedef typename Traits::GridPartType                               GridPartType;
+    typedef typename Traits::JacobianRangeType                          JacobianRangeType;
+    typedef typename GridPartType::IntersectionType                     IntersectionType;
     typedef typename GridPartType::template Codim<0>::EntityType        EntityType;
     typedef typename GridType::template Codim<0>::Entity                GridEntityType;
 
     enum { dimGrid = GridType :: dimension };
 
     // type of surface domain type
-    typedef FieldVector< DomainFieldType, dimGrid - 1 > FaceLocalDomainType;
+    typedef FieldVector< DomainFieldType, dimGrid - 1 >                 FaceLocalDomainType;
 
-    typedef typename Traits::DestinationType DestinationType;
+    typedef typename Traits::DestinationType                            DestinationType;
 
     // Indicator for Limiter
-    typedef typename Traits::IndicatorType   IndicatorType ;
+    typedef typename Traits::IndicatorType                              IndicatorType ;
 
     enum { dimRange = Traits :: dimRange };
     enum { evaluateJacobian = false };
diff --git a/dune/fem-dg/operator/limiter/limitpass.hh b/dune/fem-dg/operator/limiter/limitpass.hh
index be42ab4283bcea2176e6ab075b191370e26c1a48..c5231796317e5ab7041070185367f680d695a76e 100644
--- a/dune/fem-dg/operator/limiter/limitpass.hh
+++ b/dune/fem-dg/operator/limiter/limitpass.hh
@@ -298,21 +298,21 @@ namespace Fem
   public:
     typedef LimiterDefaultTraits<GlobalTraitsImp,Model,passId> Traits;
 
-    typedef typename Traits::DomainType DomainType;
-    typedef typename Traits::LocalDomainType LocalDomainType;
-    typedef typename Traits::FaceLocalDomainType FaceLocalDomainType;
-    typedef typename Traits::RangeType RangeType;
-    typedef typename Traits::GridType GridType;
-    typedef typename Traits::GridPartType GridPartType;
-    typedef typename Traits::JacobianRangeType JacobianRangeType;
+    typedef typename Traits::DomainType                     DomainType;
+    typedef typename Traits::LocalDomainType                LocalDomainType;
+    typedef typename Traits::FaceLocalDomainType            FaceLocalDomainType;
+    typedef typename Traits::RangeType                      RangeType;
+    typedef typename Traits::GridType                       GridType;
+    typedef typename Traits::GridPartType                   GridPartType;
+    typedef typename Traits::JacobianRangeType              JacobianRangeType;
     typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
-    typedef typename IntersectionIteratorType :: Intersection IntersectionType;
-    typedef typename GridPartType::template Codim<0>::EntityType        EntityType;
-    typedef typename DomainType :: field_type DomainFieldType;
+    typedef typename GridPartType::IntersectionType         IntersectionType;
+    typedef typename GridPartType::template Codim<0>::EntityType EntityType;
+    typedef typename DomainType::field_type                 DomainFieldType;
 
-    typedef typename Traits :: LimiterFunctionType  LimiterFunctionType;
+    typedef typename Traits::LimiterFunctionType            LimiterFunctionType;
 
-    enum { dimRange = RangeType :: dimension };
+    enum { dimRange = RangeType::dimension };
 
   public:
     /** \brief default limiter discrete model */
@@ -608,92 +608,90 @@ namespace Fem
     //- Typedefs and enums
 
     //! Repetition of template arguments
-    typedef DiscreteModelImp DiscreteModelType;
+    typedef DiscreteModelImp                                             DiscreteModelType;
     //! Repetition of template arguments
-    typedef PreviousPassImp PreviousPassType;
+    typedef PreviousPassImp                                              PreviousPassType;
 
-    typedef typename BaseType::PassIds PassIds;
+    typedef typename BaseType::PassIds                                   PassIds;
 
     // Types from the base class
-    typedef typename BaseType::EntityType EntityType;
+    typedef typename BaseType::EntityType                                EntityType;
 
-    typedef typename BaseType::ArgumentType ArgumentType;
+    typedef typename BaseType::ArgumentType                              ArgumentType;
 
   private:
-   typedef typename DiscreteModelType::Selector Selector;
-   typedef typename Dune::tuple_element< 0, Selector >::type ArgumentIdType;
+   typedef typename DiscreteModelType::Selector                          Selector;
+   typedef typename Dune::tuple_element< 0, Selector >::type             ArgumentIdType;
    static const std::size_t argumentPosition
      = Dune::FirstTypeIndex< PassIds, ArgumentIdType >::type::value;
    typedef typename Dune::tuple_element< argumentPosition, ArgumentType >::type ArgumentFunctionPtrType;
 
   public:
-    typedef typename PreviousPassType::GlobalArgumentType ArgumentFunctionType;
-    typedef typename ArgumentFunctionType :: LocalFunctionType LocalFunctionType;
+    typedef typename PreviousPassType::GlobalArgumentType                 ArgumentFunctionType;
+    typedef typename ArgumentFunctionType::LocalFunctionType              LocalFunctionType;
 
     // Types from the traits
-    typedef typename DiscreteModelType::Traits::DestinationType DestinationType;
-    typedef typename DiscreteModelType::Traits::VolumeQuadratureType VolumeQuadratureType;
-    typedef typename DiscreteModelType::Traits::FaceQuadratureType FaceQuadratureType;
-    typedef typename DiscreteModelType::Traits::DiscreteFunctionSpaceType
-    DiscreteFunctionSpaceType;
-    typedef typename DiscreteFunctionSpaceType :: IteratorType IteratorType;
+    typedef typename DiscreteModelType::Traits::DestinationType           DestinationType;
+    typedef typename DiscreteModelType::Traits::VolumeQuadratureType      VolumeQuadratureType;
+    typedef typename DiscreteModelType::Traits::FaceQuadratureType        FaceQuadratureType;
+    typedef typename DiscreteModelType::Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
+    typedef typename DiscreteFunctionSpaceType::IteratorType              IteratorType;
 
     // Types extracted from the discrete function space type
-    typedef typename DiscreteFunctionSpaceType::GridType GridType;
-    typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
-    typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
-    typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
+    typedef typename DiscreteFunctionSpaceType::GridType                  GridType;
+    typedef typename DiscreteFunctionSpaceType::GridPartType              GridPartType;
+    typedef typename DiscreteFunctionSpaceType::DomainType                DomainType;
+    typedef typename DiscreteFunctionSpaceType::DomainFieldType           DomainFieldType;
     typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
-    typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
-    typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
+    typedef typename DiscreteFunctionSpaceType::RangeFieldType            RangeFieldType;
+    typedef typename DiscreteFunctionSpaceType::JacobianRangeType         JacobianRangeType;
 
-    typedef typename GridType :: Traits :: LocalIdSet LocalIdSetType;
-    typedef typename LocalIdSetType :: IdType IdType;
+    typedef typename GridType::Traits::LocalIdSet                         LocalIdSetType;
+    typedef typename LocalIdSetType::IdType                               IdType;
 
     // Types extracted from the underlying grids
-    typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
-    typedef typename IntersectionIteratorType :: Intersection IntersectionType;
-    typedef typename GridPartType::template Codim<0>::GeometryType       Geometry;
-    typedef typename Geometry::LocalCoordinate LocalDomainType;
+    typedef typename GridPartType::IntersectionIteratorType               IntersectionIteratorType;
+    typedef typename GridPartType::IntersectionType                       IntersectionType;
+    typedef typename GridPartType::template Codim<0>::GeometryType        Geometry;
+    typedef typename Geometry::LocalCoordinate                            LocalDomainType;
 
     // Various other types
-    typedef typename DestinationType::LocalFunctionType DestLocalFunctionType;
+    typedef typename DestinationType::LocalFunctionType                   DestLocalFunctionType;
 
     typedef LimiterDiscreteModelCaller< DiscreteModelType, ArgumentType, PassIds > DiscreteModelCallerType;
 
     // type of Communication Manager
-    typedef CommunicationManager< DiscreteFunctionSpaceType > CommunicationManagerType;
+    typedef CommunicationManager< DiscreteFunctionSpaceType >             CommunicationManagerType;
 
     // Range of the destination
     enum { dimRange = DiscreteFunctionSpaceType::dimRange,
-     dimDomain = DiscreteFunctionSpaceType::dimDomain};
+           dimDomain = DiscreteFunctionSpaceType::dimDomain};
     enum { dimGrid = GridType :: dimension };
-    typedef typename GridType :: ctype ctype;
-    typedef FieldVector<ctype, dimGrid-1> FaceLocalDomainType;
+    typedef typename GridType::ctype                                       ctype;
+    typedef FieldVector<ctype, dimGrid-1>                                  FaceLocalDomainType;
 
-    typedef PointBasedDofConversionUtility< dimRange > DofConversionUtilityType;
+    typedef PointBasedDofConversionUtility< dimRange >                     DofConversionUtilityType;
 
     static const bool StructuredGrid     = GridPartCapabilities::isCartesian< GridPartType >::v;
     static const bool conformingGridPart = GridPartCapabilities::isConforming< GridPartType >::v;
 
-    typedef FieldVector< DomainType , dimRange > DeoModType;
-    typedef FieldMatrix< DomainFieldType, dimDomain , dimDomain > MatrixType;
+    typedef FieldVector< DomainType , dimRange >                           DeoModType;
+    typedef FieldMatrix< DomainFieldType, dimDomain , dimDomain >          MatrixType;
 
-    typedef typename GridPartType :: IndexSetType IndexSetType;
-    typedef AllGeomTypes< IndexSetType, GridType> GeometryInformationType;
+    typedef typename GridPartType :: IndexSetType                          IndexSetType;
+    typedef AllGeomTypes< IndexSetType, GridType>                          GeometryInformationType;
 
-    typedef GeometryInformation< GridType, 1 > FaceGeometryInformationType;
+    typedef GeometryInformation< GridType, 1 >                             FaceGeometryInformationType;
 
     // get LagrangePointSet of pol order 1
-    typedef LagrangePointSet< GridPartType, 1 > LagrangePointSetType;
+    typedef LagrangePointSet< GridPartType, 1 >                            LagrangePointSetType;
     // get lagrange point set of order 1
-    typedef CompiledLocalKeyContainer< LagrangePointSetType, 1 , 1 >
-      LagrangePointSetContainerType;
+    typedef CompiledLocalKeyContainer< LagrangePointSetType, 1 , 1 >       LagrangePointSetContainerType;
 
-    typedef DGFEntityKey<int> KeyType;
-    typedef std::vector<int> CheckType;
-    typedef std::pair< KeyType, CheckType > VectorCompType;
-    typedef std::set< VectorCompType > ComboSetType;
+    typedef DGFEntityKey<int>                                              KeyType;
+    typedef std::vector<int>                                               CheckType;
+    typedef std::pair< KeyType, CheckType >                                VectorCompType;
+    typedef std::set< VectorCompType >                                     ComboSetType;
 
     struct MatrixIF
     {
@@ -1100,11 +1098,10 @@ namespace Fem
 
         elementCounter_ = 0;
         // dod limitation
-        IteratorType endit = spc_.end();
-        for (IteratorType it = spc_.begin(); it != endit; ++it)
+        for( const auto& en : elements( spc_.gridPart() ) )
         {
           Dune::Timer localTime;
-          applyLocalImp(*it);
+          applyLocalImp(en);
           stepTime_[2] += localTime.elapsed();
           ++elementCounter_;
         }
@@ -1252,10 +1249,8 @@ namespace Fem
       {
         // check whether on of the intersections is with ghost element
         // and if so, skip the computation of the limited solution for now
-        const IntersectionIteratorType endnit = gridPart_.iend(entity);
-        for (IntersectionIteratorType nit = gridPart_.ibegin(entity); nit != endnit; ++nit)
+        for (const auto& intersection : intersections(gridPart_, entity) )
         {
-          const IntersectionType& intersection = *nit;
           if( intersection.neighbor() )
           {
             // get neighbor
@@ -1904,11 +1899,8 @@ namespace Fem
           if( ! checkPhysicalQuad(volQuad, uEn) ) return false;
         }
 
-        const IntersectionIteratorType endnit = gridPart_.iend(en);
-        for (IntersectionIteratorType nit = gridPart_.ibegin(en);
-             nit != endnit; ++nit)
+        for (const auto& intersection : intersections(gridPart_, en) )
         {
-          const IntersectionType& intersection = *nit;
           if( intersection.neighbor() && ! intersection.conforming() )
           {
             typedef typename FaceQuadratureType :: NonConformingQuadratureType NonConformingQuadratureType;
@@ -2374,12 +2366,8 @@ namespace Fem
       shockIndicator = 0;
       adaptIndicator = 0;
 
-      const IntersectionIteratorType endnit = gridPart_.iend(en);
-      for (IntersectionIteratorType niter = gridPart_.ibegin(en);
-           niter != endnit; ++niter)
+      for (const auto& intersection : intersections(gridPart_, en) )
       {
-        const IntersectionType& intersection = *niter ;
-
         typedef typename IntersectionType :: Geometry LocalGeometryType;
         const LocalGeometryType& interGeo = intersection.geometry();
 
diff --git a/dune/fem-dg/pass/dgmasspass.hh b/dune/fem-dg/pass/dgmasspass.hh
index ad702ffb486cc317ec2bebed198699917e3daaa9..8a07ac9613d7ec5ac83d37429bafac7d58df28e1 100644
--- a/dune/fem-dg/pass/dgmasspass.hh
+++ b/dune/fem-dg/pass/dgmasspass.hh
@@ -213,12 +213,9 @@ namespace Dune
         // set pointer
         prepare( argument, destination, true );
 
-        typedef typename GridPartType::template Codim< 0 >::template Partition< All_Partition >::IteratorType IteratorType;
-        const GridPartType &gridPart = scalarSpace_.gridPart();
-        const IteratorType end = gridPart.template end< 0, All_Partition >();
-        for( IteratorType it = gridPart.template begin< 0, All_Partition >(); it != end; ++it )
+        for( const auto& en : elements( scalarSpace_.gridPart(), Dune::Partitions::all ) )
         {
-          applyLocal( *it, argument, destination );
+          applyLocal( en, argument, destination );
         }
 
         // remove pointer
diff --git a/dune/fem-dg/pass/dgpass.hh b/dune/fem-dg/pass/dgpass.hh
index ebb655edd8418f902372f17a08557338186dc387..a5e7cb3ceb1dd323111ebfdab6154ffc51f10223 100644
--- a/dune/fem-dg/pass/dgpass.hh
+++ b/dune/fem-dg/pass/dgpass.hh
@@ -59,64 +59,63 @@ namespace Fem
     //- Typedefs and enums
     //! Base class
     typedef LocalPass< DiscreteModelImp , PreviousPassImp , passIdImp > BaseType;
-    typedef typename BaseType::PassIds PassIds;
+    typedef typename BaseType::PassIds                                  PassIds;
 
     //! Repetition of template arguments
-    typedef DiscreteModelImp DiscreteModelType;
-    typedef PreviousPassImp PreviousPassType;
+    typedef DiscreteModelImp                                            DiscreteModelType;
+    typedef PreviousPassImp                                             PreviousPassType;
 
     // Types from the base class
-    typedef typename BaseType::EntityType  EntityType;
-    typedef typename BaseType::ArgumentType ArgumentType;
+    typedef typename BaseType::EntityType                               EntityType;
+    typedef typename BaseType::ArgumentType                             ArgumentType;
 
     // Types from the traits
-    typedef typename DiscreteModelType::Traits::DestinationType DestinationType;
-    typedef typename DiscreteModelType::Traits::VolumeQuadratureType VolumeQuadratureType;
-    typedef typename DiscreteModelType::Traits::FaceQuadratureType FaceQuadratureType;
+    typedef typename DiscreteModelType::Traits::DestinationType         DestinationType;
+    typedef typename DiscreteModelType::Traits::VolumeQuadratureType    VolumeQuadratureType;
+    typedef typename DiscreteModelType::Traits::FaceQuadratureType      FaceQuadratureType;
     typedef typename DiscreteModelType::Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
     //! Iterator over the space
-    typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
+    typedef typename DiscreteFunctionSpaceType::IteratorType            IteratorType;
 
     // Types extracted from the discrete function space type
-    typedef typename DiscreteFunctionSpaceType::GridType GridType;
-    typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
-    typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
-    typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
-    typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
-    typedef typename DiscreteFunctionSpaceType::RangeFieldType  RangeFieldType;
-    typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
-    typedef typename DiscreteFunctionSpaceType:: BasisFunctionSetType
-      BasisFunctionSetType;
+    typedef typename DiscreteFunctionSpaceType::GridType                GridType;
+    typedef typename DiscreteFunctionSpaceType::GridPartType            GridPartType;
+    typedef typename DiscreteFunctionSpaceType::DomainType              DomainType;
+    typedef typename DiscreteFunctionSpaceType::RangeType               RangeType;
+    typedef typename DiscreteFunctionSpaceType::DomainFieldType         DomainFieldType;
+    typedef typename DiscreteFunctionSpaceType::RangeFieldType          RangeFieldType;
+    typedef typename DiscreteFunctionSpaceType::JacobianRangeType       JacobianRangeType;
+    typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType    BasisFunctionSetType;
 
     // Types extracted from the underlying grid part types
-    typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
-    typedef typename IntersectionIteratorType::Intersection IntersectionType;
-    typedef typename GridPartType :: template Codim<0>:: GeometryType  Geometry;
+    typedef typename GridPartType::IntersectionIteratorType             IntersectionIteratorType;
+    typedef typename GridPartType::IntersectionType                     IntersectionType;
+    typedef typename GridPartType::template Codim<0>::GeometryType      Geometry;
 
 
     // Various other types
-    typedef typename DestinationType::LocalFunctionType LocalFunctionType;
+    typedef typename DestinationType::LocalFunctionType                 LocalFunctionType;
 
     typedef CDGDiscreteModelCaller< DiscreteModelType, ArgumentType, PassIds > DiscreteModelCallerType;
 
-    typedef typename DestinationType :: DofBlockPtrType DofBlockPtrType;
+    typedef typename DestinationType::DofBlockPtrType                   DofBlockPtrType;
 
     // type of local id set
-    typedef typename GridPartType::IndexSetType IndexSetType;
-    typedef Fem::TemporaryLocalFunction< DiscreteFunctionSpaceType > TemporaryLocalFunctionType;
+    typedef typename GridPartType::IndexSetType                         IndexSetType;
+    typedef Fem::TemporaryLocalFunction< DiscreteFunctionSpaceType >    TemporaryLocalFunctionType;
 
 #ifdef USE_CACHED_INVERSE_MASSMATRIX
     // type of communication manager object which does communication
     typedef typename DiscreteFunctionSpaceType::template ToNewDimRange< 1 >::Type ScalarDiscreteFunctionSpaceType;
     typedef Fem::DGMassInverseMassImplementation< ScalarDiscreteFunctionSpaceType, true > MassInverseMassType ;
-    typedef typename MassInverseMassType :: KeyType MassKeyType;
-    typedef Fem::SingletonList< MassKeyType, MassInverseMassType >  InverseMassProviderType;
+    typedef typename MassInverseMassType::KeyType                        MassKeyType;
+    typedef Fem::SingletonList< MassKeyType, MassInverseMassType >       InverseMassProviderType;
     // store a reference to the mass matrix implementation
-    typedef MassInverseMassType&  LocalMassMatrixStorageType ;
+    typedef MassInverseMassType&                                         LocalMassMatrixStorageType;
 #else
     //! type of local mass matrix
     typedef Fem::LocalMassMatrix< DiscreteFunctionSpaceType, VolumeQuadratureType > LocalMassMatrixType;
-    typedef LocalMassMatrixType  LocalMassMatrixStorageType;
+    typedef LocalMassMatrixType                                          LocalMassMatrixStorageType;
 #endif
 
     // true if all intersections between element in this grid part are conforming
@@ -221,10 +220,9 @@ namespace Fem
 
       if( reallyCompute_ )
       {
-        const IteratorType endit = spc_.end();
-        for (IteratorType it = spc_.begin(); it != endit; ++it)
+        for( const auto& en : elements( spc_.gridPart() ) )
         {
-          applyLocal(*it);
+          applyLocal(en);
         }
 
         finalize(arg, dest);
@@ -557,11 +555,8 @@ namespace Fem
         // get volume of element divided by the DG polynomial factor
         const double envol = entity.geometry().volume() / ( 2.0 * spc_.order( entity ) + 1.0 ) ;
 
-        const IntersectionIteratorType endnit = gridPart_.iend(entity);
-        for (IntersectionIteratorType nit = gridPart_.ibegin(entity); nit != endnit; ++nit)
+        for (const auto& intersection : intersections(gridPart_, entity) )
         {
-          const IntersectionType& intersection = *nit;
-
           double nbvol = envol;
           double wspeedS = 0.0;
 
diff --git a/dune/fem-dg/pass/threadpass.hh b/dune/fem-dg/pass/threadpass.hh
index 5078e00285e45bfe5eccb7266c21481047f4610f..c04129bef2f7b9062f832f10f30cdd34a36a6795 100644
--- a/dune/fem-dg/pass/threadpass.hh
+++ b/dune/fem-dg/pass/threadpass.hh
@@ -151,51 +151,51 @@ namespace Fem
   {
     typedef ThreadPass< InnerPass, ThreadIterator, nonblockingcomm > ThisType;
   public:
-    typedef InnerPass InnerPassType;
-    typedef typename InnerPass :: DiscreteModelType  DiscreteModelType;
-    typedef typename InnerPass :: PreviousPassType   PreviousPassType;
+    typedef InnerPass                                                     InnerPassType;
+    typedef typename InnerPass::DiscreteModelType                         DiscreteModelType;
+    typedef typename InnerPass::PreviousPassType                          PreviousPassType;
 
     //- Typedefs and enums
     //! Base class
     typedef LocalPass< DiscreteModelType, PreviousPassType, InnerPass :: passId>  BaseType;
 
     // Types from the base class
-    typedef typename BaseType::EntityType  EntityType;
-    typedef typename BaseType::ArgumentType ArgumentType;
+    typedef typename BaseType::EntityType                                 EntityType;
+    typedef typename BaseType::ArgumentType                               ArgumentType;
 
     // Types from the traits
-    typedef typename DiscreteModelType::Traits::DestinationType DestinationType;
-    typedef typename DiscreteModelType::Traits::VolumeQuadratureType VolumeQuadratureType;
-    typedef typename DiscreteModelType::Traits::FaceQuadratureType FaceQuadratureType;
+    typedef typename DiscreteModelType::Traits::DestinationType           DestinationType;
+    typedef typename DiscreteModelType::Traits::VolumeQuadratureType      VolumeQuadratureType;
+    typedef typename DiscreteModelType::Traits::FaceQuadratureType        FaceQuadratureType;
     typedef typename DiscreteModelType::Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
     //! Iterator over the space
-    typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
+    typedef typename DiscreteFunctionSpaceType::IteratorType              IteratorType;
 
     // Types extracted from the discrete function space type
-    typedef typename DiscreteFunctionSpaceType::GridType GridType;
-    typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
-    typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
-    typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
-    typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
+    typedef typename DiscreteFunctionSpaceType::GridType                  GridType;
+    typedef typename DiscreteFunctionSpaceType::GridPartType              GridPartType;
+    typedef typename DiscreteFunctionSpaceType::DomainType                DomainType;
+    typedef typename DiscreteFunctionSpaceType::RangeType                 RangeType;
+    typedef typename DiscreteFunctionSpaceType::JacobianRangeType         JacobianRangeType;
 
     // Types extracted from the underlying grids
-    typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
-    typedef typename IntersectionIteratorType::Intersection IntersectionType;
-    typedef typename GridType::template Codim<0>::Geometry Geometry;
+    typedef typename GridPartType::IntersectionIteratorType               IntersectionIteratorType;
+    typedef typename GridPartType::IntersectionType                       IntersectionType;
+    typedef typename GridType::template Codim<0>::Geometry                Geometry;
 
     // Various other types
-    typedef typename DestinationType::LocalFunctionType LocalFunctionType;
+    typedef typename DestinationType::LocalFunctionType                   LocalFunctionType;
 
-    typedef NonBlockingCommHandle< DestinationType > NonBlockingCommHandleType ;
+    typedef NonBlockingCommHandle< DestinationType >                      NonBlockingCommHandleType ;
 
     // Range of the destination
     enum { dimRange = DiscreteFunctionSpaceType::dimRange };
 
     // type of local id set
-    typedef typename GridPartType::IndexSetType IndexSetType;
+    typedef typename GridPartType::IndexSetType                           IndexSetType;
 
     // type of thread iterators (e.g. Fem::DomainDecomposedIteratorStorage or Fem::ThreadIterator)
-    typedef ThreadIterator  ThreadIteratorType;
+    typedef ThreadIterator                                                ThreadIteratorType;
 
   protected:
     using BaseType :: spc_;