/*
 * cielothermal_core.c:
 */

#include "../include/cielo.h"
#include "../modules.h"
#include "./parallel.h"

#undef __FUNCT__ 
#define __FUNCT__ "cielothermal_core"

#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)

int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel){

	
	
	/*Femmodel: */
	BatchParams* batchparams=NULL;
	Mat* G_mn=NULL;
	Vec* y_s=NULL;
	int flag_y_s;
	DataSet* bgpdt=NULL;
	DataSet* bgpdtb=NULL;
	DataSet* est=NULL;
	DataSet* lst=NULL;
	DataSet* new_lst=NULL;
	DataSet* ept=NULL;
	DataSet* mpt=NULL;
	DataSet* geom3=NULL;
	
	int  noerr=1;

	int  dummy;
	Vec* t_f=NULL;
	Vec* t_g=NULL;

	Mat* K_gg=NULL;
	Mat* K_gg_initial=NULL;
	Mat* KT_gg=NULL;
	Vec* p_g=NULL;
	Vec* p_g_initial=NULL;
	Mat* K_ff=NULL;
	Mat* K_fs=NULL;
	Vec* p_f=NULL;

	int kflag,ktflag,pflag;
	int converged=0;
	int count;

	/*some parameters:*/
	double sparsity;
	int    connectivity;
	int    analysis_type_enum;
	char*  solverstring=NULL;
	int    debug;
	int    min_thermal_constraints;
	
	/*intermediary data: */
	int    num_unstable_constraints;

	/*Recover model parameters: */
	batchparams=femmodel->batchparams;
	G_mn=femmodel->G_mn;
	y_s=femmodel->y_s;
	flag_y_s=femmodel->flag_y_s;
	bgpdt=femmodel->bgpdt;
	bgpdtb=femmodel->bgpdtb;
	est=femmodel->est;
	lst=femmodel->lst;
	ept=femmodel->ept;
	mpt=femmodel->mpt;
	geom3=femmodel->geom3;
	uset=femmodel->uset; //external variable

	sparsity=batchparams->sparsity;
	connectivity=batchparams->connectivity;
	analysis_type_enum=batchparams->analysis_type_enum;
	solverstring=batchparams->solverstring;
	debug=batchparams->debug;
	min_thermal_constraints=batchparams->min_thermal_constraints;

	
	//Start iteration on non-linearity
	kflag=1; pflag=1; ktflag=0; //stiffness and load generation only:
	converged=0;  //flag to break loop on non-linearity
	count=1;
		
	for(;;){
		if(count==1){
			//*Generate system matrices
			Emgx(&K_gg, &p_g, &KT_gg,kflag,pflag,ktflag,sparsity,connectivity,bgpdt, bgpdtb, est, lst, ept, mpt, geom3, inputs, analysis_type_enum);
		
			//Keep initial K_gg, so that penalties don't keep piling up onto K_gg: */
			K_gg_initial=xmalloc(sizeof(Mat));
			MatDuplicate(*K_gg,MAT_COPY_VALUES,K_gg_initial);
			
			p_g_initial=xmalloc(sizeof(Vec));
			VecDuplicate(*p_g,p_g_initial);VecCopy(*p_g,*p_g_initial);
		}
		else{
			/*Copy K_gg_initial into K_gg, same for p_g: */
			K_gg=xmalloc(sizeof(Mat));
			MatDuplicate(*K_gg_initial,MAT_COPY_VALUES,K_gg);
	
			p_g=xmalloc(sizeof(Vec));
			VecDuplicate(*p_g_initial,p_g);VecCopy(*p_g_initial,*p_g);
		}

		//Add penalties
		PenaltyEmgx(K_gg, p_g, kflag,pflag,bgpdt, bgpdtb, est, lst, ept, mpt, geom3, inputs, analysis_type_enum);

		//Reduce tangent matrix from g size to f size
		Reducematrixfromgtof(&K_ff,&dummy,&dummy, &K_fs,&dummy,&dummy, K_gg, G_mn, flag_y_s, 
				uset->pv_m,uset->msize, uset->pv_n,uset->nsize, uset->pv_f,uset->fsize, uset->pv_s,uset->ssize,uset->gsize,uset->msize,uset->ssize);

		//no need for K_gg and KT_gg anymore
		MatFree(&K_gg); MatFree(&KT_gg);

		//Reduce load from g size to f size
		Reducerightside(&p_f, p_g, G_mn, K_fs, y_s, flag_y_s,
				     uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, 
					 uset->msize,uset->ssize,uset->fsize);

		//no need for p_g and K_fs anymore 
		VecFree(&p_g); MatFree(&K_fs);

		//solve
		if(uset->fsize>0){
			t_f=xmalloc(sizeof(Vec));
			Solverx(t_f,&dummy,K_ff,dummy,dummy,p_f,dummy,NULL,0,solverstring);
		}
		else{
			_printf_("All dof are constrained, f_set is empty...\n");
			goto cleanup_and_return;
		}

		//no need for K_ff and p_f anymore
		MatFree(&K_ff);VecFree(&p_f);

		//Merge back to g set
		Mergesolvec( &t_g, t_f, G_mn, y_s, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->pv_s, uset->ssize,
				uset->gsize,uset->msize,uset->nsize);

		//t_f not  needed anymore
		VecFree(&t_f);

		//Deal with penalty loads
		ParameterInputsAddFromVec(inputs,t_g,"temperature");
		PenaltyConstraintsx(&new_lst, &converged,&num_unstable_constraints, bgpdt, bgpdtb, est, lst, mpt, inputs,analysis_type_enum);
		
		//Figure out if convergence is reached.
		if(!converged){
			if (debug){
				_printf_("%s %i\n","   #unstable constraints ",num_unstable_constraints);
				if (num_unstable_constraints<=min_thermal_constraints){
					converged=1;
				}
			}
		}

		/*Reset lst to new_lst:*/
		DeleteDataSet(&lst); lst=new_lst;

		count++;
		if(converged==1)break;
	}

	cleanup_and_return:

	/*Assign output pointers: */
	*pt_g=t_g;
	return noerr;
}
#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)

