/*!\file: gia_core.cpp
 * \brief: core of the GIA solution 
 */ 

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../classes/Inputs/TriaInput.h"
#include "../classes/Inputs/DatasetInput.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"
void gia_core(FemModel* femmodel){

	Vector<IssmDouble> *wg    = NULL;
	Vector<IssmDouble> *dwdtg = NULL;
	IssmDouble          *x    = NULL;
	IssmDouble          *y    = NULL;
	int nummodels,giamodel;
	IssmDouble modelid; 

	/*parameters: */
	bool save_results;
	int  gsize;
	
	/*Start profiler*/
	femmodel->profiler->Start(GIACORE);

	/*Recover some parameters: */
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);

	if(VerboseSolution()) _printf0_("   computing GIA\n");

	/*Call on core computations: */
	femmodel->SetCurrentConfiguration(GiaAnalysisEnum);

	/*Figure out type of gia model:*/
	femmodel->parameters->FindParam(&giamodel,GiaModelEnum);

	if(giamodel==1){ //GiaIvins
		/*Figure out size of g-set deflection vector and allocate solution vector: */
		gsize      = femmodel->nodes->NumberOfDofs(GsetEnum);
		wg = new Vector<IssmDouble>(gsize);
		dwdtg = new Vector<IssmDouble>(gsize);

		/*first, recover x and y vectors from vertices: */
		VertexCoordinatesx(&x,&y,NULL,femmodel->vertices); //no need for z coordinate

		/*call the main module: */
		femmodel->Deflection(wg,dwdtg,x,y);

		/*assemble vector: */
		wg->Assemble();
		dwdtg->Assemble();

		InputUpdateFromVectorx(femmodel,wg,UGiaEnum,VertexSIdEnum);
		InputUpdateFromVectorx(femmodel,dwdtg,UGiaRateEnum,VertexSIdEnum);

		if(save_results){
			int outputs[2] = {UGiaEnum,UGiaRateEnum};
			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
		}

		xDelete<IssmDouble>(x);
		xDelete<IssmDouble>(y);
	}
	else if(giamodel==2){ //GiaCaron
		/*not implemneted yet*/
	}
	else if(giamodel==3){ //GiaMme
		
		/*Go grab the values for Ngia and Ugia in the offline GiaMme class:*/ 
		femmodel->parameters->FindParam(&modelid,GiaModelidEnum); 
		modelid--; //from matlab

		/*find the Ngia and Ugia dataset:*/
		DatasetInput* dataset_input_ngia=femmodel->inputs->GetDatasetInput(GiaMmeNgiaEnum);
		DatasetInput* dataset_input_ugia=femmodel->inputs->GetDatasetInput(GiaMmeUgiaEnum);
		/*Go find the modelid'th input:*/
		TriaInput* tria_input_ngia=dataset_input_ngia->GetTriaInputByOffset(modelid);
		TriaInput* tria_input_ugia=dataset_input_ugia->GetTriaInputByOffset(modelid);

		/*Plug into SealevelUGiaRate and SealevelNGiaRate inputs:*/
		Input* tria_input_copy_ngia=tria_input_ngia->copy();
		Input* tria_input_copy_ugia=tria_input_ugia->copy();
		tria_input_copy_ngia->ChangeEnum(NGiaRateEnum);
		tria_input_copy_ugia->ChangeEnum(UGiaRateEnum);
		femmodel->inputs->AddInput(tria_input_copy_ngia);
		femmodel->inputs->AddInput(tria_input_copy_ugia);
	}
	
	/*End profiler*/
	femmodel->profiler->Stop(GIACORE);
}
