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

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

void solutionsequence_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  configuration_type;

	/*Recover parameters: */
	femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
	femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);

	count=1;
	converged=false;

	InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
	InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
	femmodel->UpdateConstraintsx();

	for(;;){

		delete tf_old; tf_old=tf;
		SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
		Reduceloadx(pf, Kfs, ys); delete Kfs;
		Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
		delete Kff;delete pf;delete tg; delete df;
		Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys;
		InputUpdateFromSolutionx(femmodel,tg);

		ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);

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

		InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);

		if(converged)break;
	}

	InputUpdateFromSolutionx(femmodel,tg);
	femmodel->parameters->SetParam(melting_offset,MeltingOffsetEnum);

	/*Free ressources: */
	delete tg;
	delete tf;
	delete tf_old;
}
