PentaInput2.cpp

Go to the documentation of this file.

 
#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 "./PentaInput2.h"
 
/*PentaInput2 constructors and destructor*/
PentaInput2::PentaInput2(void){/*{{{*/
 
   this->numberofelements_local = -1;
   this->numberofvertices_local = -1;
   this->isserved       = false;
   this->isserved_collapsed= 0;
   this->M = -1;
   this->N = -1;
   this->values         = NULL;
   this->element_values = NULL;
 
}/*}}}*/
PentaInput2::PentaInput2(int nbe_in,int nbv_in,int interp_in){/*{{{*/
 
   _assert_(nbe_in>0);
   _assert_(nbe_in<1e11);
   _assert_(nbv_in>0);
   _assert_(nbv_in<1e11);
   this->numberofelements_local = nbe_in;
   this->numberofvertices_local = nbv_in;
   this->isserved           = false;
   this->isserved_collapsed = 0;
 
   /*Reset takes care of the rest*/
   this->Reset(interp_in);
 
}/*}}}*/
PentaInput2::~PentaInput2(){/*{{{*/
   if(this->element_values) xDelete<IssmDouble>(this->element_values);
   if(this->values)         xDelete<IssmDouble>(this->values);
}
/*}}}*/
void PentaInput2::Reset(int interp_in){/*{{{*/
 
   /*Clean up*/
   if(this->values)         xDelete<IssmDouble>(this->values);
   if(this->element_values) xDelete<IssmDouble>(this->element_values);
 
   /*Set interpolation*/
   this->interpolation  = interp_in;
 
   /*Create Sizes*/
   if(this->interpolation==P1Enum){
      this->M = this->numberofvertices_local;
      this->N = 1;
   }
   else{
      this->M = this->numberofelements_local;
      this->N = PentaRef::NumberofNodes(interp_in);
   }
 
   /*Allocate Pointers*/
   this->values         = xNewZeroInit<IssmDouble>(this->M*this->N);
   this->element_values = xNewZeroInit<IssmDouble>(PentaRef::NumberofNodes(interp_in));
}/*}}}*/
 
/*Object virtual functions definitions:*/
Input2* PentaInput2::copy() {/*{{{*/
 
   /*Create output*/
   PentaInput2* output = new PentaInput2(this->numberofelements_local,this->numberofvertices_local,this->interpolation);
 
   /*Copy values*/
   xMemCpy<IssmDouble>(output->values,this->values,this->M*this->N);
 
   /*Return output*/
   return output;
 
}
/*}}}*/
void PentaInput2::DeepEcho(void){/*{{{*/
   _printf_("PentaInput2 Echo:\n");
   _printf_("   interpolation:      "<<EnumToStringx(this->interpolation)<<"\n");
   _printf_("   nbe_local:          "<<this->numberofvertices_local<<"\n");
   _printf_("   nbv_local:          "<<this->numberofelements_local<<"\n");
   _printf_("   Size:               "<<M<<"x"<<N<<"\n");
   _printf_("   isserved:           "<<(isserved?"true":"false") << "\n");
   _printf_("   isserved_collapsed: "<<isserved_collapsed << "\n");
   if(isserved){
      _printf_("   current values:      ");
      if(isserved_collapsed){
         _printf_("[ ");
         for(int i=0;i<3;i++) _printf_(" "<<this->element_values[i]);
         _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
      }
      else{
         _printf_("[ ");
         for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) _printf_(" "<<this->element_values[i]);
         _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
      }
   }
}
/*}}}*/
void PentaInput2::Echo(void){/*{{{*/
   _printf_(setw(15)<<"   PentaInput "<<setw(25)<<left<<EnumToStringx(-1));
   if(isserved){
      if(isserved_collapsed){
         _printf_("[ ");
         for(int i=0;i<3;i++) _printf_(" "<<this->element_values[i]);
         _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
      }
      else{
         _printf_("[ ");
         for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) _printf_(" "<<this->element_values[i]);
         _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
      }
   }
}
/*}}}*/
int  PentaInput2::Id(void){/*{{{*/
   return -1;
}/*}}}*/
void PentaInput2::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   MARSHALLING_ENUM(PentaInput2Enum);
   MARSHALLING(this->numberofelements_local);
   MARSHALLING(this->numberofvertices_local);
   MARSHALLING(this->interpolation);
   MARSHALLING(this->M);
   MARSHALLING(this->N);
   this->isserved = false;
   this->isserved_collapsed = 0;
   if(this->M*this->N){
      MARSHALLING_DYNAMIC(this->values,IssmDouble,this->M*this->N);
   }
   else this->values = NULL;
 
   if(marshall_direction == MARSHALLING_BACKWARD){
      this->element_values = xNewZeroInit<IssmDouble>(PentaRef::NumberofNodes(this->interpolation));
   }
}
/*}}}*/
int  PentaInput2::ObjectEnum(void){/*{{{*/
   return PentaInput2Enum;
}
/*}}}*/
 
/*PentaInput2 management*/
void PentaInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
 
   _assert_(this);
   _assert_(row>=0);
   _assert_(row<this->M);
   _assert_(this->N==1);
 
   this->values[row] = value_in;
   this->isserved    = false;
}
/*}}}*/
void PentaInput2::SetInput(int interp_in,int numindices,int* indices,IssmDouble* values_in){/*{{{*/
 
   _assert_(this);
   if(interp_in==P1Enum && this->interpolation==P1Enum){
      _assert_(this->N==1);
      for(int i=0;i<numindices;i++){
         int row = indices[i];
         _assert_(row>=0);
         _assert_(row<this->M);
         this->values[row] = values_in[i];
      }
   }
   else if(this->interpolation!=P1Enum && interp_in==P1Enum){
      this->Reset(interp_in);
      for(int i=0;i<numindices;i++){
         int row = indices[i];
         _assert_(row>=0);
         _assert_(row<this->M);
         this->values[row] = values_in[i];
      }
   }
   else{
      _error_("not supported");
   }
 
   this->isserved    = false;
}
/*}}}*/
void PentaInput2::SetInput(int interp_in,int row,int numindices,IssmDouble* values_in){/*{{{*/
 
   _assert_(this);
   if(interp_in==this->interpolation){
      _assert_(this->N==numindices);
   }
   else{
      this->Reset(interp_in);
      _assert_(this->N==numindices);
   }
   for(int i=0;i<numindices;i++) this->values[row*this->N+i] = values_in[i];
 
   this->isserved    = false;
}
/*}}}*/
void PentaInput2::Serve(int numindices,int* indices){/*{{{*/
 
   _assert_(this);
   _assert_(this->N==1);
 
   for(int i=0;i<numindices;i++){
      int row = indices[i];
      _assert_(row>=0);
      _assert_(row<this->M);
      this->element_values[i] = this->values[row];
   }
 
   /*Set input as served*/
   this->isserved = true;
   this->isserved_collapsed = 0;
}
/*}}}*/
void PentaInput2::Serve(int row,int numindices){/*{{{*/
 
   _assert_(this);
   _assert_(this->N==numindices);
   _assert_(row<this->M);
   _assert_(row>=0);
 
   for(int i=0;i<numindices;i++){
      this->element_values[i] = this->values[row*this->N+i];
   }
 
   /*Set input as served*/
   this->isserved = true;
   this->isserved_collapsed = 0;
}/*}}}*/
void PentaInput2::ServeCollapsed(int row,int state){/*{{{*/
 
   _assert_(this);
   _assert_(this->N>=3);
   _assert_(row<this->M);
   _assert_(row>=0);
 
   if(state==1){
      for(int i=0;i<3;i++) this->element_values[i] = this->values[row*this->N+i];
      for(int i=3;i<6;i++) this->element_values[i] = 0.;
   }
   else if(state==2){
      for(int i=0;i<3;i++) this->element_values[i] = this->values[row*this->N+3+i];
      for(int i=3;i<6;i++) this->element_values[i] = 0.;
   }
   else{
      _error_("not supported");
   }
 
   /*Set input as served*/
   this->isserved = true;
   this->isserved_collapsed = state;
}/*}}}*/
void PentaInput2::SetServeCollapsed(int state){/*{{{*/
   this->isserved_collapsed = state;
}/*}}}*/
int  PentaInput2::GetInterpolation(){/*{{{*/
   return this->interpolation;
}/*}}}*/
void PentaInput2::GetInputAverage(IssmDouble* pvalue){/*{{{*/
   _assert_(this);
   _assert_(this->isserved);
 
   /*Output*/
   IssmDouble value = 0.;
 
   if(this->isserved_collapsed){
      if(this->interpolation==P0Enum){
         value = this->element_values[0];
      }
      else{
         /*Assume P1...*/
         value = 1./3.*(this->element_values[0] +  this->element_values[1] +  this->element_values[2]);
      }
   }
   else{
      int        numnodes  = this->NumberofNodes(this->interpolation);
      IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
 
      for(int i=0;i<numnodes;i++) value+=this->element_values[i];
      value = value/numnodesd;
   }
 
   *pvalue=value;
}/*}}}*/
IssmDouble PentaInput2::GetInputMin(void){/*{{{*/
   _assert_(this);
   _assert_(this->isserved);
 
   int  numnodes  = this->NumberofNodes(this->interpolation);
   if(this->isserved_collapsed) numnodes = 3;
   IssmDouble min=this->element_values[0];
 
   for(int i=1;i<numnodes;i++){
      if(this->element_values[i]<min) min=this->element_values[i];
   }
   return min;
}/*}}}*/
IssmDouble PentaInput2::GetInputMax(void){/*{{{*/
   _assert_(this);
   _assert_(this->isserved);
 
   int numnodes  = this->NumberofNodes(this->interpolation);
   if(this->isserved_collapsed) numnodes = 3;
   IssmDouble max=this->element_values[0];
 
   for(int i=1;i<numnodes;i++){
      if(this->element_values[i]>max) max=this->element_values[i];
   }
   return max;
}/*}}}*/
IssmDouble PentaInput2::GetInputMaxAbs(void){/*{{{*/
   _assert_(this);
   _assert_(this->isserved);
 
   int numnodes  = this->NumberofNodes(this->interpolation);
   if(this->isserved_collapsed) numnodes = 3;
   IssmDouble maxabs=fabs(this->element_values[0]);
 
   for(int i=1;i<numnodes;i++){
      if(fabs(this->element_values[i])>maxabs) maxabs=fabs(this->element_values[i]);
   }
   return maxabs;
}/*}}}*/
void PentaInput2::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
   _assert_(this);
   _assert_(this->isserved);
   if(this->isserved_collapsed){
      _assert_(gauss->Enum()==GaussTriaEnum);
      if(this->interpolation==P0Enum){
         derivativevalues[0] = 0.;
         derivativevalues[1] = 0.;
      }
      else{
         TriaRef temp;
         temp.GetInputDerivativeValue(derivativevalues,this->element_values,xyz_list,(GaussTria*)gauss,P1Enum);
      }
   }
   else{
      _assert_(gauss->Enum()==GaussPentaEnum);
      PentaRef::GetInputDerivativeValue(derivativevalues,this->element_values,xyz_list,(GaussPenta*)gauss,this->interpolation);
   }
}/*}}}*/
void PentaInput2::GetInputValue(IssmDouble* pvalue,Gauss* gauss){/*{{{*/
   _assert_(this);
   _assert_(this->isserved);
   if(this->isserved_collapsed){
      _assert_(gauss->Enum()==GaussTriaEnum);
      if(this->interpolation==P0Enum){
         *pvalue = this->element_values[0];
      }
      else{
         TriaRef temp;
         temp.GetInputValue(pvalue,this->element_values,(GaussTria*)gauss,P1Enum);
      }
   }
   else{
      _assert_(gauss->Enum()==GaussPentaEnum);
      PentaRef::GetInputValue(pvalue,this->element_values,(GaussPenta*)gauss,this->interpolation);
   }
}/*}}}*/
int  PentaInput2::GetResultArraySize(void){/*{{{*/
   return 1;
}
/*}}}*/
int  PentaInput2::GetResultInterpolation(void){/*{{{*/
   if(this->interpolation==P0Enum || this->interpolation==P0DGEnum){
      return P0Enum;
   }
   return P1Enum;
}/*}}}*/
int  PentaInput2::GetResultNumberOfNodes(void){/*{{{*/
   return this->N;
}
/*}}}*/
void PentaInput2::Scale(IssmDouble alpha){/*{{{*/
 
   for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*this->values[i];
   for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
}
/*}}}*/
void PentaInput2::Pow(IssmDouble alpha){/*{{{*/
 
   for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
   for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
}
/*}}}*/
void PentaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
 
   /*xinput is of the same type, so cast it: */
   if(xinput->ObjectEnum()!=PentaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
   PentaInput2* xpentainput=xDynamicCast<PentaInput2*>(xinput);
   if(xpentainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 
   /*Carry out the AXPY operation depending on type:*/
   for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xpentainput->values[i] + this->values[i];
   for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xpentainput->element_values[i] + this->element_values[i];
}
/*}}}*/
void PentaInput2::PointWiseMult(Input2* xinput){/*{{{*/
 
   /*xinput is of the same type, so cast it: */
   if(xinput->ObjectEnum()!=PentaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
   PentaInput2* xpentainput=xDynamicCast<PentaInput2*>(xinput);
   if(xpentainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 
   /* we need to check that the vector sizes are identical*/
   if(xpentainput->M!=this->M||xpentainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
 
   /*Carry out the pointwise operation depending on type:*/
   for(int i=0;i<this->M*this->N;i++) this->values[i] = xpentainput->values[i] * this->values[i];
   for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xpentainput->element_values[i] * this->element_values[i];
}
/*}}}*/
 
/*Object functions*/