ellipticestimator.hh

#ifndef __ELLIPTIC_ESTIMATOR_HH__
#define __ELLIPTIC_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 DiscreteFunctionSpace, class Model, enum FunctionSpaceNorm Norm = H1Norm>
    class EllipticEstimator
      : public DefaultEstimator<EllipticEstimator<DiscreteFunctionSpace, Model, Norm> >
    {
      typedef EllipticEstimator<DiscreteFunctionSpace, Model> ThisType;
      typedef DefaultEstimator<ThisType> BaseType;
     public:
      // Interface types
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType; // aka double
      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 ModelInterface<Model> ModelType;
      typedef ModelConstituents<ModelType> Constituents;
      typedef typename ModelType::OperatorPartsType OperatorPartsType;
 
      typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
      typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
      typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;
 
      // Various data functions ...
      typedef typename ModelType::BulkForcesFunctionType         BulkForcesFunctionType;
      typedef typename BulkForcesFunctionType::LocalFunctionType BulkForcesLocalFunctionType;
      typedef typename ModelType::DirichletBoundaryFunctionType  DirichletFunctionType;
      typedef typename DirichletFunctionType::LocalFunctionType  DirichletLocalFunctionType;
      typedef typename ModelType::NeumannBoundaryFunctionType    NeumannFunctionType;
      typedef typename NeumannFunctionType::LocalFunctionType    NeumannLocalFunctionType;
 
      // Boundary classification
      typedef typename ModelType::DirichletIndicatorType DirichletIndicatorType;
      typedef typename ModelType::NeumannIndicatorType   NeumannIndicatorType;
 
      //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 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 ModelType& model_;              //< The underlying model.
      OperatorPartsType operatorParts_;
      OperatorPartsType operatorPartsNeighbour_;
 
      const DiscreteFunctionSpaceType &dfSpace_;
      GridPartType &gridPart_;
      const IndexSetType &indexSet_;
      GridType &grid_;
 
     private:
      BulkForcesFunctionType bulkForcesFunction_;
      NeumannFunctionType    neumannFunction_;
      NeumannIndicatorType   neumannBndry_;
      DirichletIndicatorType dirichletBndry_;
      ErrorIndicatorType     localIndicators_;
      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 EllipticEstimator(const ModelType& model, const DiscreteFunctionSpaceType& dfSpace)
        : model_(model),
          operatorParts_(model.operatorParts()),
          operatorPartsNeighbour_(model.operatorParts()),
          dfSpace_(dfSpace),
          gridPart_(dfSpace_.gridPart()),
          indexSet_(gridPart_.indexSet()),
          grid_(gridPart_.grid()),
          bulkForcesFunction_(model.bulkForcesFunction(gridPart_)),
          neumannFunction_(model.neumannBoundaryFunction(gridPart_)),
          neumannBndry_(model.neumannIndicator()),
          dirichletBndry_(model.dirichletIndicator())
      {
        localIndicators_.resize(indexSet_.size(0));
      }
 
      template<class DiscreteFunctionType>
      RangeFieldType estimate(const DiscreteFunctionType& uh)
      {
        assert(!ModelConstituents<ModelType>::hasForcesFunctional);
 
        // clear all local estimators
        clear();
 
        RangeFieldType hMax = 0.0;
 
        // calculate local estimator
        const GridIteratorType end = dfSpace_.end();
        for (GridIteratorType it = dfSpace_.begin(); it != end; ++it) {
          estimateLocal(*it, uh);
          hMax = std::max(hMax, hEstimate(it->geometry()));
        }
 
        // do the comm-stuff at inter-process boundaries
        communicateJumpContributions();
 
        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
        error = grid_.comm().sum(error);
 
        grid_.comm().max(errorMax, numContributions);
        grid_.comm().min(errorMin, numContributions);
 
        if (Dune::Fem::Parameter::verbose()) {
          // print error estimator if on node 0
          std::cout << "Estimated " << (Norm == L2Norm ? "L2" : "H1") << "-Error: "
                    << std::sqrt(error) << " (squared: " << error << ")" << std::endl;
          const char *errorNames[numContributions] = { "bulk", "jump", "bdry", "comm" };
          for (int i = 0; i < numContributions; ++i) {
            if (errorMax[i] == 0.0) {
              continue; // no contribution, just skip
            }
            std::cout << "  " << errorNames[i] << "^2 min: " << errorMin[i] << std::endl;
            std::cout << "  " << errorNames[i] << "^2 max: " << errorMax[i] << std::endl;
          }
          std::cout << "hMax: " << hMax << std::endl;
        }
        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();
      }
 
      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;
        }
 
        // 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;
 
        const typename ElementType::Geometry &geometry = entity.geometry();
        const RangeFieldType h2 = h2Estimate(geometry);
        const RangeFieldType hFactor = Norm == L2Norm ? h2*h2 : h2;
        const int index = indexSet_.index(entity);
        const LocalFunctionType           uLocal = uh.localFunction(entity);
        const BulkForcesLocalFunctionType fLocal = bulkForcesFunction_.localFunction(entity);
 
        ElementQuadratureType quad(entity, 2*(dfSpace_.order() + 1));
 
        operatorParts_.setEntity(entity);
 
        RangeFieldType localEstimate = 0.0;
        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);
 
          typename LocalFunctionType::RangeType values;
          uLocal.evaluate(x, values);
 
          typename LocalFunctionType::JacobianRangeType jacobian;
          uLocal.jacobian(x, jacobian);
 
          typename LocalFunctionType::HessianRangeType hessian;
          uLocal.hessian(x, hessian);
 
          // Now compute the element residual
          RangeType el_res(0.);
          if (Constituents::hasFlux) {
            RangeType tmp;
            operatorParts_.fluxDivergence(entity, quad[qp], values, jacobian, hessian, tmp);
            el_res += tmp;
          }
          if (Constituents::hasSources) {
            RangeType tmp;
            operatorParts_.source(entity, quad[qp], values, jacobian, tmp);
            el_res += tmp;
          }
          if (Constituents::hasBulkForces) {
            RangeType tmp;
            fLocal.evaluate(x, tmp);
            el_res -= tmp;
          }
 
          localEstimate += hFactor * weight * (el_res * el_res);
        }
        localIndicators_[index] += localEstimate;
 
        errorMin[bulk] = std::min(errorMin[bulk], localEstimate);
        errorMax[bulk] = std::max(errorMax[bulk], localEstimate);
 
        // calculate face contribution
        IntersectionIteratorType end = gridPart_.iend(entity);
        for (IntersectionIteratorType it = gridPart_.ibegin(entity); it != end; ++it) {
 
          const IntersectionType &intersection = *it;
 
          if (isProcessorBoundary(intersection)) {
            // interior face accross a process boundary
            estimateProcessorBoundary(intersection, entity, uLocal);
          } else if (intersection.neighbor()) {
            // interior face
            estimateIntersection(intersection, entity, uh, uLocal);
          } else if (intersection.boundary() && !dirichletBndry_.applies(intersection)) {
            // boundary face
            bool isRobin = operatorParts_.setIntersection(intersection);
 
            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 h = 0.5 * (hEstimate(inside.geometry()) + hEstimate(outside.geometry()));
            const RangeFieldType hFactor = Norm == L2Norm ? h * h * h : h;
            error *= hFactor;
            localIndicators_[insideIndex] += error;
            if (isOutsideInterior) {
              localIndicators_[outsideIndex] += error;
            }
          }
 
          errorMin[jump] = std::min(errorMin[jump], error);
          errorMax[jump] = std::max(errorMax[jump], 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();
 
        operatorPartsNeighbour_.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, valueOutside;
          JacobianRangeType jacobianInside, jacobianOutside;
          JacobianRangeType fluxInside, fluxOutside;
          uInside.evaluate(quadInside[qp], valueInside);
          uInside.jacobian(quadInside[qp], jacobianInside);
          operatorParts_.flux(inside, quadInside[qp], valueInside, jacobianInside,
                              fluxInside);
          uOutside.evaluate(quadOutside[qp], valueOutside);
          uOutside.jacobian(quadOutside[qp], jacobianOutside);
          operatorPartsNeighbour_.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 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
          // Storing the h-estimate in the last component of commData is a hack, PLEASE FIX ME ;)
          // Note that after communication commData[#QP][0] already
          // contains the sum of the h-estimate from both side.
          RangeFieldType h = 0.5 * 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 == localJumpStorage_.end());
 
        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 estimateProcessorBoundary(const IntersectionType &intersection,
                                     const ElementType &inside,
                                     const LocalFunctionType &uInside)
      {
        if (intersection.conforming()) {
          estimateProcessorBoundary<true>(intersection, inside, uInside);
        } else {
          estimateProcessorBoundary<false>(intersection, inside, uInside);
        }
      }
 
      template<bool conforming, class LocalFunctionType>
      void estimateProcessorBoundary(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);
          operatorParts_.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;
          operatorParts_.flux(inside, bndryQuad[qp], value, jacobian, flux);
 
          RangeType residual;
          flux.mv(integrationNormal, residual);
 
          if (isRobin) {
            RangeType alphau;
            operatorParts_.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;
          error *= hFactor;
          localIndicators_[insideIndex] += error;
        }
 
        errorMin[bdry] = std::min(errorMin[bdry], error);
        errorMax[bdry] = std::max(errorMax[bdry], error);
      }
    };
 
    template<class DiscreteFunctionSpace, class Model, enum FunctionSpaceNorm Norm>
    struct EstimatorTraits<EllipticEstimator<DiscreteFunctionSpace, Model, Norm> >
    {
      typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
      typedef typename
      std::vector<typename DiscreteFunctionSpaceType::RangeFieldType>::const_iterator
      IteratorType;
    };
 
 
 
  } // ACFem:.
 
} // Dune::
 
#endif // __ELLIPTIC_ESTIMATOR_HH__