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

#ifndef _TRIA_H_
#define _TRIA_H_

/*Headers:*/
/*{{{*/
#include "./Element.h"
#include "./ElementHook.h"
#include "./TriaRef.h"
class Parameters;
class Inputs;
class IoModel;
class Results;
class Node;
class Material;
class Matpar;
class ElementMatrix;
class ElementVector;
class Vertex;

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

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

	public:

		int          id;
		int          sid;

		Node       **nodes;                       // nodes
		Vertex     **vertices;                    // 3 vertices
		Material    *material;                    // 1 material ice
		Matpar      *matpar;                      // 1 material parameter

		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);
		#ifdef _HAVE_DAKOTA_
		void  InputUpdateFromVectorDakota(IssmDouble* 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        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
		void        ComputeStrainRate(Vector<IssmDouble>* eps);
		void        ComputeStressTensor();
		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
		void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
		void        CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
		void        CreateDVector(Vector<IssmDouble>* df);
		void        CreatePVector(Vector<IssmDouble>* pf);
		void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
		void        Delta18oParameterization(void);
		int         GetNodeIndex(Node* node);
		void        GetNodesSidList(int* sidlist);
		void        GetNodesLidList(int* lidlist);
		int         GetNumberOfNodes(void);
		int         Sid();
		bool        IsOnBed();
		bool        IsFloating();
		bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
		bool        NoIceInElement();
		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);
		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");};
		int         NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
		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        ResetCoordinateSystem(void){_error_("not implemented yet");};
		void	      SmbGradients();
		IssmDouble  SurfaceArea(void);
		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
		IssmDouble  TimeAdapt();

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

		#ifdef _HAVE_GIA_
		void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
		#endif

		#ifdef _HAVE_CONTROL_
		IssmDouble DragCoefficientAbsGradient(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       GradjBSSA(Vector<IssmDouble>* gradient,int control_index);
		void       GradjZSSA(Vector<IssmDouble>* gradient,int control_index);
		void       GradjDragSSA(Vector<IssmDouble>* gradient,int control_index);
		void       GradjDragFS(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       GradjThicknessBalancethicknessSoft(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(int weight_index);
		IssmDouble ThicknessAbsMisfit(     int weight_index);
		IssmDouble SurfaceAbsVelMisfit(    int weight_index);
		IssmDouble ThicknessAbsGradient(   int weight_index);
		IssmDouble ThicknessAlongGradient( int weight_index);
		IssmDouble ThicknessAcrossGradient(int weight_index);
		IssmDouble BalancethicknessMisfit(     int weight_index);
		IssmDouble SurfaceRelVelMisfit(    int weight_index);
		IssmDouble SurfaceLogVelMisfit(    int weight_index);
		IssmDouble SurfaceLogVxVyMisfit(   int weight_index);
		IssmDouble SurfaceAverageVelMisfit(int weight_index);
		void       InputControlUpdate(IssmDouble scalar,bool save_parameter);
		#endif

		#ifdef _HAVE_GROUNDINGLINE_
		void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
		void   MigrateGroundingLine(IssmDouble* sheet_ungrounding);
		int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
		#endif

		/*}}}*/
		/*Tria specific routines:{{{*/
		ElementMatrix* CreateKMatrix(void);
		ElementMatrix* CreateKMatrixBalancethickness(void);
		ElementMatrix* CreateKMatrixBalancethickness_DG(void);
		ElementMatrix* CreateKMatrixBalancethickness_CG(void);
		ElementMatrix* CreateKMatrixBalancevelocity(void);
		ElementMatrix* CreateKMatrixSmoothedSlope(void);
		ElementMatrix* CreateKMatrixMelting(void);
		ElementMatrix* CreateKMatrixMasstransport(void);
		ElementMatrix* CreateKMatrixMasstransport_CG(void);
		ElementMatrix* CreateKMatrixMasstransport_DG(void);
		ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
		ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
		ElementMatrix* CreateMassMatrix(void);
		ElementVector* CreatePVector(void);
		ElementVector* CreatePVectorBalancethickness(void);
		ElementVector* CreatePVectorBalancethickness_DG(void);
		ElementVector* CreatePVectorBalancethickness_CG(void);
		ElementVector* CreatePVectorBalancevelocity(void);
		ElementVector* CreatePVectorSmoothedSlopeX(void);
		ElementVector* CreatePVectorSmoothedSlopeY(void);
		ElementVector* CreatePVectorMasstransport(void);
		ElementVector* CreatePVectorMasstransport_CG(void);
		ElementVector* CreatePVectorMasstransport_DG(void);
		ElementVector* CreatePVectorFreeSurfaceTop(void);
		ElementVector* CreatePVectorFreeSurfaceBase(void);
		ElementVector* CreatePVectorSlope(void);
		IssmDouble     GetArea(void);
		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
		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           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
		void           GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]);
		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
		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           GetInputListOnNodes(IssmDouble* pvalue,int enumtype);
		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
		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	         InputUpdateFromSolutionMasstransport(IssmDouble* solution);
		bool	         IsInput(int name);
		void	         SetClone(int* minranks);
		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);

		#ifdef _HAVE_STRESSBALANCE_
		ElementMatrix* CreateKMatrixStressbalanceSSA(void);
		ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void);
		ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void);
		ElementMatrix* CreateKMatrixStressbalanceSIA(void);
		ElementVector* CreatePVectorStressbalanceSSA(void);
		ElementVector* CreatePVectorStressbalanceSSADrivingStress(void);
		ElementVector* CreatePVectorStressbalanceSSAFront(void);
		ElementVector* CreatePVectorStressbalanceSIA(void);
		ElementMatrix* CreateJacobianStressbalanceSSA(void);
		void	  GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
		void	  GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
		void	  InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
		void	  InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
		#endif

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

		#ifdef _HAVE_THERMAL_
		void UpdateThermalBasalConstraints(void){_error_("not implemented yet");};
		void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
	        void DrainWaterfraction(void){_error_("not implemented yet");};
		#endif

		#ifdef _HAVE_HYDROLOGY_
		ElementMatrix* CreateKMatrixHydrologyShreve(void);
		ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
		ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
		ElementVector* CreatePVectorHydrologyShreve(void);
		ElementVector* CreatePVectorHydrologyDCInefficient(void);
		ElementVector* CreatePVectorHydrologyDCEfficient(void);
		void    GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
		void    CreateHydrologyWaterVelocityInput(void);
		void	  InputUpdateFromSolutionHydrology(IssmDouble* solution);
		void	  InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
		void    InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
		void	  InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
		void	  InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
		void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
		void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
		void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
		void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
		bool    AllActive(void);
		bool    AnyActive(void);
		#endif
		#ifdef _HAVE_BALANCED_
		#endif

		/*}}}*/

};
#endif  /* _TRIA_H */
