/*
 * UpdateElementsDiagnosticHoriz:
 */

#include "../../../Container/Container.h"
#include "../../../modules/modules.h"
#include "../../../toolkits/toolkits.h"
#include "../../../io/io.h"
#include "../../../EnumDefinitions/EnumDefinitions.h"
#include "../../../objects/objects.h"
#include "../../../shared/shared.h"
#include "../../MeshPartitionx/MeshPartitionx.h"
#include "../../../include/include.h"
#include "../ModelProcessorx.h"

void	UpdateElementsDiagnosticHoriz(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){

	int    dim;
	int    numberofelements;
	bool   ismacayealpattyn;
	bool   isstokes;
	bool   control_analysis;
	bool   dakota_analysis;

	/*Fetch constants needed: */
	iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
	iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);

	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
	if(!ismacayealpattyn & !isstokes) return;

	/*Fetch data needed: */
	iomodel->FetchData(2,MeshElementsEnum,FlowequationElementEquationEnum);

	/*Update elements: */
	int counter=0;
	for(int i=0;i<numberofelements;i++){
		if(iomodel->my_elements[i]){
			Element* element=(Element*)elements->GetObjectByOffset(counter);
			element->Update(i,iomodel,analysis_counter,analysis_type); //we need i to index into elements.
			counter++;
		}
	}

	/*Create inputs: */
	iomodel->FetchDataToInput(elements,ThicknessEnum);
	iomodel->FetchDataToInput(elements,GeometryHydrostaticRatioEnum);
	iomodel->FetchDataToInput(elements,SurfaceEnum);
	iomodel->FetchDataToInput(elements,BedEnum);
	iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
	iomodel->FetchDataToInput(elements,FrictionPEnum);
	iomodel->FetchDataToInput(elements,FrictionQEnum);
	iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum);
	iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
	iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
	iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
	iomodel->FetchDataToInput(elements,VxEnum);
	iomodel->FetchDataToInput(elements,VyEnum);

	if (dim==3){
		iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
		iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
		iomodel->FetchDataToInput(elements,VzEnum);
		iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
		iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
		iomodel->FetchDataToInput(elements,FlowequationBorderstokesEnum);
	}

	if(control_analysis){
		iomodel->FetchDataToInput(elements,InversionVxObsEnum);
		iomodel->FetchDataToInput(elements,InversionVyObsEnum);
	}

	/*Post-process: */
	elements->InputDuplicate(VxEnum,VxPicardEnum);
	elements->InputDuplicate(VxEnum,InversionVxObsEnum);
	if(dakota_analysis)elements->InputDuplicate(VxEnum,QmuVxEnum);
	
	elements->InputDuplicate(VyEnum,VyPicardEnum);
	elements->InputDuplicate(VyEnum,InversionVyObsEnum);
	if(dakota_analysis)elements->InputDuplicate(VyEnum,QmuVyEnum);
	
	if(dim==3){
		elements->InputDuplicate(VzEnum,VzPicardEnum);
		elements->InputDuplicate(VzEnum,InversionVzObsEnum);
		if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
	}
	
	/*Free data: */
	iomodel->DeleteData(2,MeshElementsEnum,FlowequationElementEquationEnum);
}
