/*
 * CreateNodesStressbalance.c:
 */

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

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

	/*Intermediary*/
	bool isSSA,isL1L2,isHO,isFS,iscoupling;
	int  finiteelement=-1,approximation=-1;

	/*Fetch parameters: */
	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
	iomodel->Constant(&isHO,FlowequationIsHOEnum);
	iomodel->Constant(&isFS,FlowequationIsFSEnum);

	/*Now, check that we have non SIA elements */
	if(!isSSA & !isL1L2 & !isHO & !isFS) return;

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

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

		/*Get finite element type*/
		if(isSSA){
			approximation=SSAApproximationEnum;
			iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
		}
		else if(isL1L2){
			approximation = L1L2ApproximationEnum;
			finiteelement = P1Enum;
		}
		else if(isHO){
			approximation = HOApproximationEnum;
			iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
		}
		else if(isFS){
			approximation = FSApproximationEnum;
			iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
		}
		iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
		CreateNodes(pnodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation);
		iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
	}
	else{
		/*Coupling: we are going to create P1 Elements only*/

		/*Recover nodes*/
		Nodes* nodes = *pnodes;
		Node*  node  = NULL;
		int    lid=0;
		if(!nodes) nodes = new Nodes();

		iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
		if(isFS){
			/*P1+ velocity*/
			for(int i=0;i<iomodel->numberofvertices;i++){
				if(iomodel->my_vertices[i]){
					approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
					if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation));
				}
			}
			for(int i=0;i<iomodel->numberofelements;i++){
				if(iomodel->my_elements[i]){
					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
					node->Deactivate();
					nodes->AddObject(node);
				}
			}
			/*P1 pressure*/
			for(int i=0;i<iomodel->numberofvertices;i++){
				if(iomodel->my_vertices[i]){
					approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
					if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
						node->Deactivate();
					}
					nodes->AddObject(node);
				}
			}
		}
		else{
			for(int i=0;i<iomodel->numberofvertices;i++){
				if(iomodel->my_vertices[i]){
					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])));
				}
			}
		}
		iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);

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