/*
 * CreateConstraintsStressbalance.c:
 */

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

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

	/*Intermediary*/
	int        i,j;
	int        count,finiteelement;
	IssmDouble g;
	IssmDouble rho_ice;
	IssmDouble FSreconditioning;
	bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
	bool       spcpresent = false;
	int        Mx,Nx;
	int        My,Ny;
	int        Mz,Nz;
	IssmDouble *spcvx          = NULL;
	IssmDouble *spcvy          = NULL;
	IssmDouble *spcvz          = NULL;
	IssmDouble *nodeonSSA = NULL;
	IssmDouble *nodeonHO   = NULL;
	IssmDouble *nodeonFS   = NULL;
	IssmDouble *nodeonbed      = NULL;
	IssmDouble *groundedice_ls = NULL;
	IssmDouble *vertices_type  = NULL;
	IssmDouble *surface        = NULL;
	IssmDouble *z              = NULL;
	IssmDouble *timesx=NULL;
	IssmDouble *timesy=NULL;
	IssmDouble *timesz=NULL;
   IssmDouble* values=NULL;

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

	/*Fetch parameters: */
	iomodel->Constant(&g,ConstantsGEnum);
	iomodel->Constant(&rho_ice,MaterialsRhoIceEnum);
	iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum);
	iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
	iomodel->Constant(&isHO,FlowequationIsHOEnum);
	iomodel->Constant(&isFS,FlowequationIsFSEnum);

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

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

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

	/*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
	if(!iscoupling){

		/*Get finite element type*/
		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);
			/*Deduce velocity interpolation from finite element*/
			switch(finiteelement){
				case P1P1Enum          : finiteelement = P1Enum;       break;
				case P1P1GLSEnum       : finiteelement = P1Enum;       break;
				case MINIcondensedEnum : finiteelement = P1bubbleEnum; break;
				case MINIEnum          : finiteelement = P1bubbleEnum; break;
				case TaylorHoodEnum    : finiteelement = P2Enum;       break;
				default: _error_("finite element "<<finiteelement<<" not supported");
			}
		}

		IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
		IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);

		if(isFS){

			/*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
			iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
			iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
			iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
			iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
			iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
			for(i=0;i<iomodel->numberofvertices;i++){
				if(iomodel->my_vertices[i]){
					if(nodeonbed[i]>0. && groundedice_ls[i]>0. && nodeonFS[i]>0.){
						if(vertices_type[i] == FSApproximationEnum){
							for(j=0;j<Nz;j++) spcvz[i*Nz+j] = 0.;
						}
						else{
							_error_("not supported");
						}
					}
				}
			}
			IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,3);
			iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
			iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
			iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
			iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
			iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);

			/*Pressure spc*/
			count = constraints->Size();
			iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
			iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
			iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
			switch(finiteelement){
				case P1bubbleEnum:
					for(i=0;i<iomodel->numberofvertices;i++){
						if(iomodel->my_vertices[i]){
							if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
								constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
								count++;
							}
						}
					}
					break;
				case P2Enum:
					for(i=0;i<iomodel->numberofvertices;i++){
						if(iomodel->my_vertices[i]){
							if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
								constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
								count++;
							}
						}
					}
					break;
				default:
					_error_("not implemented yet");
			}
			iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
			iomodel->DeleteData(surface,SurfaceEnum);
			iomodel->DeleteData(z,MeshZEnum);
		}

		*pconstraints=constraints;
		return;
	}

	/*Constraints: fetch data: */
	iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum);
	iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum);
	iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
	iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum);
	if(iomodel->dim==3)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum);
	if(iomodel->dim==3)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
	if(iomodel->dim==3)iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
	if(iomodel->dim==3)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
	iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
	iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
	iomodel->FetchData(&z,NULL,NULL,MeshZEnum);

	/*Initialize counter: */
	count=0;

	/*figure out times: */
	timesx=xNew<IssmDouble>(Nx);
	for(j=0;j<Nx;j++){
		timesx[j]=spcvx[(Mx-1)*Nx+j];
	}
	/*figure out times: */
	timesy=xNew<IssmDouble>(Ny);
	for(j=0;j<Ny;j++){
		timesy[j]=spcvy[(My-1)*Ny+j];
	}
	/*figure out times: */
	timesz=xNew<IssmDouble>(Nz);
	for(j=0;j<Nz;j++){
		timesz[j]=spcvz[(Mz-1)*Nz+j];
	}

	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
	for(i=0;i<iomodel->numberofvertices;i++){
		if(iomodel->my_vertices[i]){

			/*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/
			if(reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
				/*If grionSSA, spc HO dofs: 3 & 4*/
					if (reCast<int,IssmDouble>(nodeonHO[i])){
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						if (!xIsNan<IssmDouble>(spcvx[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvy[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}

					}
					else if (reCast<int,IssmDouble>(nodeonSSA[i])){
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						if (!xIsNan<IssmDouble>(spcvx[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvy[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}

					}
					else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA");
			}
			/*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
			else if (reCast<int,IssmDouble>(vertices_type[i])==HOFSApproximationEnum){
				/*If grion,HO spc FS dofs: 3 4 & 5*/
					if (reCast<int,IssmDouble>(nodeonHO[i])){
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						if (!xIsNan<IssmDouble>(spcvx[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvy[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}

					}
					else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc HO nodes: 1 & 2
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						if (!xIsNan<IssmDouble>(spcvx[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvy[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvz[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
					}
					else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS");
			}
			/*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
			else if (reCast<int,IssmDouble>(vertices_type[i])==SSAFSApproximationEnum){
				/*If grion,HO spc FS dofs: 3 4 & 5*/
					if (reCast<int,IssmDouble>(nodeonSSA[i])){
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						if (!xIsNan<IssmDouble>(spcvx[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvy[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}

					}
					else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc SSA nodes: 1 & 2
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
						count++;
						if (!xIsNan<IssmDouble>(spcvx[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvy[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
						if (!xIsNan<IssmDouble>(spcvz[i])){
							constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
							count++;
						}
					}
					else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS");
			}
			/*Now add the regular spcs*/
			else{
				if (Mx==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvx[i])){
					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
					count++;

				}
				else if (Mx==iomodel->numberofvertices+1) {
					/*figure out times and values: */
					values=xNew<IssmDouble>(Nx);
					spcpresent=false;
					for(j=0;j<Nx;j++){
						values[j]=spcvx[i*Nx+j];
						if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
					}

					if(spcpresent){
						constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Nx,timesx,values,StressbalanceAnalysisEnum));
						count++;
					}
					xDelete<IssmDouble>(values);
				}
				else if (vertices_type[i]==SIAApproximationEnum){
					constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum));
					count++;
				}

				if (My==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvy[i])){
					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy.
					count++;
				}
				else if (My==iomodel->numberofvertices+1){
					/*figure out times and values: */
					values=xNew<IssmDouble>(Ny);
					spcpresent=false;
					for(j=0;j<Ny;j++){
						values[j]=spcvy[i*Ny+j];
						if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
					}
					if(spcpresent){
						constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Ny,timesy,values,StressbalanceAnalysisEnum));
						count++;
					}
					xDelete<IssmDouble>(values);
				}
				else if (vertices_type[i]==SIAApproximationEnum){
					constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,StressbalanceAnalysisEnum));
					count++;
				}

				if (reCast<int,IssmDouble>(vertices_type[i])==FSApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum)){
					if (Mz==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvz[i])){
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
						count++;
					}
					else if (Mz==iomodel->numberofvertices+1){
						/*figure out times and values: */
						values=xNew<IssmDouble>(Nz);
						spcpresent=false;
						for(j=0;j<Nz;j++){
							values[j]=spcvz[i*Nz+j];
							if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
						}
						if(spcpresent){
							constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,Nz,timesz,values,StressbalanceAnalysisEnum));
							count++;
						}
						xDelete<IssmDouble>(values);
					}

				}
				if (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
					count++;
				}
			}

			/*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
			if (iomodel->dim==3) if(nodeonbed[i]>0. && groundedice_ls[i]>=0. && nodeonFS[i]>0.){
				 switch(reCast<int,IssmDouble>(vertices_type[i])){
					case SSAFSApproximationEnum:
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0.,StressbalanceAnalysisEnum));
						count++;
						break;
					case HOFSApproximationEnum:
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,0.,StressbalanceAnalysisEnum));
						count++;
						break;
					case FSApproximationEnum:
						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0.,StressbalanceAnalysisEnum));
						count++;
						break;
					default: _error_("Vertex approximation " << EnumToStringx(reCast<int,IssmDouble>(vertices_type[i])) << " not supported");
				}
			}
		}
	}

	/*Free data: */
	iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum);
	iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum);
	iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
	iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum);
	if(iomodel->dim==3)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum);
	if(iomodel->dim==3)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
	if(iomodel->dim==3)iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
	if(iomodel->dim==3)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
	iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
	iomodel->DeleteData(surface,SurfaceEnum);
	iomodel->DeleteData(z,MeshZEnum);

	/*Free resources:*/
	xDelete<IssmDouble>(timesx);
	xDelete<IssmDouble>(timesy);
	xDelete<IssmDouble>(timesz);
	xDelete<IssmDouble>(values);

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