ellipticestimator.hh

#ifndef __ELLIPTIC_ESTIMATOR_HH__
#define __ELLIPTIC_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 "../common/entitystorage.hh"
#include "../models/modeltraits.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"
#include "../algorithms/operatorselector.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunctionSpace, class Model, enum FunctionSpaceNorm Norm = H1Norm>
    class EllipticEstimator
      : public EntityStorage<typename DiscreteFunctionSpace::GridPartType,
                             typename DiscreteFunctionSpace::GridPartType::ctype>
    {
      using ThisType = EllipticEstimator;
      using BaseType = EntityStorage<typename DiscreteFunctionSpace::GridPartType,
                                     typename DiscreteFunctionSpace::GridPartType::ctype>;
     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:
      using ModelType = EffectiveModel<Model>;
      using ModelTraits = EffectiveModelTraits<Model>;
      using ModelClosure = ModelFacade<Model>;
 
      typedef typename ModelType::BoundaryConditionsType BoundaryConditionsType;
 
      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::IndexSetType IndexSetType;
      typedef typename IndexSetType::IndexType IndexType;
 
      //intersection type
      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;
 
      using DGSelectorType = OperatorTypeSelector<DiscreteFunctionSpaceType, ModelType>;
      static constexpr bool useDG = DGSelectorType::useDG;
      static constexpr DGFlavour dgFlavour = DGSelectorType::dgFlavour;
      static constexpr int sFormDG = DGSelectorType::sFormDG;
 
     protected:
      ModelClosure model_;              //< The underlying model.
      ModelClosure modelNeighbour_;
 
      const DiscreteFunctionSpaceType &dfSpace_;
      GridPartType &gridPart_;
      const IndexSetType &indexSet_;
      GridType &grid_;
 
     private:
      mutable std::vector<JacobianRangeType> outsideFaceFluxCache;
 
      enum { bulkPart, jumpPart, bdryPart, commPart, numContributions };
      RangeFieldType errorMin[numContributions];
      RangeFieldType errorMax[numContributions];
 
      class IntersectionStorage
        : public std::vector<RangeFieldType>
      {
        typedef std::vector<RangeFieldType> BaseType;
        enum {
          dimRange = RangeType::dimension
        };
        typedef Fem::StaticOffsetSubMapper<dimRange> MapperType;
        typedef Fem::SubVector<BaseType, MapperType> RangeBlockType;
        typedef Fem::Envelope<RangeBlockType> RangeBlockPtrType;
        typedef Fem::SubVector<const BaseType, MapperType> ConstRangeBlockType;
       public:
        using typename BaseType::vector;
        using BaseType::size;
        IntersectionStorage(size_t numQuadraturePoints)
          : BaseType(numQuadraturePoints*dimRange + 1)
        {}
 
        IntersectionStorage()
          : BaseType()
        {}
 
        size_t nop() const
        {
          assert(dimRange == 1 || size() % dimRange == 1);
          return size() / dimRange;
        }
 
        using BaseType::back;
        RangeFieldType& hEstimate() { return back(); }
        const RangeFieldType& hEstimate() const { return back(); }
 
 
        RangeBlockPtrType operator()(size_t i)
        {
          assert(i < nop());
          return RangeBlockPtrType(RangeBlockType(*this, MapperType(i*dimRange)));
        }
        ConstRangeBlockType operator()(size_t i) const
        {
          assert(i < nop());
          return ConstRangeBlockType(*this, MapperType(i*dimRange));
        }
      };
 
      typedef IntersectionStorage 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;
 
      typedef std::vector<RangeType> LocalIntersectionDGStorage;
      typedef std::map<IndexType, LocalIntersectionDGStorage> LocalDGStorageType;
 
      // Store for every intersection the element index, h and all
      // quad weights. This buffer is local for each computation node
      // and is not communicated.
      LocalStorageType localJumpStorage_;
 
      // Store for every intersection the contribution to the DG-jump
      // from the interior for each each quadrature point.  This
      // buffer is local for each computation node and is not
      // commmunicated.
      LocalDGStorageType localDGJumpStorage_;
 
      // Store for every intersection the contribution to all jump
      // terms (jump in the values for DG and the jumps of the flux)
      // at all quadrature points. During communication the
      // communicated values from the "outside" are added For the
      // jumps of the flux this results in the jump term as the
      // normals in the outside have the opposite direction. For the
      // jumps in the values (DG case) the jump finally is computed by
      // subtracting the values of localDGJumpStorage_ twice.
      CommunicationStorageType commJumpStorage_;
 
     public:
      explicit EllipticEstimator(const ModelType& model, const DiscreteFunctionSpaceType& dfSpace)
        : BaseType(dfSpace.gridPart()),
          model_(model),
          modelNeighbour_(model),
          dfSpace_(dfSpace),
          gridPart_(dfSpace_.gridPart()),
          indexSet_(gridPart_.indexSet()),
          grid_(gridPart_.grid())
      {
      }
 
      template<class DiscreteFunctionType>
      RangeFieldType estimate(const DiscreteFunctionType& uh)
      {
        // clear all local estimators
        clear();
 
        RangeFieldType hMax = 0.0;
 
        // calculate local estimator
        const auto end = dfSpace_.end();
        for (auto 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. NB: accumulate does not overwrite
        // "error", we need the result of accumulate.
        error = std::accumulate(this->storage_.begin(), this->storage_.end(), error);
 
        // 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 DiscreteFunctionSpaceType& space() const
      {
        return dfSpace_;
      }
 
     protected:
      void clear()
      {
        // resize and clear
        BaseType::clear();
 
        // the communication hack for parallel runs. Gnah.
        localJumpStorage_.clear();
        commJumpStorage_.clear();
        if (useDG) {
          localDGJumpStorage_.clear();
        }
 
        // debugging and stuff
        errorMin[bulkPart] =
          errorMin[jumpPart] =
          errorMin[bdryPart] =
          errorMin[commPart] = std::numeric_limits<RangeFieldType>::max();
 
        errorMax[bulkPart] =
          errorMax[jumpPart] =
          errorMax[bdryPart] =
          errorMax[commPart] = 0.0;
      }
 
      template<class DiscreteFunction>
      void estimateLocal(const ElementType &entity, const DiscreteFunction& uh)
      {
        const auto& geometry = entity.geometry();
        const auto h2 = h2Estimate(geometry);
        const auto hFactor = Norm == L2Norm ? h2*h2 : h2;
        const auto index = indexSet_.index(entity);
        const auto uLocal = uh.localFunction(entity);
 
        ElementQuadratureType quad(entity, 2*(dfSpace_.order() + 1));
 
        model_.bind(entity);
 
        RangeFieldType localEstimate = 0.0;
        const auto numQuadraturePoints = quad.nop();
        for (size_t qp = 0; qp < numQuadraturePoints; ++qp) {
 
          const auto& x = quad.point(qp); //
          const auto 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 (ModelTraits::template Exists<PDEModel::flux>::value) {
            el_res += model_.fluxDivergence(quad[qp], values, jacobian, hessian);
          }
          if (ModelTraits::template Exists<PDEModel::source>::value) {
            el_res += model_.source(quad[qp], values, jacobian);
          }
 
          localEstimate += hFactor * weight * (el_res * el_res);
        }
        this->storage_[index] += localEstimate;
 
        errorMin[bulkPart] = std::min(errorMin[bulkPart], localEstimate);
        errorMax[bulkPart] = std::max(errorMax[bulkPart], localEstimate);
 
        // calculate face contribution
        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
            estimateProcessorBoundary(intersection, entity, uLocal);
          } else if (intersection.neighbor()) {
            // interior face
            estimateIntersection(intersection, entity, uh, uLocal);
          } else if (intersection.boundary()) {
            const auto active = model_.classifyBoundary(intersection);
 
            assert(!active.first || active.second.all() || active.second.none());
 
            if (!active.second.all()) {
              estimateBoundary(intersection, entity, uLocal, active);
            }
          }
        }
 
        model_.unbind();
      }
 
      template<class DiscreteFunctionType>
      void estimateIntersection(const IntersectionType &intersection,
                                const ElementType &inside,
                                const DiscreteFunctionType& uh,
                                const typename DiscreteFunctionType::LocalFunctionType &uInside)
      {
        if (!ModelTraits::template Exists<PDEModel::flux>::value) {
          return; // 0 is a very smooth function ...
        }
 
        // get outside entity pointer
        const auto outside = intersection.outside();
 
        const auto insideIndex = indexSet_.index(inside);
        const auto outsideIndex = indexSet_.index(outside);
 
        const bool isOutsideInterior = (outside.partitionType() == Dune::InteriorEntity);
        if (!isOutsideInterior || (insideIndex < outsideIndex)) {
          const auto 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) {
            this->storage_[insideIndex] += error;
            if (isOutsideInterior) {
              this->storage_[outsideIndex] += error;
            }
          }
 
          errorMin[jumpPart] = std::min(errorMin[jumpPart], error);
          errorMax[jumpPart] = std::max(errorMax[jumpPart], 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);
 
        const auto h = 0.5 * (hEstimate(inside.geometry()) + hEstimate(outside.geometry()));
 
        // use IntersectionQuadrature to create appropriate face quadratures
        typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, conforming> IntersectionQuadratureType;
 
        const int quadOrder = 2 * (dfSpace_.order() - 1);
        // create intersection quadrature
        IntersectionQuadratureType intersectionQuad(gridPart_, intersection, quadOrder);
 
        // get appropriate quadrature references
        const auto& quadInside  = intersectionQuad.inside();
        const auto& quadOutside = intersectionQuad.outside();
 
        const auto numQuadraturePoints = quadInside.nop();
 
        modelNeighbour_.bind(outside); // go to neighbour
        RangeFieldType error = 0;
        RangeFieldType dgError = 0;
        for (size_t qp = 0; qp < numQuadraturePoints; ++qp) {
          const auto integrationNormal
            = intersection.integrationOuterNormal(quadInside.localPoint(qp));
          const auto integrationElement = integrationNormal.two_norm();
 
          // new stuff for dg estim
          DomainType unitNormal = integrationNormal;
          unitNormal/=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, valueOutside;
          JacobianRangeType jacobianInside, jacobianOutside;
          uInside.evaluate(quadInside[qp], valueInside);
          uInside.jacobian(quadInside[qp], jacobianInside);
          auto fluxInside = model_.flux(quadInside[qp], valueInside, jacobianInside);
          uOutside.evaluate(quadOutside[qp], valueOutside);
          uOutside.jacobian(quadOutside[qp], jacobianOutside);
          auto fluxOutside = modelNeighbour_.flux(quadOutside[qp], valueOutside, jacobianOutside);
 
          // Compute difference and multiply with normal
          auto jump = contractInner(fluxInside - fluxOutside, integrationNormal);
 
          error += quadInside.weight(qp) * (jump * jump) / integrationElement;
 
          if (useDG) {
            auto jump = valueInside - valueOutside;
            dgError += quadInside.weight( qp ) * (jump * jump) * integrationElement;
          }
 
        }
        modelNeighbour_.unbind();
 
        const auto hFactor = Norm == L2Norm ? h * h * h : h;
        error *= hFactor;
 
        if (useDG) {
          assert(Norm == H1Norm);
          dgError /= h;
          error += dgError;
        }
 
        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 localIt = localJumpStorage_.begin();
        auto localDGIt = localDGJumpStorage_.begin();
        for (const auto& remote : commJumpStorage_) {
          const auto& local = *localIt;
          assert(remote.first == local.first);
 
          const auto& localData = local.second;
          const auto& commData = remote.second;
 
          const auto numQuadraturePoints = localData.size();
 
          // form the square sum
          RangeFieldType error = 0;
          RangeFieldType dgError = 0;
          for (size_t qp = 0; qp < numQuadraturePoints; ++qp) {
            error += localData[qp] * (commData(qp) * commData(qp));
            if (useDG) {
              // comm uses "add", which is transformed into a "sub" by
              // subtracting one the values twice.
              RangeType jump = commData(numQuadraturePoints + qp);
              jump.axpy(-2, localDGIt->second[qp]);
              dgError += localData[qp] * (jump * jump);
            }
          }
 
          // 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.
          const auto h = 0.5 * commData.hEstimate();
          const auto hFactor = Norm == L2Norm ? h * h * h : h;
 
          // add contribution to the estimator array
          error *= hFactor;
 
          if (useDG) {
            assert(Norm == H1Norm);
            dgError /= h;
            error += dgError;
          }
 
          this->storage_[localData.bulkEntityIndex()] += error;
 
          errorMin[commPart] = std::min(errorMin[commPart], error);
          errorMax[commPart] = std::max(errorMax[commPart], error);
 
          ++localIt;
          if (useDG) {
            ++localDGIt;
          }
        }
 
        assert(localIt == localJumpStorage_.end());
        assert(localDGIt == localDGJumpStorage_.end());
 
        if (verbose) {
          auto 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 (!ModelTraits::template Exists<PDEModel::flux>::value) {
          return; // 0 is a very smooth function ...
        }
 
        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;
 
        const auto quadOrder = 2 * (dfSpace_.order() - 1);
        // create intersection quadrature
        IntersectionQuadratureType intersectionQuad(gridPart_, intersection, quadOrder);
 
        // get appropriate quadrature references
        const auto& quadInside  = intersectionQuad.inside();
        const auto numQuadraturePoints = quadInside.nop();
 
        // FIXME: is there any easier way to get hold of the face-index?
        const auto faceIndex = indexSet_.index(inside.template subEntity<1>(intersection.indexInInside()));
        const auto bulkIndex = indexSet_.index(inside);
 
        // commJumpStorage_ and localJumpStorage_ are maps, hence we
        // first accumulate the results in the following two local
        // variables and then add those variables to the communication
        // arrays.
        auto commSize = useDG ? 2*numQuadraturePoints : numQuadraturePoints;
        IntersectionStorageType commData(commSize);
        LocalIntersectionStorage localData(bulkIndex, numQuadraturePoints);
        LocalIntersectionDGStorage localDGData(useDG ? numQuadraturePoints : 0);
 
        for (size_t qp = 0; qp < numQuadraturePoints; ++qp) {
          auto normal = intersection.integrationOuterNormal(quadInside.localPoint(qp));
          const RangeFieldType 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);
          const auto fluxInside = model_.flux(quadInside[qp], valueInside, jacobianInside);
 
          // Multiply with normal and store in the communication array.
          *commData(qp) = contractInner(fluxInside, normal);
          //fluxInside.mv(normal, *commData(qp));
 
          if (useDG) {
            localDGData[qp] = valueInside;
            for (unsigned i = 0; i < valueInside.size(); ++i) {
              (*commData(numQuadraturePoints + qp))[i] = valueInside[i];
            }
          }
 
          // Remember weight and integration element for this point
          localData[qp] = quadInside.weight(qp) * integrationElement;
        }
        localData.hEstimate() = hEstimate(inside.geometry());
        // add for communication with the neighbour
        commData.hEstimate() = localData.hEstimate();
 
        commJumpStorage_[faceIndex] = commData;
        localJumpStorage_[faceIndex] = localData;
        if (useDG) {
          localDGJumpStorage_[faceIndex] = localDGData;
        }
      }
 
      template<class LocalFunctionType>
      void estimateBoundary(const IntersectionType &intersection,
                            const ElementType &inside,
                            const LocalFunctionType &uLocal,
                            const BoundaryConditionsType& active)
      {
        if (!ModelTraits::template Exists<PDEModel::flux>::value
            &&
            !ModelTraits::template Exists<PDEModel::robinFlux>::value) {
          return;
        }
 
        assert(intersection.conforming() == true);
 
        typedef Dune::Fem::IntersectionQuadrature<FaceQuadratureType, true /* conforming */> IntersectionQuadratureType;
 
        // 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 auto& bndryQuad = intersectionQuad.inside();
        const auto numQuadraturePoints = bndryQuad.nop();
 
        const auto insideIndex = indexSet_.index(inside);
 
        const auto h = hEstimate(inside.geometry());
 
        RangeFieldType error = 0.0;
 
        for (size_t qp = 0; qp < numQuadraturePoints; ++qp) {
          const auto weight = bndryQuad.weight(qp);
 
          // Get the outer normal and integration element
          auto normal = intersection.integrationOuterNormal(bndryQuad.localPoint(qp));
          const auto integrationElement = normal.two_norm();
          normal /= integrationElement;
 
          RangeType value;
          uLocal.evaluate(bndryQuad[qp], value);
 
          JacobianRangeType jacobian;
          if (ModelTraits::template Exists<PDEModel::flux>::value
              ||
              ModelTraits::template Exists<PDEModel::robinFlux>::value) {
            uLocal.jacobian(bndryQuad[qp], jacobian);
          }
 
          RangeType residual(0);
 
          if (ModelTraits::template Exists<PDEModel::flux>::value) {
            const auto flux = model_.flux(bndryQuad[qp], value, jacobian);
            residual += contractInner(flux, normal);
            //flux.umv(normal, residual);
          }
 
          if (ModelTraits::template Exists<PDEModel::robinFlux>::value) {
            if (active.first) { // TODO: component-wise
              residual += model_.robinFlux(bndryQuad[qp], normal, value, jacobian);
            }
          }
 
          error += weight * integrationElement * (residual * residual);
        }
 
        if (error > 0.0) {
          const auto hFactor = Norm == L2Norm ? h * h * h : h;
 
          error *= hFactor;
 
          this->storage_[insideIndex] += error;
        }
 
        errorMin[bdryPart] = std::min(errorMin[bdryPart], error);
        errorMax[bdryPart] = std::max(errorMax[bdryPart], error);
      }
    };
 
 
 
  } // ACFem:.
 
} // Dune::
 
#endif // __ELLIPTIC_ESTIMATOR_HH__