parabolicestimator.hh

#ifndef __PARABOLIC_ESTIMATOR_HH__
#define __PARABOLIC_ESTIMATOR_HH__
 
//- Dune-fem includes
#include <dune/fem/quadrature/caching/twistutility.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem/quadrature/intersectionquadrature.hh>
#include <dune/fem/operator/common/spaceoperatorif.hh>
#include <dune/fem/operator/matrix/blockmatrix.hh>
#include <dune/fem/space/discontinuousgalerkin.hh>
 
// Local includes
#include "../estimators/estimatorinterface.hh"
#include "../models/modelinterface.hh"
#include "../common/geometry.hh"
#include "../common/intersectiondatahandle.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class OldSolutionFunction, class TimeProvider, class ImplicitModel, class ExplicitModel>
    class ParabolicEulerEstimator
      : public DefaultEstimator<ParabolicEulerEstimator<OldSolutionFunction, TimeProvider,
                                                        ImplicitModel, ExplicitModel> >
    {
      enum FunctionSpaceNorm
      {
        L2Norm,
        H1Norm,
        Norm = L2Norm
      };
      typedef ParabolicEulerEstimator    ThisType;
      typedef DefaultEstimator<ThisType> BaseType;
     public:
      // Interface types
      typedef typename OldSolutionFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
      typedef typename DiscreteFunctionSpaceType::GridPartType   GridPartType;
      typedef typename GridPartType::GridType                    GridType;
      typedef typename GridType::template Codim<0>::Entity       ElementType;
     protected:
      typedef std::vector<RangeFieldType> ErrorIndicatorType;
     public:
      typedef typename ErrorIndicatorType::const_iterator IteratorType;
 
     protected:
      // Everything else why might need
      typedef TimeProvider  TimeProviderType;
      typedef ModelInterface<ImplicitModel> ImplicitModelType;
      typedef ModelInterface<ExplicitModel> ExplicitModelType;
      typedef typename ImplicitModelType::OperatorPartsType ImplicitOperatorPartsType;
      typedef typename ExplicitModelType::OperatorPartsType ExplicitOperatorPartsType;
      typedef ModelConstituents<ImplicitModel> ImplicitConstituents;
      typedef ModelConstituents<ExplicitModel> ExplicitConstituents;
 
      typedef OldSolutionFunction OldSolutionFunctionType;
 
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
      typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
      typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename DiscreteFunctionSpaceType::HessianRangeType HessianRangeType;
 
      // Various data functions ...
      typedef typename ImplicitModelType::BulkForcesFunctionType        BulkForcesFunctionType;
      typedef typename BulkForcesFunctionType::LocalFunctionType        BulkForcesLocalFunctionType;
      typedef typename ImplicitModelType::DirichletBoundaryFunctionType DirichletFunctionType;
      typedef typename DirichletFunctionType::LocalFunctionType         DirichletLocalFunctionType;
      typedef typename ImplicitModelType::NeumannBoundaryFunctionType   NeumannFunctionType;
      typedef typename NeumannFunctionType::LocalFunctionType           NeumannLocalFunctionType;
 
      // Boundary classification
      typedef typename ImplicitModelType::DirichletIndicatorType DirichletIndicatorType;
      typedef typename ImplicitModelType::NeumannIndicatorType   NeumannIndicatorType;
 
      // Just to have a draft what might be needed to take contributions
      // from the explicit model into account
      typedef typename ExplicitModelType::BulkForcesFunctionType         ExplicitBulkForcesFunctionType;
      typedef typename ExplicitBulkForcesFunctionType::LocalFunctionType ExplicitBulkForcesLocalFunctionType;
      // and so on, but even the above is not really needed here
 
      //typedef typename GridPartType::GridType GridType;
      typedef typename DiscreteFunctionSpaceType::IteratorType GridIteratorType;
      typedef typename GridPartType::IndexSetType IndexSetType;
      typedef typename IndexSetType::IndexType IndexType;
      typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType;
 
      typedef typename IntersectionIteratorType::Intersection IntersectionType;
 
      //typedef typename GridType::template Codim<0>::Entity ElementType; see above
      typedef typename ElementType::Geometry GeometryType;
      static const int dimension = GridType::dimension;
 
      typedef Dune::Fem::CachingQuadrature<GridPartType, 0> ElementQuadratureType;
      typedef Dune::Fem::CachingQuadrature<GridPartType, 1> FaceQuadratureType;
 
      //typedef Dune::Fem::ElementQuadrature<GridPartType, 0> ElementQuadratureType;
      //typedef Dune::Fem::ElementQuadrature<GridPartType, 1> FaceQuadratureType;
 
      typedef Dune::FieldMatrix<DomainFieldType,dimension,dimension> JacobianInverseType;
 
     protected:
      const TimeProviderType&  timeProvider_;
      const ImplicitModelType& implicitModel_;
      const ImplicitOperatorPartsType implicitParts_;
      const ImplicitOperatorPartsType implicitPartsNeighbour_;
      const ExplicitModelType& explicitModel_;
      const ExplicitOperatorPartsType explicitParts_;
      const OldSolutionFunctionType& oldSolution_;
 
      const DiscreteFunctionSpaceType &dfSpace_;
      GridPartType &gridPart_;
      const IndexSetType &indexSet_;
      GridType &grid_;
 
     private:
      BulkForcesFunctionType         bulkForcesFunction_;
      ExplicitBulkForcesFunctionType explicitBulkForcesFunction_;
      NeumannFunctionType            neumannFunction_;
      NeumannIndicatorType           neumannBndry_;
      ErrorIndicatorType             localIndicators_;
      RangeFieldType                 timeIndicator_;
 
      mutable std::vector<JacobianRangeType> outsideFaceFluxCache;
 
      enum { bulk, jump, bdry, comm, numContributions };
      RangeFieldType errorMin[numContributions];
      RangeFieldType errorMax[numContributions];
 
      typedef RangeType DataItemType;
      typedef std::vector<DataItemType> IntersectionStorageType;
      typedef std::map<IndexType, IntersectionStorageType> CommunicationStorageType;
      typedef typename Fem::DFCommunicationOperation::Add OperationType;
      typedef ACFem::IntersectionDataHandle<IndexSetType, CommunicationStorageType, OperationType> DataHandleType;
 
      class LocalIntersectionStorage
      {
        typedef std::vector<RangeFieldType> WeightStorageType;
       public:
        LocalIntersectionStorage()
        {}
        LocalIntersectionStorage(IndexType bulkEntityIndex, size_t nop)
          : bulkEntityIndex_(bulkEntityIndex), weights_(nop), h_(0)
        {}
        IndexType bulkEntityIndex() const { return bulkEntityIndex_; }
        IndexType& bulkEntityIndex() { return bulkEntityIndex_; }
        const RangeFieldType& operator[](size_t idx) const { return weights_[idx]; }
        RangeFieldType& operator[](size_t idx) { return weights_[idx]; }
        size_t size() const { return weights_.size(); }
        RangeFieldType& hEstimate() { return h_; }
        const RangeFieldType& hEstimate() const { return h_; }
       private:
        IndexType bulkEntityIndex_;
        WeightStorageType weights_;
        RangeFieldType h_;
      };
      typedef std::map<IndexType, LocalIntersectionStorage> LocalStorageType;
 
      LocalStorageType localJumpStorage_;
      CommunicationStorageType commJumpStorage_;
 
     public:
      explicit ParabolicEulerEstimator(const TimeProviderType& timeProvider,
                                       const ImplicitModelType& implicitModel,
                                       const ExplicitModelType& explicitModel,
                                       const OldSolutionFunctionType& oldSolution)
        : timeProvider_(timeProvider),
          implicitModel_(implicitModel),
          implicitParts_(implicitModel_.operatorParts()),
          implicitPartsNeighbour_(implicitModel_.operatorParts()),
          explicitModel_(explicitModel),
          explicitParts_(explicitModel_.operatorParts()),
          oldSolution_(oldSolution),
          dfSpace_(oldSolution_.space()),
          gridPart_(dfSpace_.gridPart()),
        indexSet_(gridPart_.indexSet()),
        grid_(gridPart_.grid()),
        bulkForcesFunction_(implicitModel_.bulkForcesFunction(gridPart_)),
        explicitBulkForcesFunction_(explicitModel_.bulkForcesFunction(gridPart_)),
        neumannFunction_(implicitModel_.neumannBoundaryFunction(gridPart_)),
        neumannBndry_(implicitModel_.neumannIndicator()),
        timeIndicator_(0)
      {
        localIndicators_.resize(indexSet_.size(0));
      }
 
      template<class DiscreteFunctionType>
      RangeFieldType estimate(const DiscreteFunctionType& uh)
      {
        assert(!(ModelConstituents<ImplicitModelType>::hasForcesFunctional ||
                 ModelConstituents<ExplicitModelType>::hasForcesFunctional));
 
        // clear all local estimators
        clear();
 
        // calculate local estimator
        const GridIteratorType end = dfSpace_.end();
        for (GridIteratorType it = dfSpace_.begin(); it != end; ++it) {
          estimateLocal(*it, uh);
        }
 
        RangeFieldType error = 0.0;
 
        // sum up local estimators
        const IteratorType endEstimator = localIndicators_.end();
        for (IteratorType it = localIndicators_.begin(); it != endEstimator; ++it) {
          error += *it;
        }
 
        // obtain global sum
        RangeFieldType commarray[] = { error, timeIndicator_ };
        grid_.comm().sum(commarray, 2);
        error = commarray[0];
        timeIndicator_ = commarray[1];
 
        if (Dune::Fem::Parameter::verbose()) {
          // print error estimator if on node 0
          std::cout << "Estimated " << (Norm == L2Norm ? "L2" : "H1") << "-Error: " << std::sqrt(error) << std::endl;
        }
 
        timeIndicator_ = std::sqrt(timeIndicator_);
 
        return std::sqrt(error);
      }
 
      const RangeFieldType& operator[](const size_t& idx) const
      {
        return localIndicators_[idx];
      }
 
      const RangeFieldType& operator[](const ElementType& entity) const
      {
        return (*this)[indexSet_.index(entity)];
      }
 
      size_t size() const
      {
        return indexSet_.size(0);
      }
 
      IteratorType begin() const {
        return localIndicators_.begin();
      }
 
      IteratorType end() const {
        return localIndicators_.end();
      }
 
      RangeFieldType timeEstimate() const {
        return timeIndicator_;
      }
 
      GridType& grid() const
      {
        return grid_;
      }
 
      const DiscreteFunctionSpaceType& space() const
      {
        return dfSpace_;
      }
 
     protected:
      void clear()
      {
        // resize and clear
        localIndicators_.resize(indexSet_.size(0));
        const typename ErrorIndicatorType::iterator end = localIndicators_.end();
        for (typename ErrorIndicatorType::iterator  it = localIndicators_.begin();
             it != end; ++it) {
          *it = 0.0;
        }
        timeIndicator_ = 0;
 
        // the communication hack for parallel runs. Gnah.
        localJumpStorage_.clear();
        commJumpStorage_.clear();
 
        // debugging and stuff
        errorMin[bulk] =
          errorMin[jump] =
          errorMin[bdry] =
          errorMin[comm] = std::numeric_limits<RangeFieldType>::max();
 
        errorMax[bulk] =
          errorMax[jump] =
          errorMax[bdry] =
          errorMax[comm] = 0.0;
      }
 
      template<class DiscreteFunctionType>
      void estimateLocal(const ElementType &entity, const DiscreteFunctionType& uh)
      {
        typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
        typedef typename OldSolutionFunctionType::LocalFunctionType ExplicitLocalFunctionType;
 
        const typename ElementType::Geometry &geometry = entity.geometry();
 
        const RangeFieldType volume = geometry.volume();
        const RangeFieldType h2 = (dimension == 2 ? volume : std::pow(volume, 2.0 / (RangeFieldType)dimension));
        const RangeFieldType hFactor = Norm == L2Norm ? h2 * h2 : h2;
        const int index = indexSet_.index(entity);
        const LocalFunctionType           uLocal = uh.localFunction(entity);
        const ExplicitLocalFunctionType   uOldLocal = oldSolution_.localFunction(entity);
        const BulkForcesLocalFunctionType fLocal = bulkForcesFunction_.localFunction(entity);
        const ExplicitBulkForcesLocalFunctionType fOldLocal = explicitBulkForcesFunction_.localFunction(entity);
 
        implicitParts_.setEntity(entity);
        explicitParts_.setEntity(entity);
 
        ElementQuadratureType quad(entity, 2*(dfSpace_.order() + 1));
 
        const int numQuadraturePoints = quad.nop();
        for (int qp = 0; qp < numQuadraturePoints; ++qp) {
          const typename ElementQuadratureType::CoordinateType &x = quad.point(qp); //
          const RangeFieldType weight = quad.weight(qp) * geometry.integrationElement(x);
 
          RangeType values;
          uLocal.evaluate(x, values);
 
          JacobianRangeType jacobian;
          uLocal.jacobian(x, jacobian);
 
          HessianRangeType hessian;
          uLocal.hessian(x, hessian);
 
          // Now compute the element residual
          RangeType el_res = 0;
          if (ImplicitConstituents::hasFlux) {
            RangeType tmp;
            implicitParts_.fluxDivergence(entity, quad[qp], values, jacobian, hessian, tmp);
            el_res += tmp;
          }
          if (ImplicitConstituents::hasSources) {
            RangeType tmp;
            implicitParts_.source(entity, quad[qp], values, jacobian, tmp);
            el_res += tmp;
          }
          if (ImplicitConstituents::hasBulkForces) {
            RangeType tmp;
            fLocal.evaluate(x, tmp);
            el_res -= tmp;
          }
 
          // Compute the contribution from the explicit part
          RangeType oldValues;
          uOldLocal.evaluate(x, oldValues);
 
          JacobianRangeType oldJacobian;
          uOldLocal.jacobian(x, oldJacobian);
 
          HessianRangeType oldHessian;
          uOldLocal.hessian(x, oldHessian);
 
          // add the "explicit" part. In principle overkill here
          if (ExplicitConstituents::hasFlux) {
            RangeType tmp;
            explicitParts_.fluxDivergence(entity, quad[qp], values, jacobian, hessian, tmp);
            el_res += tmp;
          }
          if (ExplicitConstituents::hasSources) {
            RangeType tmp;
            explicitParts_.source(entity, quad[qp], values, jacobian, tmp);
            el_res += tmp;
          }
          if (ExplicitConstituents::hasBulkForces) {
            // should not happen here, in principle ...
            RangeType tmp;
            fOldLocal.evaluate(x, tmp);
            el_res -= tmp;
          }
 
          localIndicators_[index] += hFactor * weight * (el_res * el_res);
 
          // Now el_res contains the space error. We also need the error
          // from the time discretization
          RangeType deltaU = values;
          deltaU -= oldValues;
 
          timeIndicator_ += weight * (deltaU * deltaU);
        }
 
        // Calculate face contribution. The explicit part is completely ignored here.
        IntersectionIteratorType end = gridPart_.iend(entity);
        for (IntersectionIteratorType it = gridPart_.ibegin(entity); it != end; ++it) {
          const IntersectionType &intersection = *it;
 
          bool isRobin = implicitParts_.setIntersection(intersection);
 
          if (isProcessorBoundary(intersection)) {
            // interior face accross a process boundary
            estimateProcessBoundary(intersection, entity, uLocal);
          } else if (intersection.neighbor()) {
            // interior face
            estimateIntersection(intersection, entity, uh, uLocal);
          } else if (intersection.boundary()) {
            // boundary face
            estimateBoundary(intersection, entity, uLocal,
                             neumannBndry_.applies(intersection),
                             isRobin);
          }
        }
      }
 
      template<class DiscreteFunctionType>
      void estimateIntersection(const IntersectionType &intersection,
                                const ElementType &inside,
                                const DiscreteFunctionType& uh,
                                const typename DiscreteFunctionType::LocalFunctionType &uInside)
      {
        typedef typename DiscreteFunctionType::LocalFunctionType LocalFunctionType;
 
        // get outside entity pointer
        const auto outside = intersection.outside();
 
        const int insideIndex = indexSet_.index(inside);
        const int outsideIndex = indexSet_.index(outside);
 
        const bool isOutsideInterior = (outside.partitionType() == Dune::InteriorEntity);
        if (!isOutsideInterior || (insideIndex < outsideIndex)) {
          const LocalFunctionType uOutside = uh.localFunction(outside);
 
          RangeFieldType error;
 
          if (intersection.conforming())
            error = estimateIntersection<true>(intersection, inside, uInside, outside, uOutside);
          else
            error = estimateIntersection<false>(intersection, inside, uInside, outside, uOutside);
 
          if (error > 0.0) {
            const RangeFieldType volume
              = 0.5 * (inside.geometry().volume() + outside.geometry().volume());
            const RangeFieldType h = (dimension == 1 ? volume : std::pow(volume, 1.0 / (RangeFieldType)dimension));
            const RangeFieldType hFactor = Norm == L2Norm ? h * h * h : h;
            localIndicators_[insideIndex] += hFactor * error;
            if (isOutsideInterior) {
              localIndicators_[outsideIndex] += hFactor * error;
            }
          }
        }
      }
 
      template<bool conforming, class LocalFunctionType>
      RangeFieldType estimateIntersection(const IntersectionType &intersection,
                                  const ElementType &inside,
                                  const LocalFunctionType &uInside,
                                  const ElementType &outside,
                                  const LocalFunctionType &uOutside) const
      {
        // make sure correct method is called
        assert(intersection.conforming() == conforming);
 
        // use IntersectionQuadrature to create appropriate face quadratures
        typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, conforming> IntersectionQuadratureType;
        typedef typename IntersectionQuadratureType::FaceQuadratureType Quadrature ;
 
        const int quadOrder = 2 * (dfSpace_.order() - 1);
        // create intersection quadrature
        IntersectionQuadratureType intersectionQuad(gridPart_, intersection, quadOrder);
 
        // get appropriate quadrature references
        const Quadrature &quadInside  = intersectionQuad.inside();
        const Quadrature &quadOutside = intersectionQuad.outside();
 
        const int numQuadraturePoints = quadInside.nop();
 
        implicitPartsNeighbour_.setEntity(outside); // go to neighbour
 
        RangeFieldType error = 0.0;
        for (int qp = 0; qp < numQuadraturePoints; ++qp) {
          DomainType integrationNormal
            = intersection.integrationOuterNormal(quadInside.localPoint(qp));
          const RangeFieldType integrationElement = integrationNormal.two_norm();
 
          // evaluate | (d u_l * n_l) + (d u_r * n_r) | = | (d u_l - d u_r) * n_l |
          // evaluate jacobian left and right of the interface and apply
          // the diffusion tensor
 
          RangeType valueInside;
          JacobianRangeType jacobianInside;
          JacobianRangeType fluxInside;
          uInside.evaluate(quadInside[qp], valueInside);
          uInside.jacobian(quadInside[qp], jacobianInside);
          implicitParts_.flux(inside, quadInside[qp], valueInside, jacobianInside,
                              fluxInside);
 
          RangeType valueOutside;
          JacobianRangeType jacobianOutside;
          JacobianRangeType fluxOutside;
          uOutside.evaluate(quadOutside[qp], valueOutside);
          uOutside.jacobian(quadOutside[qp], jacobianOutside);
          implicitPartsNeighbour_.flux(outside, quadOutside[qp], valueOutside, jacobianOutside,
                                       fluxOutside);
 
          // Compute difference
          JacobianRangeType& fluxJump(fluxInside);
          fluxJump -= fluxOutside;
 
          // Multiply with normal
          RangeType jump;
          fluxJump.mv(integrationNormal, jump);
 
          error += quadInside.weight(qp) * (jump * jump) / integrationElement;
        }
        return error;
      }
 
      void communicateJumpContributions()
      {
        if (commJumpStorage_.size() == 0) {
          // not a parallel run.
          return;
        }
 
        const bool verbose = Dune::Fem::Parameter::getValue<int>("fem.verboserank", -1) != -1;
 
        DataHandleType dataHandle(indexSet_, commJumpStorage_);
 
        gridPart_.communicate(dataHandle,
                              InteriorBorder_InteriorBorder_Interface,
                              ForwardCommunication);
 
        // At this point commJumpStorage_ should contain the proper
        // differences of the boundary flux contributions. It remains
        // to performs the actual numerical quadrature summation.
 
        auto commEnd = commJumpStorage_.end();
        auto localEnd = localJumpStorage_.end();
        auto commIt = commJumpStorage_.begin();
        auto localIt = localJumpStorage_.begin();
 
        while (commIt != commEnd) {
          assert(commIt->first == localIt->first);
 
          const LocalIntersectionStorage& localData = localIt->second;
          const IntersectionStorageType& commData = commIt->second;
 
          const size_t numQuadraturePoints = localData.size();
          assert(numQuadraturePoints+1 == commData.size());
 
          // form the square sum
          RangeFieldType error = 0.0;
          for (unsigned qp = 0; qp < numQuadraturePoints; ++qp) {
            error += localData[qp] * (commData[qp] * commData[qp]);
          }
 
          // add the appropriate weights
          // HACK PLEASE FIXME;
          RangeFieldType h = localData.hEstimate() + commData[numQuadraturePoints][0];
          RangeFieldType hFactor = Norm == L2Norm ? h * h * h : h;
 
          // add contribution to the estimator array
          error *= hFactor;
          localIndicators_[localData.bulkEntityIndex()] += error;
 
          errorMin[comm] = std::min(errorMin[comm], error);
          errorMax[comm] = std::max(errorMax[comm], error);
 
          ++commIt;
          ++localIt;
        }
 
        assert(localIt == localEnd);
 
        if (verbose) {
          size_t nFaces = commJumpStorage_.size();
          nFaces = grid_.comm().sum(nFaces);
 
          if (Dune::Fem::Parameter::verbose()) {
            dwarn << "*** Processor Boundaries: "
                  << (RangeFieldType)nFaces/2.0
                  << " Faces ***" << std::endl;
          }
        }
      }
 
      template<class LocalFunctionType>
      void estimateProcessBoundary(const IntersectionType &intersection,
                                   const ElementType &inside,
                                   const LocalFunctionType &uInside)
      {
        if (intersection.conforming()) {
          estimateProcessBoundary<true>(intersection, inside, uInside);
        } else {
          estimateProcessBoundary<false>(intersection, inside, uInside);
        }
      }
 
      template<bool conforming, class LocalFunctionType>
      void estimateProcessBoundary(const IntersectionType &intersection,
                                   const ElementType &inside,
                                   const LocalFunctionType &uInside)
      {
        // make sure correct method is called
        assert(intersection.conforming() == conforming);
 
        // use IntersectionQuadrature to create appropriate face quadratures
        typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, conforming> IntersectionQuadratureType;
        typedef typename IntersectionQuadratureType::FaceQuadratureType Quadrature ;
 
        const int quadOrder = 2 * (dfSpace_.order() - 1);
        // create intersection quadrature
        IntersectionQuadratureType intersectionQuad(gridPart_, intersection, quadOrder);
 
        // get appropriate quadrature references
        const Quadrature &quadInside  = intersectionQuad.inside();
        const int numQuadraturePoints = quadInside.nop();
 
        // FIXME: is there any easier way to get hold of the face-index?
        IndexType faceIndex = indexSet_.index(inside.template subEntity<1>(intersection.indexInInside()));
        IndexType bulkIndex = indexSet_.index(inside);
 
        // IntersectionStorageType& commData =
        //   (commJumpStorage_[faceIndex] = IntersectionStorageType(numQuadraturePoints+1));
        // LocalIntersectionStorage& localData =
        //   (localJumpStorage_[faceIndex] = LocalIntersectionStorage(bulkIndex, numQuadraturePoints));
        IntersectionStorageType commData = IntersectionStorageType(numQuadraturePoints+1);
        LocalIntersectionStorage localData = LocalIntersectionStorage(bulkIndex, numQuadraturePoints);
 
        for (int qp = 0; qp < numQuadraturePoints; ++qp) {
          DomainType integrationNormal
            = intersection.integrationOuterNormal(quadInside.localPoint(qp));
          const RangeFieldType integrationElement = integrationNormal.two_norm();
          integrationNormal /= integrationElement;
 
          // evaluate | (d u_l * n_l) + (d u_r * n_r) | = | (d u_l - d u_r) * n_l |
          // evaluate jacobian left and right of the interface and apply
          // the diffusion tensor
 
          RangeType valueInside;
          JacobianRangeType jacobianInside;
          JacobianRangeType fluxInside;
          uInside.evaluate(quadInside[qp], valueInside);
          uInside.jacobian(quadInside[qp], jacobianInside);
          implicitParts_.flux(inside, quadInside[qp], valueInside, jacobianInside,
                              fluxInside);
 
          // Multiply with normal and store in the communication array
          fluxInside.mv(integrationNormal, commData[qp]);
 
          // Remember weight and integration element for this point
          localData[qp] = quadInside.weight(qp) * integrationElement;
        }
        localData.hEstimate() = hEstimate(inside.geometry());
        // THIS IS A HACK. PLEASE FIXME
        commData[numQuadraturePoints][0] = localData.hEstimate();
 
        commJumpStorage_[faceIndex] = commData;
        localJumpStorage_[faceIndex] = localData;
      }
 
      template<class LocalFunctionType>
      void estimateBoundary(const IntersectionType &intersection,
                            const ElementType &inside,
                            const LocalFunctionType &uLocal,
                            bool isNeumann, bool isRobin)
      {
        assert(intersection.conforming() == true);
 
        if (!isNeumann && !isRobin) {
          return;
        }
 
        // Maybe optimize away if not Neumann.
        const NeumannLocalFunctionType neumannLocal = neumannFunction_.localFunction(inside, intersection);
 
        typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, true /* conforming */> IntersectionQuadratureType;
        typedef typename IntersectionQuadratureType::FaceQuadratureType Quadrature ;
 
        // TODO: make more realistic assumptions on the
        // quad-order. HOWTO?
        int quadOrder = isRobin ? 2*dfSpace_.order() : 2*(dfSpace_.order() - 1);
 
        // create intersection quadrature
        IntersectionQuadratureType intersectionQuad(gridPart_, intersection, quadOrder);
        // get appropriate quadrature references
        const Quadrature &bndryQuad  = intersectionQuad.inside();
        const int numQuadraturePoints = bndryQuad.nop();
 
        RangeFieldType error = 0.0;
 
        for (int qp = 0; qp < numQuadraturePoints; ++qp) {
          const RangeFieldType weight = bndryQuad.weight(qp);
 
          // Get the un-normalized outer normal
          DomainType integrationNormal
            = intersection.integrationOuterNormal(bndryQuad.localPoint(qp));
 
          const RangeFieldType integrationElement = integrationNormal.two_norm();
          integrationNormal /= integrationElement;
 
          RangeType value;
          uLocal.evaluate(bndryQuad[qp], value);
          JacobianRangeType jacobian;
          uLocal.jacobian(bndryQuad[qp], jacobian);
 
          JacobianRangeType flux;
          implicitParts_.flux(inside, bndryQuad[qp], value, jacobian, flux);
 
          RangeType residual;
          flux.mv(integrationNormal, residual);
 
          if (isRobin) {
            RangeType alphau;
            implicitParts_.robinFlux(intersection, bndryQuad[qp], integrationNormal, value, alphau);
            residual += alphau;
          }
 
          if (isNeumann) {
            RangeType gN;
            neumannLocal.evaluate(bndryQuad[qp], gN);
            residual -= gN;
          }
 
          error += weight * integrationElement * (residual * residual);
        }
 
        if (error > 0.0) {
          // Get index into estimator vector
          const int insideIndex = indexSet_.index(inside);
 
          const RangeFieldType volume = inside.geometry().volume();
 
          const RangeFieldType h = (dimension == 1 ? volume : std::pow(volume, 1.0 / (RangeFieldType)dimension));
          const RangeFieldType hFactor = Norm == L2Norm ? h * h * h : h;
          localIndicators_[insideIndex] += hFactor * error;
        }
      }
    };
 
    template<class OldSolutionFunction, class TimeProvider, class ImplicitModel, class ExplicitModel>
    struct EstimatorTraits<ParabolicEulerEstimator<OldSolutionFunction, TimeProvider, ImplicitModel, ExplicitModel> >
    {
      typedef typename OldSolutionFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename
      std::vector<typename DiscreteFunctionSpaceType::RangeFieldType>::const_iterator
      IteratorType;
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
#endif // __PARABOLIC_ESTIMATOR_HH__