/*
 * CreateConstraintsStressbalanceSIA.c:
 */

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

void	CreateConstraintsStressbalanceSIA(Constraints** pconstraints, IoModel* iomodel){

	/*Intermediary*/
	int        count;
	IssmDouble yts;
	bool       isSIA;

	/*Output*/
	Constraints* constraints = NULL;
	SpcStatic*    spcstatic  = NULL;

	/*Recover pointer: */
	constraints=*pconstraints;

	/*Fetch parameters: */
	iomodel->Constant(&yts,ConstantsYtsEnum);
	iomodel->Constant(&isSIA,FlowequationIsSIAEnum);

	/*Now, is the flag isSIA on? otherwise, do nothing: */
	if (!isSIA) return;

	/*Fetch data: */
	iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);

	/*Initialize conunter*/
	count=0;

	/*vx and vy are spc'd if we are not on nodeonSIA: */
	for(int i=0;i<iomodel->numberofvertices;i++){
		/*keep only this partition's nodes:*/
		if((iomodel->my_vertices[i])){
			if (!reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])==SIAApproximationEnum){

				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
				count++;

				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceSIAAnalysisEnum));
				count++;
			}
			else{
				if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvxEnum)[i]/yts,StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
					count++;
				}

				if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,iomodel->Data(StressbalanceSpcvyEnum)[i]/yts,StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
					count++;
				}
			}
		}
	}

	/*Free data: */
	iomodel->DeleteData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);

	/*Assign output pointer: */
	*pconstraints=constraints;
}
