/*!\file: diagnostic_core.cpp
 * \brief: core of the diagnostic solution 
 */ 

#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 "../include/include.h"
#include "../solvers/solvers.h"

void diagnostic_core(FemModel* femmodel){

	/*parameters: */
	bool  dakota_analysis   = false;
	int   dim               = -1;
	bool  ishutter          = false;
	bool  ismacayealpattyn  = false;
	bool  isstokes          = false;
	bool  isnewton          = false;
	bool  conserve_loads    = true;
	bool  modify_loads      = true;
	bool  save_results;
	int   solution_type;
	int   numoutputs        = 0;
	int  *requested_outputs = NULL;

	/* recover parameters:*/
	femmodel->parameters->FindParam(&dim,MeshDimensionEnum);
	femmodel->parameters->FindParam(&ishutter,FlowequationIshutterEnum);
	femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
	femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum);
	femmodel->parameters->FindParam(&isnewton,DiagnosticIsnewtonEnum);
	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
	femmodel->parameters->FindParam(&numoutputs,DiagnosticNumRequestedOutputsEnum);
	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DiagnosticRequestedOutputsEnum);

	/*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==DiagnosticSolutionEnum){
		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(ishutter) surfaceslope_core(femmodel);
	if(isstokes){
		bedslope_core(femmodel);
		femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
		ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
	}

	if(ishutter){
			
		if(VerboseSolution()) _pprintLine_("   computing hutter velocities");

		//Take the last velocity into account so that the velocity on the MacAyeal domain is not zero
		if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum);

		femmodel->SetCurrentConfiguration(DiagnosticHutterAnalysisEnum);
		solver_linear(femmodel);
		
		if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum);
	}

	if (ismacayealpattyn ^ isstokes){ // ^ = xor
		
		if(VerboseSolution()) _pprintLine_("   computing velocities");
		femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
		if(isnewton)
		 solver_newton(femmodel); 
		else
		 solver_nonlinear(femmodel,conserve_loads); 
	}

	if (ismacayealpattyn && isstokes){

		if(VerboseSolution()) _pprintLine_("   computing coupling macayealpattyn and stokes velocities and pressure ");
		solver_stokescoupling_nonlinear(femmodel,conserve_loads);
	}

	if (dim==3 & (ishutter || ismacayealpattyn)){

		if(VerboseSolution()) _pprintLine_("   computing vertical velocities");
		femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
		solver_linear(femmodel);
	}


	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,VelEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum);
		if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum);
		RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs);
	}

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