/*!\file:  adjointdiagnostic_core.cpp
 * \brief compute inverse method adjoint state
 */ 

#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 adjointdiagnostic_core(FemModel* femmodel){

	/*parameters: */
	bool isstokes;
	bool save_results;
	bool conserve_loads   = true;

	/*retrieve parameters:*/
	femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum);
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);

	/*Compute velocities*/
	if(VerboseSolution()) _pprintLine_("   computing velocities");
	femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
	solver_nonlinear(femmodel,conserve_loads); 

	/*Call SurfaceAreax, because some it might be needed by PVector*/
	SurfaceAreax(NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);

	/*Compute adjoint*/
	if(VerboseSolution()) _pprintLine_("   computing adjoint");
	femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
	solver_adjoint_linear(femmodel);

	/*Save results*/
	if(save_results){
		if(VerboseSolution()) _pprintLine_("   saving results");
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointxEnum);
		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointyEnum);
		if (isstokes){
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointpEnum);
		}
	}
}
