#ifndef ADAPTIVEMESHREFINEMENT
#define ADAPTIVEMESHREFINEMENT

/*Includes*/
/*{{{*/
/*Common includes*/
#include <iostream>
#include <fstream>
#include <string>
#include <climits>
#include <cfloat>

/*NeoPZ includes*/
#include <pz_config.h>
#include <pzreal.h>
#include <pzgmesh.h>
#include <pzvec.h>
#include <pzeltype.h>

#include <TPZRefPatternTools.h>
#include <TPZRefPatternDataBase.h>
#include <TPZRefPattern.h>

#include <tpzchangeel.h>
#include <TPZGeoElement.h>
#include <pzreftriangle.h>
#include <pzgeotriangle.h>
#include <tpzgeoelrefpattern.h>
#include <pzgraphmesh.h>
#include <TPZVTKGeoMesh.h>

#include "../shared/shared.h"

/*}}}*/

class AdaptiveMeshRefinement{

public:
	/*Public attributes{{{*/
	/* 
	 * to refine:
	 * distance_h = initial_distance * gradation ^ (level_max-h)
	 * to unrefine:
	 * distance_h = lag * initial_distance * gradation ^ (level_max-h)
	 */
	int refinement_type;						//0 uniform (faster); 1 refpattern  
	int level_max;								//max level of refinement
	double gradation;							//geometric progression ratio to calculate radius of level h
	double lag;									//lag used in the unrefine process
	/*Target and estimators*/
	double groundingline_distance;		//all elements with distance from grounding line <= groundingline_distance will be refined
	double icefront_distance;				//all elements with distance from ice front <= icefront_distance will be refined
	double thicknesserror_threshold;		//if ==0, it will not be used
	double thicknesserror_groupthreshold;//group threshold
	double thicknesserror_maximum;		//max value of the error estimator; in general, it is defined in the first time step. Attention with restart
	double deviatoricerror_threshold;	//if ==0, it will not be used
	double deviatoricerror_groupthreshold;//group threshold
	double deviatoricerror_maximum;		//max value of the error estimator; in general, it is defined in the first time step. Attention with restart
	/*}}}*/
	/*Public methods{{{*/
	/* Constructor, destructor etc*/
	AdaptiveMeshRefinement();															
	AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp); 					
	AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp);	
	virtual ~AdaptiveMeshRefinement();												
	/*General methods*/
	void CleanUp();
	void Initialize();
	void ExecuteRefinement(int numberofpoints,double* xylist,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
	void ExecuteRefinement(int numberofpoints,double* xylist,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
   void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist);
	void CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements);
	void CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements);
	/*}}}*/
private:
	/*Private attributes{{{*/
	std::vector<int> sid2index;					// Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) 
	std::vector<int> index2sid;					// Vector that keeps sid of issm mesh elements used in the neopz mesh (index) 
	std::vector<int> specialelementsindex;		// Vector that keeps index of the special elements (created to avoid haning nodes) 
	TPZGeoMesh *fathermesh;							// Entire mesh without refinement if refinement_type==1; refined with hanging nodes if efinement_type==0
	TPZGeoMesh *previousmesh;						// Refined mesh without hanging nodes (it is always refpattern type), used to generate ISSM mesh
	/*}}}*/
	/*Private methods{{{*/
   void RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror);
	void RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh);
	void RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh);
	void DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh);
	void GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py,int** pelements);
	TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh);
   inline int GetElemMaterialID(){return 1;} 
	inline int GetNumberOfNodes(){return 3;}
	void PrintGMeshVTK(TPZGeoMesh *gmesh,std::ofstream &file,bool matColor=true);
	int GetVTK_ElType(TPZGeoEl* gel);
	int VerifyRefinementType(TPZGeoEl* geoel);
	/*}}}*/
};

#endif