/*
 * CreateNodesBalancethickness.c:
 */

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

void	CreateNodesBalancethickness(Nodes** pnodes, IoModel* iomodel){

	/*Intermediary*/
	int  i,j;
	int  node_id;
	int  vertex_id;
	int  io_index;
	bool continuous_galerkin=true;
	int    dim;
	int    numberofelements;
	int    numberofvertices;
	int    stabilization;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum);

	/*Recover pointer: */
	Nodes* nodes=*pnodes;

	/*Create nodes if they do not exist yet*/
	if(!nodes) nodes = new Nodes();

	/*Continuous Galerkin partition of nodes: */
	if(stabilization==3) continuous_galerkin=false;
	NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin);

	/*Check in 3d*/
	if(stabilization==3 && dim==3) _error_("DG 3d not implemented yet");

	/*First fetch data: */
	iomodel->FetchData(7,MeshElementsEnum,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum);

	if(continuous_galerkin){

		/*Build Nodes dataset (Continuous Galerkin)*/
		for (i=0;i<numberofvertices;i++){
			if(iomodel->my_vertices[i]){

				/*Add node to nodes dataset: */
				nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,BalancethicknessAnalysisEnum));

			}
		}
	}
	else{

		/*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/
		for (i=0;i<numberofelements;i++){
			for (j=0;j<3;j++){

				if(iomodel->my_nodes[3*i+j]){ 

					//Get index of the vertex on which the current node is located
					vertex_id=reCast<int,IssmDouble>(*(iomodel->Data(MeshElementsEnum)+3*i+j)); //(Matlab indexing)
					io_index=vertex_id-1;                      //(C indexing)
					_assert_(vertex_id>0 && vertex_id<=numberofvertices);

					//Compute Node id
					node_id=iomodel->nodecounter+3*i+j+1;

					/*Add node to nodes dataset: */
					nodes->AddObject(new Node(node_id,node_id-1,vertex_id,io_index,iomodel,BalancethicknessAnalysisEnum));

				}
			}
		}
	}

	/*Clean fetched data: */
	iomodel->DeleteData(7,MeshElementsEnum,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum);

	/*Assign output pointer: */
	*pnodes=nodes;
}
