/*
 * UpdateElementsStressbalance:
 */
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../../../toolkits/toolkits.h"
#include "../../../classes/classes.h"
#include "../../../shared/shared.h"
#include "../../MeshPartitionx/MeshPartitionx.h"
#include "../ModelProcessorx.h"

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

	int    materials_type,finiteelement;
	int    approximation;
	int*   finiteelement_list=NULL;
	bool   isSSA,isL1L2,isHO,isFS,iscoupling;
	bool   control_analysis;
	bool   dakota_analysis;

	/*Fetch constants needed: */
	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
	iomodel->Constant(&isHO,FlowequationIsHOEnum);
	iomodel->Constant(&isFS,FlowequationIsFSEnum);
	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
	iomodel->Constant(&materials_type,MaterialsEnum);

	/*Now, is the flag macayaealHO on? otherwise, do nothing: */
	if(!isSSA & !isL1L2 & !isHO & !isFS) return;

	/*Fetch data needed and allocate vectors: */
	iomodel->FetchData(1,FlowequationElementEquationEnum);
	finiteelement_list=xNewZeroInit<int>(iomodel->numberofelements);


	/*Do we have coupling*/
	if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
	 iscoupling = true;
	else
	 iscoupling = false;

	/*Get finite element type*/
	if(!iscoupling){
		if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
		else if(isL1L2) finiteelement = P1Enum;
		else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
		else if(isFS)   iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
		for(int i=0;i<iomodel->numberofelements;i++){
			finiteelement_list[i]=finiteelement;
		}
	}
	else{
		if(isFS){
			for(int i=0;i<iomodel->numberofelements;i++){
				approximation=reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[i]);
				if(approximation==FSApproximationEnum || approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
					finiteelement_list[i]=MINIcondensedEnum;
				}
				else{
					finiteelement_list[i]=P1Enum;
				}
			}
		}
		else{
			finiteelement = P1Enum;
			for(int i=0;i<iomodel->numberofelements;i++){
				finiteelement_list[i]=finiteelement;
			}
		}
	}

	/*Update elements: */
	int counter=0;
	for(int i=0;i<iomodel->numberofelements;i++){
		if(iomodel->my_elements[i]){
			Element* element=(Element*)elements->GetObjectByOffset(counter);
			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
			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,MaskIceLevelsetEnum);
	iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
	iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
	iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
	iomodel->FetchDataToInput(elements,VxEnum);
	iomodel->FetchDataToInput(elements,VyEnum);
	iomodel->FetchDataToInput(elements,LoadingforceXEnum);
	iomodel->FetchDataToInput(elements,LoadingforceYEnum);
	if(materials_type==MatdamageiceEnum){
		iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum);
	}
	if(iomodel->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,FlowequationBorderFSEnum);
		iomodel->FetchDataToInput(elements,LoadingforceZEnum);
	}

	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);

	#ifdef _HAVE_ANDROID_
	elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum);
	#endif

	elements->InputDuplicate(VyEnum,VyPicardEnum);
	elements->InputDuplicate(VyEnum,InversionVyObsEnum);
	if(dakota_analysis)elements->InputDuplicate(VyEnum,QmuVyEnum);

	if(iomodel->dim==3){
		elements->InputDuplicate(VzEnum,VzPicardEnum);
		elements->InputDuplicate(VzEnum,InversionVzObsEnum);
		if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
	}

	/*Free data: */
	iomodel->DeleteData(1,FlowequationElementEquationEnum);
	xDelete<int>(finiteelement_list);
}
