Node.cpp

Go to the documentation of this file.

 
/*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"
#include "../analyses/analyses.h"
/*}}}*/
 
/*Node constructors and destructors:*/
Node::Node(){/*{{{*/
   this->approximation  = 0;
   this->gsize          = -1;
   this->fsize          = -1;
   this->ssize          = -1;
   this->clone          = false;
   this->active         = true;
   this->freeze         = false;
   this->f_set          = NULL;
   this->s_set          = NULL;
   this->svalues        = NULL;
   this->doftype        = NULL;
   this->gdoflist       = NULL;
   this->fdoflist       = NULL;
   this->sdoflist       = NULL;
   this->gdoflist_local = NULL;
   this->fdoflist_local = NULL;
   this->sdoflist_local = NULL;
}
/*}}}*/
Node::Node(int node_id,int node_sid,int io_index,bool node_clone,IoModel* iomodel,int node_analysis,int in_approximation,bool isamr){/*{{{*/
 
   /*id: */
   this->id            = node_id;
   this->sid           = node_sid;
   this->lid           = -1; /*Assigned by Finalize*/
   this->pid           = -1; /*Assigned by Finalize*/
   this->analysis_enum = node_analysis;
   this->clone         = node_clone;
   this->active        = true;
   this->freeze        = false;
 
   /*Initialize coord_system: Identity matrix by default*/
   for(int k=0;k<3;k++) for(int l=0;l<3;l++) this->coord_system[k][l]=0.0;
   for(int k=0;k<3;k++) this->coord_system[k][k]=1.0;
 
   this->approximation=0;
   if(analysis_enum==StressbalanceAnalysisEnum) this->approximation=in_approximation;
 
   /*indexing:*/
   this->indexingupdate = true;
   this->doftype        = NULL;
   Analysis *analysis = EnumToAnalysis(analysis_enum);
   this->gsize = analysis->DofsPerNode(&this->doftype,iomodel->domaintype,in_approximation);
   delete analysis;
 
   if(this->gsize>0){
      this->f_set          = xNew<bool>(this->gsize);
      this->s_set          = xNew<bool>(this->gsize);
      this->svalues        = xNew<IssmDouble>(this->gsize);
      this->gdoflist       = xNew<int>(this->gsize);
      this->gdoflist_local = xNew<int>(this->gsize);
      this->fsize          = -1;
      this->ssize          = -1;
      this->fdoflist       = NULL;
      this->sdoflist       = NULL;
      this->fdoflist_local = NULL;
      this->sdoflist_local = NULL;
   }
   else{
      this->f_set          = NULL;
      this->s_set          = NULL;
      this->svalues        = NULL;
      this->gdoflist       = NULL;
      this->gdoflist_local = NULL;
      this->fsize          = -1;
      this->ssize          = -1;
      this->fdoflist       = NULL;
      this->sdoflist       = NULL;
      this->fdoflist_local = NULL;
      this->sdoflist_local = NULL;
   }
 
   /*Assign values assuming no Dirichlet at this point*/
   for(int i=0;i<this->gsize;i++){
      this->f_set[i]    = true;
      this->s_set[i]    = false;
      this->svalues[i]  = 0.;
      this->gdoflist[i] = -1;
   }
 
   /*Stop here if AMR*/
   if(isamr) return;
 
   /*Stressbalance Horiz*/
   if(analysis_enum==StressbalanceAnalysisEnum){
 
      /*Coordinate system provided, convert to coord_system matrix*/
      _assert_(iomodel->Data("md.stressbalance.referential")); 
      XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data("md.stressbalance.referential")[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->domaintype!=Domain2DhorizontalEnum && iomodel->domaintype!=Domain3DsurfaceEnum){
         /*We have a  3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
         _assert_(iomodel->Data("md.mesh.vertexonbase")); 
         _assert_(iomodel->Data("md.flowequation.vertex_equation"));
         if(in_approximation==SSAApproximationEnum && !reCast<int>(iomodel->Data("md.mesh.vertexonbase")[io_index])){
            this->HardDeactivate();
         }
         if(in_approximation==L1L2ApproximationEnum && !reCast<int>(iomodel->Data("md.mesh.vertexonbase")[io_index])){
            this->HardDeactivate();
         }
         if(in_approximation==SSAHOApproximationEnum && reCast<int>(iomodel->Data("md.flowequation.borderSSA")[io_index])){
            if(!reCast<int>(iomodel->Data("md.mesh.vertexonbase")[io_index])){
               this->HardDeactivate();
            }
         }
         if(in_approximation==SSAFSApproximationEnum && reCast<int>(iomodel->Data("md.flowequation.borderSSA")[io_index])){
            if(!reCast<int>(iomodel->Data("md.mesh.vertexonbase")[io_index])){
               for(int k=0;k<=1;k++) this->FreezeDof(k);
            }
         }
      }
      /*spc all nodes on SIA*/
      if(in_approximation==SIAApproximationEnum){
         this->HardDeactivate();
      }
   }
 
   /*2d solutions in 3d, we need to constrain all the nodes that are not on base*/
   if(
            analysis_enum==FreeSurfaceBaseAnalysisEnum || 
            analysis_enum==MasstransportAnalysisEnum || 
            analysis_enum==MeltingAnalysisEnum || 
            analysis_enum==L2ProjectionBaseAnalysisEnum || 
            analysis_enum==BalancethicknessAnalysisEnum ||
            analysis_enum==HydrologyDCInefficientAnalysisEnum ||
            analysis_enum==HydrologyDCEfficientAnalysisEnum ||
            analysis_enum==GLheightadvectionAnalysisEnum ||
            analysis_enum==LevelsetAnalysisEnum
     ){
      if(iomodel->domaintype!=Domain2DhorizontalEnum & iomodel->domaintype!=Domain3DsurfaceEnum){
         /*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
         _assert_(iomodel->Data("md.mesh.vertexonbase"));
         if(!(reCast<bool>(iomodel->Data("md.mesh.vertexonbase")[io_index]))){
            this->HardDeactivate();
         }
      }
   }
   if(
            analysis_enum==FreeSurfaceTopAnalysisEnum
     ){
      if(iomodel->domaintype!=Domain2DhorizontalEnum){
         /*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
         _assert_(iomodel->Data("md.mesh.vertexonsurface"));
         if(!(reCast<bool>(iomodel->Data("md.mesh.vertexonsurface")[io_index]))){
            this->HardDeactivate();
         }
      }
   }
 
}
/*}}}*/
Node::~Node(){/*{{{*/
 
   if(this->f_set)          xDelete<bool>(f_set);
   if(this->s_set)          xDelete<bool>(s_set);
   if(this->svalues)        xDelete<IssmDouble>(svalues);
   if(this->doftype)        xDelete<int>(doftype);
   if(this->gdoflist)       xDelete<int>(gdoflist);
   if(this->fdoflist)       xDelete<int>(fdoflist);
   if(this->sdoflist)       xDelete<int>(sdoflist);
   if(this->gdoflist_local) xDelete<int>(gdoflist_local);
   if(this->fdoflist_local) xDelete<int>(fdoflist_local);
   if(this->sdoflist_local) xDelete<int>(sdoflist_local);
   return;
}
/*}}}*/
Object* Node::copy(void){/*{{{*/
 
   /*output: */
   Node* output=NULL;
 
   /*initalize output: */
   output=new Node();
 
   /*id: */
   output->id  = this->id;
   output->sid = this->sid;
   output->lid = this->lid;
   output->pid = this->pid;
   output->analysis_enum = this->analysis_enum;
   output->approximation = this->approximation;
 
   /*Initialize coord_system: */
   for(int k=0;k<3;k++) for(int l=0;l<3;l++) output->coord_system[k][l]=this->coord_system[k][l];
 
   /*indexing:*/
   output->indexingupdate = this->indexingupdate;
   output->gsize  = this->gsize;
   output->fsize  = this->fsize;
   output->ssize  = this->ssize;
   output->clone  = this->clone;
   output->active = this->active;
   output->freeze = this->freeze;
   if(output->gsize>0){
      output->f_set=xNew<bool>(output->gsize);
      output->s_set=xNew<bool>(output->gsize);
      output->svalues=xNew<IssmDouble>(output->gsize);
      if(this->doftype) output->doftype=xNew<int>(output->gsize);
      output->gdoflist=xNew<int>(output->gsize);
      output->gdoflist_local=xNew<int>(output->gsize);
   }
   if(output->fsize>0){
      output->fdoflist=xNew<int>(output->fsize);
      output->fdoflist_local=xNew<int>(output->fsize);
   }
   if(output->ssize>0){
      output->sdoflist=xNew<int>(output->ssize);
      output->sdoflist_local=xNew<int>(output->ssize);
   }
 
   if(output->gsize>0){
      memcpy(output->f_set,this->f_set,output->gsize*sizeof(bool));
      memcpy(output->s_set,this->s_set,output->gsize*sizeof(bool));
      xMemCpy<IssmDouble>(output->svalues,this->svalues,output->gsize);
      if(output->doftype)memcpy(output->doftype,this->doftype,output->gsize*sizeof(int));
      memcpy(output->gdoflist,this->gdoflist,output->gsize*sizeof(int));
      memcpy(output->gdoflist_local,this->gdoflist_local,output->gsize*sizeof(int));
   }
   if(output->fsize>0){
      memcpy(output->fdoflist,this->fdoflist,output->fsize*sizeof(int));
      memcpy(output->fdoflist_local,this->fdoflist_local,output->fsize*sizeof(int));
   }
   if(output->ssize>0){
      memcpy(output->sdoflist,this->sdoflist,output->ssize*sizeof(int));
      memcpy(output->sdoflist_local,this->sdoflist_local,output->ssize*sizeof(int));
   }
 
   return (Object*)output; 
}
/*}}}*/
void Node::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   MARSHALLING_ENUM(NodeEnum);
   MARSHALLING(id);
   MARSHALLING(sid);
   MARSHALLING(lid);
   MARSHALLING(pid);
   MARSHALLING(indexingupdate);
   MARSHALLING(analysis_enum);
   MARSHALLING_ARRAY(coord_system,IssmDouble,9);
 
   MARSHALLING(gsize);
   MARSHALLING(fsize);
   MARSHALLING(ssize);
   MARSHALLING(clone);
   MARSHALLING(active);
   MARSHALLING(freeze);
   MARSHALLING_DYNAMIC(f_set,bool,gsize);
   MARSHALLING_DYNAMIC(s_set,bool,gsize);
   MARSHALLING_DYNAMIC(svalues,IssmDouble,gsize);
   MARSHALLING_DYNAMIC(doftype,int,gsize);
   MARSHALLING_DYNAMIC(gdoflist,int,gsize);
   MARSHALLING_DYNAMIC(fdoflist,int,fsize);
   MARSHALLING_DYNAMIC(sdoflist,int,ssize);
   MARSHALLING_DYNAMIC(gdoflist_local,int,gsize);
   MARSHALLING_DYNAMIC(fdoflist_local,int,fsize);
   MARSHALLING_DYNAMIC(sdoflist_local,int,ssize);
} /*}}}*/
 
/*Object virtual functions definitions:*/
void Node::DeepEcho(void){/*{{{*/
 
   int i;
   _printf_("Node:\n");
   _printf_("   id: " << id << "\n");
   _printf_("   sid: " << sid << "\n");
   _printf_("   analysis_enum: " << EnumToStringx(analysis_enum) << "\n");
   _printf_("   approximation: " << EnumToStringx(approximation) << "\n");
   _printf_("   indexingupdate: " << indexingupdate << "\n");
   _printf_("   gsize:  " << gsize << "\n");
   _printf_("   fsize:  " << fsize << "\n");
   _printf_("   ssize:  " << ssize << "\n");
   _printf_("   clone:  " << clone << "\n");
   _printf_("   active: " << active << "\n");
   _printf_("   freeze: " << freeze << "\n");
   _printf_("   f_set = [ ");
   for(i=0;i<gsize;i++) _printf_((f_set[i]?1:0)<< " ");
   _printf_("]\n");
   _printf_("   s_set = [ ");
   for(i=0;i<gsize;i++) _printf_((s_set[i]?1:0)<< " ");
   _printf_("]\n");
   _printf_("   svalues (" << this->ssize << "): |");
   for(i=0;i<this->gsize;i++){
      if(this->s_set[i])_printf_(" " << svalues[i] << " |");
   }
   _printf_("\n");
   if(doftype){
      _printf_("   doftype: |");
      for(i=0;i<gsize;i++){
         _printf_(" " << doftype[i] << " |");
      }
      _printf_("\n");
   }
   else _printf_("   doftype: NULL\n");
 
   _printf_("   g_doflist (" << this->gsize << "): |");
   for(i=0;i<this->gsize;i++) _printf_(" " << gdoflist[i] << " |");
   _printf_("\n");
   _printf_("   g_doflist_local (" << this->gsize << "): |");
   for(i=0;i<this->gsize;i++) _printf_(" " << gdoflist_local[i] << " |");
   _printf_("\n");
 
   _printf_("   f_doflist (" << this->fsize << "): |");
   for(i=0;i<this->fsize;i++) _printf_(" " << fdoflist[i] << " |");
   _printf_("\n");
   _printf_("   f_doflist_local (" << this->fsize << "): |");
   for(i=0;i<this->fsize;i++) _printf_(" " << fdoflist_local[i] << " |");
   _printf_("\n");
 
   _printf_("   s_doflist (" << this->ssize << "): |");
   for(i=0;i<this->ssize;i++) _printf_(" " << sdoflist[i] << " |");
   _printf_("\n");
   _printf_("   s_doflist_local (" << this->ssize << "): |");
   for(i=0;i<this->ssize;i++) _printf_(" " << sdoflist_local[i] << " |");
   _printf_("\n");
 
}
/*}}}*/
void Node::Echo(void){/*{{{*/
 
   _printf_("Node:\n");
   _printf_("   id : " << id << "\n");
   _printf_("   sid: " << sid << "\n");
   _printf_("   lid: " << lid << "\n");
   _printf_("   pid: " << pid << "\n");
   _printf_("   analysis_enum: " << EnumToStringx(analysis_enum) << "\n");
   _printf_("   approximation: " << EnumToStringx(approximation) << "\n");
   _printf_("   indexingupdate: " << indexingupdate << "\n");
   _printf_("   gsize:  " << gsize << "\n");
   _printf_("   fsize:  " << fsize << "\n");
   _printf_("   ssize:  " << ssize << "\n");
   _printf_("   clone:  " << clone << "\n");
   _printf_("   active: " << active << "\n");
   _printf_("   freeze: " << freeze << "\n");
}
/*}}}*/
int  Node::Id(void){ return id; }/*{{{*/
/*}}}*/
int  Node::ObjectEnum(void){/*{{{*/
 
   return NodeEnum;
 
}
/*}}}*/
 
/*Node management:*/
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];
 
}
/*}}}*/
int  Node::GetDof(int dofindex,int setenum){/*{{{*/
 
   _assert_(!this->indexingupdate);
   if(setenum==GsetEnum){
      _assert_(dofindex>=0 && dofindex<gsize);
      return gdoflist[dofindex];
   }
   else if(setenum==FsetEnum){
      _assert_(dofindex>=0 && dofindex<fsize);
      return fdoflist[dofindex];
   }
   else if(setenum==SsetEnum){
      _assert_(dofindex>=0 && dofindex<ssize);
      return sdoflist[dofindex];
   }
   else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
 
} /*}}}*/
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->gsize;i++) outdoflist[i]=gdoflist[i];
      if(setenum==FsetEnum)for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist[i];
      if(setenum==SsetEnum)for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist[i];
   }
   else{
 
      if(setenum==GsetEnum){
         if(doftype){
            count=0;
            for(i=0;i<this->gsize;i++){
               if(doftype[i]==approximation_enum){
                  outdoflist[count]=gdoflist[i];
                  count++;
               }
            }
            _assert_(count); //at least one dof should be the approximation requested
         }
         else for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist[i];
      }
      else if(setenum==FsetEnum){
         if(doftype){
            count=0;
            count2=0;
            for(i=0;i<this->gsize;i++){
               if(f_set[i]){
                  if(doftype[i]==approximation_enum){
                     outdoflist[count]=fdoflist[count2];
                     count++;
                  }
                  count2++;
               }
            }
         }
         else for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist[i];
      }
      else if(setenum==SsetEnum){
         if(doftype){
            count=0;
            count2=0;
            for(i=0;i<this->gsize;i++){
               if(s_set[i]){
                  if(doftype[i]==approximation_enum){
                     outdoflist[count]=sdoflist[count2];
                     count++;
                  }
                  count2++;
               }
            }
         }
         else for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist[i];
      }
      else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
}
/*}}}*/
void Node::GetDofListLocal(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->gsize;i++) outdoflist[i]=gdoflist_local[i];
      if(setenum==FsetEnum)for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist_local[i];
      if(setenum==SsetEnum)for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist_local[i];
   }
   else{
 
      if(setenum==GsetEnum){
         if(doftype){
            count=0;
            for(i=0;i<this->gsize;i++){
               if(doftype[i]==approximation_enum){
                  outdoflist[count]=gdoflist_local[i];
                  count++;
               }
            }
            _assert_(count); //at least one dof should be the approximation requested
         }
         else for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist_local[i];
      }
      else if(setenum==FsetEnum){
         if(doftype){
            count=0;
            count2=0;
            for(i=0;i<this->gsize;i++){
               if(f_set[i]){
                  if(doftype[i]==approximation_enum){
                     outdoflist[count]=fdoflist_local[count2];
                     count++;
                  }
                  count2++;
               }
            }
         }
         else for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist_local[i];
      }
      else if(setenum==SsetEnum){
         if(doftype){
            count=0;
            count2=0;
            for(i=0;i<this->gsize;i++){
               if(s_set[i]){
                  if(doftype[i]==approximation_enum){
                     outdoflist[count]=sdoflist_local[count2];
                     count++;
                  }
                  count2++;
               }
            }
         }
         else for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist_local[i];
      }
      else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
}
/*}}}*/
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->gsize;i++) outdoflist[i]=i;
      else if(setenum==FsetEnum){
         count=0;
         for(i=0;i<this->gsize;i++){
            if(f_set[i]){
               outdoflist[count]=i;
               count++;
            }
         }
      }
      else if(setenum==SsetEnum){
         count=0;
         for(i=0;i<this->gsize;i++){
            if(s_set[i]){
               outdoflist[count]=i;
               count++;
            }
         }
      }
      else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
   else{
 
      if(setenum==GsetEnum){
         if(doftype){
            count=0;
            for(i=0;i<this->gsize;i++){
               if(doftype[i]==approximation_enum){
                  outdoflist[count]=count;
                  count++;
               }
            }
            _assert_(count);
         }
         else for(i=0;i<this->gsize;i++) outdoflist[i]=i;
      }
      else if(setenum==FsetEnum){
 
         if(doftype){
            count=0;
            count2=0;
            for(i=0;i<this->gsize;i++){
               if(doftype[i]==approximation_enum){
                  if(f_set[i]){
                     outdoflist[count]=count2;
                     count++;
                  }
                  count2++;
               }
            }
            _assert_(count2);
         }
         else{
 
            count=0;
            for(i=0;i<this->gsize;i++){
               if(f_set[i]){
                  outdoflist[count]=i;
                  count++;
               }
            }
         }
      }
      else if(setenum==SsetEnum){
         if(doftype){
            count=0;
            count2=0;
            for(i=0;i<this->gsize;i++){
               if(doftype[i]==approximation_enum){
                  if(s_set[i]){
                     outdoflist[count]=count2;
                     count++;
                  }
                  count2++;
               }
            }
            _assert_(count2);
         }
         else{
            count=0;
            for(i=0;i<this->gsize;i++){
               if(s_set[i]){
                  outdoflist[count]=i;
                  count++;
               }
            }
         }
      }
      else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
}
/*}}}*/
int  Node::Lid(void){/*{{{*/
   return lid; 
}
/*}}}*/
int  Node::Sid(void){/*{{{*/
   return sid; 
}
/*}}}*/
int  Node::Pid(void){/*{{{*/
   return this->pid; 
}
/*}}}*/
 
/*Node numerics:*/
void Node::Activate(void){/*{{{*/
 
   if(!IsActive() && !this->freeze){
      this->indexingupdate = true;
      this->active = true;
      for(int i=0;i<this->gsize;i++){
         this->f_set[i]    = true;
         this->s_set[i]    = false;
         this->svalues[i]  = 0.; 
      }
   }
 
}
/*}}}*/
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);
   this->svalues[dof]=value;
}
/*}}}*/
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->ssize){
      values=xNew<IssmDouble>(this->ssize);
      count=0;
      for(i=0;i<this->gsize;i++){
         if(this->s_set[i]){
            values[count]=this->svalues[i];
            _assert_(!xIsNan<IssmDouble>(values[count]));
            count++;
         }
      }
 
      /*Add values into constraint vector: */
      ys->SetValues(this->ssize,this->sdoflist,values,INS_VAL);
   }
 
   /*Free ressources:*/
   xDelete<IssmDouble>(values);
 
}
/*}}}*/
void Node::Deactivate(void){/*{{{*/
 
   if(IsActive() && !this->freeze){
      this->indexingupdate = true;
      this->active = false;
      /*Constrain to 0. at this point*/
      for(int i=0;i<this->gsize;i++){
         this->f_set[i]    = false;
         this->s_set[i]    = true;
         this->svalues[i]  = 0.; 
      }
   }
} /*}}}*/
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. Only do this for active nodes. */
   _assert_(dof<this->gsize);
   _assert_(this->active);
 
   if(this->f_set[dof] == 0){
      if(this->freeze) _error_("Cannot change dof of frozen node");
      this->indexingupdate = true;
      this->f_set[dof]=1; 
      this->s_set[dof]=0;
   }
}
/*}}}*/
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->gsize);
 
   if(this->f_set[dof] == 1){
      //if(this->freeze) _error_("Cannot change dof of frozen node");
      this->indexingupdate = true;
      this->f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
      this->s_set[dof]=1;
   }
}
/*}}}*/
void Node::FreezeDof(int dof){/*{{{*/
 
   DofInSSet(dof); //with 0 displacement for this dof.
   //FIXME: for now we don't want this element to change so we use freeze
   this->freeze =true;
 
}
/*}}}*/
int  Node::GetApproximation(){/*{{{*/
 
   return approximation;
}
/*}}}*/
void Node::SetApproximation(int in_approximation){/*{{{*/
 
   this->approximation = in_approximation;
}
/*}}}*/
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->gsize;
      else if (setenum==FsetEnum) numdofs=this->fsize;
      else if (setenum==SsetEnum) numdofs=this->ssize;
      else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
   else{
      if(setenum==GsetEnum){
         if(this->doftype){
            numdofs=0;
            for(i=0;i<this->gsize;i++){
               if(this->doftype[i]==approximation_enum) numdofs++;
            }
         }
         else numdofs=this->gsize;
      }
      else if (setenum==FsetEnum){
         if(this->doftype){
            numdofs=0;
            for(i=0;i<this->gsize;i++){
               if((this->doftype[i]==approximation_enum) && (this->f_set[i])) numdofs++;
            }
         }
         else numdofs=this->fsize;
      }
      else if (setenum==SsetEnum){
         if(this->doftype){
         numdofs=0;
            for(i=0;i<this->gsize;i++){
               if((this->doftype[i]==approximation_enum) && (this->s_set[i])) numdofs++;
            }
         }
         else numdofs=this->ssize;
      }
      else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
   return numdofs;
}
/*}}}*/
void Node::HardDeactivate(void){/*{{{*/
 
   this->Deactivate();
   this->freeze =true;
 
}
/*}}}*/
bool Node::IsActive(void){/*{{{*/
 
   return active;
 
}
/*}}}*/
int  Node::IsClone(){/*{{{*/
 
   return clone;
 
}
/*}}}*/
void Node::ReindexingDone(void){/*{{{*/
 
   this->indexingupdate = false;
 
}
/*}}}*/
void Node::RelaxConstraint(int dof){/*{{{*/
 
   /*Dof should be added to the f-set, and taken out of the s-set:*/
   if(this->active){
      DofInFSet(dof);
      this->svalues[dof]=0.;
   }
}
/*}}}*/
bool Node::RequiresDofReindexing(void){/*{{{*/
 
   return this->indexingupdate;
 
}
/*}}}*/
void Node::VecMerge(Vector<IssmDouble>* ug,IssmDouble* local_uf,int* indices_uf,IssmDouble* local_ys,int* indices_ys){/*{{{*/
 
   /*Only perform operation if not clone*/
   if(this->IsClone()) return;
 
   if(this->fsize){
      int*        indices = xNew<int>(this->fsize);
      IssmDouble* values  = xNew<IssmDouble>(this->fsize);
 
      int count = 0;
      for(int i=0;i<this->gsize;i++){
         if(this->f_set[i]){
            _assert_(local_uf);
            _assert_(this->fdoflist[count]==indices_uf[this->fdoflist_local[count]]);
 
            values[count] =local_uf[this->fdoflist_local[count]];
            indices[count]=this->gdoflist[i];
            count++;
         }
      }
      ug->SetValues(this->fsize,indices,values,INS_VAL);
 
      xDelete<IssmDouble>(values);
      xDelete<int>(indices);
   }
   if(this->ssize){
      int*        indices = xNew<int>(this->ssize);
      IssmDouble* values  = xNew<IssmDouble>(this->ssize);
 
      int count = 0;
      for(int i=0;i<this->gsize;i++){
         if(this->s_set[i]){
            _assert_(local_ys);
            _assert_(this->sdoflist[count]==indices_ys[this->sdoflist_local[count]]);
 
            values[count] =local_ys[this->sdoflist_local[count]];
            indices[count]=this->gdoflist[i];
            count++;
         }
      }
      ug->SetValues(this->ssize,indices,values,INS_VAL);
      
      xDelete<IssmDouble>(values);
      xDelete<int>(indices);
   }
}
/*}}}*/
void Node::VecReduce(Vector<IssmDouble>* uf, IssmDouble* local_ug,int* indices_ug){/*{{{*/
 
 
   /*Only perform operation if not clone*/
   if(this->IsClone()) return;
 
   if(this->fsize){
      int*        indices = xNew<int>(this->fsize);
      IssmDouble* values  = xNew<IssmDouble>(this->fsize);
 
      int count = 0;
      for(int i=0;i<this->gsize;i++){
         if(this->f_set[i]){
            _assert_(local_ug);
            _assert_(this->gdoflist[i]==indices_ug[this->gdoflist_local[i]]);
 
            values[count] =local_ug[this->gdoflist_local[i]];
            indices[count]=this->fdoflist[count];
            count++;
         }
      }
      uf->SetValues(this->fsize,indices,values,INS_VAL);
 
      xDelete<IssmDouble>(values);
      xDelete<int>(indices);
   }
}
/*}}}*/
 
/* indexing routines:*/
void Node::AllocateDofLists(int setenum){/*{{{*/
 
   /*Initialize: */
   int size=0;
 
   if(setenum==FsetEnum){
      for(int i=0;i<this->gsize;i++) if(f_set[i])size++;
      this->fsize=size;
      xDelete<int>(this->fdoflist);
      xDelete<int>(this->fdoflist_local);
 
      if(this->fsize){
         this->fdoflist=xNew<int>(size);
         this->fdoflist_local=xNew<int>(size);
      }
      else{
         this->fdoflist=NULL;
         this->fdoflist_local=NULL;
      }
   }
   else if(setenum==SsetEnum){
      for(int i=0;i<this->gsize;i++) if(s_set[i])size++;
      this->ssize=size;
      xDelete<int>(this->sdoflist);
      xDelete<int>(this->sdoflist_local);
 
      if(this->ssize){
         this->sdoflist=xNew<int>(size);
         this->sdoflist_local=xNew<int>(size);
      }
      else{
         this->sdoflist=NULL;
         this->sdoflist_local=NULL;
      }
   }
   /*TODO, we should never be here for GSET! add an assert and track down whether this happens*/
}
/*}}}*/
void Node::DistributeLocalDofs(int* pdofcount,int setenum){/*{{{*/
 
   /*Get current count*/
   int dofcount=*pdofcount;
 
   /*This node should distribute dofs for setenum set (eg, f_set or s_set), go ahead: */
   if(setenum==GsetEnum){
      _assert_(this->gsize==0 || this->gdoflist_local);
      for(int i=0;i<this->gsize;i++) gdoflist_local[i]=dofcount+i;
      dofcount+=this->gsize;
   }
   else if(setenum==FsetEnum){
      _assert_(this->fsize==0 || this->fdoflist_local);
      for(int i=0;i<this->fsize;i++) fdoflist_local[i]=dofcount+i;
      dofcount+=this->fsize;
   }
   else if(setenum==SsetEnum){
      _assert_(this->ssize==0 || this->sdoflist_local);
      for(int i=0;i<this->ssize;i++) sdoflist_local[i]=dofcount+i;
      dofcount+=this->ssize;
   }
   else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
 
   /*Assign output pointers: */
   *pdofcount=dofcount;
}
/*}}}*/
void Node::DistributeGlobalDofsMasters(int dofcount,int setenum){/*{{{*/
 
   /*This node is a clone, don't offset the dofs!: */
   if(clone) return;
 
   /*This node should off_set the dofs, go ahead: */
   if(setenum==GsetEnum){
      _assert_(this->gsize==0 || this->gdoflist);
      for(int i=0;i<this->gsize;i++) this->gdoflist[i]=this->gdoflist_local[i]+dofcount;
   }
   else if(setenum==FsetEnum){
      _assert_(this->fsize==0 || this->fdoflist);
      for(int i=0;i<this->fsize;i++) this->fdoflist[i]=this->fdoflist_local[i]+dofcount;
   }
   else if(setenum==SsetEnum){
      _assert_(this->ssize==0 || this->sdoflist);
      for(int i=0;i<this->ssize;i++) this->sdoflist[i]=this->sdoflist_local[i]+dofcount;
   }
   else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
}
/*}}}*/
void Node::ShowMasterDofs(int* truedofs,int setenum){/*{{{*/
 
   _assert_(!this->clone);
 
   /*Ok, we are not a clone, just plug our dofs into truedofs: */
   switch(setenum){
      case GsetEnum:
         for(int j=0;j<this->gsize;j++) truedofs[j]=gdoflist[j];
         break;
      case FsetEnum:
         for(int j=0;j<this->fsize;j++) truedofs[j]=fdoflist[j];
         break;
      case SsetEnum:
         for(int j=0;j<this->ssize;j++) truedofs[j]=sdoflist[j];
         break;
      default:
         _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
 
}
/*}}}*/
void Node::UpdateCloneDofs(int* alltruedofs,int setenum){/*{{{*/
 
   _assert_(this->clone);
 
   /*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(int j=0;j<this->gsize;j++) gdoflist[j]=alltruedofs[j];
         break;
      case FsetEnum:
         for(int j=0;j<this->fsize;j++) fdoflist[j]=alltruedofs[j];
         break;
      case SsetEnum:
         for(int j=0;j<this->ssize;j++) sdoflist[j]=alltruedofs[j];
         break;
      default:
         _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
   }
}
/*}}}*/
 
/*Methods inherent to Node: */
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* 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 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;
}
/*}}}*/