/*! \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 Material;
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
		Material    *material;                    // 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     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<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
		void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
		void   ComputeStrainRate(Vector<IssmDouble>* 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   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum);
		void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df);
		void   CreatePVector(Vector<IssmDouble>* pf);
		void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
		void   Delta18oParameterization(void);
		int    GetNodeIndex(Node* node);
		int    GetNumberOfNodes(void);
		void   GetNodesSidList(int* sidlist);
		int    Sid();
		bool   IsOnBed();
		bool   IsFloating(); 
		bool   IsNodeOnShelf(); 
		bool   IsNodeOnShelfFromFlags(IssmDouble* flags);
		bool   IsOnWater(); 
		void   GetSolutionFromInputs(Vector<IssmDouble>* solution);
		void   GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
		void   GetVectorFromResults(Vector<IssmDouble>* 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){_error_("not implemented yet");};
		void   MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding);
		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units);
		void   PotentialSheetUngrounding(Vector<IssmDouble>* 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){_error_("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<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
		IssmDouble TimeAdapt();
		int*   GetHorizontalNeighboorSids(void);
		void   SmearFunction(Vector<IssmDouble>* 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<IssmDouble>* gradient,int control_type,int control_index);
		void   GradjBGradient(Vector<IssmDouble>* gradient,int weight_index,int control_index);
		void   GradjZGradient(Vector<IssmDouble>* gradient,int weight_index,int control_index);
		void   GradjBMacAyeal(Vector<IssmDouble>* gradient,int control_index);
		void   GradjZMacAyeal(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragMacAyeal(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragStokes(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragGradient(Vector<IssmDouble>* gradient,int weight_index,int control_index);
		void   GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
		void   GradjVxBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
		void   GradjVyBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
		void   GradjThicknessWeakBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
		void   GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
		void   SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
		void   ControlInputGetGradient(Vector<IssmDouble>* 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 ThicknessAlongGradient( bool process_units,int weight_index);
		IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index);
		IssmDouble BalancethicknessMisfit(     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* CreateMassMatrix(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	         GetVertexPidList(int* doflist);
		void           GetVertexSidList(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<IssmDouble>* solution);
		void	  GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* 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<IssmDouble>* solution);
		void	  InputUpdateFromSolutionHydrology(IssmDouble* solution);
		#endif
		#ifdef _HAVE_BALANCED_
		#endif

		/*}}}*/

};
#endif  /* _TRIA_H */
