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

#ifndef _PENTA_H_
#define _PENTA_H_

/*Headers:*/
/*{{{1*/
#include "./Element.h"
#include "./PentaHook.h"
#include "./PentaRef.h"
class  Object;
class Parameters;
class Inputs;
class IoModel;
class Node;
class Matice;
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
		Matice      *matice;       // 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: {{{1*/
		Penta();
		Penta(int penta_id,int penta_sid,int i, IoModel* iomodel,int nummodels);
		~Penta();
		/*}}}*/
		/*Object virtual functions definitions: {{{1*/
		Object*   copy();
		void	  DeepEcho();
		void	  Echo();
		int		  ObjectEnum();
		int		  Id(); 
		#ifdef _SERIAL_
		void	  Marshall(char** pmarshalled_dataset);
		int		  MarshallSize();
		void	  Demarshall(char** pmarshalled_dataset);
		#endif
		int		  MyRank();
		/*}}}*/
		/*Update virtual functions definitions: {{{1*/
		void  InputUpdateFromConstant(bool constant, int name);
		void  InputUpdateFromConstant(double constant, int name);
		void  InputUpdateFromConstant(int constant, int name);
		void  InputUpdateFromSolution(double* solutiong);
		void  InputUpdateFromVector(bool* vector, int name, int type);
		void  InputUpdateFromVector(double* 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(double* vector, int name, int type);
		void  InputUpdateFromVectorDakota(int* vector, int name, int type);
		void  InputUpdateFromMatrixDakota(double* matrix, int nows, int ncols, int name, int type);
		#endif
		void  InputUpdateFromIoModel(int index, IoModel* iomodel);
		/*}}}*/
		/*Element virtual functions definitions: {{{1*/
		void   AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part);
		void   BasalFrictionCreateInput(void);
		void   ComputeBasalStress(Vec sigma_b);
		void   ComputeStrainRate(Vec 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   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
		void   CreatePVector(Vec pf);
		void   CreateJacobianMatrix(Mat Jff);
		void   DeleteResults(void);
		int    GetNodeIndex(Node* node);
		void   GetSolutionFromInputs(Vec solution);
		double GetZcoord(GaussPenta* gauss);
		void   GetVectorFromInputs(Vec vector,int name_enum);
		void   GetVectorFromResults(Vec vector,int offset,int interp);
		
		int    Sid();
		void   InputArtificialNoise(int enum_type,double min, double max);
		bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
		void   InputCreate(double scalar,int name,int code);
		void   InputCreate(double* 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,double scale_factor);
		
		void   InputToResult(int enum_type,int step,double time);
		void   MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding);
		void   PotentialSheetUngrounding(Vec potential_sheet_ungrounding);
		void   RequestedOutput(int output_enum,int step,double time);
		void   ListResultsInfo(int** results_enums,int** results_size,double** results_times,int** results_steps,int* num_results);
		void   PatchFill(int* pcount, Patch* patch);
		void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
		void   PositiveDegreeDay(void);
		void   ProcessResultsUnits(void);
		void   ResetCoordinateSystem(void);
		double SurfaceArea(void);
		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
		int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf);
		double TimeAdapt();
		int*   GetHorizontalNeighboorSids(void);
		void   ViscousHeatingCreateInput(void);
		void   SmearFunction(Vec smearedvector,double (*WeightFunction)(double distance,double radius),double radius);

		 #ifdef _HAVE_RESPONSES_
		double IceVolume(void);
		void   MinVel(double* pminvel, bool process_units);
		void   MinVx(double* pminvx, bool process_units);
		void   MinVy(double* pminvy, bool process_units);
		void   MinVz(double* pminvz, bool process_units);
		int    NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units);
		double MassFlux(double* segment,bool process_units);
		void   MaxAbsVx(double* pmaxabsvx, bool process_units);
		void   MaxAbsVy(double* pmaxabsvy, bool process_units);
		void   MaxAbsVz(double* pmaxabsvz, bool process_units);
		void   MaxVel(double* pmaxvel, bool process_units);
		void   ElementResponse(double* presponse,int response_enum,bool process_units);
		void   MaxVx(double* pmaxvx, bool process_units);
		void   MaxVy(double* pmaxvy, bool process_units);
		void   MaxVz(double* pmaxvz, bool process_units);
		#endif

		#ifdef _HAVE_CONTROL_
		double DragCoefficientAbsGradient(bool process_units,int weight_index);
		void   GradientIndexing(int* indexing,int control_index);
		void   Gradj(Vec gradient,int control_type,int control_index);
		void   GradjDragMacAyeal(Vec gradient,int control_index);
		void   GradjDragPattyn(Vec gradient,int control_index);
		void   GradjDragStokes(Vec gradient,int control_index);
		void   GradjBbarMacAyeal(Vec gradient,int control_index);
		void   GradjBbarPattyn(Vec gradient,int control_index);
		void   GradjBbarStokes(Vec gradient,int control_index);
		void   GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index,const char* data);
		void   SetControlInputsFromVector(double* vector,int control_enum,int control_index);
		void   ControlInputGetGradient(Vec gradient,int enum_type,int control_index);
		void   ControlInputScaleGradient(int enum_type,double scale);
		void   ControlInputSetGradient(double* gradient,int enum_type,int control_index);
		double RheologyBbarAbsGradient(bool process_units,int weight_index);
		double ThicknessAbsMisfit(     bool process_units,int weight_index);
		double SurfaceAbsVelMisfit(    bool process_units,int weight_index);
		double SurfaceRelVelMisfit(    bool process_units,int weight_index);
		double SurfaceLogVelMisfit(    bool process_units,int weight_index);
		double SurfaceLogVxVyMisfit(   bool process_units,int weight_index);
		double SurfaceAverageVelMisfit(bool process_units,int weight_index);
		double ThicknessAbsGradient(bool process_units,int weight_index);
		void   InputControlUpdate(double scalar,bool save_parameter);
		#endif
		/*}}}*/
		/*Penta specific routines:{{{1*/
		void	  BedNormal(double* bed_normal, double xyz_list[3][3]);
		ElementMatrix* CreateKMatrixPrognostic(void);
		ElementMatrix* CreateKMatrixSlope(void);
		ElementVector* CreatePVectorPrognostic(void);
		ElementVector* CreatePVectorSlope(void);
		void	  GetDofList(int** pdoflist,int approximation_enum,int setenum);
		void	  GetDofList1(int* doflist);
		void    GetSidList(int* sidlist);
		void    GetConnectivityList(int* connectivity);
		int     GetElementType(void);
		void    GetElementSizes(double* hx,double* hy,double* hz);
		void    GetInputListOnVertices(double* pvalue,int enumtype);
		void    GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue);
		void    GetInputValue(double* pvalue,Node* node,int enumtype);
		void	  GetPhi(double* phi, double*  epsilon, double viscosity);
		void	  GetSolutionFromInputsEnthalpy(Vec solutiong);
		double  GetStabilizationParameter(double u, double v, double w, double diameter, double kappa);
		void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
		void    GetStrainRate3d(double* epsilon,double* 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(double* solutiong);
		void    InputUpdateFromSolutionOneDof(double* solutiong,int enum_type);
		void    InputUpdateFromSolutionOneDofCollapsed(double* solutiong,int enum_type);
		bool	  IsInput(int name);
		bool	  IsOnSurface(void);
		bool	  IsOnBed(void);
		bool    IsFloating(void); 
		bool    IsNodeOnShelf(); 
		bool    IsNodeOnShelfFromFlags(double* flags);
		bool    IsOnWater(void); 
		double  MinEdgeLength(double xyz_list[6][3]);
		void	  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
		void	  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
		void	  SetClone(int* minranks);
		Tria*	  SpawnTria(int g0, int g1, int g2);
		void	  SurfaceNormal(double* surface_normal, double 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* 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( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticMacAyealStokes( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticPattynStokes( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticHutter( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticVert( double* solutiong);
		void           InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
		void	         GetSolutionFromInputsDiagnosticHoriz(Vec solutiong);
		void	         GetSolutionFromInputsDiagnosticHutter(Vec solutiong);
		void	         GetSolutionFromInputsDiagnosticStokes(Vec solutiong);
		void	         GetSolutionFromInputsDiagnosticVert(Vec 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* CreatePVectorDiagnosticPattyn(void);
		ElementVector* CreatePVectorDiagnosticPattynStokes(void);
		ElementVector* CreatePVectorDiagnosticStokes(void);
		ElementVector* CreatePVectorDiagnosticStokesViscous(void);
		ElementVector* CreatePVectorDiagnosticStokesShelf(void);
		ElementVector* CreatePVectorDiagnosticVert(void);
		ElementVector* CreatePVectorDiagnosticVertVolume(void);
		ElementVector* CreatePVectorDiagnosticVertBase(void);
		#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( double* solutiong);
		void    InputUpdateFromSolutionAdjointStokes( double* 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(Vec solutiong);
		void           InputUpdateFromSolutionThermal( double* solutiong);
		void           InputUpdateFromSolutionEnthalpy( double* solutiong);
		#endif
		#ifdef _HAVE_BALANCED_
		ElementMatrix* CreateKMatrixBalancethickness(void);
		ElementVector* CreatePVectorBalancethickness(void);
		#endif
		/*}}}*/
};
#endif  /* _PENTA_H */
