/** \file Views.h */

#ifndef MY_VIEWS_H
#define MY_VIEWS_H

#include <stdexcept>
#include <boost/type_traits/is_base_of.hpp>
#include <boost/utility/enable_if.hpp>

#include "AMDiS.h"
#ifdef USE_BACKGROUNDMESH
  #include "BackgroundMesh.h"
#endif
#ifdef USE_KD_TREE
  #include "kdtree_nanoflann_dof.h"
#endif

using namespace AMDiS; 

/// Abstract-Function-Views
/// ________________________________________________________________________________________________

/// result = 0.5*(f(x)+1)
template<typename T>
struct PhaseView : AbstractFunction<T, WorldVector<double> > {
  PhaseView(AbstractFunction<T, WorldVector<double> > &fct_) : fct(fct_) {}
  PhaseView(AbstractFunction<T, WorldVector<double> > *fct_) : fct(*fct_) {}
  T operator()(const WorldVector<double>& x) const {
	  return std::max(0.0, std::min(1.0, 0.5*(fct(x)+1.0)));
  }
private:
  AbstractFunction<T, WorldVector<double> > fct;
};

/// result = 2*f(x)-1
template<typename T>
struct ChView : AbstractFunction<T, WorldVector<double> > {
  ChView(AbstractFunction<T, WorldVector<double> > &fct_) : fct(fct_) {}
  ChView(AbstractFunction<T, WorldVector<double> > *fct_) : fct(*fct_) {}
  T operator()(const WorldVector<double>& x) const {
	  return std::max(-1.0, std::min(1.0, 2.0*fct(x)-1.0));
  }
private:
  AbstractFunction<T, WorldVector<double> > fct;
};

template<typename T>
struct MultView : AbstractFunction<T, WorldVector<double> > {
  MultView(AbstractFunction<T, WorldVector<double> > &fct_, double factor_) : fct(fct_), factor(factor_) {}
  MultView(AbstractFunction<T, WorldVector<double> > *fct_, double factor_) : fct(*fct_), factor(factor_) {}
  T operator()(const WorldVector<double>& x) const {
    return factor*fct(x);
  }
private:
  AbstractFunction<T, WorldVector<double> > fct;
  double factor;
};

/// DOFVector-Views
/// ________________________________________________________________________________________________

template<typename T>
struct DOFView : public DOFVector<T> {
  typedef T result_type;
  
  DOFView(DOFVector<T> *vec_)
  : DOFVector<T>(vec_->getFeSpace(), "dof_view"), 
    vector(vec_)
  {
    nullify(value0);
  }
    
  DOFView(DOFVector<T> &vec_)
  : DOFVector<T>(vec_.getFeSpace(), "dof_view"), 
    vector(&vec_)
  {
    nullify(value0);
  }
  
  /// Returns \ref vec[i]
  inline const T& operator[](WorldVector<double> &x) const 
  {
    FUNCNAME("DOFVector<T>::operator[x]");
    DegreeOfFreedom idx = -1;
    WorldVector<double> x_ = coordsView(x);
    bool inside = vector->getDofIdxAtPoint(x_, idx);
    
    return (inside ? operator[](idx) : value0);
  } 

  /// Returns \ref vec[i]
  inline T& operator[](WorldVector<double> &x) 
  {
    FUNCNAME("DOFVector<T>::operator[x]");
    DegreeOfFreedom idx = -1;
    WorldVector<double> x_ = coordsView(x);
    bool inside = vector->getDofIdxAtPoint(x_, idx);

    return (inside ? operator[](idx) : value0);
  }
  
  inline T& operator[](DegreeOfFreedom idx)
  {
    DOFVector<T>::operator[](idx)= valueView(idx);
    return DOFVector<T>::operator[](idx);
  }
  
  inline const T& operator[](DegreeOfFreedom idx) const 
  {
    const T& swap(valueView(idx));
    T& swap2 = const_cast<T&>(DOFVector<T>::operator[](idx));
    swap2 = swap;
    return swap2;
  }
  
  inline const T evalAtPoint(WorldVector<double> &p, 
			      ElInfo *oldElInfo = NULL, 
			      T* valueReturn = NULL
			    ) const 
  {
    FUNCNAME("DOFVector<T>::operator[x]");
    
    WorldVector<double> x_ = coordsView(p);
    const FiniteElemSpace* feSpace = vector->getFeSpace();
    
    if (oldElInfo && feSpace->getMesh() == oldElInfo->getMesh()) {
      Mesh* mesh = feSpace->getMesh();
      DimVec<double> lambda(mesh->getDim(), NO_INIT);
      int k = oldElInfo->worldToCoord(x_, &lambda);
      if (k < 0) {
	std::vector<DegreeOfFreedom> localIndices;	// DOF-indices of all DOFs in trinangle
	mtl::dense_vector<double> uh; 		// evaluation of DOFVectors at QPs

	localIndices.resize(feSpace->getBasisFcts()->getNumber());
	uh.change_dim(mesh->getDim()+1);
	
	feSpace->getBasisFcts()->getLocalIndices(oldElInfo->getElement(), feSpace->getAdmin(), localIndices);
	for (int l = 0; l < feSpace->getBasisFcts()->getNumber(); l++)
	  uh[l] = (*vector)[localIndices[l]];
	T value = feSpace->getBasisFcts()->evalUh(lambda, uh);
	return value;
      }
    }
    
#ifdef USE_BACKGROUNDMESH
    using namespace experimental;
    bool use_backgroundmesh = false;
    Parameters::get("backgroundMesh->enabled",use_backgroundmesh);
    if (use_backgroundmesh) {
      Box* box = Box::provideBackgroundMesh(feSpace);
      int evaluation_accuracy = 0;
      Parameters::get("backgroundMesh->evaluation accuracy",evaluation_accuracy);
      T value = box->evalAtPoint(*vector, x_, evaluation_accuracy);
      return value;
    }
#endif
#ifdef USE_KD_TREE
    using namespace experimental;
    bool use_kdtree = false;
    Parameters::get("KD-Tree->enabled",use_kdtree);
    if (use_kdtree) {
      KD_Tree_Dof* tree = provideKDTree(feSpace);
      T value = tree->evalAtPoint(*vector, x_);
      return value;
    }
#endif
    DegreeOfFreedom idx = -1;
    T value = value0;    
    Mesh* mesh = feSpace->getMesh();
    
    ElInfo *elInfo = mesh->createNewElInfo();
    DimVec<double> lambda(mesh->getDim(), NO_INIT);

    bool inside = mesh->findElInfoAtPoint(x_, elInfo, lambda, (oldElInfo ? oldElInfo->getMacroElement() : NULL), NULL, NULL);
    
    if (inside) {      
      std::vector<DegreeOfFreedom> localIndices;
      mtl::dense_vector<T> uh; 	

      localIndices.resize(feSpace->getBasisFcts()->getNumber());
      uh.change_dim(feSpace->getBasisFcts()->getNumber());
    
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
      for (int l = 0; l < feSpace->getBasisFcts()->getNumber(); l++)
	uh[l] = valueView(localIndices[l]);
      value = feSpace->getBasisFcts()->evalUh(lambda, uh);
    }
    if (valueReturn)
      *valueReturn = value;
    delete elInfo;
    
    return value;
  } 
  
protected:
  
  virtual inline T valueView(DegreeOfFreedom idx) const
  {
    return (*vector)[idx];
  }
  
  virtual inline WorldVector<double> coordsView(WorldVector<double> &x) const
  {
    return x;
  }
  
  void setDefaultValue(T& value0_)
  {
    value0 = value0_;
  };
  
  DOFVector<T> *vector;
  T value0;
};

/// DOFVector-Views ... Value-Views
/// ________________________________________________________________________________________________

template<typename T> 
struct PhaseDOFView : public DOFView<T> {
  typedef T result_type;
  
  PhaseDOFView(DOFVector<T> &vec_) : DOFView<T>(&vec_) {}
  PhaseDOFView(DOFVector<T> *vec_) : DOFView<T>(vec_) {}
	
protected:
  
  inline T valueView(DegreeOfFreedom idx) const
  {
    return std::max(0.0, std::min(1.0, 0.5*((*DOFView<T>::vector)[idx]+1.0)));
  }
};


template<typename T> 
struct InvPhaseView : public DOFView<T> {
  typedef T result_type;

  InvPhaseView(DOFVector<T> &vec_) : DOFView<T>(&vec_) {}
  InvPhaseView(DOFVector<T> *vec_) : DOFView<T>(vec_) {}
 
protected:
  
  inline T valueView(DegreeOfFreedom idx) const
  {
    return -(*DOFView<T>::vector)[idx];
  }
};

/// DOFVector-Views ... Coord-Views
/// ________________________________________________________________________________________________

/// vec[x] := vec.evalAtPoint(R^(-1) * x)
template<typename T> 
struct RotateView : public DOFView<T> {
  typedef T result_type;

  RotateView(DOFVector<T> &vec_, WorldMatrix<double> R_) : DOFView<T>(&vec_), R(R_) {}
  RotateView(DOFVector<T> *vec_, WorldMatrix<double> R_) : DOFView<T>(vec_), R(R_) {}
  
  RotateView(DOFVector<T> &vec_, double alpha) : DOFView<T>(&vec_) {
    R.setDiag(1.0);
    R[0][0] = cos(alpha);  R[0][1] = -sin(alpha);
    R[1][0] = sin(alpha);  R[1][1] = cos(alpha);
  }
  RotateView(DOFVector<T> *vec_, double alpha) : DOFView<T>(vec_) {
    R.setDiag(1.0);
    R[0][0] = cos(alpha);  R[0][1] = -sin(alpha);
    R[1][0] = sin(alpha);  R[1][1] = cos(alpha);
  }
 
protected:
  
  inline WorldVector<double> coordsView(WorldVector<double> &x) const
  {
    WorldVector<double> x_;
    WorldMatrix<double> R_ = R;
    for (int i = 0; i < R_.getNumRows(); i++)
      for (int j = 0; j < i; j++)
	std::swap(R_[i][j], R_[j][i]);
    x_.multMatrixVec(R_, x);
    return x_;
  }
  
private:
  
  WorldMatrix<double> R;
};


/// vec[x] := vec.evalAtPoint(x ./ scale)
template<typename T> 
struct ScaleView : public DOFView<T> {
  typedef T result_type;

  ScaleView(DOFVector<T> &vec_, WorldVector<double> scale_) : DOFView<T>(&vec_), scale(scale_) {}
  ScaleView(DOFVector<T> *vec_, WorldVector<double> scale_) : DOFView<T>(vec_), scale(scale_) {}
 
  ScaleView(DOFVector<T> &vec_, double scale_) : DOFView<T>(&vec_) { scale.set(scale_); }
  ScaleView(DOFVector<T> *vec_, double scale_) : DOFView<T>(vec_) { scale.set(scale_); }
  
protected:
  
  inline WorldVector<double> coordsView(WorldVector<double> &x) const
  {
    WorldVector<double> x_;
    for (int i = 0; i < x_.getSize(); i++)
      x_[i] = x[i]/scale[i];
    return x_;
  }
  
private:
  
  WorldVector<double> scale;
};

/// vec[x] := vec.evalAtPoint(x - shift)
template<typename T> 
struct ShiftView : public DOFView<T> {
  typedef T result_type;

  ShiftView(DOFVector<T> &vec_, WorldVector<double> shift_) : DOFView<T>(&vec_), shift(shift_) {}
  ShiftView(DOFVector<T> *vec_, WorldVector<double> shift_) : DOFView<T>(vec_), shift(shift_) {}
  
  void setShift(WorldVector<double> &shift_) { shift = shift_; }
  
protected:
  
  inline WorldVector<double> coordsView(WorldVector<double> &x) const
  {
    WorldVector<double> x_;
    for (int i = 0; i < x_.getSize(); i++)
      x_[i] = x[i]-shift[i];
    return x_;
  }
  
private:
  
  WorldVector<double> shift;
};

/// evalAtPoint for several Data-Structures
/// ________________________________________________________________________________________________

// DOFVector can be accessed by locating the elInfo and than using barycentric coordinates
// The default value, if p outside of geometry, is 0.0
template<typename T>
inline T evalAtPoint(const DOFVector<T> &obj, WorldVector<double> &p, ElInfo* oldElInfo = NULL)
{
  const FiniteElemSpace* feSpace = obj.getFeSpace();
  Mesh* mesh = feSpace->getMesh();
  
  if (oldElInfo && feSpace->getMesh() == oldElInfo->getMesh()) {
    DimVec<double> lambda(mesh->getDim(), NO_INIT);
    oldElInfo->worldToCoord(p, &lambda);
    
    std::vector<DegreeOfFreedom> localIndices;	// DOF-indices of all DOFs in trinangle
    mtl::dense_vector<T> uh; 		// evaluation of DOFVectors at QPs

    localIndices.resize(feSpace->getBasisFcts()->getNumber());
    uh.change_dim(mesh->getDim()+1);
    
    feSpace->getBasisFcts()->getLocalIndices(oldElInfo->getElement(), feSpace->getAdmin(), localIndices);
    for (int l = 0; l < feSpace->getBasisFcts()->getNumber(); l++)
      uh[l] = obj[localIndices[l]];
    T value = feSpace->getBasisFcts()->evalUh(lambda, uh);
    return value;
  }
    
#ifdef USE_BACKGROUNDMESH
  using namespace experimental;
  bool use_backgroundmesh = false;
  Parameters::get("backgroundMesh->enabled",use_backgroundmesh);
  if (use_backgroundmesh) {
    Box* box = Box::provideBackgroundMesh(feSpace);
    int evaluation_accuracy = 0;
    Parameters::get("backgroundMesh->evaluation accuracy",evaluation_accuracy);
    T value = box->evalAtPoint(obj, p, evaluation_accuracy);
    return value;
  }
#endif
#ifdef USE_KD_TREE
    using namespace experimental;
    bool use_kdtree = false;
    Parameters::get("KD-Tree->enabled",use_kdtree);
    if (use_kdtree) {
      KD_Tree_Dof* tree = provideKDTree(feSpace);
      T value = tree->evalAtPoint(obj, p);
      return value;
    }
#endif   
  T value; nullify(value);
  
  ElInfo *elInfo = mesh->createNewElInfo();
  DimVec<double> lambda(mesh->getDim(), NO_INIT);

  std::vector<DegreeOfFreedom> localIndices;	// DOF-indices of all DOFs in trinangle
  mtl::dense_vector<T> uh; 	

  localIndices.resize(feSpace->getBasisFcts()->getNumber());
  uh.change_dim(feSpace->getBasisFcts()->getNumber());

  bool inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
  
  if (inside) {
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), obj.getFeSpace()->getAdmin(), localIndices);
    for (int l = 0; l < feSpace->getBasisFcts()->getNumber(); l++)
      uh[l] = obj[localIndices[l]];
    value = feSpace->getBasisFcts()->evalUh(lambda, uh);
  }
  delete elInfo;
  return value;
};

// DOFView can be accessed by own method
template<typename T >
inline T evalAtPoint(const DOFView<T> &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
{
  return obj.evalAtPoint(p, elInfo);
};

// intrinsic types can return values directly
template<typename T>
inline T evalAtPoint(T &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
{
  return obj;
};

// AbstractFunctions with argument-type = WorldVector can be accessed directly
template<typename T>
inline T evalAtPoint(AbstractFunction<T, WorldVector<double> > &obj, WorldVector<double> &p, ElInfo* elInfo = NULL)
{
  return obj(p);
};

/// Type-Traits for value_type of Data-Structures
/// ________________________________________________________________________________________________

template<typename T, typename enable=void> struct ValueType { typedef T type; };
template<typename T> struct ValueType<DOFVector<T> > { typedef T type; };
template<typename T> struct ValueType<std::vector<T> > { typedef T type; };
template<typename T> struct ValueType<std::list<T> > { typedef T type; };
template<typename T> struct ValueType<std::set<T> > { typedef T type; };
template<typename T> struct ValueType<mtl::dense_vector<T> > { typedef T type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<DOFView<double>, Derived > >::type> { typedef double type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<DOFView<WorldVector<double> >, Derived > >::type> { typedef WorldVector<double> type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<AbstractFunction<double, WorldVector<double> >, Derived > >::type> { typedef double type; };
template<typename Derived> struct ValueType<Derived, typename boost::enable_if<boost::is_base_of<AbstractFunction<WorldVector<double>, WorldVector<double> >, Derived > >::type> { typedef WorldVector<double> type; };
  
/// transformDOF by coords (independent of mesh and feSpace)
/// ________________________________________________________________________________________________

// result = vec
template<typename TOut, typename Vector>
inline void interpol_coords(Vector &vec, DOFVector<TOut> &result)
{
  typedef typename ValueType<Vector>::type T;
  DOFVector<WorldVector<double> > coords(result.getFeSpace(), "coords");
  result.getFeSpace()->getMesh()->getDofIndexCoords(result.getFeSpace(), coords);

  DOFIterator<WorldVector<double> > it_coords(&coords, USED_DOFS);
  DOFIterator<TOut> it_result(&result, USED_DOFS);
  for (it_coords.reset(), it_result.reset(); !it_coords.end(); ++it_coords, ++it_result) {
    *it_result = evalAtPoint<T>(vec, *it_coords);
  }
};

// result = op(vec)
template<typename TOut, typename Vector>
inline void transformDOF_coords(Vector &vec,
  DOFVector<TOut> &result,
  AbstractFunction<TOut, typename ValueType<Vector>::type> *op)
{
  typedef typename ValueType<Vector>::type T;
  DOFVector<WorldVector<double> > coords(result.getFeSpace(), "coords");
  result.getFeSpace()->getMesh()->getDofIndexCoords(result.getFeSpace(), coords);

  DOFIterator<WorldVector<double> > it_coords(&coords, USED_DOFS);
  DOFIterator<TOut> it_result(&result, USED_DOFS);
  for (it_coords.reset(), it_result.reset(); !it_coords.end(); ++it_coords, ++it_result) {
    T value = evalAtPoint<T>(vec, *it_coords);
    *it_result = (*op)(value);
  }
};

// result = binary_op(vec1, vec2)
template<typename TOut, typename Vector1, typename Vector2>
inline void transformDOF_coords(Vector1 &vec1, Vector2 &vec2,
  DOFVector<TOut> &result,
  BinaryAbstractFunction<TOut, typename ValueType<Vector1>::type, 
			       typename ValueType<Vector2>::type> *binary_op)
{
  typedef typename ValueType<Vector1>::type T1;
  typedef typename ValueType<Vector2>::type T2;
  DOFVector<WorldVector<double> > coords(result.getFeSpace(), "coords");
  result.getFeSpace()->getMesh()->getDofIndexCoords(result.getFeSpace(), coords);

  DOFIterator<WorldVector<double> > it_coords(&coords, USED_DOFS);
  DOFIterator<TOut> it_result(&result, USED_DOFS);
  for (it_coords.reset(), it_result.reset(); !it_coords.end(); ++it_coords, ++it_result) {
    T1 value1 = evalAtPoint<T1>(vec1, *it_coords);
    T2 value2 = evalAtPoint<T2>(vec2, *it_coords);
    
    *it_result = (*binary_op)(value1, value2);
  }
};

// result = tertiary_op(vec1, vec2, vec3)
template<typename TOut, typename Vector1, typename Vector2, typename Vector3>
inline void transformDOF_coords(Vector1 &vec1, Vector2 &vec2, Vector3 &vec3,
  DOFVector<TOut> &result,
  TertiaryAbstractFunction<TOut, typename ValueType<Vector1>::type, 
				 typename ValueType<Vector2>::type, 
				 typename ValueType<Vector3>::type> *tertiary_op)
{
  typedef typename ValueType<Vector1>::type T1;
  typedef typename ValueType<Vector2>::type T2;
  typedef typename ValueType<Vector3>::type T3;
  DOFVector<WorldVector<double> > coords(result.getFeSpace(), "coords");
  result.getFeSpace()->getMesh()->getDofIndexCoords(result.getFeSpace(), coords);

  DOFIterator<WorldVector<double> > it_coords(&coords, USED_DOFS);
  DOFIterator<TOut> it_result(&result, USED_DOFS);
  for (it_coords.reset(), it_result.reset(); !it_coords.end(); ++it_coords, ++it_result) {
    T1 value1 = evalAtPoint<T1>(vec1, *it_coords);
    T2 value2 = evalAtPoint<T2>(vec2, *it_coords);
    T3 value3 = evalAtPoint<T3>(vec3, *it_coords);

    *it_result = (*tertiary_op)(value1, value2, value3);
  }
};


// result = quart_op(vec1, vec2, vec3, vec4)
template<typename TOut, typename Vector1, typename Vector2, typename Vector3, typename Vector4>
inline void transformDOF_coords(Vector1 &vec1, Vector2 &vec2, Vector3 &vec3, Vector4 &vec4,
  DOFVector<TOut> &result,
  QuartAbstractFunction<TOut, typename ValueType<Vector1>::type, 
			      typename ValueType<Vector2>::type, 
			      typename ValueType<Vector3>::type, 
			      typename ValueType<Vector4>::type> *quart_op)
{
  typedef typename ValueType<Vector1>::type T1;
  typedef typename ValueType<Vector2>::type T2;
  typedef typename ValueType<Vector3>::type T3;
  typedef typename ValueType<Vector4>::type T4;
  DOFVector<WorldVector<double> > coords(result.getFeSpace(), "coords");
  result.getFeSpace()->getMesh()->getDofIndexCoords(result.getFeSpace(), coords);

  DOFIterator<WorldVector<double> > it_coords(&coords, USED_DOFS);
  DOFIterator<TOut> it_result(&result, USED_DOFS);
  for (it_coords.reset(), it_result.reset(); !it_coords.end(); ++it_coords, ++it_result) {
    T1 value1 = evalAtPoint<T1>(vec1, *it_coords);
    T2 value2 = evalAtPoint<T2>(vec2, *it_coords);
    T3 value3 = evalAtPoint<T3>(vec3, *it_coords);
    T4 value4 = evalAtPoint<T3>(vec4, *it_coords);

    *it_result = (*quart_op)(value1, value2, value3, value4);
  }
};

#endif // MY_VIEWS_H