#ifndef __Viscosity_Takahashi2015_h__
#define __Viscosity_Takahashi2015_h__

#include "SPlisHSPlasH/Common.h"
#include "SPlisHSPlasH/FluidModel.h"
#include "SPlisHSPlasH/NonPressureForceBase.h"
#include "SPlisHSPlasH/Utilities/MatrixFreeSolver.h"


namespace SPH
{
	/** \brief This class implements a variant of the implicit simulation method for
	* viscous fluids introduced by Takahashi et al. [TDF+15].
	* In the original work of Takahashi et al. the second-ring neighbors are required
	* to create the matrix of the linear system. In contrast we use a meshless 
	* conjugate gradient solver which performs the required matrix-vector multiplication 
	* in two sequential loops. In this way only the one-ring neighbors are required
	* in each loop which increases the performance significantly.\n\n
	* Thanks to Anreas Peer who helped us with the implementation.
	*
	* References:
	* - [TDF+15] T. Takahashi, Y. Dobashi, I. Fujishiro, T. Nishita, and M.C. Lin. Implicit Formulation for SPH-based Viscous Fluids. Computer Graphics Forum, 34(2):493-502, 2015. URL: http://dx.doi.org/10.1111/cgf.12578
	*/
	class Viscosity_Takahashi2015 : public NonPressureForceBase
	{
	protected: 
		std::vector<Vector3r> m_accel;
		std::vector<Matrix3r> m_viscousStress;
		typedef Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> Solver;
		Solver m_solver;
		Real m_viscosity;
		unsigned int m_iterations;
		unsigned int m_maxIter;
		Real m_maxError;

		virtual void initParameters();
		static void computeViscosityAcceleration(Viscosity_Takahashi2015 *visco, const Real* v);

	public:
		static std::string METHOD_NAME;
		static int VISCOSITY_COEFFICIENT;
		static int ITERATIONS;
		static int MAX_ITERATIONS;
		static int MAX_ERROR;

		Viscosity_Takahashi2015(FluidModel *model);
		virtual ~Viscosity_Takahashi2015(void);

		static NonPressureForceBase* creator(FluidModel* model) { return new Viscosity_Takahashi2015(model); }
		virtual std::string getMethodName() { return METHOD_NAME; }

		virtual void step();
		virtual void reset();

		virtual void performNeighborhoodSearchSort();

		static void matrixVecProd(const Real* vec, Real *result, void *userData);
		FORCE_INLINE static void diagonalMatrixElement(const unsigned int row, Real &result, void *userData);

		FORCE_INLINE const Matrix3r& getViscousStress(const unsigned int i) const
		{
			return m_viscousStress[i];
		}

		FORCE_INLINE Matrix3r& getViscousStress(const unsigned int i)
		{
			return m_viscousStress[i];
		}

		FORCE_INLINE void setViscousStress(const unsigned int i, const Matrix3r &val)
		{
			m_viscousStress[i] = val;
		}

		FORCE_INLINE const Vector3r& getAccel(const unsigned int i) const
		{
			return m_accel[i];
		}

		FORCE_INLINE Vector3r& getAccel(const unsigned int i)
		{
			return m_accel[i];
		}

		FORCE_INLINE void setAccel(const unsigned int i, const Vector3r &val)
		{
			m_accel[i] = val;
		}
	};
}

#endif