TransientArrayParam.cpp

Go to the documentation of this file.

 
/*header 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"
/*}}}*/
 
/*TransientArrayParam constructors and destructor*/
TransientArrayParam::TransientArrayParam(){/*{{{*/
   return;
}
/*}}}*/
TransientArrayParam::TransientArrayParam(int in_enum_type,IssmDouble* in_values,IssmDouble* in_time,bool interpolation_on,int in_N,int in_M){/*{{{*/
 
   _assert_(in_values && in_time);
 
   this->enum_type=in_enum_type;
   this->M=in_M; //Number of rows
   this->N=in_N; //Number of timesteps
   this->interpolation=interpolation_on;
 
   this->values=xNew<IssmDouble>(M*N);
   xMemCpy<IssmDouble>(values,in_values,M*N);
 
   this->timesteps=xNew<IssmDouble>(N);
   xMemCpy<IssmDouble>(timesteps,in_time,N);
}
/*}}}*/
TransientArrayParam::~TransientArrayParam(){/*{{{*/
   xDelete<IssmDouble>(values);
   xDelete<IssmDouble>(timesteps);
}/*}}}*/
 
/*Object virtual functions definitions:*/
Param* TransientArrayParam::copy() {/*{{{*/
 
   return new TransientArrayParam(this->enum_type,this->values,this->timesteps,this->interpolation,this->M,this->N);
 
}
/*}}}*/
void TransientArrayParam::DeepEcho(void){/*{{{*/
 
   _printf_("TransientArrayParam:\n");
   _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
   _printf_("   number of time steps: " << this->N << "\n");
   _printf_("   number of rows: " << this->M << "\n");
   for(int i=0;i<this->N;i++){
      _printf_("  time: " << this->timesteps[i] << "\n");
      for(int k=0;k<this->M;k++){
         _printf_("     values: " << this->values[k*N+i] << "\n");
      }
      _printf_("\n");
   }
}
/*}}}*/
void TransientArrayParam::Echo(void){/*{{{*/
 
   _printf_("TransientArrayParam:\n");
   _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
   _printf_("   size: " << this->M << " by " << this->N << "\n");
 
}
/*}}}*/
int  TransientArrayParam::Id(void){ return -1; }/*{{{*/
/*}}}*/
void TransientArrayParam::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   MARSHALLING_ENUM(TransientArrayParamEnum);
 
   MARSHALLING(enum_type);
   MARSHALLING(interpolation);
   MARSHALLING(M);
   MARSHALLING(N);
   if(marshall_direction==MARSHALLING_BACKWARD){
      values    = xNew<IssmDouble>(M*N);
      timesteps = xNew<IssmDouble>(N);
   }
   MARSHALLING_ARRAY(values,IssmDouble,M*N);
   MARSHALLING_ARRAY(timesteps,IssmDouble,N);
 
}/*}}}*/
int  TransientArrayParam::ObjectEnum(void){/*{{{*/
 
   return TransientArrayParamEnum;
 
}/*}}}*/
 
/*TransientArrayParam virtual functions definitions: */
void  TransientArrayParam::GetParameterValue(IssmDouble* pdouble,int row,IssmDouble time){/*{{{*/
 
   IssmDouble output;
   bool       found;
   _assert_(row>=0 && row<this->M); 
 
   /*Ok, we have the time and row, go through the timesteps, and figure out which interval we 
    *fall within. Then interpolate the values on this interval: */
   if(time<this->timesteps[0]){
      /*get values for the first time: */
      output=this->values[row*this->N];
      found=true;
   }
   else if(time>this->timesteps[this->N-1]){
      /*get values for the last time: */
      output=this->values[(row+1)*this->N-1];
      found=true;
   }
   else{
      /*Find which interval we fall within: */
      for(int i=0;i<this->N;i++){
         if(time==this->timesteps[i]){
            /*We are right on one step time: */
            output = this->values[row*this->N+i];
            found=true;
            break; //we are done with the time interpolation.
         }
         else{
            if(this->timesteps[i]<time && time<this->timesteps[i+1]){
               /*ok, we have the interval [i:i+1]. Interpolate linearly for now: */
               IssmDouble deltat = this->timesteps[i+1]-this->timesteps[i];
               IssmDouble alpha  = (time-this->timesteps[i])/deltat;
               if(interpolation==true) output=(1.0-alpha)*this->values[row*this->N+i] + alpha*this->values[row*this->N+i+1];
               else output=this->values[row*this->N+i];
               found=true;
               break;
            }
            else continue; //keep looking on the next interval
         }
      }
   }
   if(!found)_error_("did not find time interval on which to interpolate values");
   *pdouble=output;
}
/*}}}*/