<commit_before>#ifndef DUNE_MSFEM_ERRORESTIMATOR_HH
#define DUNE_MSFEM_ERRORESTIMATOR_HH

#include <dune/common/unused.hh>
#include <memory>
#include <array>
#include <boost/range/adaptor/map.hpp>
#include <dune/stuff/common/fixed_map.hh>
#include <dune/stuff/common/memory.hh>
#include <dune/stuff/aliases.hh>
#include <dune/multiscale/tools/misc.hh>

// where the quadratures are defined
#include <dune/fem/quadrature/cachingquadrature.hh>

#include <dune/multiscale/tools/errorestimation/MsFEM/conservative_flux_solver.hh>
#include "estimator_utils.hh"

namespace Dune {
template< class DiscreteFunctionImp,
          class DiffusionImp,
          class SourceImp,
          class MacroMicroGridSpecifierImp,
          class SubGridListImp >
class MsFEMErrorEstimator
{
  typedef MsFEMErrorEstimator< DiscreteFunctionImp,
              DiffusionImp,
              SourceImp,
              MacroMicroGridSpecifierImp,
              SubGridListImp >
    ThisType;
  typedef DiffusionImp               DiffusionOperatorType;
  typedef SourceImp                  SourceType;
  typedef MacroMicroGridSpecifierImp MacroMicroGridSpecifierType;
  typedef SubGridListImp             SubGridListType;

  //! Necessary typedefs for the DiscreteFunctionImp:

  typedef DiscreteFunctionImp                              DiscreteFunctionType;
  typedef typename DiscreteFunctionType::FunctionSpaceType FunctionSpaceType;

  typedef typename DiscreteFunctionType::LocalFunctionType         LocalFunctionType;
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;

  typedef typename DiscreteFunctionSpaceType::RangeFieldType    RangeFieldType;
  typedef typename DiscreteFunctionSpaceType::DomainType        DomainType;
  typedef typename DiscreteFunctionSpaceType::RangeType         RangeType;
  typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;

  typedef typename DiscreteFunctionType::DofIteratorType           DofIteratorType;
  typedef typename DiscreteFunctionSpaceType::GridPartType         GridPartType;
  typedef typename DiscreteFunctionSpaceType::GridType             GridType;
  typedef typename DiscreteFunctionSpaceType::IteratorType         IteratorType;
  typedef typename DiscreteFunctionSpaceType::LagrangePointSetType LagrangePointSetType;

  typedef typename GridType::Traits::LeafIndexSet LeafIndexSetType;

  typedef typename GridPartType::IntersectionIteratorType         IntersectionIteratorType;
  typedef typename IntersectionIteratorType::Intersection         Intersection;
  typedef typename GridType::template Codim< 0 >::Entity          EntityType;
  typedef typename GridType::template Codim< 0 >::EntityPointer   EntityPointerType;
  typedef typename GridType::template Codim< 1 >::Entity          FaceType;
  typedef typename GridType::template Codim< 1 >::EntityPointer   FacePointerType;
  typedef typename GridType::template Codim< 0 >::Geometry        EntityGeometryType;
  typedef typename GridType::template Codim< 1 >::Geometry        FaceGeometryType;
  typedef typename DiscreteFunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;

  typedef CachingQuadrature< GridPartType, 0 > EntityQuadratureType;
  typedef CachingQuadrature< GridPartType, 1 > FaceQuadratureType;

  // --------------------------- subgrid typedefs ------------------------------------

  typedef SubGrid< GridType::dimension, GridType > SubGridType;
  typedef LeafGridPart< SubGridType >              SubGridPart;

  typedef typename SubGridType::Traits::LeafIndexSet SubGridLeafIndexSet;

  typedef LagrangeDiscreteFunctionSpace< FunctionSpaceType, SubGridPart,
                                         1 /*POLORDER*/ > SubGridDiscreteFunctionSpaceType;

  typedef AdaptiveDiscreteFunction< SubGridDiscreteFunctionSpaceType > SubGridDiscreteFunctionType;

  typedef typename SubGridDiscreteFunctionSpaceType::IteratorType SubGridIteratorType;

  typedef typename SubGridIteratorType::Entity SubGridEntityType;

  typedef typename SubGridEntityType::EntityPointer SubGridEntityPointerType;

  typedef typename SubGridDiscreteFunctionType::LocalFunctionType SubGridLocalFunctionType;

  typedef typename SubGridDiscreteFunctionSpaceType::LagrangePointSetType
  SubGridLagrangePointSetType;

  enum { faceCodim = 1 };
  typedef typename SubGridLagrangePointSetType::template Codim< faceCodim >
    ::SubEntityIteratorType
  SubGridFaceDofIteratorType;

  typedef typename LagrangePointSetType::template Codim< faceCodim >
    ::SubEntityIteratorType
  FaceDofIteratorType;

  typedef typename SubGridType::template Codim< 0 >::Geometry SubGridEntityGeometryType;

  typedef CachingQuadrature< SubGridPart, 0 > LocalGridEntityQuadratureType;

  typedef std::array<RangeType, 3> JumpArray;
  typedef std::array<const Intersection*, 3> IntersectionArray;
  typedef std::unique_ptr< DiscreteFunctionType > DiscreteFunctionPointer;
  typedef std::array< DiscreteFunctionPointer, 2 > DiscreteFunctionPointerPair;
  typedef Stuff::Common::FixedMap<std::string, RangeType, 6> ErrormapType;

  typedef EstimatorUtils<ThisType> EstimatorUtilsType;
  template <class T>
  friend struct EstimatorUtils;

  //!-----------------------------------------------------------------------------------------

  enum { dimension = GridType::dimension };
  enum { spacePolOrd = DiscreteFunctionSpaceType::polynomialOrder };
  enum { maxnumOfBaseFct = 100 };

  const DiscreteFunctionSpaceType& fineDiscreteFunctionSpace_;
  MacroMicroGridSpecifierType& specifier_;
  SubGridListType& subgrid_list_;
  const DiffusionOperatorType& diffusion_;
  const SourceType& f_;
  const std::string& path_;

public:
  MsFEMErrorEstimator(const DiscreteFunctionSpaceType& fineDiscreteFunctionSpace,
                      MacroMicroGridSpecifierType& specifier,
                      SubGridListType& subgrid_list,
                      const DiffusionOperatorType& diffusion,
                      const SourceType& f,
                      const std::string& path)
    : fineDiscreteFunctionSpace_(fineDiscreteFunctionSpace)
      , specifier_(specifier)
      , subgrid_list_(subgrid_list)
      , diffusion_(diffusion)
      , f_(f)
      , path_(path)
  {}

private:

  //! method to get the local mesh size H of a coarse grid entity 'T'
  // works only for our 2D examples!!!!
  RangeType get_coarse_grid_H(const EntityType& entity) const {
    // entity_H means H (the diameter of the entity)
    RangeType entity_H = 0.0;

    const GridPartType& coarseGridPart = specifier_.coarseSpace().gridPart();

    // compute the size of the faces of the entities and selected the largest.
    IntersectionIteratorType endnit = coarseGridPart.iend(entity);

    for (IntersectionIteratorType nit = coarseGridPart.ibegin(entity); nit != endnit; ++nit)
    {
      FaceQuadratureType innerFaceQuadrature(coarseGridPart, *nit, 0, FaceQuadratureType::INSIDE);

      DomainType scaledOuterNormal
        = nit->integrationOuterNormal( innerFaceQuadrature.localPoint(0) );

      // get 'volume' of the visited face (this only works because we do not have curved faces):
      RangeType visitedFaceVolume(0.0);
      for (int k = 0; k < dimension; ++k)
        visitedFaceVolume += scaledOuterNormal[k] * scaledOuterNormal[k];
      visitedFaceVolume = sqrt(visitedFaceVolume);

      if (visitedFaceVolume > entity_H)
        entity_H = visitedFaceVolume;
    }

    return entity_H;
  } // get_coarse_grid_H

  // for a coarse grid entity T:
  // return:  H_T ||f||_{L^2(T)}
  RangeType indicator_f(const EntityType& entity) const {
    // create quadrature for given geometry type
    const CachingQuadrature< GridPartType, 0 > entityQuadrature(entity, 2 * spacePolOrd + 2);

    // get geoemetry of entity
    const EntityGeometryType& geometry = entity.geometry();

    const RangeType H_T = get_coarse_grid_H(entity);

    RangeType y(0);
    RangeType local_indicator(0);

    const int quadratureNop = entityQuadrature.nop();
    for (int quadraturePoint = 0; quadraturePoint < quadratureNop; ++quadraturePoint)
    {
      const double weight = entityQuadrature.weight(quadraturePoint)
                            * geometry.integrationElement( entityQuadrature.point(quadraturePoint) );

      f_.evaluate(geometry.global( entityQuadrature.point(quadraturePoint) ), y);
      y = y * y;

      local_indicator += weight * y;
    }

    local_indicator *= std::pow(H_T, 2.0);

    return sqrt(local_indicator);
  } // indicator_f


  // jump in conservative flux
  void getFluxes(const EntityType& coarse_entity,
                 const DiscreteFunctionType& msfem_coarse_part,
                 RangeType& jump_conservative_flux,
                 RangeType& jump_coarse_flux) const {
    const DiscreteFunctionSpaceType& coarseDiscreteFunctionSpace = specifier_.coarseSpace();
    const LeafIndexSetType& coarseGridLeafIndexSet = coarseDiscreteFunctionSpace.gridPart().grid().leafIndexSet();
    const int index_coarse_entity = coarseGridLeafIndexSet.index(coarse_entity);

    //! ---- get sub grid that for coarse entity
    // the sub grid U(T) that belongs to the coarse_grid_entity T
    SubGridType& sub_grid_U_T = subgrid_list_.getSubGrid(index_coarse_entity);
    SubGridPart subGridPart(sub_grid_U_T);

    SubGridDiscreteFunctionSpaceType localDiscreteFunctionSpace(subGridPart);
    std::array<SubGridDiscreteFunctionType, 2> conservative_flux_coarse_ent = {{
        SubGridDiscreteFunctionType ("Conservative Flux on coarse entity for e_0",
                                                                    localDiscreteFunctionSpace),
        SubGridDiscreteFunctionType ("Conservative Flux on coarse entity for e_1",
                                                                    localDiscreteFunctionSpace) }};
    // --------- load local solutions -------
    boost::format flux_location("%s/cf_problems/_conservativeFlux_e_%d_sg_%d");
    for( int i : {0,1}) {
      conservative_flux_coarse_ent[i].clear();
      const std::string cf_solution_location = (flux_location % path_ % i % index_coarse_entity).str();
      DiscreteFunctionReader(cf_solution_location).read(0, conservative_flux_coarse_ent[i]);
    }


    DiscreteFunctionPointerPair cflux_coarse_ent_host = {{
    DiscreteFunctionPointer( new DiscreteFunctionType ("Conservative Flux on coarse entity for e_0",
                                                  fineDiscreteFunctionSpace_)),
    DiscreteFunctionPointer( new DiscreteFunctionType ("Conservative Flux on coarse entity for e_1",
                                                  fineDiscreteFunctionSpace_)) }};

    EstimatorUtilsType::subgrid_to_hostrid_function(conservative_flux_coarse_ent,
                                cflux_coarse_ent_host);

    // flux for each neighbor entity
    std::array< DiscreteFunctionPointerPair, 3 > cflux_neighbor_ent_host;

    //!TODO save Intersection(s) instead
    IntersectionArray coarse_face;
    std::array<RangeType, 3> coarse_face_volume;

    const GridPartType& coarseGridPart = specifier_.coarseSpace().gridPart();

    int local_face_index = 0;

    IntersectionIteratorType endnit = coarseGridPart.iend(coarse_entity);
    for (IntersectionIteratorType face_it = coarseGridPart.ibegin(coarse_entity); face_it != endnit; ++face_it)
    {
      coarse_face[local_face_index] = face_it.operator->();
      coarse_face_volume[local_face_index] = face_it->geometry().volume();

      if ( face_it->neighbor() )
      {
        const EntityPointerType outside_it = face_it->outside();

        const int index_coarse_neighbor_entity = coarseGridLeafIndexSet.index(*outside_it);

        // --- get subgrids and load fluxes ---

        SubGridType& sub_grid_neighbor_U_T = subgrid_list_.getSubGrid(index_coarse_neighbor_entity);
        SubGridPart subGridPart_neighbor(sub_grid_neighbor_U_T);
        SubGridDiscreteFunctionSpaceType localDiscreteFunctionSpace_neighbor(subGridPart_neighbor);

        std::array<SubGridDiscreteFunctionType, 2> conservative_flux_coarse_ent_neighbor = {{
          SubGridDiscreteFunctionType("Conservative Flux on neighbor coarse entity for e_0",
                    localDiscreteFunctionSpace_neighbor),
          SubGridDiscreteFunctionType("Conservative Flux on neighbor coarse entity for e_1",
                    localDiscreteFunctionSpace_neighbor) }};

        // --------- load local solutions -------
        // the file/place, where we saved the solutions conservative flux problems problems
        for( int i : {0,1}) {
          conservative_flux_coarse_ent_neighbor[i].clear();
          const std::string cf_solution_location_neighbor = (flux_location
                                            % path_ % i % index_coarse_neighbor_entity).str();
          // reader for data file:
          DiscreteFunctionReader(cf_solution_location_neighbor).read(0, conservative_flux_coarse_ent_neighbor[i]);
          cflux_neighbor_ent_host[local_face_index][i] = DSC::make_unique<DiscreteFunctionType>(
            "Conservative Flux on neighbor coarse entity for e_" + Stuff::Common::toString(i),
            fineDiscreteFunctionSpace_);
        }
        EstimatorUtilsType::subgrid_to_hostrid_function(conservative_flux_coarse_ent_neighbor,
                                    cflux_neighbor_ent_host[local_face_index]);
      }

      local_face_index += 1;
    }

    if (local_face_index != 3)
    {
      DUNE_THROW(Dune::InvalidStateException,"Error! Implementation only for triangular mesh in 2d!");
    }

    const auto contributions = EstimatorUtilsType::flux_contributions(localDiscreteFunctionSpace,
                                                  sub_grid_U_T,
                                                  coarseGridLeafIndexSet,
                                                  cflux_coarse_ent_host,
                                                  msfem_coarse_part,
                                                  coarse_face,
                                                  cflux_neighbor_ent_host,
                                                  index_coarse_entity,
                                                  coarse_face_volume,
                                                  specifier_.getLevelDifference(),
                                                  fineDiscreteFunctionSpace_);
    const auto jump = contributions.first;
    const auto coarse_jump = contributions.second;

    jump_conservative_flux = ( sqrt(jump[0]) + sqrt(jump[1]) + sqrt(jump[2]) );
    jump_coarse_flux = sqrt(coarse_jump[0]) + sqrt(coarse_jump[1]) + sqrt(coarse_jump[2]);
  } // getFluxes

  void fine_contribution(const LeafIndexSetType& coarseGridLeafIndexSet,
                         const DiscreteFunctionType& msfem_solution,
                         std::vector< RangeType >& loc_approximation_error,
                         std::vector< RangeType >& loc_fine_grid_jumps,
                         ErrormapType& errors) const
  {
    // fine-grid iterator:
    const IteratorType fine_grid_end = fineDiscreteFunctionSpace_.end();
    for (IteratorType fine_grid_it = fineDiscreteFunctionSpace_.begin(); fine_grid_it != fine_grid_end; ++fine_grid_it)
    {
      EntityType& entity = *fine_grid_it;

      // identify coarse grid father entity
      EntityPointerType coarse_father = Stuff::Grid::make_father(coarseGridLeafIndexSet,
                                                                 EntityPointerType(*fine_grid_it),
                                                                 specifier_.getLevelDifference());
      const int coarse_father_index = coarseGridLeafIndexSet.index(*coarse_father);

      const EntityGeometryType& entityGeometry = entity.geometry();

      const EntityQuadratureType entityQuadrature(entity, 0);    // 0 = polynomial order
      const DomainType& x = entityGeometry.global( entityQuadrature.point(0) );

      const LocalFunctionType local_msfem_sol = msfem_solution.localFunction(entity);
      JacobianRangeType gradient_msfem_sol(0.);
      local_msfem_sol.jacobian(entityQuadrature[0], gradient_msfem_sol);

      JacobianRangeType diffusive_flux_x;
      diffusion_.diffusiveFlux(x, gradient_msfem_sol, diffusive_flux_x);

      const EntityQuadratureType highOrder_entityQuadrature(entity, 2 * spacePolOrd + 2);

      const int quadratureNop = highOrder_entityQuadrature.nop();
      for (int quadraturePoint = 0; quadraturePoint < quadratureNop; ++quadraturePoint)
      {
        const double weight = highOrder_entityQuadrature.weight(quadraturePoint)
                              * entityGeometry.integrationElement( highOrder_entityQuadrature.point(quadraturePoint) );

        DomainType point = entityGeometry.global( highOrder_entityQuadrature.point(quadraturePoint) );

        JacobianRangeType diffusive_flux_high_order;
        diffusion_.diffusiveFlux(point, gradient_msfem_sol, diffusive_flux_high_order);

        RangeType value = 0.0;
        for (int i = 0; i < dimension; ++i)
          value += pow(diffusive_flux_x[0][i] - diffusive_flux_high_order[0][i], 2.0);

        loc_approximation_error[coarse_father_index] += weight * value;
        errors["total_approximation_error"] += weight * value;
      }

      const GridPartType& fineGridPart = fineDiscreteFunctionSpace_.gridPart();

      IntersectionIteratorType endnit = fineGridPart.iend(*fine_grid_it);
      for (IntersectionIteratorType nit = fineGridPart.ibegin(*fine_grid_it); nit != endnit; ++nit)
      {
        FaceQuadratureType innerFaceQuadrature(fineGridPart, *nit, 0, FaceQuadratureType::INSIDE);
        const FaceGeometryType& faceGeometry = nit->geometry();

        if (!nit->neighbor())
        { continue; }

        EntityPointerType outer_fine_grid_it = nit->outside();
        EntityType& outer_entity = *outer_fine_grid_it;

        EntityQuadratureType outer_entityQuadrature(outer_entity, 0);    // 0 = polynomial order
        const EntityGeometryType& outer_entityGeometry = outer_entity.geometry();
        const DomainType& outer_x = outer_entityGeometry.global( outer_entityQuadrature.point(0) );

        LocalFunctionType outer_local_msfem_sol = msfem_solution.localFunction(outer_entity);
        JacobianRangeType outer_gradient_msfem_sol(0.);
        outer_local_msfem_sol.jacobian(outer_entityQuadrature[0], outer_gradient_msfem_sol);

        JacobianRangeType diffusive_flux_outside;
        diffusion_.diffusiveFlux(outer_x, outer_gradient_msfem_sol, diffusive_flux_outside);

        DomainType unitOuterNormal
          = nit->unitOuterNormal( innerFaceQuadrature.localPoint(0) );

        const RangeType edge_length = faceGeometry.volume();

        RangeType int_value = 0.0;
        for (int i = 0; i < dimension; ++i)
          int_value[i] += (diffusive_flux_x[0][i] - diffusive_flux_outside[0][i]) * unitOuterNormal[i];
        int_value = pow(int_value, 2.0);

        loc_fine_grid_jumps[coarse_father_index] += edge_length * edge_length * int_value;
        errors["total_fine_grid_jumps"] += edge_length * edge_length * int_value;
      }
    }

  }

  void coarse_contribution(const DiscreteFunctionSpaceType& coarseDiscreteFunctionSpace,
                           const DiscreteFunctionType& msfem_coarse_part,
                           const DiscreteFunctionType& msfem_fine_part,<commit_msg>fix scalar vs vector access<commit_after>#ifndef DUNE_MSFEM_ERRORESTIMATOR_HH
#define DUNE_MSFEM_ERRORESTIMATOR_HH

#include <dune/common/unused.hh>
#include <memory>
#include <array>
#include <boost/range/adaptor/map.hpp>
#include <dune/stuff/common/fixed_map.hh>
#include <dune/stuff/common/memory.hh>
#include <dune/stuff/aliases.hh>
#include <dune/multiscale/tools/misc.hh>

// where the quadratures are defined
#include <dune/fem/quadrature/cachingquadrature.hh>

#include <dune/multiscale/tools/errorestimation/MsFEM/conservative_flux_solver.hh>
#include "estimator_utils.hh"

namespace Dune {
template< class DiscreteFunctionImp,
          class DiffusionImp,
          class SourceImp,
          class MacroMicroGridSpecifierImp,
          class SubGridListImp >
class MsFEMErrorEstimator
{
  typedef MsFEMErrorEstimator< DiscreteFunctionImp,
              DiffusionImp,
              SourceImp,
              MacroMicroGridSpecifierImp,
              SubGridListImp >
    ThisType;
  typedef DiffusionImp               DiffusionOperatorType;
  typedef SourceImp                  SourceType;
  typedef MacroMicroGridSpecifierImp MacroMicroGridSpecifierType;
  typedef SubGridListImp             SubGridListType;

  //! Necessary typedefs for the DiscreteFunctionImp:

  typedef DiscreteFunctionImp                              DiscreteFunctionType;
  typedef typename DiscreteFunctionType::FunctionSpaceType FunctionSpaceType;

  typedef typename DiscreteFunctionType::LocalFunctionType         LocalFunctionType;
  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;

  typedef typename DiscreteFunctionSpaceType::RangeFieldType    RangeFieldType;
  typedef typename DiscreteFunctionSpaceType::DomainType        DomainType;
  typedef typename DiscreteFunctionSpaceType::RangeType         RangeType;
  typedef typename DiscreteFunctionSpaceType::JacobianRangeType JacobianRangeType;

  typedef typename DiscreteFunctionType::DofIteratorType           DofIteratorType;
  typedef typename DiscreteFunctionSpaceType::GridPartType         GridPartType;
  typedef typename DiscreteFunctionSpaceType::GridType             GridType;
  typedef typename DiscreteFunctionSpaceType::IteratorType         IteratorType;
  typedef typename DiscreteFunctionSpaceType::LagrangePointSetType LagrangePointSetType;

  typedef typename GridType::Traits::LeafIndexSet LeafIndexSetType;

  typedef typename GridPartType::IntersectionIteratorType         IntersectionIteratorType;
  typedef typename IntersectionIteratorType::Intersection         Intersection;
  typedef typename GridType::template Codim< 0 >::Entity          EntityType;
  typedef typename GridType::template Codim< 0 >::EntityPointer   EntityPointerType;
  typedef typename GridType::template Codim< 1 >::Entity          FaceType;
  typedef typename GridType::template Codim< 1 >::EntityPointer   FacePointerType;
  typedef typename GridType::template Codim< 0 >::Geometry        EntityGeometryType;
  typedef typename GridType::template Codim< 1 >::Geometry        FaceGeometryType;
  typedef typename DiscreteFunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;

  typedef CachingQuadrature< GridPartType, 0 > EntityQuadratureType;
  typedef CachingQuadrature< GridPartType, 1 > FaceQuadratureType;

  // --------------------------- subgrid typedefs ------------------------------------

  typedef SubGrid< GridType::dimension, GridType > SubGridType;
  typedef LeafGridPart< SubGridType >              SubGridPart;

  typedef typename SubGridType::Traits::LeafIndexSet SubGridLeafIndexSet;

  typedef LagrangeDiscreteFunctionSpace< FunctionSpaceType, SubGridPart,
                                         1 /*POLORDER*/ > SubGridDiscreteFunctionSpaceType;

  typedef AdaptiveDiscreteFunction< SubGridDiscreteFunctionSpaceType > SubGridDiscreteFunctionType;

  typedef typename SubGridDiscreteFunctionSpaceType::IteratorType SubGridIteratorType;

  typedef typename SubGridIteratorType::Entity SubGridEntityType;

  typedef typename SubGridEntityType::EntityPointer SubGridEntityPointerType;

  typedef typename SubGridDiscreteFunctionType::LocalFunctionType SubGridLocalFunctionType;

  typedef typename SubGridDiscreteFunctionSpaceType::LagrangePointSetType
  SubGridLagrangePointSetType;

  enum { faceCodim = 1 };
  typedef typename SubGridLagrangePointSetType::template Codim< faceCodim >
    ::SubEntityIteratorType
  SubGridFaceDofIteratorType;

  typedef typename LagrangePointSetType::template Codim< faceCodim >
    ::SubEntityIteratorType
  FaceDofIteratorType;

  typedef typename SubGridType::template Codim< 0 >::Geometry SubGridEntityGeometryType;

  typedef CachingQuadrature< SubGridPart, 0 > LocalGridEntityQuadratureType;

  typedef std::array<RangeType, 3> JumpArray;
  typedef std::array<const Intersection*, 3> IntersectionArray;
  typedef std::unique_ptr< DiscreteFunctionType > DiscreteFunctionPointer;
  typedef std::array< DiscreteFunctionPointer, 2 > DiscreteFunctionPointerPair;
  typedef Stuff::Common::FixedMap<std::string, RangeType, 6> ErrormapType;

  typedef EstimatorUtils<ThisType> EstimatorUtilsType;
  template <class T>
  friend struct EstimatorUtils;

  //!-----------------------------------------------------------------------------------------

  enum { dimension = GridType::dimension };
  enum { spacePolOrd = DiscreteFunctionSpaceType::polynomialOrder };
  enum { maxnumOfBaseFct = 100 };

  const DiscreteFunctionSpaceType& fineDiscreteFunctionSpace_;
  MacroMicroGridSpecifierType& specifier_;
  SubGridListType& subgrid_list_;
  const DiffusionOperatorType& diffusion_;
  const SourceType& f_;
  const std::string& path_;

public:
  MsFEMErrorEstimator(const DiscreteFunctionSpaceType& fineDiscreteFunctionSpace,
                      MacroMicroGridSpecifierType& specifier,
                      SubGridListType& subgrid_list,
                      const DiffusionOperatorType& diffusion,
                      const SourceType& f,
                      const std::string& path)
    : fineDiscreteFunctionSpace_(fineDiscreteFunctionSpace)
      , specifier_(specifier)
      , subgrid_list_(subgrid_list)
      , diffusion_(diffusion)
      , f_(f)
      , path_(path)
  {}

private:

  //! method to get the local mesh size H of a coarse grid entity 'T'
  // works only for our 2D examples!!!!
  RangeType get_coarse_grid_H(const EntityType& entity) const {
    // entity_H means H (the diameter of the entity)
    RangeType entity_H = 0.0;

    const GridPartType& coarseGridPart = specifier_.coarseSpace().gridPart();

    // compute the size of the faces of the entities and selected the largest.
    IntersectionIteratorType endnit = coarseGridPart.iend(entity);

    for (IntersectionIteratorType nit = coarseGridPart.ibegin(entity); nit != endnit; ++nit)
    {
      FaceQuadratureType innerFaceQuadrature(coarseGridPart, *nit, 0, FaceQuadratureType::INSIDE);

      DomainType scaledOuterNormal
        = nit->integrationOuterNormal( innerFaceQuadrature.localPoint(0) );

      // get 'volume' of the visited face (this only works because we do not have curved faces):
      RangeType visitedFaceVolume(0.0);
      for (int k = 0; k < dimension; ++k)
        visitedFaceVolume += scaledOuterNormal[k] * scaledOuterNormal[k];
      visitedFaceVolume = sqrt(visitedFaceVolume);

      if (visitedFaceVolume > entity_H)
        entity_H = visitedFaceVolume;
    }

    return entity_H;
  } // get_coarse_grid_H

  // for a coarse grid entity T:
  // return:  H_T ||f||_{L^2(T)}
  RangeType indicator_f(const EntityType& entity) const {
    // create quadrature for given geometry type
    const CachingQuadrature< GridPartType, 0 > entityQuadrature(entity, 2 * spacePolOrd + 2);

    // get geoemetry of entity
    const EntityGeometryType& geometry = entity.geometry();

    const RangeType H_T = get_coarse_grid_H(entity);

    RangeType y(0);
    RangeType local_indicator(0);

    const int quadratureNop = entityQuadrature.nop();
    for (int quadraturePoint = 0; quadraturePoint < quadratureNop; ++quadraturePoint)
    {
      const double weight = entityQuadrature.weight(quadraturePoint)
                            * geometry.integrationElement( entityQuadrature.point(quadraturePoint) );

      f_.evaluate(geometry.global( entityQuadrature.point(quadraturePoint) ), y);
      y = y * y;

      local_indicator += weight * y;
    }

    local_indicator *= std::pow(H_T, 2.0);

    return sqrt(local_indicator);
  } // indicator_f


  // jump in conservative flux
  void getFluxes(const EntityType& coarse_entity,
                 const DiscreteFunctionType& msfem_coarse_part,
                 RangeType& jump_conservative_flux,
                 RangeType& jump_coarse_flux) const {
    const DiscreteFunctionSpaceType& coarseDiscreteFunctionSpace = specifier_.coarseSpace();
    const LeafIndexSetType& coarseGridLeafIndexSet = coarseDiscreteFunctionSpace.gridPart().grid().leafIndexSet();
    const int index_coarse_entity = coarseGridLeafIndexSet.index(coarse_entity);

    //! ---- get sub grid that for coarse entity
    // the sub grid U(T) that belongs to the coarse_grid_entity T
    SubGridType& sub_grid_U_T = subgrid_list_.getSubGrid(index_coarse_entity);
    SubGridPart subGridPart(sub_grid_U_T);

    SubGridDiscreteFunctionSpaceType localDiscreteFunctionSpace(subGridPart);
    std::array<SubGridDiscreteFunctionType, 2> conservative_flux_coarse_ent = {{
        SubGridDiscreteFunctionType ("Conservative Flux on coarse entity for e_0",
                                                                    localDiscreteFunctionSpace),
        SubGridDiscreteFunctionType ("Conservative Flux on coarse entity for e_1",
                                                                    localDiscreteFunctionSpace) }};
    // --------- load local solutions -------
    boost::format flux_location("%s/cf_problems/_conservativeFlux_e_%d_sg_%d");
    for( int i : {0,1}) {
      conservative_flux_coarse_ent[i].clear();
      const std::string cf_solution_location = (flux_location % path_ % i % index_coarse_entity).str();
      DiscreteFunctionReader(cf_solution_location).read(0, conservative_flux_coarse_ent[i]);
    }


    DiscreteFunctionPointerPair cflux_coarse_ent_host = {{
    DiscreteFunctionPointer( new DiscreteFunctionType ("Conservative Flux on coarse entity for e_0",
                                                  fineDiscreteFunctionSpace_)),
    DiscreteFunctionPointer( new DiscreteFunctionType ("Conservative Flux on coarse entity for e_1",
                                                  fineDiscreteFunctionSpace_)) }};

    EstimatorUtilsType::subgrid_to_hostrid_function(conservative_flux_coarse_ent,
                                cflux_coarse_ent_host);

    // flux for each neighbor entity
    std::array< DiscreteFunctionPointerPair, 3 > cflux_neighbor_ent_host;

    //!TODO save Intersection(s) instead
    IntersectionArray coarse_face;
    std::array<RangeType, 3> coarse_face_volume;

    const GridPartType& coarseGridPart = specifier_.coarseSpace().gridPart();

    int local_face_index = 0;

    IntersectionIteratorType endnit = coarseGridPart.iend(coarse_entity);
    for (IntersectionIteratorType face_it = coarseGridPart.ibegin(coarse_entity); face_it != endnit; ++face_it)
    {
      coarse_face[local_face_index] = face_it.operator->();
      coarse_face_volume[local_face_index] = face_it->geometry().volume();

      if ( face_it->neighbor() )
      {
        const EntityPointerType outside_it = face_it->outside();

        const int index_coarse_neighbor_entity = coarseGridLeafIndexSet.index(*outside_it);

        // --- get subgrids and load fluxes ---

        SubGridType& sub_grid_neighbor_U_T = subgrid_list_.getSubGrid(index_coarse_neighbor_entity);
        SubGridPart subGridPart_neighbor(sub_grid_neighbor_U_T);
        SubGridDiscreteFunctionSpaceType localDiscreteFunctionSpace_neighbor(subGridPart_neighbor);

        std::array<SubGridDiscreteFunctionType, 2> conservative_flux_coarse_ent_neighbor = {{
          SubGridDiscreteFunctionType("Conservative Flux on neighbor coarse entity for e_0",
                    localDiscreteFunctionSpace_neighbor),
          SubGridDiscreteFunctionType("Conservative Flux on neighbor coarse entity for e_1",
                    localDiscreteFunctionSpace_neighbor) }};

        // --------- load local solutions -------
        // the file/place, where we saved the solutions conservative flux problems problems
        for( int i : {0,1}) {
          conservative_flux_coarse_ent_neighbor[i].clear();
          const std::string cf_solution_location_neighbor = (flux_location
                                            % path_ % i % index_coarse_neighbor_entity).str();
          // reader for data file:
          DiscreteFunctionReader(cf_solution_location_neighbor).read(0, conservative_flux_coarse_ent_neighbor[i]);
          cflux_neighbor_ent_host[local_face_index][i] = DSC::make_unique<DiscreteFunctionType>(
            "Conservative Flux on neighbor coarse entity for e_" + Stuff::Common::toString(i),
            fineDiscreteFunctionSpace_);
        }
        EstimatorUtilsType::subgrid_to_hostrid_function(conservative_flux_coarse_ent_neighbor,
                                    cflux_neighbor_ent_host[local_face_index]);
      }

      local_face_index += 1;
    }

    if (local_face_index != 3)
    {
      DUNE_THROW(Dune::InvalidStateException,"Error! Implementation only for triangular mesh in 2d!");
    }

    const auto contributions = EstimatorUtilsType::flux_contributions(localDiscreteFunctionSpace,
                                                  sub_grid_U_T,
                                                  coarseGridLeafIndexSet,
                                                  cflux_coarse_ent_host,
                                                  msfem_coarse_part,
                                                  coarse_face,
                                                  cflux_neighbor_ent_host,
                                                  index_coarse_entity,
                                                  coarse_face_volume,
                                                  specifier_.getLevelDifference(),
                                                  fineDiscreteFunctionSpace_);
    const auto jump = contributions.first;
    const auto coarse_jump = contributions.second;

    jump_conservative_flux = ( sqrt(jump[0]) + sqrt(jump[1]) + sqrt(jump[2]) );
    jump_coarse_flux = sqrt(coarse_jump[0]) + sqrt(coarse_jump[1]) + sqrt(coarse_jump[2]);
  } // getFluxes

  void fine_contribution(const LeafIndexSetType& coarseGridLeafIndexSet,
                         const DiscreteFunctionType& msfem_solution,
                         std::vector< RangeType >& loc_approximation_error,
                         std::vector< RangeType >& loc_fine_grid_jumps,
                         ErrormapType& errors) const
  {
    // fine-grid iterator:
    const IteratorType fine_grid_end = fineDiscreteFunctionSpace_.end();
    for (IteratorType fine_grid_it = fineDiscreteFunctionSpace_.begin(); fine_grid_it != fine_grid_end; ++fine_grid_it)
    {
      EntityType& entity = *fine_grid_it;

      // identify coarse grid father entity
      EntityPointerType coarse_father = Stuff::Grid::make_father(coarseGridLeafIndexSet,
                                                                 EntityPointerType(*fine_grid_it),
                                                                 specifier_.getLevelDifference());
      const int coarse_father_index = coarseGridLeafIndexSet.index(*coarse_father);

      const EntityGeometryType& entityGeometry = entity.geometry();

      const EntityQuadratureType entityQuadrature(entity, 0);    // 0 = polynomial order
      const DomainType& x = entityGeometry.global( entityQuadrature.point(0) );

      const LocalFunctionType local_msfem_sol = msfem_solution.localFunction(entity);
      JacobianRangeType gradient_msfem_sol(0.);
      local_msfem_sol.jacobian(entityQuadrature[0], gradient_msfem_sol);

      JacobianRangeType diffusive_flux_x;
      diffusion_.diffusiveFlux(x, gradient_msfem_sol, diffusive_flux_x);

      const EntityQuadratureType highOrder_entityQuadrature(entity, 2 * spacePolOrd + 2);

      const int quadratureNop = highOrder_entityQuadrature.nop();
      for (int quadraturePoint = 0; quadraturePoint < quadratureNop; ++quadraturePoint)
      {
        const double weight = highOrder_entityQuadrature.weight(quadraturePoint)
                              * entityGeometry.integrationElement( highOrder_entityQuadrature.point(quadraturePoint) );

        DomainType point = entityGeometry.global( highOrder_entityQuadrature.point(quadraturePoint) );

        JacobianRangeType diffusive_flux_high_order;
        diffusion_.diffusiveFlux(point, gradient_msfem_sol, diffusive_flux_high_order);

        RangeType value = 0.0;
        for (int i = 0; i < dimension; ++i)
          value += pow(diffusive_flux_x[0][i] - diffusive_flux_high_order[0][i], 2.0);

        loc_approximation_error[coarse_father_index] += weight * value;
        errors["total_approximation_error"] += weight * value;
      }

      const GridPartType& fineGridPart = fineDiscreteFunctionSpace_.gridPart();

      IntersectionIteratorType endnit = fineGridPart.iend(*fine_grid_it);
      for (IntersectionIteratorType nit = fineGridPart.ibegin(*fine_grid_it); nit != endnit; ++nit)
      {
        FaceQuadratureType innerFaceQuadrature(fineGridPart, *nit, 0, FaceQuadratureType::INSIDE);
        const FaceGeometryType& faceGeometry = nit->geometry();

        if (!nit->neighbor())
        { continue; }

        EntityPointerType outer_fine_grid_it = nit->outside();
        EntityType& outer_entity = *outer_fine_grid_it;

        EntityQuadratureType outer_entityQuadrature(outer_entity, 0);    // 0 = polynomial order
        const EntityGeometryType& outer_entityGeometry = outer_entity.geometry();
        const DomainType& outer_x = outer_entityGeometry.global( outer_entityQuadrature.point(0) );

        LocalFunctionType outer_local_msfem_sol = msfem_solution.localFunction(outer_entity);
        JacobianRangeType outer_gradient_msfem_sol(0.);
        outer_local_msfem_sol.jacobian(outer_entityQuadrature[0], outer_gradient_msfem_sol);

        JacobianRangeType diffusive_flux_outside;
        diffusion_.diffusiveFlux(outer_x, outer_gradient_msfem_sol, diffusive_flux_outside);

        DomainType unitOuterNormal
          = nit->unitOuterNormal( innerFaceQuadrature.localPoint(0) );

        const RangeType edge_length = faceGeometry.volume();

        RangeType int_value = 0.0;
        for (int i = 0; i < dimension; ++i)
          int_value += (diffusive_flux_x[0][i] - diffusive_flux_outside[0][i]) * unitOuterNormal[i];
        int_value = pow(int_value, 2.0);

        loc_fine_grid_jumps[coarse_father_index] += edge_length * edge_length * int_value;
        errors["total_fine_grid_jumps"] += edge_length * edge_length * int_value;
      }
    }

  }

  void coarse_contribution(const DiscreteFunctionSpaceType& coarseDiscreteFunctionSpace,
                           const DiscreteFunctionType& msfem_coarse_part,
                           const DiscreteFunctionType& msfem_fine_part,
{'fraction_non_alphanumeric': 0.08348222377586162, 'fraction_numerical': 0.004712706180895414, 'mean_word_length': 3.08634006983388, 'pattern_counts': {'":': 0, '<': 99, '<?xml version=': 0, '>': 109, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 74, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1568276', 'n_tokens_mistral': 11044, 'n_tokens_neox': 10317, 'n_words': 2105}