/*
 * CreateConstraintsStressbalanceHoriz.c:
 */

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

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

	/*Intermediary*/
	int count;
	IssmDouble yts;

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

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

	/*return if 2d mesh*/
	if(iomodel->dim==2) return;

	/*Fetch data: */
	iomodel->FetchData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);

	/*Initialize counter*/
	count=0;

	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
	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(FlowequationBorderFSEnum)[i])){
				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
				count++;
			}
			else if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvzEnum)[i])){
				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,
								iomodel->Data(StressbalanceSpcvzEnum)[i]/yts,StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
				count++;

			}
		} 
	}

	/*Free data: */
	iomodel->DeleteData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);

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