/*! \file Penta.h 
 *  \brief: header file for penta object
 */

#ifndef _PENTA_H_
#define _PENTA_H_

/*Headers:*/
/*{{{*/
#include "./Element.h"
#include "./ElementHook.h"
#include "./PentaRef.h"
class Object;
class Parameters;
class Results;
class Inputs;
class Input;
class IoModel;
class Node;
class Material;
class Matpar;
class Tria;
class ElementMatrix;
class ElementVector;

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

class Penta: public Element,public ElementHook,public PentaRef{

	public:

		int          id;
		int          sid;

		Node       **nodes;                       // set of nodes
		Vertex     **vertices;                    // 6 vertices
		Material    *material;                    // 1 material ice
		Matpar      *matpar;                      // 1 material parameter
		Penta      **verticalneighbors;           // 2 neighbors: first one under, second one above

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

		/*Penta constructors and destructor: {{{*/
		Penta();
		Penta(int penta_id,int penta_sid,int i, IoModel* iomodel,int nummodels);
		~Penta();
		/*}}}*/
		/*Object virtual functions definitions: {{{*/
		Object *copy();
		void    DeepEcho();
		void    Echo();
		int     ObjectEnum();
		int     Id();
		/*}}}*/
		/*Update virtual functions definitions: {{{*/
		void  InputUpdateFromConstant(bool constant, int name);
		void  InputUpdateFromConstant(IssmDouble constant, int name);
		void  InputUpdateFromConstant(int constant, int name);
		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  InputUpdateFromIoModel(int index, IoModel* iomodel);
		/*}}}*/
		/*Element virtual functions definitions: {{{*/
		void   BasalFrictionCreateInput(void);
		void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
		void   ComputeStrainRate(Vector<IssmDouble>* eps);
		void   ComputeStressTensor();
		void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,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);
		void   DeleteResults(void);
		int    GetNodeIndex(Node* node);
		void   GetNodesSidList(int* sidlist);
		void   GetNodesLidList(int* lidlist);
		int    GetNumberOfNodes(void);
		void   GetSolutionFromInputs(Vector<IssmDouble>* solution);
		IssmDouble GetZcoord(GaussPenta* gauss);
		void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
		void   GetVectorFromResults(Vector<IssmDouble>* vector,int offset,int name_enum,int interp);

		int    Sid();
		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   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   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
		void   ResetCoordinateSystem(void);
		void   SmbGradients();
		IssmDouble SurfaceArea(void);
		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
		IssmDouble TimeAdapt();
		void   ViscousHeatingCreateInput(void);

		 #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   MaxVel(IssmDouble* pmaxvel);
		void   ElementResponse(IssmDouble* presponse,int response_enum);
		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   GradjDragSSA(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragHO(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
		void   GradjBbarSSA(Vector<IssmDouble>* gradient,int control_index);
		void   GradjBbarHO(Vector<IssmDouble>* gradient,int control_index);
		void   GradjBbarFS(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 SurfaceRelVelMisfit(    int weight_index);
		IssmDouble SurfaceLogVelMisfit(    int weight_index);
		IssmDouble SurfaceLogVxVyMisfit(   int weight_index);
		IssmDouble SurfaceAverageVelMisfit(int weight_index);
		IssmDouble ThicknessAbsGradient(int weight_index);
		IssmDouble ThicknessAlongGradient( int weight_index){_error_("not supported");};
		IssmDouble ThicknessAcrossGradient(int weight_index){_error_("not supported");};
		IssmDouble BalancethicknessMisfit(int weight_index){_error_("not supported");};
		void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
		#endif

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

		/*}}}*/
		/*Penta specific routines:{{{*/
		void	         BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
		ElementMatrix* CreateBasalMassMatrix(void);
		ElementMatrix* CreateKMatrix(void);
		ElementMatrix* CreateKMatrixMasstransport(void);
		ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
		ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
		ElementVector* CreatePVector(void);
		ElementVector* CreatePVectorMasstransport(void);
		ElementVector* CreatePVectorFreeSurfaceTop(void);
		ElementVector* CreatePVectorFreeSurfaceBase(void);
		ElementVector* CreatePVectorSlope(void);
		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
		void	         GetDofList(int** pdoflist,int approximation_enum,int setenum);
		void	         GetDofListVelocity(int** pdoflist,int setenum);
		void	         GetDofListPressure(int** pdoflist,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);
		int            GetElementType(void);
		void           GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype);
		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
		void	         GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
		void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
		void	         GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
		IssmDouble     GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
		void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
		void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
		Penta*         GetUpperElement(void);
		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
		Penta*         GetLowerElement(void);
		Penta*         GetBasalElement(void);
		void	         InputExtrude(int enum_type,int object_type);
		void           InputUpdateFromSolutionMasstransport(IssmDouble* solutiong);
		void           InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solutiong);
		void           InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solutiong);
		void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
		void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
		bool	         IsInput(int name);
		bool	         IsOnSurface(void);
		bool	         IsOnBed(void);
		bool           IsFloating(void); 
		bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
		bool           NoIceInElement(void); 
		IssmDouble     MinEdgeLength(IssmDouble xyz_list[6][3]);
		void	         ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
		void	         ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
		void	         SetClone(int* minranks);
		Tria*	         SpawnTria(int location);
		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);

		#ifdef _HAVE_STRESSBALANCE_
		ElementMatrix* CreateKMatrixCouplingSSAHO(void);
		ElementMatrix* CreateKMatrixCouplingSSAHOViscous(void);
		ElementMatrix* CreateKMatrixCouplingSSAHOFriction(void);
		ElementMatrix* CreateKMatrixCouplingSSAFS(void);
		ElementMatrix* CreateKMatrixCouplingSSAFSViscous(void);
		ElementMatrix* CreateKMatrixCouplingSSAFSFriction(void);
		ElementMatrix* CreateKMatrixCouplingHOFS(void);
		ElementMatrix* CreateKMatrixStressbalanceHoriz(void);
		ElementMatrix* CreateKMatrixAdjointHoriz(void);
		ElementVector* CreateDVectorStressbalanceHoriz(void);
		ElementVector* CreateDVectorStressbalanceFS(void);
		ElementMatrix* CreateKMatrixStressbalanceSIA(void);
		ElementMatrix* CreateKMatrixStressbalanceSSA2d(void);
		ElementMatrix* CreateKMatrixStressbalanceSSA3d(void);
		ElementMatrix* CreateKMatrixStressbalanceSSA3dViscous(void);
		ElementMatrix* CreateKMatrixStressbalanceSSA3dFriction(void);
		ElementMatrix* CreateKMatrixStressbalanceSSAHO(void);
		ElementMatrix* CreateKMatrixStressbalanceSSAFS(void);
		ElementMatrix* CreateKMatrixStressbalanceL1L2(void);
		ElementMatrix* CreateKMatrixStressbalanceL1L2Viscous(void);
		ElementMatrix* CreateKMatrixStressbalanceL1L2Friction(void);
		ElementMatrix* CreateKMatrixStressbalanceHO(void);
		ElementMatrix* CreateKMatrixStressbalanceHOViscous(void);
		ElementMatrix* CreateKMatrixStressbalanceHOFriction(void);
		ElementMatrix* CreateKMatrixStressbalanceHOFS(void);
		ElementMatrix* CreateKMatrixStressbalanceFS(void);
		ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
		void           KMatrixGLSstabilization(ElementMatrix* Ke);
		ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
		ElementMatrix* CreateKMatrixStressbalanceVert(void);
		ElementMatrix* CreateKMatrixStressbalanceVertVolume(void);
		ElementMatrix* CreateKMatrixStressbalanceVertSurface(void);
		ElementMatrix* CreateJacobianStressbalanceHoriz(void);
		ElementMatrix* CreateJacobianStressbalanceSSA2d(void);
		ElementMatrix* CreateJacobianStressbalanceHO(void);
		ElementMatrix* CreateJacobianStressbalanceFS(void);
		void           InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceSSA( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceSSAHO( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceSSAFS( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceL1L2( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceHO( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceHOFS( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceVert( IssmDouble* solutiong);
		void           InputUpdateFromSolutionStressbalanceFS( IssmDouble* solutiong);
		void	         GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solutiong);
		void	         GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solutiong);
		void	         GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solutiong);
		void	         GetSolutionFromInputsStressbalanceVert(Vector<IssmDouble>* solutiong);
		ElementVector* CreatePVectorCouplingSSAFS(void);
		ElementVector* CreatePVectorCouplingSSAFSViscous(void);
		ElementVector* CreatePVectorCouplingSSAFSFriction(void);
		ElementVector* CreatePVectorCouplingHOFS(void);
		ElementVector* CreatePVectorCouplingHOFSViscous(void);
		ElementVector* CreatePVectorCouplingHOFSFriction(void);
		ElementVector* CreatePVectorStressbalanceHoriz(void);
		ElementVector* CreatePVectorStressbalanceSIA(void);
		ElementVector* CreatePVectorStressbalanceSSA(void);
		ElementVector* CreatePVectorStressbalanceSSAHO(void);
		ElementVector* CreatePVectorStressbalanceSSAFS(void);
		ElementVector* CreatePVectorStressbalanceL1L2(void);
		ElementVector* CreatePVectorStressbalanceHO(void);
		ElementVector* CreatePVectorStressbalanceHODrivingStress(void);
		ElementVector* CreatePVectorStressbalanceHOFront(void);
		ElementVector* CreatePVectorStressbalanceHOFS(void);
		ElementVector* CreatePVectorStressbalanceFS(void);
		ElementVector* CreatePVectorStressbalanceFSFront(void);
		ElementVector* CreatePVectorStressbalanceFSViscous(void);
		void           PVectorGLSstabilization(ElementVector* pe);
		ElementVector* CreatePVectorStressbalanceFSShelf(void);
		ElementVector* CreatePVectorStressbalanceVert(void);
		ElementVector* CreatePVectorStressbalanceVertVolume(void);
		ElementVector* CreatePVectorStressbalanceVertBase(void);
		void           GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*);
		#endif

		#ifdef _HAVE_CONTROL_
		ElementVector* CreatePVectorAdjointHoriz(void);
		ElementMatrix* CreateKMatrixAdjointSSA2d(void);
		ElementMatrix* CreateKMatrixAdjointHO(void);
		ElementMatrix* CreateKMatrixAdjointFS(void);
		ElementVector* CreatePVectorAdjointSSA(void);
		ElementVector* CreatePVectorAdjointHO(void);
		ElementVector* CreatePVectorAdjointFS(void);
		void           InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong);
		void           InputUpdateFromSolutionAdjointFS( IssmDouble* solutiong);
		#endif

		#ifdef _HAVE_HYDROLOGY_
		ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
		ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
		ElementVector* CreatePVectorHydrologyDCInefficient(void);
		ElementVector* CreatePVectorHydrologyDCEfficient(void);
		void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
		void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
		void    GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type);
		void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
		void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
		void    InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
		#endif

		#ifdef _HAVE_THERMAL_
		ElementMatrix* CreateKMatrixEnthalpy(void);
		ElementMatrix* CreateKMatrixEnthalpyVolume(void);
		ElementMatrix* CreateKMatrixEnthalpyShelf(void);
		ElementMatrix* CreateKMatrixThermal(void);
		ElementMatrix* CreateKMatrixMelting(void);
		ElementMatrix* CreateKMatrixThermalVolume(void);
		ElementMatrix* CreateKMatrixThermalShelf(void);
		ElementVector* CreatePVectorEnthalpy(void);
		ElementVector* CreatePVectorEnthalpyVolume(void);
		ElementVector* CreatePVectorEnthalpyShelf(void);
		ElementVector* CreatePVectorEnthalpySheet(void);
		ElementVector* CreatePVectorMelting(void);
		ElementVector* CreatePVectorThermal(void);
		ElementVector* CreatePVectorThermalVolume(void);
		ElementVector* CreatePVectorThermalShelf(void);
		ElementVector* CreatePVectorThermalSheet(void);
		void	         GetSolutionFromInputsThermal(Vector<IssmDouble>* solutiong);
		void           InputUpdateFromSolutionThermal(IssmDouble* solutiong);
		void           InputUpdateFromSolutionEnthalpy(IssmDouble* solutiong);
		void           UpdateThermalBasalConstraints(void);
		void           ComputeBasalMeltingrate(void);
		void           DrainWaterfraction(void);
		#endif

		#ifdef _HAVE_BALANCED_
		ElementMatrix* CreateKMatrixBalancethickness(void);
		ElementVector* CreatePVectorBalancethickness(void);
		#endif
		/*}}}*/
};
#endif  /* _PENTA_H */
