/*! \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;
class GaussPenta;
#include "../../shared/Exceptions/exceptions.h"
#include "../../shared/Enum/Enum.h"
/*}}}*/

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

	public:

		Penta      **verticalneighbors;           // 2 neighbors: first one under, second one above

		/*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();
		int     ObjectEnum();
		/*}}}*/
		/*Update virtual functions definitions: {{{*/
		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);
		IssmDouble CharacteristicLength(void){_error_("not implemented yet");};
		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   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
		int    FiniteElement(void);
		void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
		void   Delta18oParameterization(void);
		Penta* GetUpperPenta(void);
		Penta* GetLowerPenta(void);
		Penta* GetSurfacePenta(void);
		Penta* GetBasalPenta(void);
		Element* GetUpperElement(void);
		Element* GetLowerElement(void);
		Element* GetSurfaceElement(void);
		Element* GetBasalElement(void);
		void   GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
		IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
		int    GetNodeIndex(Node* node);
		int    GetNumberOfNodes(void);
		int    GetNumberOfVertices(void);
		void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
		IssmDouble GetXcoord(Gauss* gauss);
		IssmDouble GetYcoord(Gauss* gauss);
		IssmDouble GetZcoord(Gauss* gauss);
		void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
		void   GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
		void   GetVerticesCoordinatesTop(IssmDouble** pxyz_list);

		void   InputDepthAverageAtBase(int enum_type,int average_enum_type);
		void   InputDuplicate(int original_enum,int new_enum);
		void   InputScale(int enum_type,IssmDouble scale_factor);
		int    NumberofNodesVelocity(void);
		int    NumberofNodesPressure(void);
		int    VelocityInterpolation();
		int    PressureInterpolation();
		bool   IsZeroLevelset(int levelset_enum);
		bool   IsIcefront(void);
		void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
		void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
		void   GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");};
		void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
		void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
		void   ResetFSBasalBoundaryCondition(void);
		void   SmbGradients();
		void   SetTemporaryElementType(int element_type_in);
	   Element* SpawnBasalElement(void);
		Element* SpawnTopElement(void);
		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   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
		void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
		int    VertexConnectivity(int vertexindex);
		void   VerticalSegmentIndices(int** pindices,int* pnumseg);
		void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);

		void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
		IssmDouble IceVolume(void);
		IssmDouble IceVolumeAboveFloatation(void);
		IssmDouble TotalSmb(void);
		void   MinVel(IssmDouble* pminvel);
		void   MinVx(IssmDouble* pminvx);
		void   MinVy(IssmDouble* pminvy);
		void   MinVz(IssmDouble* pminvz);
		IssmDouble MassFlux(IssmDouble* segment);
		IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
		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);

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

		IssmDouble DragCoefficientAbsGradient(void);
		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);
		void   ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
		IssmDouble RheologyBbarAbsGradient(void);
		IssmDouble ThicknessAbsMisfit(void);
		IssmDouble SurfaceAbsVelMisfit(void);
		IssmDouble SurfaceRelVelMisfit(void);
		IssmDouble SurfaceLogVelMisfit(void);
		IssmDouble SurfaceLogVxVyMisfit(void);
		IssmDouble SurfaceAverageVelMisfit(void);
		IssmDouble ThicknessAbsGradient(void);
		IssmDouble ThicknessAlongGradient(void){_error_("not supported");};
		IssmDouble ThicknessAcrossGradient(void){_error_("not supported");};
		IssmDouble BalancethicknessMisfit(void){_error_("not supported");};
		void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
		IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
		IssmDouble MisfitArea(int weightsenum){_error_("not implemented yet");};

		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);
		/*}}}*/
		/*Penta specific routines:{{{*/
		void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
		void	         NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list);
		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
		void	         NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
		ElementMatrix* CreateBasalMassMatrix(void);
		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);

		void	         GetVertexPidList(int* doflist);
		int            GetElementType(void);
		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
		Node*          GetNode(int node_number);
		void	         InputExtrude(int enum_type);
		void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
		void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
		bool	         IsOnSurface(void);
		bool	         IsOnBed(void);
		bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
		void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
		void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
		void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
		void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
		Gauss*         NewGauss(void);
		Gauss*         NewGauss(int order);
		Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
		Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
		Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
		Gauss*         NewGaussBase(int order);
		Gauss*         NewGaussLine(int vertex1,int vertex2,int order);
		Gauss*         NewGaussTop(int order);
		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
		void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
		void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
		void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
		void           NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
		void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
		void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
		void           NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
		IssmDouble     MinEdgeLength(IssmDouble* xyz_list);
		void	         SetClone(int* minranks);
		Tria*	         SpawnTria(int index1,int index2,int index3);
		IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);

		void           UpdateConstraintsExtrudeFromBase(void);
		void           UpdateConstraintsExtrudeFromTop(void);
		/*}}}*/
};
#endif  /* _PENTA_H */
