/*!\file: steadystate_core.cpp
 * \brief: core of the steadystate solution 
 */ 

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

#include "../include/include.h"
#include "../toolkits/toolkits.h"
#include "../classes/objects/objects.h"
#include "../shared/shared.h"
#include "../io/io.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "./solutions.h"
#include "../modules/modules.h"
#include "../solvers/solvers.h"

void steadystate_core(FemModel* femmodel){

	/*intermediary: */
	int step; 

	/*parameters: */
	bool save_results,isenthalpy;
	int  maxiter;
	int  numoutputs         = 0;
	int  *requested_outputs = NULL;

	/* recover parameters:*/
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
	femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum);
	femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
	femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
	femmodel->parameters->SetParam(false,SaveResultsEnum);
	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);

	/*intialize counters: */
	step=1;

	for(;;){

		if(VerboseSolution()) _pprintLine_("   computing temperature and velocity for step: " << step);
		#ifdef _HAVE_THERMAL_
		if(isenthalpy==0){
			thermal_core(femmodel);
		}
		else{
			enthalpy_core(femmodel);
		}
		#else
		_error_("ISSM was not compiled with thermal capabilities. Exiting");
		#endif

		if(VerboseSolution()) _pprintLine_("   computing new velocity");
		diagnostic_core(femmodel);

		if (step>1){
			if(VerboseSolution()) _pprintLine_("   checking velocity, temperature and pressure convergence");
			if(steadystateconvergence(femmodel)) break;
		}
		if(step>maxiter){
			if(VerboseSolution()) _pprintLine_("   maximum number steadystate iterations " << maxiter << " reached");
			break;
		}

		if(VerboseSolution()) _pprintLine_("   saving velocity, temperature and pressure to check for convergence at next step");
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxPicardEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyPicardEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzPicardEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressurePicardEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum);

		//increase counter
		step++;
	}

	if(save_results){
		if(VerboseSolution()) _pprintLine_("   saving results");
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum);
		if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum);
		if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum);
		if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum);
		femmodel->RequestedOutputsx(requested_outputs,numoutputs);
	}

	/*Free ressources:*/
	xDelete<int>(requested_outputs);
}
