Vertex.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 <string.h>
#include "classes.h"
#include "shared/shared.h"
/*}}}*/
 
/*Vertex constructors and destructor:*/
Vertex::Vertex(){/*{{{*/
   return;
}
/*}}}*/
Vertex::Vertex(int vertex_id, int vertex_sid,bool vertex_clone, IoModel* iomodel,bool isamr){/*{{{*/
 
   /*Checks in debugging mode*/
   _assert_(vertex_sid>=0 && vertex_sid<iomodel->numberofvertices);
 
   /*IDs*/
   this->id    = vertex_id;
   this->sid   = vertex_sid;
   this->pid   = -1; /*Assigned later*/
   this->lid   = -1; /*Assigned later*/
   this->clone = vertex_clone;
 
   /*Properties from iomodel*/
   _assert_(iomodel->numbernodetoelementconnectivity);
   this->connectivity = iomodel->numbernodetoelementconnectivity[vertex_sid];
   this->domaintype   = iomodel->domaintype;
 
   /*Coordinates, only if not AMR*/
   if(!isamr){
      _assert_(iomodel->Data("md.mesh.x") && iomodel->Data("md.mesh.y") && iomodel->Data("md.mesh.z"));
      this->x            = iomodel->Data("md.mesh.x")[vertex_sid];
      this->y            = iomodel->Data("md.mesh.y")[vertex_sid];
      this->z            = iomodel->Data("md.mesh.z")[vertex_sid];
      if(iomodel->Data("md.mesh.lat") && iomodel->Data("md.mesh.long")){
         this->latitute     = iomodel->Data("md.mesh.lat")[vertex_sid];
         this->longitude    = iomodel->Data("md.mesh.long")[vertex_sid];
      }
 
      switch(iomodel->domaintype){
         case Domain3DEnum:
            _assert_(iomodel->Data("md.geometry.base") && iomodel->Data("md.geometry.thickness"));
            this->sigma = (iomodel->Data("md.mesh.z")[vertex_sid]-iomodel->Data("md.geometry.base")[vertex_sid])/(iomodel->Data("md.geometry.thickness")[vertex_sid]);
            break;
         case Domain3DsurfaceEnum:
            _assert_(iomodel->Data("md.mesh.lat") && iomodel->Data("md.mesh.long") && iomodel->Data("md.mesh.r"));
            this->latitute     = iomodel->Data("md.mesh.lat")[vertex_sid];
            this->longitude    = iomodel->Data("md.mesh.long")[vertex_sid];
            this->R            = iomodel->Data("md.mesh.r")[vertex_sid];
            break;
         case Domain2DhorizontalEnum:
            this->sigma = 0.;
            break;
         case Domain2DverticalEnum:
            _assert_(iomodel->Data("md.geometry.base") && iomodel->Data("md.geometry.thickness"));
            this->sigma = (iomodel->Data("md.mesh.y")[vertex_sid]-iomodel->Data("md.geometry.base")[vertex_sid])/(iomodel->Data("md.geometry.thickness")[vertex_sid]);
            break;
      }
   }
   else{
      this->x         = 0.;
      this->y         = 0.;
      this->z         = 0.;
      this->latitute  = 0.;
      this->longitude = 0.;
      this->R         = 0.;
      this->sigma     = 0.;
   }
 
}/*}}}*/
Vertex::~Vertex(){/*{{{*/
   return;
}
/*}}}*/
 
/*Object virtual functions definitions:*/
Object* Vertex::copy() {/*{{{*/
 
   return new Vertex(*this); 
 
}
/*}}}*/
void Vertex::DeepEcho(void){/*{{{*/
   this->Echo();
}
/*}}}*/
void Vertex::Echo(void){/*{{{*/
 
   _printf_("Vertex:\n");
   _printf_("   id: " << id << "\n");
   _printf_("   sid: " << sid << "\n");
   _printf_("   pid: " << pid << "\n");
   _printf_("   lid: " << lid << "\n");
   _printf_("   x: " << x << "\n");
   _printf_("   y: " << y << "\n");
   _printf_("   z: " << z << "\n");
   _printf_("   sigma: " << sigma << "\n");
   _printf_("   connectivity: " << connectivity << "\n");
   _printf_("   clone: " << clone << "\n");
 
   return;
}
/*}}}*/
int Vertex::Id(void){ return id; }/*{{{*/
/*}}}*/
void Vertex::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   MARSHALLING_ENUM(VertexEnum);
   MARSHALLING(clone);
   MARSHALLING(domaintype);
   MARSHALLING(id);
   MARSHALLING(sid);
   MARSHALLING(pid);
   MARSHALLING(lid);
   MARSHALLING(x);
   MARSHALLING(y);
   MARSHALLING(z);
   MARSHALLING(sigma);
   MARSHALLING(connectivity);
 
}
/*}}}*/
int Vertex::ObjectEnum(void){/*{{{*/
   return VertexEnum;
}/*}}}*/
 
/*Vertex management: */
int        Vertex::Connectivity(void){return connectivity;}/*{{{*/
/*}}}*/
IssmDouble Vertex::GetLatitude(){/*{{{*/
   return this->latitute;
}
/*}}}*/
IssmDouble Vertex::GetLongitude(){/*{{{*/
   return this->longitude;
}
/*}}}*/
IssmDouble Vertex::GetRadius(){/*{{{*/
   return this->R;
}
/*}}}*/
IssmDouble Vertex::GetX(){/*{{{*/
   return this->x;
}
/*}}}*/
IssmDouble Vertex::GetY(){/*{{{*/
   return this->y;
}
/*}}}*/
IssmDouble Vertex::GetZ(){/*{{{*/
   return this->z;
}
/*}}}*/
int        Vertex::Pid(void){ return pid; }/*{{{*/
/*}}}*/
int        Vertex::Lid(void){ return lid; }/*{{{*/
/*}}}*/
int        Vertex::Sid(void){ return sid; }/*{{{*/
/*}}}*/
void       Vertex::UpdatePosition(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* surface,IssmDouble* bed){/*{{{*/
 
   IssmDouble oldy,newy,vely;
   IssmDouble oldz,newz,velz;
   IssmDouble dt;
 
   /*Get time stepping*/
   parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 
   /*sigma remains constant. z=bed+sigma*thickness*/
   switch(this->domaintype){
      case Domain2DhorizontalEnum:
         /*Nothing*/
         return;
      case Domain2DverticalEnum:
         oldy = this->y;
         newy = bed[this->pid]+sigma*(surface[this->pid] - bed[this->pid]);
         vely = (newy-oldy)/dt;
         this->y = newy;
         vy->SetValue(this->pid,vely,INS_VAL);
         _assert_(!xIsNan<IssmDouble>(vely));
         return;
      case Domain3DEnum:
         oldz = this->z;
         newz = bed[this->pid]+sigma*(surface[this->pid] - bed[this->pid]);
         velz = (newz-oldz)/dt;
         this->z = newz;
         vz->SetValue(this->pid,velz,INS_VAL);
         _assert_(!xIsNan<IssmDouble>(velz));
         return;
      default:
         _error_("not implemented");
   }
}
/*}}}*/
void       Vertex::VertexCoordinates(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz, bool spherical){/*{{{*/
 
   if(this->clone==true) return;
 
   if(!spherical){
      vx->SetValue(this->sid,this->x,INS_VAL);
      vy->SetValue(this->sid,this->y,INS_VAL);
      vz->SetValue(this->sid,this->z,INS_VAL);
   }
   else{
      vx->SetValue(this->sid,this->latitute,INS_VAL);
      vy->SetValue(this->sid,this->longitude,INS_VAL);
      vz->SetValue(this->sid,this->R,INS_VAL);
   }
 
   return;
}
/*}}}*/
 
/*Methods relating to Vertex, but not internal methods: */
void GetVerticesCoordinates(IssmDouble* xyz,Vertex** vertices, int numvertices,bool spherical){ /*{{{*/
 
   _assert_(vertices);
   _assert_(xyz);
 
   if(!spherical){
      for(int i=0;i<numvertices;i++) {
         xyz[i*3+0]=vertices[i]->GetX();
         xyz[i*3+1]=vertices[i]->GetY();
         xyz[i*3+2]=vertices[i]->GetZ();
      }
   }
   else{
      for(int i=0;i<numvertices;i++) {
         xyz[i*3+0]=vertices[i]->GetLatitude();
         xyz[i*3+1]=vertices[i]->GetLongitude();
         xyz[i*3+2]=vertices[i]->GetRadius();
      }
   }
}/*}}}*/