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

#ifndef _PENTA_H_
#define _PENTA_H_

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

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

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

	public:

		int          id;
		int          sid;

		Node       **nodes;                       // 6 nodes
		Material    *material;                    // 1 material ice
		Matpar      *matpar;                      // 1 material parameter
		Penta      **verticalneighbors;           // 2 neighbors: first one under, second one above
		int          horizontalneighborsids[3];

		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(bool* vector, int name, int type);
		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
		void  InputUpdateFromVector(int* vector, int name, int type);
		#ifdef _HAVE_DAKOTA_
		void  InputUpdateFromVectorDakota(bool* vector, int name, int type);
		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
		void  InputUpdateFromVectorDakota(int* 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   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
		void   BasalFrictionCreateInput(void);
		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);
		void   DeleteResults(void);
		int    GetNodeIndex(Node* node);
		void   GetNodesSidList(int* sidlist);
		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);
		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   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding);
		void   PotentialSheetUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
		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   ProcessResultsUnits(void);
		void   ResetCoordinateSystem(void);
		void   SmbGradients();
		IssmDouble SurfaceArea(void);
		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
		int    UpdatePotentialSheetUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum,bool process_units);
		IssmDouble TimeAdapt();
		int*   GetHorizontalNeighboorSids(void);
		void   ViscousHeatingCreateInput(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   MaxVel(IssmDouble* pmaxvel, bool process_units);
		void   ElementResponse(IssmDouble* presponse,int response_enum,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   GradjDragMacAyeal(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragPattyn(Vector<IssmDouble>* gradient,int control_index);
		void   GradjDragStokes(Vector<IssmDouble>* gradient,int control_index);
		void   GradjBbarMacAyeal(Vector<IssmDouble>* gradient,int control_index);
		void   GradjBbarPattyn(Vector<IssmDouble>* gradient,int control_index);
		void   GradjBbarStokes(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 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);
		IssmDouble ThicknessAbsGradient(bool process_units,int weight_index);
		IssmDouble ThicknessAlongGradient( bool process_units,int weight_index){_error_("not supported");};
		IssmDouble ThicknessAcrossGradient(bool process_units,int weight_index){_error_("not supported");};
		void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
		#endif
		/*}}}*/
		/*Penta specific routines:{{{*/
		void	  BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
		ElementMatrix* CreateBasalMassMatrix(void);
		ElementMatrix* CreateKMatrixPrognostic(void);
		ElementVector* CreatePVectorPrognostic(void);
		ElementVector* CreatePVectorSlope(void);
		void	  GetDofList(int** pdoflist,int approximation_enum,int setenum);
		void	  GetVertexPidList(int* doflist);
		void    GetVertexSidList(int* sidlist);
		void    GetConnectivityList(int* connectivity);
		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    GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
		void	  GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
		void	  GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
		IssmDouble  GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
		void    GetStrainRate3dPattyn(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);
		Penta*  GetLowerElement(void);
		Penta*  GetBasalElement(void);
		void	  InputExtrude(int enum_type,int object_type);
		void    InputUpdateFromSolutionPrognostic(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    IsNodeOnShelf(); 
		bool    IsNodeOnShelfFromFlags(IssmDouble* flags);
		bool    IsOnWater(void); 
		IssmDouble  MinEdgeLength(IssmDouble xyz_list[6][3]);
		void	  ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
		void	  ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
		void	  SetClone(int* minranks);
		Tria*	  SpawnTria(int g0, int g1, int g2);
		void	  SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);

		#ifdef _HAVE_DIAGNOSTIC_
		ElementMatrix* CreateKMatrixCouplingMacAyealPattyn(void);
		ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void);
		ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
		ElementMatrix* CreateKMatrixCouplingMacAyealStokes(void);
		ElementMatrix* CreateKMatrixCouplingMacAyealStokesViscous(void);
		ElementMatrix* CreateKMatrixCouplingMacAyealStokesFriction(void);
		ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
		ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
		ElementMatrix* CreateKMatrixAdjointHoriz(void);
		ElementVector* CreateDVectorDiagnosticHoriz(void);
		ElementVector* CreateDVectorDiagnosticStokes(void);
		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dViscous(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void);
		ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void);
		ElementMatrix* CreateKMatrixDiagnosticL1L2(void);
		ElementMatrix* CreateKMatrixDiagnosticL1L2Viscous(void);
		ElementMatrix* CreateKMatrixDiagnosticL1L2Friction(void);
		ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
		ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
		ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
		ElementMatrix* CreateKMatrixDiagnosticPattynStokes(void);
		ElementMatrix* CreateKMatrixDiagnosticStokes(void);
		ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
		ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
		ElementMatrix* CreateKMatrixDiagnosticVert(void);
		ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
		ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
		ElementMatrix* CreateJacobianDiagnosticHoriz(void);
		ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void);
		ElementMatrix* CreateJacobianDiagnosticPattyn(void);
		ElementMatrix* CreateJacobianDiagnosticStokes(void);
		void           InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticMacAyeal( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticL1L2( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong);
		void           InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong);
		void	         GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solutiong);
		void	         GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solutiong);
		void	         GetSolutionFromInputsDiagnosticStokes(Vector<IssmDouble>* solutiong);
		void	         GetSolutionFromInputsDiagnosticVert(Vector<IssmDouble>* solutiong);
		ElementVector* CreatePVectorCouplingMacAyealStokes(void);
		ElementVector* CreatePVectorCouplingMacAyealStokesViscous(void);
		ElementVector* CreatePVectorCouplingMacAyealStokesFriction(void);
		ElementVector* CreatePVectorCouplingPattynStokes(void);
		ElementVector* CreatePVectorCouplingPattynStokesViscous(void);
		ElementVector* CreatePVectorCouplingPattynStokesFriction(void);
		ElementVector* CreatePVectorDiagnosticHoriz(void);
		ElementVector* CreatePVectorDiagnosticHutter(void);
		ElementVector* CreatePVectorDiagnosticMacAyeal(void);
		ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
		ElementVector* CreatePVectorDiagnosticMacAyealStokes(void);
		ElementVector* CreatePVectorDiagnosticL1L2(void);
		ElementVector* CreatePVectorDiagnosticPattyn(void);
		ElementVector* CreatePVectorDiagnosticPattynStokes(void);
		ElementVector* CreatePVectorDiagnosticStokes(void);
		ElementVector* CreatePVectorDiagnosticStokesViscous(void);
		ElementVector* CreatePVectorDiagnosticStokesShelf(void);
		ElementVector* CreatePVectorDiagnosticVert(void);
		ElementVector* CreatePVectorDiagnosticVertVolume(void);
		ElementVector* CreatePVectorDiagnosticVertBase(void);
		void GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*);
		#endif

		#ifdef _HAVE_CONTROL_
		ElementVector* CreatePVectorAdjointHoriz(void);
		ElementMatrix* CreateKMatrixAdjointMacAyeal2d(void);
		ElementMatrix* CreateKMatrixAdjointPattyn(void);
		ElementMatrix* CreateKMatrixAdjointStokes(void);
		ElementVector* CreatePVectorAdjointMacAyeal(void);
		ElementVector* CreatePVectorAdjointPattyn(void);
		ElementVector* CreatePVectorAdjointStokes(void);
		void    InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong);
		void    InputUpdateFromSolutionAdjointStokes( IssmDouble* solutiong);
		#endif

		#ifdef _HAVE_HYDROLOGY_
		void    CreateHydrologyWaterVelocityInput(void);
		#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);
		#endif
		#ifdef _HAVE_BALANCED_
		ElementMatrix* CreateKMatrixBalancethickness(void);
		ElementVector* CreatePVectorBalancethickness(void);
		#endif
		/*}}}*/
};
#endif  /* _PENTA_H */
