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

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

#include "./classes.h"
#include "shared/shared.h"
#include "modules/ModelProcessorx/ModelProcessorx.h"
/*}}}*/

/*Node constructors and destructors:*/
/*FUNCTION Node::Node() default constructor {{{*/
Node::Node(){
	this->approximation=0;
}
/*}}}*/
/*FUNCTION Node::Node(int node_id,int node_sid,int node_lid,int io_index, IoModel* iomodel,int analysis_type,int approximation) {{{*/
Node::Node(int node_id,int node_sid,int node_lid,int io_index, IoModel* iomodel,int analysis_type,int in_approximation){

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

	/*id: */
	this->id            = node_id;
	this->sid           = node_sid;
	this->lid           = node_lid;
	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:*/
	this->indexingupdate = true;
	DistributeNumDofs(&this->indexing,analysis_type,in_approximation); //number of dofs per node
	gsize=this->indexing.gsize;

	if(analysis_type==StressbalanceAnalysisEnum)
	 this->approximation=in_approximation;
	else
	 this->approximation=0;

	/*Stressbalance Horiz*/
	#ifdef _HAVE_STRESSBALANCE_
	if(analysis_type==StressbalanceAnalysisEnum){

		/*Coordinate system provided, convert to coord_system matrix*/
		_assert_(iomodel->Data(StressbalanceReferentialEnum)); 
		XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[io_index*6]);
		_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4);

		if(iomodel->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(in_approximation==SSAApproximationEnum && !reCast<int>(iomodel->Data(MeshVertexonbedEnum)[io_index])){
				this->Deactivate();
			}
			if(in_approximation==L1L2ApproximationEnum && !reCast<int>(iomodel->Data(MeshVertexonbedEnum)[io_index])){
				this->Deactivate();
			}
			if(in_approximation==SSAHOApproximationEnum && reCast<int>(iomodel->Data(FlowequationBorderSSAEnum)[io_index])){
				if(!reCast<int>(iomodel->Data(MeshVertexonbedEnum)[io_index])){
					this->Deactivate();
				}
			}
			if(in_approximation==SSAFSApproximationEnum && reCast<int>(iomodel->Data(FlowequationBorderSSAEnum)[io_index])){
				if(!reCast<int>(iomodel->Data(MeshVertexonbedEnum)[io_index])){
					for(k=1;k<=2;k++) this->FreezeDof(k);
				}
			}
		}
		/*spc all nodes on SIA*/
		if(in_approximation==SIAApproximationEnum){
			this->Deactivate();
		}
	}
	#endif

	/*2d solutions in 3d, we need to constrain all the nodes that are not on base*/
	if(
				analysis_type==FreeSurfaceBaseAnalysisEnum || 
				analysis_type==MasstransportAnalysisEnum || 
				analysis_type==MeltingAnalysisEnum || 
				analysis_type==BedSlopeAnalysisEnum || 
				analysis_type==SurfaceSlopeAnalysisEnum || 
				analysis_type==BalancethicknessAnalysisEnum ||
				analysis_type==HydrologyDCInefficientAnalysisEnum ||
				analysis_type==HydrologyDCEfficientAnalysisEnum
				){
		if(iomodel->dim==3){
			/*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
			_assert_(iomodel->Data(MeshVertexonbedEnum));
			if(!(reCast<bool>(iomodel->Data(MeshVertexonbedEnum)[io_index]))){
				this->Deactivate();
			}
		}
	}
	if(analysis_type==FreeSurfaceTopAnalysisEnum){
		if(iomodel->dim==3){
			/*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
			_assert_(iomodel->Data(MeshVertexonsurfaceEnum));
			if(!(reCast<bool>(iomodel->Data(MeshVertexonsurfaceEnum)[io_index]))){
				this->Deactivate();
			}
		}
	}

}
/*}}}*/
/*FUNCTION Node::~Node(){{{*/
Node::~Node(){
	return;
}
/*}}}*/

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

	_printf_("Node:\n");
	_printf_("   id : " << id << "\n");
	_printf_("   sid: " << sid << "\n");
	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
	_printf_("   approximation: " << EnumToStringx(approximation) << "\n");
	_printf_("   indexingupdate: " << indexingupdate << "\n");
	indexing.Echo();

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

	_printf_("Node:\n");
	_printf_("   id: " << id << "\n");
	_printf_("   sid: " << sid << "\n");
	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
	_printf_("   approximation: " << EnumToStringx(approximation) << "\n");
	_printf_("   indexingupdate: " << indexingupdate << "\n");
	indexing.DeepEcho();

}
/*}}}*/
/*FUNCTION Node::Id{{{*/
int    Node::Id(void){ return id; }
/*}}}*/
/*FUNCTION Node::ObjectEnum{{{*/
int Node::ObjectEnum(void){

	return NodeEnum;

}
/*}}}*/

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

	_assert_(!this->indexingupdate);
	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_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");

} /*}}}*/
/*FUNCTION Node::GetDofList{{{*/
void  Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){
	int i;
	int count=0;
	int count2=0;

	_assert_(!this->indexingupdate);

	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_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
	}
}
/*}}}*/
/*FUNCTION Node::GetLocalDofList{{{*/
void  Node::GetLocalDofList(int* outdoflist,int approximation_enum,int setenum){
	int i;
	int count=0;
	int count2=0;

	_assert_(!this->indexingupdate);

	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_("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_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
	}
}
/*}}}*/
/*FUNCTION Node::Sid{{{*/
int Node::Sid(void){
	return sid; 
}
/*}}}*/
/*FUNCTION Node::Lid{{{*/
int  Node::Lid(void){
	return lid; 
}
/*}}}*/
#ifdef _HAVE_STRESSBALANCE_
/*FUNCTION Node::GetCoordinateSystem{{{*/
void Node::GetCoordinateSystem(IssmDouble* 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::InAnalysis{{{*/
bool Node::InAnalysis(int in_analysis_type){
	if (in_analysis_type==this->analysis_type) return true;
	else return false;
}
/*}}}*/

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

	/*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::RequiresDofReindexing{{{*/
bool Node::RequiresDofReindexing(void){

	return this->indexingupdate;

}
/*}}}*/
/*FUNCTION Node::ReindexingDone{{{*/
void Node::ReindexingDone(void){

	this->indexingupdate = false;

}
/*}}}*/
/*FUNCTION Node::RelaxConstraint{{{*/
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]=0.;
}
/*}}}*/
/*FUNCTION Node::CreateVecSets {{{*/
void  Node::CreateVecSets(Vector<IssmDouble>* pv_g,Vector<IssmDouble>* pv_f,Vector<IssmDouble>* pv_s){

	IssmDouble gvalue=1.; //all nodes are in the g set;
	IssmDouble value;

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

		/*g set: */
		pv_g->SetValue(indexing.gdoflist[i],gvalue,INS_VAL);

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

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

	}

}
/*}}}*/
/*FUNCTION Node::CreateNodalConstraints{{{*/
void  Node::CreateNodalConstraints(Vector<IssmDouble>* ys){

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

	/*Recover values for s set and plug them in constraints vector: */
	if(this->indexing.ssize){
		values=xNew<IssmDouble>(this->indexing.ssize);
		count=0;
		for(i=0;i<this->indexing.gsize;i++){
			if(this->indexing.s_set[i]){
				values[count]=this->indexing.svalues[i];
				_assert_(!xIsNan<IssmDouble>(values[count]));
				count++;
			}
		}

		/*Add values into constraint vector: */
		ys->SetValues(this->indexing.ssize,this->indexing.sdoflist,values,INS_VAL);
	}

	/*Free ressources:*/
	xDelete<IssmDouble>(values);

}
/*}}}*/
/*FUNCTION Node::DofInSSet {{{*/
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. */
	_assert_(dof<this->indexing.gsize);

	if(this->indexing.f_set[dof] == 1){
		this->indexingupdate = true;
		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 {{{*/
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. */

	if(this->indexing.f_set[dof] == 0){
		this->indexingupdate = true;
		this->indexing.f_set[dof]=1; 
		this->indexing.s_set[dof]=0;
	}
}
/*}}}*/
/*FUNCTION Node::FreezeDof{{{*/
void  Node::FreezeDof(int dof){

	DofInSSet(dof-1); //with 0 displacement for this dof.

}
/*}}}*/
/*FUNCTION Node::Deactivate{{{*/
void  Node::Deactivate(void){

	if(IsActive()){
		this->indexingupdate = true;
		indexing.Deactivate();
	}

}
/*}}}*/
/*FUNCTION Node::Activate{{{*/
void  Node::Activate(void){

	if(!IsActive()){
		this->indexingupdate = true;
		indexing.Activate();
	}

}
/*}}}*/
/*FUNCTION Node::IsActive{{{*/
bool Node::IsActive(void){

	return indexing.active;

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

	return approximation;
}
/*}}}*/
/*FUNCTION Node::GetNumberOfDofs{{{*/
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_("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_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
	}
	return numdofs;
}
/*}}}*/
/*FUNCTION Node::IsClone {{{*/
int   Node::IsClone(){

	return indexing.clone;

}
/*}}}*/
/*FUNCTION Node::UpdateSpcs {{{*/
void   Node::UpdateSpcs(IssmDouble* ys){

	int count = 0;
	for(int 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 {{{*/
void Node::VecMerge(Vector<IssmDouble>* ug, IssmDouble* vector_serial,int setenum){

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

	if(setenum==FsetEnum){
		if(this->indexing.fsize){
			indices=xNew<int>(this->indexing.fsize);
 			values=xNew<IssmDouble>(this->indexing.fsize);

			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=xNew<int>(this->indexing.ssize);
			values=xNew<IssmDouble>(this->indexing.ssize);

			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:*/
	xDelete<IssmDouble>(values);
	xDelete<int>(indices);
}
/*}}}*/
/*FUNCTION Node::VecReduce {{{*/
void   Node::VecReduce(Vector<IssmDouble>* vector, IssmDouble* ug_serial,int setenum){

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

	if(setenum==FsetEnum){
		if(this->indexing.fsize){
 			values=xNew<IssmDouble>(this->indexing.fsize);

			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=xNew<IssmDouble>(this->indexing.ssize);

			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:*/
	xDelete<IssmDouble>(values);
}
/*}}}*/

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

	int i;
	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_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");

	/*Assign output pointers: */
	*pdofcount=dofcount;
}
/*}}}*/
/*FUNCTION Node::OffsetDofs{{{*/
void  Node::OffsetDofs(int dofcount,int setenum){

	int i;

	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_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
}
/*}}}*/
/*FUNCTION Node::ShowTrueDofs{{{*/
void  Node::ShowTrueDofs(int* truedofs, int ncols,int setenum){

	int j;

	/*Are we a clone? : */
	if(indexing.clone) return;

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

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

	int j;

	/*If we are not a clone, don't update, we already have dofs!: */
	if(!indexing.clone)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: */
	switch(setenum){
		case GsetEnum:
			for(j=0;j<this->indexing.gsize;j++) indexing.gdoflist[j]=alltruedofs[ncols*sid+j];
			break;
		case FsetEnum:
			for(j=0;j<this->indexing.fsize;j++) indexing.fdoflist[j]=alltruedofs[ncols*sid+j];
			break;
		case SsetEnum:
			for(j=0;j<this->indexing.ssize;j++) indexing.sdoflist[j]=alltruedofs[ncols*sid+j];
			break;
		default:
			_error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
	}
}
/*}}}*/
/*FUNCTION Node::SetClone {{{*/
void  Node::SetClone(int* minranks){

	int my_rank;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

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

/*Methods inherent to Node: */
int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation){ /*{{{*/

	int  i,j,count,numdof,numgdof;
	int* ndof_list=NULL;
	int* ngdof_list_cumulative=NULL;
	int *doflist = NULL;

	if(numnodes){
		/*allocate: */
		ndof_list=xNew<int>(numnodes);
		ngdof_list_cumulative=xNew<int>(numnodes);

		/*Get number of dofs per node, and total for this given set*/
		numdof=0;
		numgdof=0;
		for(i=0;i<numnodes;i++){

			/*Cumulative list= number of dofs before node i*/
			ngdof_list_cumulative[i]=numgdof;

			/*Number of dofs for node i for given set and for the g set*/
			ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,setenum);
			numgdof    +=nodes[i]->GetNumberOfDofs(approximation,GsetEnum);
			numdof     +=ndof_list[i];
		}

		if(numdof){
			/*Allocate: */
			doflist=xNew<int>(numdof);

			/*Populate: */
			count=0;
			for(i=0;i<numnodes;i++){
				nodes[i]->GetLocalDofList(&doflist[count],approximation,setenum);
				count+=ndof_list[i];
			}

			/*We now have something like: [0 1 0 2 1 2]. Offset by gsize, to get something like: [0 1 2 4 6 7]:*/
			count=0;
			for(i=0;i<numnodes;i++){
				for(j=0;j<ndof_list[i];j++){
					doflist[count+j]+=ngdof_list_cumulative[i];
				}
				count+=ndof_list[i];
			}
		}
		else doflist=NULL;
	}

	/*Free ressources:*/
	xDelete<int>(ndof_list);
	xDelete<int>(ngdof_list_cumulative);

	/*CLean-up and return*/
	return doflist;
}
/*}}}*/
int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/

	int  i,numdof,count;
	int* ndof_list=NULL;
	int *doflist = NULL;

	if(numnodes){

		/*Allocate:*/
		ndof_list=xNew<int>(numnodes);

		/*First, figure out size of doflist: */
		numdof=0;
		for(i=0;i<numnodes;i++){
			ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,setenum);
			numdof+=ndof_list[i];
		}

		if(numdof){
			/*Allocate: */
			doflist=xNew<int>(numdof);

			/*Populate: */
			count=0;
			for(i=0;i<numnodes;i++){
				nodes[i]->GetDofList(&doflist[count],approximation,setenum);
				count+=ndof_list[i];
			}
		}
		else doflist=NULL;
	}
	/*Free ressources:*/
	xDelete<int>(ndof_list);

	return doflist;
}
/*}}}*/
int GetNumberOfDofs(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/

	/*output: */
	int numberofdofs=0;

	for(int i=0;i<numnodes;i++){
		numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum);
	}

	return numberofdofs;
}
/*}}}*/
#ifdef _HAVE_STRESSBALANCE_
void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum){/*{{{*/

	int* cs_array=NULL;

	/*All nodes have the same Coordinate System*/
	cs_array=xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;

	/*Call core*/
	TransformInvStiffnessMatrixCoord(Ke,nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}
/*}}}*/
void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array){/*{{{*/

	int     i,j;
	int     numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);

	/*Transform matrix: R*Ke*R^T */
	TripleMultiply(transform,numdofs,numdofs,0,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,1,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}
/*}}}*/
void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum){/*{{{*/

	int* cs_array=NULL;

	/*All nodes have the same Coordinate System*/
	cs_array=xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;

	/*Call core*/
	TransformLoadVectorCoord(pe,nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}
/*}}}*/
void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array){/*{{{*/

	int     i;
	int     numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current load vector*/
	values=xNew<IssmDouble>(pe->nrows);
	for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);

	/*Transform matrix: R^T*pe */
	MatrixMultiply(transform,numdofs,numdofs,1,
				values,pe->nrows,1,0,
				&pe->values[0],0);

	/*Free Matrices*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}
/*}}}*/
void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int cs_enum){/*{{{*/

	int* cs_array=NULL;

	/*All nodes have the same Coordinate System*/
	cs_array=xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;

	/*Call core*/
	TransformSolutionCoord(solution,nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}
/*}}}*/
void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array){/*{{{*/

	int     i;
	int     numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current solution vector*/
	values=xNew<IssmDouble>(numdofs);
	for(i=0;i<numdofs;i++) values[i]=solution[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);

	/*Transform matrix: R*U */
	MatrixMultiply(transform,numdofs,numdofs,0,
				values,numdofs,1,0,
				&solution[0],0);

	/*Free Matrices*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}
/*}}}*/
void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum){/*{{{*/

	int* cs_array=NULL;

	/*All nodes have the same Coordinate System*/
	cs_array=xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;

	/*Call core*/
	TransformStiffnessMatrixCoord(Ke,nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}
/*}}}*/
void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array){/*{{{*/

	int         numdofs = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(int i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);

	/*Transform matrix: R^T*Ke*R */
	TripleMultiply(transform,numdofs,numdofs,1,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,0,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}
/*}}}*/
void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array){/*{{{*/

	int     i,counter;
	int     numdofs           = 0;
	IssmDouble  norm;
	IssmDouble *transform         = NULL;
	IssmDouble *values            = NULL;
	IssmDouble  coord_system[3][3];

	/*Some checks in debugging mode*/
	_assert_(numnodes && nodes);

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Allocate and initialize transform matrix*/
	transform=xNew<IssmDouble>(numdofs*numdofs);
	for(i=0;i<numdofs*numdofs;i++) transform[i]=0.0;

	/*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
	 *for 3 nodes:

	 *     | T1 0  0 |
	 * Q = | 0  T2 0 |
	 *     | 0  0  T3|
	 *
	 * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
	 * associated to this node*/
	counter=0;
	for(i=0;i<numnodes;i++){
		nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
		switch(cs_array[i]){
			case PressureEnum:
				/*DO NOT change anything*/
				transform[(numdofs)*(counter) + counter] = 1.;
				counter+=1;
				break;
			case XYEnum:
				/*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
				norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(norm>1.e-4);
				transform[(numdofs)*(counter+0) + counter+0] =   coord_system[0][0]/norm;
				transform[(numdofs)*(counter+0) + counter+1] = - coord_system[1][0]/norm;
				transform[(numdofs)*(counter+1) + counter+0] =   coord_system[1][0]/norm;
				transform[(numdofs)*(counter+1) + counter+1] =   coord_system[0][0]/norm;
				counter+=2;
				break;
			case XYZEnum:
				/*The 3 coordinates are changed (x,y,z)*/
				transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
				transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
				transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
				transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
				transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
				transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
				transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
				transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
				transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
				counter+=3;
				break;
			default:
				_error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Assign output pointer*/
	*ptransform=transform;
}
/*}}}*/
#endif
