parabolicestimator.hh

#ifndef __PARABOLIC_ESTIMATOR_HH__
#define __PARABOLIC_ESTIMATOR_HH__
 
#include <numeric>
 
//- 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 "../tensors/tensor.hh"
#include "../functions/functions.hh"
#include "../common/entitystorage.hh"
#include "../common/geometry.hh"
#include "../common/intersectionclassification.hh"
#include "../common/intersectiondatahandle.hh"
#include "../common/functionspacenorm.hh"
#include "../models/expressions.hh"
#include "../models/modelfacade.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class OldSolutionFunction, class TimeProvider, class ImplicitModel, class ExplicitModel,
             enum FunctionSpaceNorm Norm = L2Norm>
    class ParabolicEulerEstimator
      : public EntityStorage<typename OldSolutionFunction::DiscreteFunctionSpaceType::GridPartType,
                             typename OldSolutionFunction::DiscreteFunctionSpaceType::GridPartType::ctype>
    {
      using ThisType = ParabolicEulerEstimator;
      using BaseType = EntityStorage<typename OldSolutionFunction::DiscreteFunctionSpaceType::GridPartType,
                                     typename OldSolutionFunction::DiscreteFunctionSpaceType::GridPartType::ctype>;
     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:
      using TimeProviderType = TimeProvider;
 
      using ImplicitModelType = EffectiveModel<ImplicitModel>;
      using ImplicitTraits = EffectiveModelTraits<ImplicitModel>;
      using ImplicitModelClosure = ModelFacade<ImplicitModel>;
      using ExplicitModelType = EffectiveModel<ExplicitModel>;
      using ExplicitTraits = EffectiveModelTraits<ExplicitModel>;
      using ExplicitModelClosure = ModelFacade<ExplicitModel>;
 
      typedef typename ImplicitModelType::BoundaryConditionsType BoundaryConditionsType;
 
      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;
 
      //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_;
      ImplicitModelClosure implicitModel_;
      ImplicitModelClosure implicitModelNeighbour_;
      ExplicitModelClosure explicitModel_;
      const OldSolutionFunctionType& oldSolution_;
 
      const DiscreteFunctionSpaceType &dfSpace_;
      GridPartType &gridPart_;
      const IndexSetType &indexSet_;
      GridType &grid_;
 
     private:
      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)
        : BaseType(oldSolution.space().gridPart()),
          timeProvider_(timeProvider),
          implicitModel_(implicitModel),
          implicitModelNeighbour_(implicitModel_),
          explicitModel_(explicitModel),
          oldSolution_(oldSolution),
          dfSpace_(oldSolution_.space()),
          gridPart_(dfSpace_.gridPart()),
        indexSet_(gridPart_.indexSet()),
        grid_(gridPart_.grid()),
        timeIndicator_(0)
      {
      }
 
      template<class DiscreteFunctionType>
      RangeFieldType estimate(const DiscreteFunctionType& uh)
      {
        // 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
        error = std::accumulate(this->storage_.begin(), this->storage_.end(), error);
 
        // 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);
      }
 
      RangeFieldType timeEstimate() const {
        return timeIndicator_;
      }
 
      const DiscreteFunctionSpaceType& space() const
      {
        return dfSpace_;
      }
 
     protected:
      void clear()
      {
        // resize and clear
        BaseType::clear();
        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);
 
        implicitModel_.bind(entity);
        explicitModel_.bind(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 (ImplicitTraits::template Exists<PDEModel::flux>::value) {
            el_res += implicitModel_.fluxDivergence(quad[qp], values, jacobian, hessian);
          }
          if (ImplicitTraits::template Exists<PDEModel::source>::value) {
            el_res += implicitModel_.source(quad[qp], values, jacobian);
          }
 
          // 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 (ExplicitTraits::template Exists<PDEModel::flux>::value) {
            el_res += explicitModel_.fluxDivergence(quad[qp], values, jacobian, hessian);
          }
          if (ExplicitTraits::template Exists<PDEModel::source>::value) {
            el_res += explicitModel_.source(quad[qp], values, jacobian);
          }
 
          this->storage_[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.
        const auto end = gridPart_.iend(entity);
        for (auto it = gridPart_.ibegin(entity); it != end; ++it) {
          const auto& intersection = *it;
 
          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()) {
            const auto active = implicitModel_.classifyBoundary(intersection);
 
            assert(!active.first || active.second.all() || active.second.none());
 
            if (!active.second.all()) {
              estimateBoundary(intersection, entity, uLocal, active);
            }
          }
        }
      }
 
      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;
            this->storage_[insideIndex] += hFactor * error;
            if (isOutsideInterior) {
              this->storage_[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)
      {
        // 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();
 
        implicitModelNeighbour_.bind(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;
          uInside.evaluate(quadInside[qp], valueInside);
          uInside.jacobian(quadInside[qp], jacobianInside);
          decltype(auto) fluxInside = implicitModel_.flux(quadInside[qp], valueInside, jacobianInside);
 
          RangeType valueOutside;
          JacobianRangeType jacobianOutside;
          uOutside.evaluate(quadOutside[qp], valueOutside);
          uOutside.jacobian(quadOutside[qp], jacobianOutside);
          decltype(auto) fluxOutside = implicitModelNeighbour_.flux(quadOutside[qp], valueOutside, jacobianOutside);
 
          // Multiply with normal
          auto jump = contractInner(fluxInside - fluxOutside, integrationNormal);
 
          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", 0) != -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;
          this->storage_[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) {
          auto normal = intersection.integrationOuterNormal(quadInside.localPoint(qp));
          const auto integrationElement = normal.two_norm();
          normal /= 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;
          uInside.evaluate(quadInside[qp], valueInside);
          uInside.jacobian(quadInside[qp], jacobianInside);
          auto fluxInside = implicitModel_.flux(quadInside[qp], valueInside, jacobianInside);
 
          // Multiply with normal and store in the communication array
          commData[qp] = contractInner(fluxInside, normal);
 
          // 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,
                            const BoundaryConditionsType& active)
      {
        if (!ImplicitTraits::template Exists<PDEModel::flux>::value
            &&
            !ImplicitTraits::template Exists<PDEModel::robinFlux>::value) {
          return;
        }
 
        assert(intersection.conforming() == true);
 
        typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, true /* conforming */> IntersectionQuadratureType;
        typedef typename IntersectionQuadratureType::FaceQuadratureType Quadrature ;
 
        // TODO: make more realistic assumptions on the
        // quad-order. HOWTO?
        const auto quadOrder = 2*dfSpace_.order();
 
        // 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 normal = intersection.integrationOuterNormal(bndryQuad.localPoint(qp));
          const auto integrationElement = normal.two_norm();
          normal /= integrationElement;
 
          RangeType value;
          uLocal.evaluate(bndryQuad[qp], value);
 
          JacobianRangeType jacobian;
          if (ImplicitTraits::template Exists<PDEModel::flux>::value
              ||
              ImplicitTraits::template Exists<PDEModel::robinFlux>::value) {
            uLocal.jacobian(bndryQuad[qp], jacobian);
          }
 
          RangeType residual(0);
          if (ImplicitTraits::template Exists<PDEModel::flux>::value) {
            residual += contractInner(implicitModel_.flux(bndryQuad[qp], value, jacobian),
                              normal);
          }
 
          if (ImplicitTraits::template Exists<PDEModel::robinFlux>::value) {
            if (active.first) { // TODO: component-wise
              residual += implicitModel_.robinFlux(bndryQuad[qp], normal, value, jacobian);
            }
          }
 
          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;
          this->storage_[insideIndex] += hFactor * error;
        }
      }
    };
 
 
 
  } // ACFem::
 
} // Dune::
 
#endif // __PARABOLIC_ESTIMATOR_HH__