/*!\file: stressbalance_core.cpp
 * \brief: core of the stressbalance solution 
 */ 

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

void stressbalance_core(FemModel* femmodel){

	/*parameters: */
	bool  dakota_analysis;
	int   dim;
	bool  isSIA,isSSA,isL1L2,isHO,isFS;
	bool  conserve_loads    = true;
	bool  save_results;
	int   newton;
	int   solution_type;
	int   numoutputs        = 0;
	int  *requested_outputs = NULL;

	/* recover parameters:*/
	femmodel->parameters->FindParam(&dim,MeshDimensionEnum);
	femmodel->parameters->FindParam(&isSIA,FlowequationIsSIAEnum);
	femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
	femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum);
	femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
	femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
	femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum);
	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum);

	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
	if(dakota_analysis && solution_type==StressbalanceSolutionEnum){
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVxEnum,VxEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVyEnum,VyEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVzEnum,VzEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuPressureEnum,PressureEnum);
	}

	/*Compute slopes: */
	if(isSIA) surfaceslope_core(femmodel);
	if(isFS){
		bedslope_core(femmodel);
		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
		ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
	}

	if(isSIA){
		if(VerboseSolution()) _printf0_("   computing SIA velocities\n");

		/*Take the last velocity into account so that the velocity on the SSA domain is not zero*/
		if(isSSA || isL1L2 || isHO ) ResetBoundaryConditions(femmodel,StressbalanceSIAAnalysisEnum);
		femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
		solutionsequence_linear(femmodel);
		if(isSSA || isL1L2 || isHO) ResetBoundaryConditions(femmodel,StressbalanceAnalysisEnum);
	}

	if ((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor
		if(VerboseSolution()) _printf0_("   computing velocities\n");

		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
		if(newton>0)
		 solutionsequence_newton(femmodel);
		else
		 solutionsequence_nonlinear(femmodel,conserve_loads); 
	}

	if ((isSSA || isL1L2 || isHO) && isFS){
		if(VerboseSolution()) _printf0_("   computing coupling betweem lower order models and full-FS\n");
		solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads);
	}

	if (dim==3 & (isSIA || isSSA || isL1L2 || isHO)){
		if(VerboseSolution()) _printf0_("   computing vertical velocities\n");
		femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
		solutionsequence_linear(femmodel);
	}

	if(save_results){
		if(VerboseSolution()) _printf0_("   saving results\n");
		InputToResultx(femmodel,VxEnum);
		InputToResultx(femmodel,VyEnum);
		InputToResultx(femmodel,VelEnum);
		InputToResultx(femmodel,PressureEnum);
		if(dim==3) InputToResultx(femmodel,VzEnum);
		femmodel->RequestedOutputsx(requested_outputs,numoutputs);
	}

	if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx();

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