#include "./HydrologyShreveAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

/*Model processing*/
int  HydrologyShreveAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
	return 1;
}/*}}}*/
void HydrologyShreveAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/

	/*retrieve some parameters: */
	int  hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);

	/*Now, do we really want Shreve?*/
	if(hydrology_model!=HydrologyshreveEnum) return;

	parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
	parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum));

}/*}}}*/
void HydrologyShreveAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

	/*Fetch data needed: */
	int    hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);

	/*Now, do we really want Shreve?*/
	if(hydrology_model!=HydrologyshreveEnum) return;

	/*Update elements: */
	int counter=0;
	for(int i=0;i<iomodel->numberofelements;i++){
		if(iomodel->my_elements[i]){
			Element* element=(Element*)elements->GetObjectByOffset(counter);
			element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
			counter++;
		}
	}

	iomodel->FetchDataToInput(elements,ThicknessEnum);
	iomodel->FetchDataToInput(elements,SurfaceEnum);
	iomodel->FetchDataToInput(elements,BedEnum);
	iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
	iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
	iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
	iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
	iomodel->FetchDataToInput(elements,WatercolumnEnum);

	elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum);
}/*}}}*/
void HydrologyShreveAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/

	/*Fetch parameters: */
	int  hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);

	/*Now, do we really want Shreve?*/
	if(hydrology_model!=HydrologyshreveEnum) return;

	if(iomodel->meshtype==Mesh3DEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
	::CreateNodes(nodes,iomodel,HydrologyShreveAnalysisEnum,P1Enum);
	iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
}/*}}}*/
void HydrologyShreveAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/

	/*retrieve some parameters: */
	int          hydrology_model;
	iomodel->Constant(&hydrology_model,HydrologyModelEnum);

	if(hydrology_model!=HydrologyshreveEnum) return;

	IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum);

}/*}}}*/
void HydrologyShreveAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
	/*No loads*/
}/*}}}*/

/*Finite Element Analysis*/
ElementMatrix* HydrologyShreveAnalysis::CreateKMatrix(Element* element){/*{{{*/
	_error_("not implemented yet");
}/*}}}*/
ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
_error_("not implemented yet");
}/*}}}*/
void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
	element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
}/*}}}*/
void HydrologyShreveAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

	/*Intermediary*/
	int* doflist = NULL;

	/*Fetch number of nodes for this finite element*/
	int numnodes = element->GetNumberOfNodes();

	/*Fetch dof list and allocate solution vector*/
	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
	IssmDouble* values = xNew<IssmDouble>(numnodes);

	/*Use the dof list to index into the solution vector: */
	for(int i=0;i<numnodes;i++){
		values[i]=solution[doflist[i]];
		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
		if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
	}

	/*Add input to the element: */
	element->AddInput(WatercolumnEnum,values,P1Enum);

	/*Free ressources:*/
	xDelete<IssmDouble>(values);
	xDelete<int>(doflist);
}/*}}}*/