/*!\file Pengrid.c
 * \brief: implementation of the Pengrid object
 */

/*Headers*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../classes.h"
#include "shared/shared.h"
#include "../../analyses/analyses.h"
/*}}}*/

/*Element macros*/
#define NUMVERTICES   1

/*Pengrid constructors and destructor*/
Pengrid::Pengrid(){/*{{{*/
	this->parameters=NULL;
	this->hnode=NULL;
	this->node=NULL;
	this->helement=NULL;
	this->element=NULL;
	this->hmatpar=NULL;
	this->matpar=NULL;

	/*not active, not zigzagging: */
	active=0;
	zigzag_counter=0;

}
/*}}}*/
Pengrid::Pengrid(int id, int index, IoModel* iomodel, int in_analysis_type){ //i is the element index/*{{{*/

	int pengrid_node_id;
	int pengrid_matpar_id;
	int pengrid_element_id;

	/*Some checks if debugging activated*/
	_assert_(iomodel->singlenodetoelementconnectivity);
	_assert_(index>=0 && index<iomodel->numberofvertices);
	_assert_(id);

	/*id: */
	this->id=id;
	this->analysis_type=in_analysis_type;

	/*hooks: */
	pengrid_node_id=iomodel->nodecounter+index+1;
	pengrid_element_id=iomodel->singlenodetoelementconnectivity[index];
	_assert_(pengrid_element_id);
	pengrid_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object

	this->hnode=new Hook(&pengrid_node_id,1);
	this->helement=new Hook(&pengrid_element_id,1);
	this->hmatpar=new Hook(&pengrid_matpar_id,1);

	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
	this->parameters=NULL;
	this->node=NULL;
	this->element=NULL;
	this->matpar=NULL;

	//let's not forget internals
	this->active=0;
	this->zigzag_counter=0;

}
/*}}}*/
Pengrid::~Pengrid(){/*{{{*/
	delete hnode;
	delete helement;
	delete hmatpar;
	return;
}
/*}}}*/

/*Object virtual functions definitions:*/
void Pengrid::Echo(void){/*{{{*/
	this->DeepEcho();
}
/*}}}*/
void Pengrid::DeepEcho(void){/*{{{*/

	_printf_("Pengrid:\n");
	_printf_("   id: " << id << "\n");
	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
	hnode->DeepEcho();
	helement->DeepEcho();
	hmatpar->DeepEcho();
	_printf_("   active " << this->active << "\n");
	_printf_("   zigzag_counter " << this->zigzag_counter << "\n");
	_printf_("   parameters\n");
	parameters->DeepEcho();
}
/*}}}*/
int    Pengrid::Id(void){ return id; }/*{{{*/
/*}}}*/
int Pengrid::ObjectEnum(void){/*{{{*/

	return PengridEnum;
}
/*}}}*/
Object* Pengrid::copy() {/*{{{*/

	Pengrid* pengrid=NULL;

	pengrid=new Pengrid();

	/*copy fields: */
	pengrid->id=this->id;
	pengrid->analysis_type=this->analysis_type;

	/*point parameters: */
	pengrid->parameters=this->parameters;

	/*now deal with hooks and objects: */
	pengrid->hnode=(Hook*)this->hnode->copy();
	pengrid->hmatpar=(Hook*)this->hmatpar->copy();
	pengrid->helement=(Hook*)this->helement->copy();

	/*corresponding fields*/
	pengrid->node  =(Node*)pengrid->hnode->delivers();
	pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers();
	pengrid->element=(Element*)pengrid->helement->delivers();

	//let's not forget internals
	pengrid->active=this->active=0;
	pengrid->zigzag_counter=this->zigzag_counter=0;

	return pengrid;

}
/*}}}*/

/*Load virtual functions definitions:*/
void  Pengrid::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/

	/*Take care of hooking up all objects for this load, ie links the objects in the hooks to their respective 
	 * datasets, using internal ids and offsets hidden in hooks: */
	hnode->configure(nodesin);
	helement->configure(elementsin);
	hmatpar->configure(materialsin);

	/*Get corresponding fields*/
	node=(Node*)hnode->delivers();
	element=(Element*)helement->delivers();
	matpar=(Matpar*)hmatpar->delivers();

	/*point parameters to real dataset: */
	this->parameters=parametersin;
}
/*}}}*/
void  Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/

}
/*}}}*/
void  Pengrid::ResetHooks(){/*{{{*/

	this->node=NULL;
	this->element=NULL;
	this->matpar=NULL;
	this->parameters=NULL;

	/*Get Element type*/
	this->hnode->reset();
	this->helement->reset();
	this->hmatpar->reset();

}
/*}}}*/
void  Pengrid::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/

	/*No loads applied, do nothing: */
	return;

}
/*}}}*/
void  Pengrid::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/

	ElementVector* pe=NULL;
	int analysis_type;
	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	switch(analysis_type){
		
	case HydrologyDCInefficientAnalysisEnum:
		pe = CreatePVectorHydrologyDCInefficient();
		break;
	default:
		/*No loads applied, do nothing: */
		return;
	}
	if(pe){
		pe->AddToGlobal(pf);
		delete pe;
	}

}
/*}}}*/
void  Pengrid::GetNodesSidList(int* sidlist){/*{{{*/

	_assert_(sidlist);
	_assert_(node);

	sidlist[0]=node->Sid();
}
/*}}}*/
void  Pengrid::GetNodesLidList(int* lidlist){/*{{{*/

	_assert_(lidlist);
	_assert_(node);

	lidlist[0]=node->Lid();
}
/*}}}*/
int   Pengrid::GetNumberOfNodes(void){/*{{{*/

	return NUMVERTICES;
}
/*}}}*/
void  Pengrid::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/

	/*Retrieve parameters: */
	ElementMatrix* Ke=NULL;
	int analysis_type;
	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	switch(analysis_type){
		case ThermalAnalysisEnum:
			Ke=PenaltyCreateKMatrixThermal(kmax);
			break;
		case MeltingAnalysisEnum:
			Ke=PenaltyCreateKMatrixMelting(kmax);
			break;
		case HydrologyDCInefficientAnalysisEnum:
			Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax);
			break;
		default:
			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
	}

	/*Add to global matrix*/
	if(Ke){
		Ke->AddToGlobal(Kff,Kfs);
		delete Ke;
	}
}
/*}}}*/
void  Pengrid::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/

	/*Retrieve parameters: */
	ElementVector* pe=NULL;
	int analysis_type;
	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	switch(analysis_type){
		case ThermalAnalysisEnum:
			pe=PenaltyCreatePVectorThermal(kmax);
			break;
		case MeltingAnalysisEnum:
			pe=PenaltyCreatePVectorMelting(kmax);
			break;
		case StressbalanceAnalysisEnum: case AdjointHorizAnalysisEnum:
			break;
		case HydrologyDCInefficientAnalysisEnum:
			pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax);
			break;
		default:
			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
	}

	/*Add to global Vector*/
	if(pe){
		pe->AddToGlobal(pf);
		delete pe;
	}
}
/*}}}*/
bool  Pengrid::InAnalysis(int in_analysis_type){/*{{{*/
	if (in_analysis_type==this->analysis_type)return true;
	else return false;
}
/*}}}*/
bool  Pengrid::IsPenalty(void){/*{{{*/
	return true;
}
/*}}}*/
void  Pengrid::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/

	/*Output */
	int d_nz = 0;
	int o_nz = 0;

	if(!flags[this->node->Lid()]){

		/*flag current node so that no other element processes it*/
		flags[this->node->Lid()]=true;

		int counter=0;
		while(flagsindices[counter]>=0) counter++;
		flagsindices[counter]=this->node->Lid();

		/*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
		switch(set2_enum){
			case FsetEnum:
				if(node->indexing.fsize){
					if(this->node->IsClone())
					 o_nz += 1;
					else
					 d_nz += 1;
				}
				break;
			case GsetEnum:
				if(node->indexing.gsize){
					if(this->node->IsClone())
					 o_nz += 1;
					else
					 d_nz += 1;
				}
				break;
			case SsetEnum:
				if(node->indexing.ssize){
					if(this->node->IsClone())
					 o_nz += 1;
					else
					 d_nz += 1;
				}
				break;
			default: _error_("not supported");
		}
	}

	/*Assign output pointers: */
	*pd_nz=d_nz;
	*po_nz=o_nz;
}
/*}}}*/

/*Update virtual functions definitions:*/
void  Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Pengrid::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
	/*Nothing*/
}
/*}}}*/
void  Pengrid::InputUpdateFromConstant(int constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Pengrid::InputUpdateFromConstant(bool constant, int name){/*{{{*/

	switch(name){

		case ResetPenaltiesEnum:
			if (constant) zigzag_counter=0;
			return;

	}
}
/*}}}*/

/*Pengrid management:*/
void  Pengrid::ConstraintActivate(int* punstable){/*{{{*/

	int analysis_type;

	/*Retrieve parameters: */
	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	if (analysis_type==StressbalanceAnalysisEnum){
		/*No penalty to check*/
		return;
	}
	else if (analysis_type==ThermalAnalysisEnum){
		ConstraintActivateThermal(punstable);
	}
	else if (analysis_type==MeltingAnalysisEnum){
		/*No penalty to check*/
		return;
	}
	else if (analysis_type==HydrologyDCInefficientAnalysisEnum){
		ConstraintActivateHydrologyDCInefficient(punstable);
		return;
	}
	else{
		_error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
	}

}
/*}}}*/
void  Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/

	//   The penalty is stable if it doesn't change during to successive iterations.   
	IssmDouble pressure;
	IssmDouble temperature;
	IssmDouble t_pmp;
	int        new_active;
	int        unstable=0;
	int        penalty_lock;

	/*recover pointers: */
	Penta* penta=(Penta*)element;

	/*check that pengrid is not a clone (penalty to be added only once)*/
	if (node->IsClone()){
		unstable=0;
		*punstable=unstable;
		return;
	}

	//First recover pressure and temperature values, using the element: */
	penta->GetInputValue(&pressure,node,PressureEnum);
	penta->GetInputValue(&temperature,node,TemperaturePicardEnum);

	//Recover our data:
	parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum);

	//Compute pressure melting point
	t_pmp=matpar->TMeltingPoint(pressure);

	//Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.

	if (temperature>t_pmp){
		new_active=1;
	}
	else{
		new_active=0;
	}

	//Figure out stability of this penalty
	if (active==new_active){
		unstable=0;
	}
	else{
		unstable=1;
		if(penalty_lock)zigzag_counter++;
	}

	/*If penalty keeps zigzagging more than 5 times: */
	if(penalty_lock){
		if(zigzag_counter>penalty_lock){
			unstable=0;
			active=1;
		}
	}

	//Set penalty flag
	active=new_active;

	//*Assign output pointers:*/
	*punstable=unstable;
}
/*}}}*/
ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){/*{{{*/

	IssmDouble pressure,temperature,t_pmp;
	IssmDouble penalty_factor;

	Penta* penta=(Penta*)element;

	/*check that pengrid is not a clone (penalty to be added only once)*/
	if (node->IsClone()) return NULL;
	ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters);

	/*Retrieve all parameters*/
	penta->GetInputValue(&pressure,node,PressureEnum);
	penta->GetInputValue(&temperature,node,TemperatureEnum);
	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);

	/*Compute pressure melting point*/
	t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;

	/*Add penalty load*/
	if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
		Ke->values[0]=kmax*pow(10.,penalty_factor);
	}

	/*Clean up and return*/
	return Ke;
}
/*}}}*/
ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(IssmDouble kmax){/*{{{*/

	IssmDouble    penalty_factor;

	/*Initialize Element matrix and return if necessary*/
	if(!this->active) return NULL;
	ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);

	/*recover parameters: */
	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);

	Ke->values[0]=kmax*pow(10.,penalty_factor);

	/*Clean up and return*/
	return Ke;
}
/*}}}*/
ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){/*{{{*/

	IssmDouble pressure;
	IssmDouble temperature;
	IssmDouble melting_offset;
	IssmDouble t_pmp;
	IssmDouble dt,penalty_factor;

	/*recover pointers: */
	Penta* penta=(Penta*)element;

	/*check that pengrid is not a clone (penalty to be added only once)*/
	if (node->IsClone()) return NULL;
	ElementVector* pe=new ElementVector(&node,NUMVERTICES,this->parameters);

	/*Retrieve all parameters*/
	penta->GetInputValue(&pressure,node,PressureEnum);
	penta->GetInputValue(&temperature,node,TemperatureEnum);
	parameters->FindParam(&melting_offset,MeltingOffsetEnum);
	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);

	/*Compute pressure melting point*/
	t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;

	/*Add penalty load
	  This time, the penalty must have the same value as the one used for the thermal computation
	  so that the corresponding melting can be computed correctly
	  In the thermal computation, we used kmax=melting_offset, and the same penalty_factor*/
	if (temperature<t_pmp){ //%no melting
		pe->values[0]=0;
	}
	else{
		if (reCast<bool>(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt;
		else    pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp);
	}

	/*Clean up and return*/
	return pe;
}
/*}}}*/
ElementVector* Pengrid::PenaltyCreatePVectorThermal(IssmDouble kmax){/*{{{*/

	IssmDouble pressure;
	IssmDouble t_pmp;
	IssmDouble penalty_factor;

	Penta* penta=(Penta*)element;

	/*Initialize Element matrix and return if necessary*/
	if(!this->active) return NULL;
	ElementVector* pe=new ElementVector(&node,1,this->parameters);

	/*Retrieve all parameters*/
	penta->GetInputValue(&pressure,node,PressureEnum);
	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);

	/*Compute pressure melting point*/
	t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;

	pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp;

	/*Clean up and return*/
	return pe;
}
/*}}}*/
void  Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/

	//   The penalty is stable if it doesn't change during two consecutive iterations.   
	int        unstable=0;
	int        new_active;
	int        penalty_lock;
	IssmDouble pressure;
	IssmDouble h;
	IssmDouble h_max;	
	HydrologyDCInefficientAnalysis* inefanalysis = NULL;

	/*check that pengrid is not a clone (penalty to be added only once)*/
	if(node->IsClone()){
		unstable=0;
		*punstable=unstable;
		return;
	}
	if(!element->IsOnBase()){
		unstable=0;
		active=0;
		*punstable=unstable;
		return;
	}

	/*Get sediment water head h*/
	inefanalysis = new HydrologyDCInefficientAnalysis();
	element->GetInputValue(&h,node,SedimentHeadEnum);
	inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
	parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);

	if (h>h_max){
	 new_active=1;
	}
	else{
	 new_active=0;
	}

	if(this->active==new_active){
		unstable=0;
	}
	else{
		unstable=1;
		if(penalty_lock)zigzag_counter++;
	}

	/*If penalty keeps zigzagging more than penalty_lock times: */
	if(penalty_lock){
		if(zigzag_counter>penalty_lock){
			unstable=0;
			active=1;
		}
	}
	/*Set penalty flag*/
	this->active=new_active;

	/*Assign output pointers:*/
	delete inefanalysis;
	*punstable=unstable;
}
/*}}}*/
ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
	IssmDouble    penalty_factor;

	/*Retrieve parameters*/
	parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);

	/*Initialize Element matrix and return if necessary*/
	if(!this->active) return NULL;
	ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);

	Ke->values[0]=kmax*pow(10.,penalty_factor);

	/*Clean up and return*/
	return Ke;
}
/*}}}*/
ElementVector* Pengrid::PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax){/*{{{*/

	IssmDouble h_max;
	IssmDouble penalty_factor;
	HydrologyDCInefficientAnalysis* inefanalysis = NULL;

	/*Initialize Element matrix and return if necessary*/
	if(!this->active) return NULL;
	ElementVector* pe=new ElementVector(&node,1,this->parameters);
	inefanalysis = new HydrologyDCInefficientAnalysis();

	/*Retrieve parameters*/
	parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);

	/*Get h_max and compute penalty*/
	inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);

	pe->values[0]=kmax*pow(10.,penalty_factor)*h_max;

	/*Clean up and return*/
	delete inefanalysis;
	return pe;
}
/*}}}*/
ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/

	IssmDouble moulin_load,dt;

	/*Initialize Element matrix*/
	ElementVector* pe=new ElementVector(&node,1,this->parameters);

	this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
	parameters->FindParam(&dt,TimesteppingTimeStepEnum);

	if(dt!=0.0) pe->values[0]=moulin_load*dt;

	/*Clean up and return*/
	return pe;
 }
/*}}}*/
void  Pengrid::ResetConstraint(void){/*{{{*/
	active         = 0;
	zigzag_counter = 0;
}
/*}}}*/
