source: issm/branches/trunk-jpl-damage/src/c/objects/Node.cpp@ 11684

/*!\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(Vec pv_g,Vec pv_f,Vec 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: */
                VecSetValues(pv_g,1,&indexing.gdoflist[i],&gvalue,INSERT_VALUES);
                
                /*f set: */
                value=(double)this->indexing.f_set[i];
                VecSetValues(pv_f,1,&indexing.gdoflist[i],&value,INSERT_VALUES);

                /*s set: */
                value=(double)this->indexing.s_set[i];
                VecSetValues(pv_s,1,&indexing.gdoflist[i],&value,INSERT_VALUES);

        }


}
/*}}}*/
/*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,INSERT_VALUES);
        }

        /*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,INSERT_VALUES);
                }
        }
        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,INSERT_VALUES);
                }
        }
        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,INSERT_VALUES);
                }
        }
        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,INSERT_VALUES);
                }
        }
        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;       
        }

}
/*}}}*/