/*! \file Tria.h 
 *  \brief: header file for tria object
 */

#ifndef _TRIA_H_
#define _TRIA_H_

/*Headers:*/
/*{{{*/
#include "./Element.h"
#include "./TriaHook.h"
#include "./TriaRef.h"
class Parameters;
class Inputs;
class IoModel;
class Node;
class Matice;
class Matpar;
class ElementMatrix;
class ElementVector;

#include "../../../include/include.h"
#include "../../../shared/Exceptions/exceptions.h"
#include "../../../EnumDefinitions/EnumDefinitions.h"
/*}}}*/

class Tria: public Element,public TriaHook,public TriaRef{

	public:

		int  id;
		int  sid;

		Node   **nodes;    // 3 nodes
		Matice  *matice;   // 1 material ice
		Matpar  *matpar;   // 1 material parameter
		int      horizontalneighborsids[3];

		Parameters *parameters;   //pointer to solution parameters
		Inputs     *inputs;
		Results    *results;

		/*Tria constructors, destructors {{{*/
		Tria();
		Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels);
		~Tria();
		/*}}}*/
		/*Object virtual functions definitions:{{{ */
		void  Echo();
		void  DeepEcho();
		int   Id(); 
		int   MyRank();
		int   ObjectEnum();
		Object* copy();
		/*}}}*/
		/*Update virtual functions resolution: {{{*/
		void  InputUpdateFromSolution(IssmDouble* solutiong);
		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
		void  InputUpdateFromVector(int* vector, int name, int type);
		void  InputUpdateFromVector(bool* vector, int name, int type);
		#ifdef _HAVE_DAKOTA_
		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
		void  InputUpdateFromVectorDakota(int* vector, int name, int type);
		void  InputUpdateFromVectorDakota(bool* vector, int name, int type);
		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
		#endif
		void  InputUpdateFromConstant(IssmDouble constant, int name);
		void  InputUpdateFromConstant(int constant, int name);
		void  InputUpdateFromConstant(bool constant, int name);
		void  InputUpdateFromIoModel(int index, IoModel* iomodel);
		/*}}}*/
		/*Element virtual functions definitions: {{{*/
		void   AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
		void   ComputeBasalStress(Vector* sigma_b);
		void   ComputeStrainRate(Vector* eps);
		void   ComputeStressTensor();
		void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
		void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
		void   CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df);
		void   CreatePVector(Vector* pf);
		void   CreateJacobianMatrix(Matrix* Jff);
		void   Delta18oParameterization(void);
		int    GetNodeIndex(Node* node);
		int    Sid();
		bool   IsOnBed();
		bool   IsFloating(); 
		bool   IsNodeOnShelf(); 
		bool   IsNodeOnShelfFromFlags(IssmDouble* flags);
		bool   IsOnWater(); 
		void   GetSolutionFromInputs(Vector* solution);
		void   GetVectorFromInputs(Vector* vector, int name_enum);
		void   GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp);
		void   InputArtificialNoise(int enum_type,IssmDouble min, IssmDouble max);
		bool   InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums);
		void   InputCreate(IssmDouble scalar,int name,int code);
		void   InputCreate(IssmDouble* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
		void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum);
		void   InputDuplicate(int original_enum,int new_enum);
		void   InputScale(int enum_type,IssmDouble scale_factor);
		void   InputToResult(int enum_type,int step,IssmDouble time);
		void   DeleteResults(void);
		void   MaterialUpdateFromTemperature(void){_error2_("not implemented yet");};
		void   MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding);
		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units);
		void   PotentialSheetUngrounding(Vector* potential_sheet_ungrounding);
		void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
		void   RequestedOutput(int output_enum,int step,IssmDouble time);
		void   ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results);
		void   PatchFill(int* pcount, Patch* patch);
		void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
		void   ProcessResultsUnits(void);
		void   ResetCoordinateSystem(void){_error2_("not implemented yet");};
		void	 SmbGradients();
		IssmDouble SurfaceArea(void);
		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
		int    UpdatePotentialSheetUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
		IssmDouble TimeAdapt();
		int*   GetHorizontalNeighboorSids(void);
		void   SmearFunction(Vector* smearedvector,IssmDouble (*WeightFunction)(IssmDouble distance,IssmDouble radius),IssmDouble radius);

		#ifdef _HAVE_RESPONSES_
		IssmDouble IceVolume(void);
		IssmDouble TotalSmb(void);
		void   MinVel(IssmDouble* pminvel, bool process_units);
		void   MinVx(IssmDouble* pminvx, bool process_units);
		void   MinVy(IssmDouble* pminvy, bool process_units);
		void   MinVz(IssmDouble* pminvz, bool process_units);
		IssmDouble MassFlux(IssmDouble* segment,bool process_units);
		void   MaxAbsVx(IssmDouble* pmaxabsvx, bool process_units);
		void   MaxAbsVy(IssmDouble* pmaxabsvy, bool process_units);
		void   MaxAbsVz(IssmDouble* pmaxabsvz, bool process_units);
		void   ElementResponse(IssmDouble* presponse,int response_enum,bool process_units);
		void   MaxVel(IssmDouble* pmaxvel, bool process_units);
		void   MaxVx(IssmDouble* pmaxvx, bool process_units);
		void   MaxVy(IssmDouble* pmaxvy, bool process_units);
		void   MaxVz(IssmDouble* pmaxvz, bool process_units);
		#endif


		#ifdef _HAVE_CONTROL_
		IssmDouble DragCoefficientAbsGradient(bool process_units,int weight_index);
		void   GradientIndexing(int* indexing,int control_index);
		void   Gradj(Vector* gradient,int control_type,int control_index);
		void   GradjBGradient(Vector* gradient,int weight_index,int control_index);
		void   GradjBMacAyeal(Vector* gradient,int control_index);
		void   GradjDragMacAyeal(Vector* gradient,int control_index);
		void   GradjDragStokes(Vector* gradient,int control_index);
		void   GradjDragGradient(Vector* gradient,int weight_index,int control_index);
		void   GradjDhDtBalancedthickness(Vector* gradient,int control_index);
		void   GradjVxBalancedthickness(Vector* gradient,int control_index);
		void   GradjVyBalancedthickness(Vector* gradient,int control_index);
		void   GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data);
		void   SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
		void   ControlInputGetGradient(Vector* gradient,int enum_type,int control_index);
		void   ControlInputScaleGradient(int enum_type,IssmDouble scale);
		void   ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
		IssmDouble RheologyBbarAbsGradient(bool process_units,int weight_index);
		IssmDouble ThicknessAbsMisfit(     bool process_units,int weight_index);
		IssmDouble SurfaceAbsVelMisfit(    bool process_units,int weight_index);
		IssmDouble ThicknessAbsGradient(bool process_units,int weight_index);
		IssmDouble SurfaceRelVelMisfit(    bool process_units,int weight_index);
		IssmDouble SurfaceLogVelMisfit(    bool process_units,int weight_index);
		IssmDouble SurfaceLogVxVyMisfit(   bool process_units,int weight_index);
		IssmDouble SurfaceAverageVelMisfit(bool process_units,int weight_index);
		void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
		#endif

		/*}}}*/
		/*Tria specific routines:{{{*/
		ElementMatrix* CreateKMatrixBalancethickness(void);
		ElementMatrix* CreateKMatrixBalancethickness_DG(void);
		ElementMatrix* CreateKMatrixBalancethickness_CG(void);
		ElementMatrix* CreateKMatrixMelting(void);
		ElementMatrix* CreateKMatrixPrognostic(void);
		ElementMatrix* CreateKMatrixPrognostic_CG(void);
		ElementMatrix* CreateKMatrixPrognostic_DG(void);
		ElementMatrix* CreateKMatrixSlope(void);
		ElementVector* CreatePVectorBalancethickness(void);
		ElementVector* CreatePVectorBalancethickness_DG(void);
		ElementVector* CreatePVectorBalancethickness_CG(void);
		ElementVector* CreatePVectorPrognostic(void);
		ElementVector* CreatePVectorPrognostic_CG(void);
		ElementVector* CreatePVectorPrognostic_DG(void);
		ElementVector* CreatePVectorSlope(void);
		IssmDouble         GetArea(void);
		int            GetElementType(void);
		void	         GetDofList(int** pdoflist,int approximation_enum,int setenum);
		void	         GetDofList1(int* doflist);
		void           GetSidList(int* sidlist);
		void           GetConnectivityList(int* connectivity);
		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue,int index); //TO BE REMOVED
		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
		void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
		void	         InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
		void	         InputUpdateFromSolutionPrognostic(IssmDouble* solution);
		bool	         IsInput(int name);
		void	         SetClone(int* minranks);
		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
		
		#ifdef _HAVE_DIAGNOSTIC_
		ElementMatrix* CreateKMatrixDiagnosticMacAyeal(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
		ElementVector* CreatePVectorDiagnosticMacAyeal(void);
		ElementVector* CreatePVectorDiagnosticHutter(void);
		ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
		void	  GetSolutionFromInputsDiagnosticHoriz(Vector* solution);
		void	  GetSolutionFromInputsDiagnosticHutter(Vector* solution);
		void	  InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solution);
		void	  InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solution);
		#endif

		#ifdef _HAVE_CONTROL_
		ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
		ElementMatrix* CreateKMatrixAdjointMacAyeal(void);
		ElementVector* CreatePVectorAdjointHoriz(void);
		ElementVector* CreatePVectorAdjointStokes(void);
		ElementVector* CreatePVectorAdjointBalancethickness(void);
		void	  InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution);
		void	  InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution);
		#endif

		#ifdef _HAVE_HYDROLOGY_
		ElementMatrix* CreateKMatrixHydrology(void);
		ElementVector* CreatePVectorHydrology(void);
		void      CreateHydrologyWaterVelocityInput(void);
		void	  GetSolutionFromInputsHydrology(Vector* solution);
		void	  InputUpdateFromSolutionHydrology(IssmDouble* solution);
		#endif
		#ifdef _HAVE_BALANCED_
		#endif

		/*}}}*/

};
#endif  /* _TRIA_H */
