/** \file BackgroundMesh.h */

#ifndef BACKGROUND_MESH_H
#define BACKGROUND_MESH_H

// AMDiS-includes
#include "DOFVector.h"		// DOFVector
#include "FiniteElemSpace.h"	// FiniteElemSpace
#include "FixVec.h"		// WorldVector
#include "Global.h"		// DegreeOfFreedom

#ifndef USE_EXPERIMENTAL
#define USE_EXPERIMENTAL
#endif

namespace experimental {

  using namespace AMDiS;
  
  typedef WorldVector<double> PointType;
  
  inline double distance(const PointType &x, const PointType &y, unsigned int n)
  {
    double d = 0.0;
    for (unsigned int i = 0; i < n; i++)
      d += sqr(x[i]-y[i]);
    return sqrt(d);
  }
  
  inline double distance2(const PointType &x, const PointType &y, unsigned int n)
  {
    double d = 0.0;
    for (unsigned int i = 0; i < n; i++)
      d += sqr(x[i]-y[i]);
    return d;
  }
  
  typedef std::pair<DegreeOfFreedom, PointType> DataType;

  struct Box {

    Box(int DOW_, std::vector<int> N_);
    Box(int DOW_, PointType min_corner_, PointType max_corner_, std::vector<int> N_);
    void init();

    void fillBox(const FiniteElemSpace* feSpace);
    void clearBox();
    void clearBoxData();

    int getMaxBoxSize();

    bool inBox(const PointType& x);
    int getBox(const PointType& x);

    bool getNearestData(const PointType& x, DataType& data);
    bool getNearestData(const PointType& x, std::vector<DataType>& data, int nData);

    /**
    * strategies:
    * 0 .. get nearest point to x and eval DOFVector at this DOF-index
    * 1 .. get n nearest points, calc weighted sum of data at these points
    * 2 .. get n nearest points, calc regression plane and eval at x
    * 3 .. get n nearest points, calc weighted regression plane and eval at x
    **/
    template<typename T>
    T evalAtPoint(const DOFVector<T>& vec, const PointType& x, int strategy = 0, int nrOfPoints = -1);

    static Box* provideBackgroundMesh(const FiniteElemSpace* feSpace);
    static void delete_all();

    void addData(PointType x, DataType data);

  protected:

    template<typename T>
    T evalAtPoint_simple(const DOFVector<T>& vec, const PointType& x);

    template<typename T>
    T evalAtPoint_weighted_sum(const DOFVector<T>& vec, const PointType& x, int nrOfPoints);

    template<typename T>
    T evalAtPoint_regression(const DOFVector<T>& vec, const PointType& x, int nrOfPoints);

    template<typename T>
    T evalAtPoint_weighted_regression(const DOFVector<T>& vec, const PointType& x, int nrOfPoints);

    inline int idx2nr(std::vector<int>& idx);

    inline void nr2idx(int nr, std::vector<int>& idx);


    void getSurroundingBoxes(int nr, std::set<int> &surrounding_nrs);
    void getSurroundingBoxes(int center, int radius, std::set<int> &nrs);
    void getSurroundingBoxes(std::set<int> &nrs, std::set<int> &surrounding_nrs);

    /// berechne die minimale Entfernung der Ränder der Box zu x
    double getBoxBoundaryDist(int nr, const PointType& x);

    static std::map<const FiniteElemSpace*, std::pair<int, Box*> > boxMap;

  private:

    int DOW;
    std::vector<int> N;
    bool boxFilled;

    PointType min_corner;
    PointType max_corner;

    std::vector<std::vector<DataType> > boxData;
  };

} // end namespace experimental

#include "BackgroundMesh.hh"

#endif // BACKGROUND_MESH_H