/*!\file Node.c
 * \brief: implementation of the Node object
 */

/*Include files: {{{1*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#include "./objects.h"
#include "../Container/Container.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../shared/shared.h"
#include "../include/include.h"
#include "../modules/modules.h"
/*}}}*/

/*Node constructors and destructors:*/
/*FUNCTION Node::Node() default constructor {{{1*/
Node::Node(){
		 this->inputs=NULL;
		 this->hvertex=NULL;
		 return;
}
/*}}}*/
/*FUNCTION Node::Node(int node_id,int node_sid,int vertex_id,int io_index, IoModel* iomodel,int analysis_type) {{{1*/
Node::Node(int node_id,int node_sid,int vertex_id,int io_index, IoModel* iomodel,int analysis_type){

	/*Intermediary*/
	int k,l;
	int gsize;
	int dim;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);

	/*id: */
	this->id=node_id; 
	this->sid=node_sid; 
	this->analysis_type=analysis_type;

	/*Initialize coord_system: Identity matrix by default*/
	for(k=0;k<3;k++) for(l=0;l<3;l++) this->coord_system[k][l]=0.0;
	for(k=0;k<3;k++) this->coord_system[k][k]=1.0;

	/*indexing:*/
	DistributeNumDofs(&this->indexing,analysis_type,iomodel->Data(FlowequationVertexEquationEnum)+io_index); //number of dofs per node
	gsize=this->indexing.gsize;

	/*Hooks*/
	this->hvertex=new Hook(&vertex_id,1); //node id is the same as the vertex id, continuous galerkin!

	//intialize inputs, and add as many inputs per element as requested: 
	this->inputs=new Inputs();
	if (iomodel->Data(MeshVertexonbedEnum))
	 this->inputs->AddInput(new BoolInput(MeshVertexonbedEnum,(IssmBool)iomodel->Data(MeshVertexonbedEnum)[io_index]));
	if (iomodel->Data(MeshVertexonsurfaceEnum))
	 this->inputs->AddInput(new BoolInput(MeshVertexonsurfaceEnum,(IssmBool)iomodel->Data(MeshVertexonsurfaceEnum)[io_index]));
	if (iomodel->Data(MaskVertexonfloatingiceEnum))
	 this->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,(IssmBool)iomodel->Data(MaskVertexonfloatingiceEnum)[io_index]));
	if (iomodel->Data(MaskVertexongroundediceEnum))
	 this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,(IssmBool)iomodel->Data(MaskVertexongroundediceEnum)[io_index]));
	if (analysis_type==DiagnosticHorizAnalysisEnum)
	 this->inputs->AddInput(new IntInput(ApproximationEnum,(IssmInt)iomodel->Data(FlowequationVertexEquationEnum)[io_index]));
	
	/*set single point constraints: */

	/*spc all nodes on water*/
	if (!iomodel->Data(MaskVertexonwaterEnum)) _error_("iomodel->nodeonwater is NULL");
	if (iomodel->Data(MaskVertexonwaterEnum)[io_index]){
		for(k=1;k<=gsize;k++){
			this->FreezeDof(k);
		}
	}

	/*Diagnostic Horiz*/
	#ifdef _HAVE_DIAGNOSTIC_
	if (analysis_type==DiagnosticHorizAnalysisEnum){

		/*Coordinate system provided, convert to coord_system matrix*/
		_assert_(iomodel->Data(DiagnosticReferentialEnum)); 
		XZvectorsToCoordinateSystem(&this->coord_system[0][0],iomodel->Data(DiagnosticReferentialEnum)+io_index*6);

		if (dim==3){
			/*We have a  3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
			_assert_(iomodel->Data(MeshVertexonbedEnum)); 
			_assert_(iomodel->Data(FlowequationVertexEquationEnum));
			if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){
				for(k=1;k<=gsize;k++) this->FreezeDof(k);
			}
			if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealPattynApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){
				if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){
					for(k=1;k<=gsize;k++) this->FreezeDof(k);
				}
			}
			if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealStokesApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){
				if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){
					for(k=1;k<=2;k++) this->FreezeDof(k);
				}
			}
		}
		/*spc all nodes on hutter*/
		if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==HutterApproximationEnum){
			for(k=1;k<=gsize;k++){
				this->FreezeDof(k);
			}
		}
	}
	#endif

	/*Diagnostic Hutter*/
	if (analysis_type==DiagnosticHutterAnalysisEnum){
		if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL");
		/*Constrain all nodes that are not Hutter*/
		if (!iomodel->Data(FlowequationVertexEquationEnum)[io_index]==HutterApproximationEnum){
			for(k=1;k<=gsize;k++){
				this->FreezeDof(k);
			}
		}
	}

	/*Prognostic/ Melting/ Slopecompute/ Balancethickness*/
	if (
				analysis_type==PrognosticAnalysisEnum || 
				analysis_type==MeltingAnalysisEnum || 
				analysis_type==BedSlopeAnalysisEnum || 
				analysis_type==SurfaceSlopeAnalysisEnum || 
				analysis_type==BalancethicknessAnalysisEnum
				){
		if (dim==3){
			/*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
			_assert_(iomodel->Data(MeshVertexonbedEnum));
			if (!iomodel->Data(MeshVertexonbedEnum)[io_index]){
				for(k=1;k<=gsize;k++){
					this->FreezeDof(k);
				}
			}
		}
	}

}
/*}}}*/
/*FUNCTION Node::~Node(){{{1*/
Node::~Node(){
	delete inputs;
	delete hvertex;
	return;
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION Node::Echo{{{1*/
void Node::Echo(void){

	printf("Node:\n");
	printf("   id: %i\n",id);
	printf("   sid: %i\n",sid);
	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
	indexing.Echo();
	printf("   hvertex:     not displayed\n");
	printf("   inputs:      %p\n",inputs);


}
/*}}}*/
/*FUNCTION Node::DeepEcho{{{1*/
void Node::DeepEcho(void){

	printf("Node:\n");
	printf("   id: %i\n",id);
	printf("   sid: %i\n",sid);
	printf("   analysis_type: %s\n",EnumToStringx(analysis_type));
	indexing.DeepEcho();
	printf("Vertex:\n");
	hvertex->DeepEcho();
	printf("   inputs\n");


}
/*}}}*/
/*FUNCTION Node::Id{{{1*/
int    Node::Id(void){ return id; }
/*}}}*/
/*FUNCTION Node::MyRank{{{1*/
int    Node::MyRank(void){ 
	extern int my_rank;

	return my_rank; 
}
/*}}}*/
#ifdef _SERIAL_
/*FUNCTION Node::Marshall{{{1*/
void  Node::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;
	char* marshalled_inputs=NULL;
	int   marshalled_inputssize;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*get enum type of Node: */
	enum_type=NodeEnum;
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	
	/*marshall Node data: */
	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
	memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
	memcpy(marshalled_dataset,&coord_system,9*sizeof(double));marshalled_dataset+=9*sizeof(double);  

	/*marshall objects: */
	indexing.Marshall(&marshalled_dataset);
	hvertex->Marshall(&marshalled_dataset);

	/*Marshall inputs: */
	marshalled_inputssize=inputs->MarshallSize();
	marshalled_inputs=inputs->Marshall();
	memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputssize*sizeof(char));
	marshalled_dataset+=marshalled_inputssize;

	/*Free ressources:*/
	xfree((void**)&marshalled_inputs);

	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}*/
/*FUNCTION Node::MarshallSize{{{1*/
int   Node::MarshallSize(){

	return sizeof(id)+
		sizeof(sid)+
		indexing.MarshallSize()+
		hvertex->MarshallSize()+
		inputs->MarshallSize()+
		sizeof(analysis_type)+
		9*sizeof(double)+
		sizeof(int); //sizeof(int) for enum type
}
/*}}}*/
/*FUNCTION Node::Demarshall{{{1*/
void  Node::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
	 *object data (thanks to DataSet::Demarshall):*/

	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
	memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
	memcpy(&coord_system,marshalled_dataset,9*sizeof(double));marshalled_dataset+=9*sizeof(double);
	
	/*demarshall objects: */
	indexing.Demarshall(&marshalled_dataset);
	hvertex=new Hook(); hvertex->Demarshall(&marshalled_dataset);

	/*demarshall inputs: */
	inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}*/
#endif
/*FUNCTION Node::ObjectEnum{{{1*/
int Node::ObjectEnum(void){

	return NodeEnum;

}
/*}}}*/

/*Node management:*/
/*FUNCTION Node::Configure {{{1*/
void  Node::Configure(DataSet* nodesin,Vertices* verticesin){

	int i;

	/*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
	 * datasets, using internal ids and off_sets hidden in hooks: */
	hvertex->configure(verticesin);

}
/*FUNCTION Node::SetCurrentConfiguration {{{1*/
void  Node::SetCurrentConfiguration(DataSet* nodesin,Vertices* verticesin){

}
/*FUNCTION Node::GetDof {{{1*/
int   Node::GetDof(int dofindex,int setenum){

	if(setenum==GsetEnum){
		_assert_(dofindex>=0 && dofindex<indexing.gsize);
		return indexing.gdoflist[dofindex];
	}
	else if(setenum==FsetEnum){
		_assert_(dofindex>=0 && dofindex<indexing.fsize);
		return indexing.fdoflist[dofindex];
	}
	else if(setenum==SsetEnum){
		_assert_(dofindex>=0 && dofindex<indexing.ssize);
		return indexing.sdoflist[dofindex];
	}
	else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");

}
/*}}}*/
/*FUNCTION Node::GetDofList1{{{1*/
int  Node::GetDofList1(void){

	Vertex* vertex=NULL;

	vertex=(Vertex*)this->hvertex->delivers();

	return vertex->dof;
}
/*}}}*/
/*FUNCTION Node::GetDofList{{{1*/
void  Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){
	int i;
	int count=0;
	int count2=0;
		
	if(approximation_enum==NoneApproximationEnum){
		if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i];
		if(setenum==FsetEnum)for(i=0;i<this->indexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i];
		if(setenum==SsetEnum)for(i=0;i<this->indexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i];
	}
	else{

		if(setenum==GsetEnum){
			if(indexing.doftype){
				count=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.doftype[i]==approximation_enum){
						outdoflist[count]=indexing.gdoflist[i];
						count++;
					}
				}
				_assert_(count); //at least one dof should be the approximation requested
			}
			else for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i];
		}
		else if(setenum==FsetEnum){
			if(indexing.doftype){
				count=0;
				count2=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.f_set[i]){
						if(indexing.doftype[i]==approximation_enum){
							outdoflist[count]=indexing.fdoflist[count2];
							count++;
						}
						count2++;
					}
				}
			}
			else for(i=0;i<this->indexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i];
		}
		else if(setenum==SsetEnum){
			if(indexing.doftype){
				count=0;
				count2=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.s_set[i]){
						if(indexing.doftype[i]==approximation_enum){
							outdoflist[count]=indexing.sdoflist[count2];
							count++;
						}
						count2++;
					}
				}
			}
			else for(i=0;i<this->indexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i];
		}
		else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
	}
}
/*}}}*/
/*FUNCTION Node::GetSidList{{{1*/
int  Node::GetSidList(void){

	Vertex* vertex=NULL;

	vertex=(Vertex*)this->hvertex->delivers();

	return vertex->sid;
}
/*}}}*/
/*FUNCTION Node::GetLocalDofList{{{1*/
void  Node::GetLocalDofList(int* outdoflist,int approximation_enum,int setenum){
	int i;
	int count=0;
	int count2=0;
		
	if(approximation_enum==NoneApproximationEnum){
		if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i;
		else if(setenum==FsetEnum){
			count=0;
			for(i=0;i<this->indexing.gsize;i++){
				if(indexing.f_set[i]){
					outdoflist[count]=i;
					count++;
				}
			}
		}
		else if(setenum==SsetEnum){
			count=0;
			for(i=0;i<this->indexing.gsize;i++){
				if(indexing.s_set[i]){
					outdoflist[count]=i;
					count++;
				}
			}
		}
		else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
	}
	else{

		if(setenum==GsetEnum){
			if(indexing.doftype){
				count=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.doftype[i]==approximation_enum){
						outdoflist[count]=count;
						count++;
					}
				}
				_assert_(count);
			}
			else for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i;
		}
		else if(setenum==FsetEnum){

			if(indexing.doftype){
				count=0;
				count2=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.doftype[i]==approximation_enum){
						if(indexing.f_set[i]){
							outdoflist[count]=count2;
							count++;
						}
						count2++;
					}
				}
				_assert_(count2);
			}
			else{

				count=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.f_set[i]){
						outdoflist[count]=i;
						count++;
					}
				}
			}
		}
		else if(setenum==SsetEnum){
			if(indexing.doftype){
				count=0;
				count2=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.doftype[i]==approximation_enum){
						if(indexing.s_set[i]){
							outdoflist[count]=count2;
							count++;
						}
						count2++;
					}
				}
				_assert_(count2);
			}
			else{
				count=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(indexing.s_set[i]){
						outdoflist[count]=i;
						count++;
					}
				}
			}
		}
		else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
	}
}
/*}}}*/
/*FUNCTION Node::Sid{{{1*/
int    Node::Sid(void){ return sid; }
/*}}}*/
/*FUNCTION Node::GetVertexId {{{1*/
int   Node::GetVertexId(void){

	Vertex*  vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	return vertex->id;
}
/*}}}*/
/*FUNCTION Node::GetVertexDof {{{1*/
int   Node::GetVertexDof(void){

	Vertex*  vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	return vertex->dof;
}
/*}}}*/
#ifdef _HAVE_DIAGNOSTIC_
/*FUNCTION Node::GetCoordinateSystem{{{1*/
void Node::GetCoordinateSystem(double* coord_system_out){

	/*Copy coord_system*/
	for(int k=0;k<3;k++) for(int l=0;l<3;l++) coord_system_out[3*k+l]=this->coord_system[k][l];

}
/*}}}*/
#endif
/*FUNCTION Node::SetVertexDof {{{1*/
void   Node::SetVertexDof(int in_dof){

	Vertex*  vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	vertex->dof=in_dof;

}
/*}}}*/
/*FUNCTION Node::InAnalysis{{{1*/
bool Node::InAnalysis(int in_analysis_type){
	if (in_analysis_type==this->analysis_type) return true;
	else return false;
}
/*}}}*/

/*Node numerics:*/
/*FUNCTION Node::ApplyConstraints{{{1*/
void  Node::ApplyConstraint(int dof,double value){

	int index;

	/*Dof should be added in the s set, describing which 
	 * dofs are constrained to a certain value (dirichlet boundary condition*/
	DofInSSet(dof-1);
	this->indexing.svalues[dof-1]=value;
}
/*}}}*/
/*FUNCTION Node::RelaxConstraint{{{1*/
void  Node::RelaxConstraint(int dof){

	/*Dof should be added to the f-set, and taken out of the s-set:*/
	DofInFSet(dof-1);
	this->indexing.svalues[dof-1]=NAN;
}
/*}}}*/
/*FUNCTION Node::CreateVecSets {{{1*/
void  Node::CreateVecSets(Vector* pv_g,Vector* pv_f,Vector* pv_s){

	double gvalue=1.0; //all nodes are in the g set;
	double value;

	int i;

	for(i=0;i<this->indexing.gsize;i++){

		/*g set: */
		pv_g->SetValue(indexing.gdoflist[i],gvalue,INS_VAL);
		
		/*f set: */
		value=(double)this->indexing.f_set[i];
		pv_f->SetValue(indexing.gdoflist[i],value,INS_VAL);

		/*s set: */
		value=(double)this->indexing.s_set[i];
		pv_s->SetValue(indexing.gdoflist[i],value,INS_VAL);

	}


}
/*}}}*/
/*FUNCTION Node::CreateNodalConstraints{{{1*/
void  Node::CreateNodalConstraints(Vector* ys){

	int i;
	double* values=NULL;
	int count;

	/*Recover values for s set and plug them in constraints vector: */
	if(this->indexing.ssize){
		values=(double*)xmalloc(this->indexing.ssize*sizeof(double));
		count=0;
		for(i=0;i<this->indexing.gsize;i++){
			if(this->indexing.s_set[i]){
				values[count]=this->indexing.svalues[i];
				_assert_(!isnan(values[count]));
				count++;
			}
		}
		
		/*Add values into constraint vector: */
		ys->SetValues(this->indexing.ssize,this->indexing.sdoflist,values,INS_VAL);
	}

	/*Free ressources:*/
	xfree((void**)&values);


}
/*}}}*/
/*FUNCTION Node::DofInSSet {{{1*/
void  Node::DofInSSet(int dof){

	/*Put dof for this node into the s set (ie, this dof will be constrained 
	 * to a fixed value during computations. */

	this->indexing.f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
	this->indexing.s_set[dof]=1;
}
/*}}}*/
/*FUNCTION Node::DofInFSet {{{1*/
void  Node::DofInFSet(int dof){

	/*Put dof for this node into the f set (ie, this dof will NOT be constrained 
	 * to a fixed value during computations. */

	this->indexing.f_set[dof]=1; 
	this->indexing.s_set[dof]=0;
}
/*}}}*/
/*FUNCTION Node::FreezeDof{{{1*/
void  Node::FreezeDof(int dof){
	
	DofInSSet(dof-1); //with 0 displacement for this dof.

}
/*}}}*/
/*FUNCTION Node::GetApproximation {{{1*/
int   Node::GetApproximation(){

	int approximation;

	/*recover parameters: */
	inputs->GetInputValue(&approximation,ApproximationEnum);

	return approximation;
}
/*}}}*/
/*FUNCTION Node::GetConnectivity {{{1*/
int Node::GetConnectivity(){

	Vertex*  vertex=NULL;
	vertex=(Vertex*)hvertex->delivers();
	return vertex->connectivity;
}
/*}}}*/
/*FUNCTION Node::GetNumberOfDofs{{{1*/
int   Node::GetNumberOfDofs(int approximation_enum,int setenum){

	/*Get number of degrees of freedom in a node, for a certain set (g,f or s-set)
	 *and for a certain approximation type: */
	
	int i;
	int numdofs=0;

	if(approximation_enum==NoneApproximationEnum){
		if (setenum==GsetEnum) numdofs=this->indexing.gsize;
		else if (setenum==FsetEnum) numdofs=this->indexing.fsize;
		else if (setenum==SsetEnum) numdofs=this->indexing.ssize;
		else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
	}
	else{
		if(setenum==GsetEnum){
			if(this->indexing.doftype){
				numdofs=0;
				for(i=0;i<this->indexing.gsize;i++){
					if(this->indexing.doftype[i]==approximation_enum) numdofs++;
				}
			}
			else numdofs=this->indexing.gsize;
		}
		else if (setenum==FsetEnum){
			if(this->indexing.doftype){
				numdofs=0;
				for(i=0;i<this->indexing.gsize;i++){
					if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.f_set[i])) numdofs++;
				}
			}
			else numdofs=this->indexing.fsize;
		}
		else if (setenum==SsetEnum){
			if(this->indexing.doftype){
				numdofs=0;
				for(i=0;i<this->indexing.gsize;i++){
					if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.s_set[i])) numdofs++;
				}
			}
			else numdofs=this->indexing.ssize;
		}
		else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
	}
	return numdofs;
}
/*}}}*/
/*FUNCTION Node::GetSigma {{{1*/
double Node::GetSigma(){
	Vertex* vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	return vertex->sigma;
}
/*}}}*/
/*FUNCTION Node::GetX {{{1*/
double Node::GetX(){
	Vertex* vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	return vertex->x;
}
/*}}}*/
/*FUNCTION Node::GetY {{{1*/
double Node::GetY(){
	Vertex* vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	return vertex->y;
}
/*}}}*/
/*FUNCTION Node::GetZ {{{1*/
double Node::GetZ(){
	Vertex* vertex=NULL;

	vertex=(Vertex*)hvertex->delivers();
	return vertex->z;
}
/*}}}*/
/*FUNCTION Node::IsClone {{{1*/
int   Node::IsClone(){
	
	return indexing.clone;

}
/*}}}*/
/*FUNCTION Node::IsOnBed {{{1*/
int   Node::IsOnBed(){

	bool onbed;

	/*recover parameters: */
	inputs->GetInputValue(&onbed,MeshVertexonbedEnum);

	return onbed;
}
/*}}}*/
/*FUNCTION Node::IsGrounded {{{1*/
int   Node::IsGrounded(){

	bool onsheet;

	/*recover parameters: */
	inputs->GetInputValue(&onsheet,MaskVertexongroundediceEnum);

	return onsheet;
}		
/*}}}*/
/*FUNCTION Node::IsFloating {{{1*/
int   Node::IsFloating(){
	
	bool onshelf;

	/*recover parameters: */
	inputs->GetInputValue(&onshelf,MaskVertexonfloatingiceEnum);

	return onshelf;
}
/*}}}*/
/*FUNCTION Node::IsOnSurface {{{1*/
int   Node::IsOnSurface(){

	bool onsurface;

	/*recover parameters: */
	inputs->GetInputValue(&onsurface,MeshVertexonsurfaceEnum);

	return onsurface;
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromVector(double* vector, int name, int type){{{1*/
void  Node::InputUpdateFromVector(double* vector, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromVector(int* vector, int name, int type){{{1*/
void  Node::InputUpdateFromVector(int* vector, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromVector(bool* vector, int name, int type){{{1*/
void  Node::InputUpdateFromVector(bool* vector, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromVectorDakota(double* vector, int name, int type){{{1*/
void  Node::InputUpdateFromVectorDakota(double* vector, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){{{1*/
void  Node::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromVectorDakota(int* vector, int name, int type){{{1*/
void  Node::InputUpdateFromVectorDakota(int* vector, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromVectorDakota(bool* vector, int name, int type){{{1*/
void  Node::InputUpdateFromVectorDakota(bool* vector, int name, int type){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromConstant(double constant, int name){{{1*/
void  Node::InputUpdateFromConstant(double constant, int name){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromConstant(int constant, int name){{{1*/
void  Node::InputUpdateFromConstant(int constant, int name){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::InputUpdateFromConstant(bool constant, int name){{{1*/
void  Node::InputUpdateFromConstant(bool constant, int name){

	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Node::UpdateSpcs {{{1*/
void   Node::UpdateSpcs(double* ys){

	int     count=0;
	int     i;

	count=0;
	for(i=0;i<this->indexing.gsize;i++){
		if(this->indexing.s_set[i]){
			this->indexing.svalues[i]=ys[this->indexing.sdoflist[count]];
			count++;
		}
	}
}
/*}}}*/
/*FUNCTION Node::VecMerge {{{1*/
void   Node::VecMerge(Vector* ug, double* vector_serial,int setenum){

	double* values=NULL;
	int*    indices=NULL;
	int     count=0;
	int     i;

	if(setenum==FsetEnum){
		if(this->indexing.fsize){
			indices=(int*)xmalloc(this->indexing.fsize*sizeof(int));
			values=(double*)xmalloc(this->indexing.fsize*sizeof(double));

			for(i=0;i<this->indexing.gsize;i++){
				if(this->indexing.f_set[i]){
					_assert_(vector_serial);
					values[count]=vector_serial[this->indexing.fdoflist[count]];
					indices[count]=this->indexing.gdoflist[i];
					count++;
				}
			}

			/*Add values into ug: */
			ug->SetValues(this->indexing.fsize,indices,values,INS_VAL);
		}
	}
	else if(setenum==SsetEnum){
		if(this->indexing.ssize){
			indices=(int*)xmalloc(this->indexing.ssize*sizeof(int));
			values=(double*)xmalloc(this->indexing.ssize*sizeof(double));

			for(i=0;i<this->indexing.gsize;i++){
				if(this->indexing.s_set[i]){
					_assert_(vector_serial);
					values[count]=vector_serial[this->indexing.sdoflist[count]];
					indices[count]=this->indexing.gdoflist[i];
					count++;
				}
			}

			/*Add values into ug: */
			ug->SetValues(this->indexing.ssize,indices,values,INS_VAL);
		}
	}
	else _error_("VecMerge can only merge from the s or f-set onto the g-set!");

	/*Free ressources:*/
	xfree((void**)&values);
	xfree((void**)&indices);
}
/*}}}*/
/*FUNCTION Node::VecReduce {{{1*/
void   Node::VecReduce(Vector* vector, double* ug_serial,int setenum){

	double* values=NULL;
	int     count=0;
	int     i;

	if(setenum==FsetEnum){
		if(this->indexing.fsize){
			values=(double*)xmalloc(this->indexing.fsize*sizeof(double));

			for(i=0;i<this->indexing.gsize;i++){
				if(this->indexing.f_set[i]){
					_assert_(ug_serial);
					values[count]=ug_serial[this->indexing.gdoflist[i]];
					count++;
				}
			}

			/*Add values into ug: */
			vector->SetValues(this->indexing.fsize,this->indexing.fdoflist,values,INS_VAL);
		}
	}
	else if(setenum==SsetEnum){
		if(this->indexing.ssize){
			values=(double*)xmalloc(this->indexing.ssize*sizeof(double));

			for(i=0;i<this->indexing.gsize;i++){
				if(this->indexing.s_set[i]){
					_assert_(ug_serial);
					values[count]=ug_serial[this->indexing.gdoflist[i]];
					count++;
				}
			}

			/*Add values into ug: */
			vector->SetValues(this->indexing.ssize,this->indexing.sdoflist,values,INS_VAL);
		}
	}
	else _error_("VecReduce can only merge from the s or f-set onto the g-set!");

	/*Free ressources:*/
	xfree((void**)&values);
}
/*}}}*/

/* DofObject routines:*/
/*FUNCTION Node::DistributeDofs{{{1*/
void  Node::DistributeDofs(int* pdofcount,int setenum){

	int i;
	extern int my_rank;
	int dofcount;

	dofcount=*pdofcount;

	/*Initialize: */
	if(setenum==FsetEnum) this->indexing.InitSet(setenum);
	if(setenum==SsetEnum) this->indexing.InitSet(setenum);
	
		
	/*For clone nodfs, don't distribute dofs, we will get them from another cpu in UpdateCloneDofs!*/
	if(indexing.clone){
		return;
	}

	/*This node should distribute dofs for setenum set (eg, f_set or s_set), go ahead: */

	if(setenum==GsetEnum){
		for(i=0;i<this->indexing.gsize;i++){
			indexing.gdoflist[i]=dofcount+i;
		}
		dofcount+=this->indexing.gsize;
	}
	else if(setenum==FsetEnum){
		for(i=0;i<this->indexing.fsize;i++){
			indexing.fdoflist[i]=dofcount+i;
		}
		dofcount+=this->indexing.fsize;
	}
	else if(setenum==SsetEnum){
		for(i=0;i<this->indexing.ssize;i++){
			indexing.sdoflist[i]=dofcount+i;
		}
		dofcount+=this->indexing.ssize;
	}
	else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");


	/*Assign output pointers: */
	*pdofcount=dofcount;

}
/*}}}*/
/*FUNCTION Node::Off_setDofs{{{1*/
void  Node::OffsetDofs(int dofcount,int setenum){
	
	int i;
	extern int my_rank;
	
	if(indexing.clone){
		/*This node is a clone, don't off_set the dofs!: */
		return;
	}

	/*This node should off_set the dofs, go ahead: */
	if(setenum==GsetEnum){
		for(i=0;i<this->indexing.gsize;i++) indexing.gdoflist[i]+=dofcount;
	}
	else if(setenum==FsetEnum){
		for(i=0;i<this->indexing.fsize;i++) indexing.fdoflist[i]+=dofcount;
	}
	else if(setenum==SsetEnum){
		for(i=0;i<this->indexing.ssize;i++) indexing.sdoflist[i]+=dofcount;
	}
	else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");
}
/*}}}*/
/*FUNCTION Node::ShowTrueDofs{{{1*/
void  Node::ShowTrueDofs(int* truedofs, int ncols,int setenum){

	int j;
	extern int my_rank;
	
	/*Are we a clone? : */
	if(indexing.clone)return;

	/*Ok, we are not a clone, just plug our dofs into truedofs: */
	if(setenum==GsetEnum)for(j=0;j<this->indexing.gsize;j++)  *(truedofs+ncols*sid+j)=indexing.gdoflist[j];
	else if(setenum==FsetEnum)for(j=0;j<this->indexing.fsize;j++)  *(truedofs+ncols*sid+j)=indexing.fdoflist[j];
	else if(setenum==SsetEnum)for(j=0;j<this->indexing.ssize;j++)  *(truedofs+ncols*sid+j)=indexing.sdoflist[j];
	else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");

}
/*}}}*/
/*FUNCTION Node::UpdateCloneDofs{{{1*/
void  Node::UpdateCloneDofs(int* alltruedofs,int ncols,int setenum){

	int j;
	extern int my_rank;
	
	/*If we are not a clone, don't update, we already have dofs!: */
	if(indexing.clone==0)return;


	/*Ok, we are a clone node, but we did not create the dofs for this node.
	 *      * Therefore, our doflist is garbage right now. Go pick it up in the alltruedofs: */
	if(setenum==GsetEnum)for(j=0;j<this->indexing.gsize;j++) indexing.gdoflist[j]=*(alltruedofs+ncols*sid+j);
	else if(setenum==FsetEnum)for(j=0;j<this->indexing.fsize;j++) indexing.fdoflist[j]=*(alltruedofs+ncols*sid+j);
	else if(setenum==SsetEnum)for(j=0;j<this->indexing.ssize;j++) indexing.sdoflist[j]=*(alltruedofs+ncols*sid+j);
	else _error_("%s%s%s"," set of enum type ",EnumToStringx(setenum)," not supported yet!");

}
/*}}}*/
/*FUNCTION Node::SetClone {{{1*/
void  Node::SetClone(int* minranks){

	extern int my_rank;

	if (minranks[sid]==my_rank){
		indexing.clone=0;
	}
	else{
		/*!there is a cpu with lower rank that has the same node, 
		therefore, I am a clone*/
		indexing.clone=1; 	
	}

}
/*}}}*/