/*!\file: solver_thermal_nonlinear.cpp
 * \brief: core of the thermal solution 
 */ 

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

void solver_thermal_nonlinear(FemModel* femmodel){

	/*solution : */
	Vector<IssmDouble>* tg=NULL; 
	Vector<IssmDouble>* tf=NULL; 
	Vector<IssmDouble>* tf_old=NULL; 
	Vector<IssmDouble>* ys=NULL; 
	IssmDouble melting_offset;

	/*intermediary: */
	Matrix<IssmDouble>* Kff=NULL;
	Matrix<IssmDouble>* Kfs=NULL;
	Vector<IssmDouble>* pf=NULL;
	Vector<IssmDouble>* df=NULL;

	bool converged;
	int constraints_converged;
	int num_unstable_constraints;
	int count;
	int thermal_penalty_threshold;
	int thermal_maxiter;

	/*parameters:*/
	int kflag,pflag;
	bool lowmem=0;
	int  configuration_type;

	/*Recover parameters: */
	kflag=1; pflag=1;
	femmodel->parameters->FindParam(&lowmem,SettingsLowmemEnum);
	femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
	femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);

	count=1;
	converged=false;

	if(VerboseSolution()) _pprintLine_("starting direct shooting method");
	InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum);
	InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum);
	femmodel->UpdateConstraintsx();

	for(;;){

		femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset);
		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
		Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf);
		Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
		xdelete(&tf_old); tf_old=tf->Duplicate();
		xdelete(&Kff);xdelete(&pf);xdelete(&tg); xdelete(&df);
		Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys);
		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,tg);

		ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);

		if (!converged){
			if(VerboseConvergence()) _pprintLine_("   #unstable constraints = " << num_unstable_constraints);
			if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
			if (count>=thermal_maxiter){
				converged=true;
				_pprintLine_("   maximum number of iterations (" << thermal_maxiter << ") exceeded"); 
			}
		}
		count++;

		InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);

		if(converged)break;
	}

	InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,tg);
	InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,melting_offset,MeltingOffsetEnum);

	/*Free ressources: */
	xdelete(&tg);
	xdelete(&tf);
	xdelete(&tf_old);
	xdelete(&ys);
}
