Penta.cpp

Go to the documentation of this file.

 
/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include "../classes.h"
#include "../Inputs2/PentaInput2.h"
#include "../Inputs2/ControlInput2.h"
#include "../Inputs2/TransientInput2.h"
#include "../Inputs2/DatasetInput2.h"
#include "../../shared/shared.h"
/*}}}*/
 
/*Element macros*/
#define NUMVERTICES   6
#define NUMVERTICES2D 3
 
/*Constructors/destructor/copy*/
Penta::~Penta(){/*{{{*/
   this->parameters=NULL;
}
/*}}}*/
Penta::Penta(int penta_id, int penta_sid,int penta_lid,IoModel* iomodel,int nummodels)/*{{{*/
   :ElementHook(nummodels,penta_id,NUMVERTICES,iomodel){
 
   int penta_elements_ids[2];
 
   /*Checks in debugging mode*/
   _assert_(iomodel->Data("md.mesh.upperelements"));
   _assert_(iomodel->Data("md.mesh.lowerelements"));
 
   /*id: */
   this->id  = penta_id;
   this->sid = penta_sid;
   this->lid = penta_lid;
 
   /*surface and base*/
   this->isonsurface = false;
   this->isonbase    = false;
 
   /*Build neighbors list*/
   if (xIsNan<IssmDouble>(iomodel->Data("md.mesh.upperelements")[penta_sid]) || iomodel->Data("md.mesh.upperelements")[penta_sid]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta
   else                                    penta_elements_ids[1]=reCast<int,IssmDouble>((iomodel->Data("md.mesh.upperelements")[penta_sid]));
   if (xIsNan<IssmDouble>(iomodel->Data("md.mesh.lowerelements")[penta_sid]) || iomodel->Data("md.mesh.lowerelements")[penta_sid]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta
   else                                    penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data("md.mesh.lowerelements")[penta_sid]));
   this->InitHookNeighbors(penta_elements_ids);
 
   //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
   this->parameters=NULL;
 
   /*initialize pointers:*/
   this->nodes             = NULL;
   this->vertices          = NULL;
   this->material          = NULL;
   this->verticalneighbors = NULL;
 
   /*Only allocate pointer*/
   this->element_type_list=xNew<int>(nummodels);
 
   /*surface and base*/
   _assert_(iomodel->Data("md.mesh.vertexonsurface"));
   _assert_(iomodel->Data("md.mesh.vertexonbase"));
   this->isonsurface = false;
   this->isonbase    = false;
   IssmDouble sum = 0.;
   for(int i=0;i<NUMVERTICES;i++) sum += iomodel->Data("md.mesh.vertexonsurface")[reCast<int>(iomodel->elements[(penta_id-1)*NUMVERTICES+i])-1];
   _assert_(sum>=0 && sum<4);
   if(sum>2.5) this->isonsurface = true;
   sum = 0.;
   for(int i=0;i<NUMVERTICES;i++) sum += iomodel->Data("md.mesh.vertexonbase")[reCast<int>(iomodel->elements[(penta_id-1)*NUMVERTICES+i])-1];
   _assert_(sum>=0 && sum<4);
   if(sum>2.5) this->isonbase = true;
}
/*}}}*/
Object* Penta::copy() {/*{{{*/
 
   int i;
   Penta* penta=NULL;
 
   penta=new Penta();
 
   //deal with PentaRef mother class
   int nanalyses = this->numanalyses;
   if(nanalyses > 0){
      penta->element_type_list=xNew<int>(nanalyses);
      for(i=0;i<nanalyses;i++) {
         if (this->element_type_list[i]) penta->element_type_list[i]=this->element_type_list[i];
         else penta->element_type_list[i] = 0;
      }
   }
   else penta->element_type_list = NULL;
   penta->element_type=this->element_type;
   penta->numanalyses=nanalyses;
 
   //deal with ElementHook mother class
   if (this->hnodes){
      penta->hnodes=xNew<Hook*>(penta->numanalyses);
      for(i=0;i<penta->numanalyses;i++){
         if (this->hnodes[i]) penta->hnodes[i] = (Hook*)(this->hnodes[i]->copy());
         else penta->hnodes[i] = NULL;
      }
   }
   else penta->hnodes = NULL;
 
   penta->hvertices = (Hook*)this->hvertices->copy();
   penta->hmaterial = (Hook*)this->hmaterial->copy();
   if (this->hneighbors) penta->hneighbors = (Hook*)(this->hneighbors->copy());
   else penta->hneighbors = NULL;
 
   /*deal with Tria fields: */
   penta->id  = this->id;
   penta->sid = this->sid;
   penta->lid = this->lid;
   penta->isonbase  = this->isonbase;
   penta->isonsurface  = this->isonsurface;
 
   /*point parameters: */
   penta->parameters=this->parameters;
 
   /*recover objects: */
   if (this->nodes) {
      unsigned int num_nodes = 6;
      penta->nodes = xNew<Node*>(num_nodes); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes.
      for(i=0;i<num_nodes;i++) if(this->nodes[i]) penta->nodes[i]=this->nodes[i]; else penta->nodes[i] = NULL;
   }
   else penta->nodes = NULL;
 
   penta->vertices = (Vertex**)this->hvertices->deliverp();
   penta->material = (Material*)this->hmaterial->delivers();
   penta->verticalneighbors = (Penta**)this->hneighbors->deliverp();
 
   return penta;
 
}
/*}}}*/
void Penta::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   MARSHALLING_ENUM(PentaEnum);
   MARSHALLING(this->isonsurface);
   MARSHALLING(this->isonbase);
 
   /*Call parent classes: */
   ElementHook::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   Element::MarshallElement(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses);
   PentaRef::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
 
   vertices = (Vertex**)this->hvertices->deliverp();
   material = (Material*)this->hmaterial->delivers();
   verticalneighbors = (Penta**)this->hneighbors->deliverp();
 
}
/*}}}*/
 
/*Other*/
void       Penta::AddBasalInput2(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
 
   _assert_(this->inputs2);
   if(!IsOnBase()) return;
   else{
      if(interpolation_enum==P1Enum || interpolation_enum==P1DGEnum){
         IssmDouble extrudedvalues[NUMVERTICES];
         for(int i=0;i<NUMVERTICES2D;i++){
            extrudedvalues[i]=values[i];
            extrudedvalues[i+NUMVERTICES2D]=values[i];
         }
         Penta* penta=this;
         for(;;){
            penta->AddInput2(input_enum,&extrudedvalues[0],interpolation_enum);
            if (penta->IsOnSurface()) break;
            penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
         }
      }
      else _error_("not implemented yet");
   }
 
}
/*}}}*/
void       Penta::AddInput2(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
 
   
   int vertexlids[NUMVERTICES];
 
   /*Call inputs method*/
   _assert_(this->inputs2);
   switch(interpolation_enum){
      case P1Enum:
         for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;
         inputs2->SetPentaInput(input_enum,interpolation_enum,NUMVERTICES,vertexlids,values);
         break;
      case P1DGEnum:
         inputs2->SetPentaInput(input_enum,interpolation_enum,this->lid,NUMVERTICES,values);
         break;
      default:
         inputs2->SetPentaInput(input_enum,interpolation_enum,this->lid,this->GetNumberOfNodes(interpolation_enum),values);
   }
 
}
/*}}}*/
void       Penta::AddControlInput(int input_enum,Inputs2* inputs2,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id){/*{{{*/
 
   /*Intermediaries*/
   int vertexlids[NUMVERTICES];
 
   _assert_(iomodel->elements);
   for(int i=0;i<NUMVERTICES;i++){
      int vertexid =reCast<int>(iomodel->elements[NUMVERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
      vertexlids[i]=iomodel->my_vertices_lids[vertexid-1];
   }
 
   /*Call inputs method*/
   switch(interpolation_enum){
      case P1Enum:
         inputs2->SetPentaControlInput(input_enum,PentaInput2Enum,interpolation_enum,id,NUMVERTICES,vertexlids,values,values_min,values_max);
         break;
      default:
         _error_("Cannot add \""<<EnumToStringx(input_enum)<<"\" interpolation "<<EnumToStringx(interpolation_enum)<<" not supported");
   }
 
}
/*}}}*/
void       Penta::DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs2* inputs2,IoModel* iomodel,int input_enum){/*{{{*/
 
   /*Intermediaries*/
   int        vertexsids[NUMVERTICES];
   int        vertexlids[NUMVERTICES];
   IssmDouble nodeinputs[NUMVERTICES];
 
   /*Some sanity checks*/
   if(num_inputs<1)                 _error_("Cannot create a DatasetInput of size <1");
   if(M!=iomodel->numberofvertices) _error_("Input size not supported yet");
   if(N!=num_inputs)                _error_("Sizes are not consistent");
 
   /*Get indices*/
   _assert_(iomodel->elements);
   for(int i=0;i<NUMVERTICES;i++){
      vertexsids[i] = reCast<int>(iomodel->elements[NUMVERTICES*this->Sid()+i])-1;
      vertexlids[i] = iomodel->my_vertices_lids[vertexsids[i]];
   }
 
   /*Create inputs and add to DataSetInput*/
   for(int i=0;i<num_inputs;i++){
      for(int j=0;j<NUMVERTICES;j++) nodeinputs[j]=array[vertexsids[j]*N+i];
      inputs2->SetPentaDatasetInput(input_enum,individual_enums[i],P1Enum,NUMVERTICES,vertexlids,nodeinputs);
   }
}
/*}}}*/
void       Penta::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){/*{{{*/
 
   PentaRef::BasalNodeIndices(pnumindices,pindices,finiteelement);
 
}
/*}}}*/
void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
   _error_("Not supported yet!");
}
/*}}}*/
void       Penta::CalvingRateVonmises(){/*{{{*/
 
   if(!this->IsOnBase()) return;
 
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
   IssmDouble  calvingratex[NUMVERTICES];
   IssmDouble  calvingratey[NUMVERTICES];
   IssmDouble  calvingrate[NUMVERTICES];
   IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
   IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
   IssmDouble  epse_2,groundedice,bed,sealevel;
   IssmDouble  sigma_vm[NUMVERTICES];
 
   /* Get node coordinates and dof list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*Depth average B for stress calculation*/
   this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
   this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 
   /*Retrieve all inputs and parameters we will need*/
   Input2* vx_input = this->GetInput2(VxAverageEnum); _assert_(vx_input);
   Input2* vy_input = this->GetInput2(VyAverageEnum); _assert_(vy_input);
   Input2* gr_input = this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
   Input2* bs_input = this->GetInput2(BaseEnum);                    _assert_(bs_input);
   Input2* B_input  = this->GetInput2(MaterialsRheologyBbarEnum);   _assert_(B_input);
   Input2* n_input  = this->GetInput2(MaterialsRheologyNEnum);   _assert_(n_input);
   Input2* smax_fl_input = this->GetInput2(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
   Input2* smax_gr_input = this->GetInput2(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
   Input2* sl_input  = this->GetInput2(SealevelEnum); _assert_(sl_input);
 
   /* Start looping on the number of vertices: */
   GaussPenta* gauss=new GaussPenta();
   for(int iv=0;iv<3;iv++){
      gauss->GaussVertex(iv);
 
      /*Get velocity components and thickness*/
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      gr_input->GetInputValue(&groundedice,gauss);
      bs_input->GetInputValue(&bed,gauss);
      smax_fl_input->GetInputValue(&sigma_max_floating,gauss);
      smax_gr_input->GetInputValue(&sigma_max_grounded,gauss);
      vel=sqrt(vx*vx+vy*vy)+1.e-14;
      sl_input->GetInputValue(&sealevel,gauss);
 
      /*Compute strain rate and viscosity: */
      this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 
      /*Get Eigen values*/
      Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
      _assert_(!xIsNan<IssmDouble>(lambda1));
      _assert_(!xIsNan<IssmDouble>(lambda2));
 
      /*Process Eigen values (only account for extension)*/
      lambda1 = max(lambda1,0.);
      lambda2 = max(lambda2,0.);
 
      /*Calculate sigma_vm*/
      epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
      sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n));
 
      /*Tensile stress threshold*/
      if(groundedice<0)
       sigma_max = sigma_max_floating;
      else
       sigma_max = sigma_max_grounded;
 
      /*Assign values*/
      if(bed>sealevel){
         calvingratex[iv]=0.;
         calvingratey[iv]=0.;
      }
      else{
         calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
         calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
      }
      calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
   }
 
   /*Add input*/
   this->AddBasalInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
   this->AddBasalInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
   this->AddBasalInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
   this->AddBasalInput2(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
 
   this->InputExtrude(CalvingratexEnum,-1);
   this->InputExtrude(CalvingrateyEnum,-1);
   this->InputExtrude(CalvingCalvingrateEnum,-1);
   this->InputExtrude(SigmaVMEnum,-1);
 
   /*Clean up and return*/
   delete gauss;
}
/*}}}*/
void       Penta::CalvingRateLevermann(){/*{{{*/
 
   IssmDouble  xyz_list[NUMVERTICES][3];
   GaussPenta* gauss=NULL;
   IssmDouble  vx,vy,vel;
   IssmDouble  strainparallel;
   IssmDouble  propcoeff;
   IssmDouble  strainperpendicular;
   IssmDouble  calvingratex[NUMVERTICES];
   IssmDouble  calvingratey[NUMVERTICES];
   IssmDouble  calvingrate[NUMVERTICES];
 
   /* Get node coordinates and dof list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*Retrieve all inputs and parameters we will need*/
   Input2* vx_input=this->GetInput2(VxEnum);                                                    _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum);                                                    _assert_(vy_input);
   Input2* strainparallel_input=this->GetInput2(StrainRateparallelEnum);                        _assert_(strainparallel_input);
   Input2* strainperpendicular_input=this->GetInput2(StrainRateperpendicularEnum);              _assert_(strainperpendicular_input);
   Input2* levermanncoeff_input=this->GetInput2(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
 
   /* Start looping on the number of vertices: */
   gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      /* Get the value we need*/
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      vel=vx*vx+vy*vy;
      strainparallel_input->GetInputValue(&strainparallel,gauss);
      strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
      levermanncoeff_input->GetInputValue(&propcoeff,gauss);
 
      /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
      calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
      if(calvingrate[iv]<0){
         calvingrate[iv]=0;
      }
      calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
      calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
   }
 
   /*Add input*/
   this->AddBasalInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
   this->AddBasalInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
   this->AddBasalInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
 
   /*Clean up and return*/
   delete gauss;
 
}
/*}}}*/
void       Penta::CalvingFluxLevelset(){/*{{{*/
 
   /*Make sure there is an ice front here*/
   if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
      IssmDouble flux_per_area=0;
      this->AddInput2(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum);
   }
   else{
      int               domaintype,index1,index2;
      const IssmPDouble epsilon = 1.e-15;
      IssmDouble        s1,s2;
      IssmDouble        gl[NUMVERTICES];
      IssmDouble        xyz_front[2][3];
 
      IssmDouble *xyz_list = NULL;
      this->GetVerticesCoordinates(&xyz_list);
 
      /*Recover parameters and values*/
      Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
      /*Be sure that values are not zero*/
      if(gl[0]==0.) gl[0]=gl[0]+epsilon;
      if(gl[1]==0.) gl[1]=gl[1]+epsilon;
      if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
      int pt1 = 0;
      int pt2 = 1;
      if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 
         /*Portion of the segments*/
         s1=gl[2]/(gl[2]-gl[1]);
         s2=gl[2]/(gl[2]-gl[0]);
         if(gl[2]<0.){
            pt1 = 1; pt2 = 0;
         }
         xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
         xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
         xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
         xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
         xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
         xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
      }
      else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
 
         /*Portion of the segments*/
         s1=gl[0]/(gl[0]-gl[1]);
         s2=gl[0]/(gl[0]-gl[2]);
         if(gl[0]<0.){
            pt1 = 1; pt2 = 0;
         }
 
         xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
         xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
         xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
         xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
         xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
         xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
      }
      else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
 
         /*Portion of the segments*/
         s1=gl[1]/(gl[1]-gl[0]);
         s2=gl[1]/(gl[1]-gl[2]);
         if(gl[1]<0.){
            pt1 = 1; pt2 = 0;
         }
 
         xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
         xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
         xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
         xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
         xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
         xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
      }
      else{
         _error_("case not possible");
      }
 
      /*Some checks in debugging mode*/
      _assert_(s1>=0 && s1<=1.);
      _assert_(s2>=0 && s2<=1.);
 
      /*Get normal vector*/
      IssmDouble normal[3];
      this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
      normal[0] = -normal[0];
      normal[1] = -normal[1];
 
      /*Get inputs*/
      IssmDouble flux = 0.;
      IssmDouble area = 0.;
      IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
      IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
      Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
      Input2* calvingratex_input=NULL;
      Input2* calvingratey_input=NULL;
      calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
      calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
 
      /*Start looping on Gaussian points*/
      Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
         thickness_input->GetInputValue(&thickness,gauss);
         calvingratex_input->GetInputValue(&calvingratex,gauss);
         calvingratey_input->GetInputValue(&calvingratey,gauss);
         this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
 
         flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
         area += Jdet*gauss->weight*thickness;
 
         flux_per_area=flux/area;
      }
 
      this->AddInput2(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum);
 
      /*Clean up and return*/
      delete gauss;
   }
}
/*}}}*/
void       Penta::CalvingMeltingFluxLevelset(){/*{{{*/
 
   /*Make sure there is an ice front here*/
   if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
      IssmDouble flux_per_area=0;
      this->AddInput2(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum);
   }
   else{
      int               domaintype,index1,index2;
      const IssmPDouble epsilon = 1.e-15;
      IssmDouble        s1,s2;
      IssmDouble        gl[NUMVERTICES];
      IssmDouble        xyz_front[2][3];
 
      IssmDouble *xyz_list = NULL;
      this->GetVerticesCoordinates(&xyz_list);
 
      /*Recover parameters and values*/
      Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
      /*Be sure that values are not zero*/
      if(gl[0]==0.) gl[0]=gl[0]+epsilon;
      if(gl[1]==0.) gl[1]=gl[1]+epsilon;
      if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
      int pt1 = 0;
      int pt2 = 1;
      if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 
         /*Portion of the segments*/
         s1=gl[2]/(gl[2]-gl[1]);
         s2=gl[2]/(gl[2]-gl[0]);
         if(gl[2]<0.){
            pt1 = 1; pt2 = 0;
         }
         xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
         xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
         xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
         xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
         xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
         xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
      }
      else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
 
         /*Portion of the segments*/
         s1=gl[0]/(gl[0]-gl[1]);
         s2=gl[0]/(gl[0]-gl[2]);
         if(gl[0]<0.){
            pt1 = 1; pt2 = 0;
         }
 
         xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
         xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
         xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
         xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
         xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
         xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
      }
      else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
 
         /*Portion of the segments*/
         s1=gl[1]/(gl[1]-gl[0]);
         s2=gl[1]/(gl[1]-gl[2]);
         if(gl[1]<0.){
            pt1 = 1; pt2 = 0;
         }
 
         xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
         xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
         xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
         xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
         xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
         xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
      }
      else{
         _error_("case not possible");
      }
 
      /*Some checks in debugging mode*/
      _assert_(s1>=0 && s1<=1.);
      _assert_(s2>=0 && s2<=1.);
 
      /*Get normal vector*/
      IssmDouble normal[3];
      this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
      normal[0] = -normal[0];
      normal[1] = -normal[1];
 
      /*Get inputs*/
      IssmDouble flux = 0.;
      IssmDouble area = 0.;
      IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
      IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
      Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
      Input2* calvingratex_input=NULL;
      Input2* calvingratey_input=NULL;
      Input2* vx_input=NULL;
      Input2* vy_input=NULL;
      Input2* meltingrate_input=NULL;
      calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
      calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
      vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
      vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
      meltingrate_input=this->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
 
      /*Start looping on Gaussian points*/
      Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
      for(int ig=gauss->begin();ig<gauss->end();ig++){
 
         gauss->GaussPoint(ig);
         thickness_input->GetInputValue(&thickness,gauss);
         calvingratex_input->GetInputValue(&calvingratex,gauss);
         calvingratey_input->GetInputValue(&calvingratey,gauss);
         vx_input->GetInputValue(&vx,gauss);
         vy_input->GetInputValue(&vy,gauss);
         vel=vx*vx+vy*vy;
         meltingrate_input->GetInputValue(&meltingrate,gauss);
         meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
         meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
         this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
 
         flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
         area += Jdet*gauss->weight*thickness;
 
         flux_per_area=flux/area;
      }
 
      this->AddInput2(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum);
 
      /*Clean up and return*/
      delete gauss;
   }
}
/*}}}*/
void       Penta::ComputeBasalStress(void){/*{{{*/
 
   _error_("not implemented (needs to be redone)");
   int         i,j;
   int         dofv[3]={0,1,2};
   int         dofp[1]={3};
   int         analysis_type,approximation;
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  xyz_list_tria[3][3];
   IssmDouble  rho_ice,gravity,FSreconditioning;
   IssmDouble  pressure,viscosity,Jdet2d;
   IssmDouble  bed_normal[3];
   IssmDouble  basalforce[3] = {0.};
   IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
   IssmDouble  stresstensor[6]={0.0};
   IssmDouble  sigma_xx,sigma_yy,sigma_zz;
   IssmDouble  sigma_xy,sigma_xz,sigma_yz;
   IssmDouble  surface=0,value=0;
   GaussPenta* gauss;
 
   /*retrive parameters: */
   parameters->FindParam(&analysis_type,AnalysisTypeEnum);
   this->GetInput2Value(&approximation,ApproximationEnum);
 
   /*Check analysis_types*/
   if (analysis_type!=StressbalanceAnalysisEnum) _error_("Not supported yet!");
   if (approximation!=FSApproximationEnum) _error_("Not supported yet!");
 
   /*retrieve some parameters: */
   this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 
   if(!IsOnBase()){
      //put zero
      //sigma_b->SetValue(id-1,0.0,INS_VAL);
      return;
   }
 
   /*recovre material parameters: */
   rho_ice=FindParam(MaterialsRhoIceEnum);
   gravity=FindParam(ConstantsGEnum);
 
   /* Get node coordinates and dof list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   for(i=0;i<3;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
 
   /*Retrieve all inputs we will be needing: */
   Input2* pressure_input=this->GetInput2(PressureEnum); _assert_(pressure_input);
   Input2* vx_input=this->GetInput2(VxEnum);             _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum);             _assert_(vy_input);
   Input2* vz_input=this->GetInput2(VzEnum);             _assert_(vz_input);
 
   /* Start  looping on the number of gaussian points: */
   gauss=new GaussPenta(0,1,2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      /*Compute strain rate viscosity and pressure: */
      this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
      this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
      pressure_input->GetInputValue(&pressure,gauss);
 
      /*Compute Stress*/
      sigma_xx=2*viscosity*epsilon[0]-pressure*FSreconditioning; // sigma = nu eps - pressure
      sigma_yy=2*viscosity*epsilon[1]-pressure*FSreconditioning;
      sigma_zz=2*viscosity*epsilon[2]-pressure*FSreconditioning;
      sigma_xy=2*viscosity*epsilon[3];
      sigma_xz=2*viscosity*epsilon[4];
      sigma_yz=2*viscosity*epsilon[5];
 
      /*Get normal vector to the bed */
      NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
 
      /*basalforce*/
      basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2];
      basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2];
      basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2];
 
      GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
      value+=sigma_zz*Jdet2d*gauss->weight;
      surface+=Jdet2d*gauss->weight;
   }
   value=value/surface;
 
   /*Add value to output*/
   //sigma_b->SetValue(id-1,value,INS_VAL);
}
/*}}}*/
void       Penta::ComputeDeviatoricStressTensor(){/*{{{*/
 
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  viscosity;
   IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
   IssmDouble  tau_xx[NUMVERTICES];
   IssmDouble  tau_yy[NUMVERTICES];
   IssmDouble  tau_zz[NUMVERTICES];
   IssmDouble  tau_xy[NUMVERTICES];
   IssmDouble  tau_xz[NUMVERTICES];
   IssmDouble  tau_yz[NUMVERTICES];
   IssmDouble  tau_eff[NUMVERTICES];
   GaussPenta* gauss=NULL;
 
   /* Get node coordinates and dof list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*Retrieve all inputs we will be needing: */
   Input2* vx_input=this->GetInput2(VxEnum);             _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum);             _assert_(vy_input);
   Input2* vz_input=this->GetInput2(VzEnum);             _assert_(vz_input);
 
   /* Start looping on the number of vertices: */
   gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      /*Compute strain rate viscosity and pressure: */
      this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
      this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 
      /*Compute Stress*/
      tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
      tau_yy[iv]=2*viscosity*epsilon[1];
      tau_zz[iv]=2*viscosity*epsilon[2];
      tau_xy[iv]=2*viscosity*epsilon[3];
      tau_xz[iv]=2*viscosity*epsilon[4];
      tau_yz[iv]=2*viscosity*epsilon[5];
 
      tau_eff[iv] = tau_xx[iv]*tau_xx[iv] + tau_yy[iv]*tau_yy[iv] + tau_zz[iv]*tau_zz[iv] +
        2*tau_xy[iv]*tau_xy[iv] + 2*tau_xz[iv]*tau_xz[iv] + 2*tau_yz[iv]*tau_yz[iv];
 
      tau_eff[iv] = sqrt(tau_eff[iv]/2.);
   }
 
   /*Add Stress tensor components into inputs*/
   this->AddInput2(DeviatoricStressxxEnum,&tau_xx[0],P1DGEnum);
   this->AddInput2(DeviatoricStressxyEnum,&tau_xy[0],P1DGEnum);
   this->AddInput2(DeviatoricStressxzEnum,&tau_xz[0],P1DGEnum);
   this->AddInput2(DeviatoricStressyyEnum,&tau_yy[0],P1DGEnum);
   this->AddInput2(DeviatoricStressyzEnum,&tau_yz[0],P1DGEnum);
   this->AddInput2(DeviatoricStresszzEnum,&tau_zz[0],P1DGEnum);
   this->AddInput2(DeviatoricStresseffectiveEnum,&tau_eff[0],P1DGEnum);
 
   /*Clean up and return*/
   delete gauss;
}
/*}}}*/
void       Penta::ComputeStressTensor(){/*{{{*/
 
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  pressure,viscosity;
   IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
   IssmDouble  sigma_xx[NUMVERTICES];
   IssmDouble  sigma_yy[NUMVERTICES];
   IssmDouble  sigma_zz[NUMVERTICES];
   IssmDouble  sigma_xy[NUMVERTICES];
   IssmDouble  sigma_xz[NUMVERTICES];
   IssmDouble  sigma_yz[NUMVERTICES];
   GaussPenta* gauss=NULL;
 
   /* Get node coordinates and dof list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*Retrieve all inputs we will be needing: */
   Input2* pressure_input=this->GetInput2(PressureEnum); _assert_(pressure_input);
   Input2* vx_input=this->GetInput2(VxEnum);             _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum);             _assert_(vy_input);
   Input2* vz_input=this->GetInput2(VzEnum);             _assert_(vz_input);
 
   /* Start looping on the number of vertices: */
   gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      /*Compute strain rate viscosity and pressure: */
      this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
      this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
      pressure_input->GetInputValue(&pressure,gauss);
 
      /*Compute Stress*/
      sigma_xx[iv]=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
      sigma_yy[iv]=2*viscosity*epsilon[1]-pressure;
      sigma_zz[iv]=2*viscosity*epsilon[2]-pressure;
      sigma_xy[iv]=2*viscosity*epsilon[3];
      sigma_xz[iv]=2*viscosity*epsilon[4];
      sigma_yz[iv]=2*viscosity*epsilon[5];
   }
 
   /*Add Stress tensor components into inputs*/
   this->AddInput2(StressTensorxxEnum,&sigma_xx[0],P1DGEnum);
   this->AddInput2(StressTensorxyEnum,&sigma_xy[0],P1DGEnum);
   this->AddInput2(StressTensorxzEnum,&sigma_xz[0],P1DGEnum);
   this->AddInput2(StressTensoryyEnum,&sigma_yy[0],P1DGEnum);
   this->AddInput2(StressTensoryzEnum,&sigma_yz[0],P1DGEnum);
   this->AddInput2(StressTensorzzEnum,&sigma_zz[0],P1DGEnum);
 
   /*Clean up and return*/
   delete gauss;
}
/*}}}*/
void       Penta::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin,Inputs2* inputs2in){/*{{{*/
 
   int analysis_counter;
 
   /*go into parameters and get the analysis_counter: */
   parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
 
   /*Get Element type*/
   this->element_type=this->element_type_list[analysis_counter];
 
   /*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 offsets hidden in hooks: */
   if (this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin);
   this->hvertices->configure(verticesin);
   this->hmaterial->configure(materialsin);
   this->hneighbors->configure(elementsin);
 
   /*Now, go pick up the objects inside the hooks: */
   if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
   else this->nodes=NULL;
   this->vertices          = (Vertex**)this->hvertices->deliverp();
   this->material          = (Material*)this->hmaterial->delivers();
   this->verticalneighbors = (Penta**)this->hneighbors->deliverp();
 
   /*point parameters to real dataset: */
   this->parameters=parametersin;
   this->inputs2=inputs2in;
}
/*}}}*/
void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){/*{{{*/
 
   if(enum_type==MaterialsRheologyBbarEnum) enum_type = MaterialsRheologyBEnum;
   if(enum_type==DamageDbarEnum)            enum_type = DamageDEnum;
 
   _error_("not implemented");
   int    vertexpidlist[NUMVERTICES];
   IssmDouble grad_list[NUMVERTICES];
   Input2* grad_input=NULL;
   Input2* input=NULL;
 
   if(enum_type==MaterialsRheologyBbarEnum){
      input=this->GetInput2(MaterialsRheologyBEnum);
   }
   else if(enum_type==DamageDbarEnum){
      input=this->GetInput2(DamageDEnum);
   }
   else{
      input=this->GetInput2(enum_type);
   }
   if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
   if(input->ObjectEnum()!=ControlInput2Enum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
 
   GradientIndexing(&vertexpidlist[0],control_index);
 
   //for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
   //grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
   //((ControlInput*)input)->SetGradient(grad_input);
   _error_("not implemented");
 
}
/*}}}*/
void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
 
   int        idlist[NUMVERTICES];
   int        vertexlids[NUMVERTICES];
   IssmDouble grad_list[NUMVERTICES];
 
   if(enum_type==MaterialsRheologyBbarEnum) enum_type = MaterialsRheologyBEnum;
   if(enum_type==DamageDbarEnum)            enum_type = DamageDEnum;
 
   GradientIndexing(&idlist[0],control_index);
   for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[idlist[i]];
   for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;
 
   this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0]);
}/*}}}*/
void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
 
   int         sidlist[NUMVERTICES];
   int         lidlist[NUMVERTICES];
   int         connectivity[NUMVERTICES];
   IssmPDouble values[NUMVERTICES];
   IssmPDouble gradients[NUMVERTICES];
   IssmDouble  value,gradient;
 
   if(control_enum==MaterialsRheologyBbarEnum) control_enum = MaterialsRheologyBEnum;
   if(control_enum==DamageDbarEnum)            control_enum = DamageDEnum;
 
   this->GetVerticesConnectivityList(&connectivity[0]);
   this->GetVerticesSidList(&sidlist[0]);
   this->GetVerticesLidList(&lidlist[0]);
 
   ElementInput2* control_value    = this->inputs2->GetControlInput2Data(control_enum,"value");    _assert_(control_value);
   ElementInput2* control_gradient = this->inputs2->GetControlInput2Data(control_enum,"gradient"); _assert_(control_gradient);
   control_value->Serve(NUMVERTICES,&lidlist[0]);
   control_gradient->Serve(NUMVERTICES,&lidlist[0]);
 
   GaussPenta* gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      control_value->GetInputValue(&value,gauss);
      control_gradient->GetInputValue(&gradient,gauss);
 
      values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
      gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
   }
   delete gauss;
 
   vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
   vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
}/*}}}*/
void       Penta::CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){/*{{{*/
 
   /*Get current field and vertex coordinates*/
   IssmDouble ls[NUMVERTICES],distance;
   Element::GetInputListOnVertices(&ls[0],distanceenum);
 
   /*Get distance from list of segments and reset ls*/
   for(int j=0;j<NUMVERTICES;j++){
      distance=distances[this->vertices[j]->Lid()];
      if(xIsNan<IssmDouble>(distance)) _error_("NaN found in vector");
      if(xIsInf<IssmDouble>(distance)) _error_("Inf found in vector");
 
      /*FIXME: do we really need this?*/
      if(distanceenum==MaskIceLevelsetEnum) if(distance>10000) distance=10000;
      if(ls[j]>0){
         ls[j] = distance;
      }
      else{
         ls[j] = - distance;
      }
   }
 
   /*Update Levelset*/
   this->AddInput2(distanceenum,&ls[0],P1Enum);
}
/*}}}*/
void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
   _assert_(end_time>start_time);
 
   /*Get transient input time steps*/
   TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
   PentaInput2* averaged_input = transient_input->GetPentaInput(start_time,end_time,averaging_method);
   Input2* averaged_copy = averaged_input->copy();
 
   averaged_copy->ChangeEnum(averagedinput_enum);
   this->inputs2->AddInput(averaged_copy);
}
/*}}}*/
void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
 
   switch(response_enum){
      case MaterialsRheologyBbarEnum:
         *presponse=this->material->GetBbar(NULL);
         break;
      case DamageDbarEnum:
         *presponse=this->material->GetDbar(NULL);
         break;
      case VelEnum:
         {
 
            /*Get input:*/
            IssmDouble vel;
            Input2* vel_input=this->GetInput2(VelEnum); _assert_(vel_input);
            vel_input->GetInputAverage(&vel);
 
            /*Assign output pointers:*/
            *presponse=vel;
         }
         break;
      default:
         _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
   }
 
}
/*}}}*/
void       Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
 
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble xmin,ymin,zmin;
   IssmDouble xmax,ymax,zmax;
 
   /*Get xyz list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
   ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
   zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
 
   for(int i=1;i<NUMVERTICES;i++){
      if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
      if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
      if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
      if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
      if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
      if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
   }
 
   *hx=xmax-xmin;
   *hy=ymax-ymin;
   *hz=zmax-zmin;
}
/*}}}*/
int        Penta::FiniteElement(void){/*{{{*/
   return this->element_type;
}
/*}}}*/
IssmDouble Penta::FloatingArea(bool scaled){/*{{{*/
 
   /*Intermediaries*/
   int         domaintype;
   IssmDouble  phi,base_area,scalefactor,floatingarea;
   IssmDouble  xyz_list[NUMVERTICES][3];
 
   if(!IsIceInElement() || !IsOnBase())return 0.;
 
   /*Get problem dimension*/
   this->FindParam(&domaintype,DomainTypeEnum);
   if(domaintype!=Domain3DEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   phi=this->GetGroundedPortion(&xyz_list[0][0]);
   base_area= 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
 
   floatingarea=(1-phi)*base_area;
 
   if(scaled==true){
      Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
      scalefactor_input->GetInputAverage(&scalefactor);
      floatingarea=floatingarea*scalefactor;
   }
 
   /*Clean up and return*/
   return floatingarea;
}
/*}}}*/
void       Penta::FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure){/*{{{*/
 
   if(!IsOnBase()) return;
 
   int approximation;
   this->GetInput2Value(&approximation,ApproximationEnum);
   if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum || approximation==HOFSApproximationEnum){
      _error_("Cannot compute contact condition for non FS elements");
   }
 
   /*Intermediaries*/
   IssmDouble* xyz_list = NULL;
   IssmDouble  bed_normal[3],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
   IssmDouble  water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
   IssmDouble  sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmazz[NUMVERTICES],sigmaxy[NUMVERTICES];
   IssmDouble  sigmayz[NUMVERTICES],sigmaxz[NUMVERTICES],sigma_nn[NUMVERTICES];
   IssmDouble  viscosity,epsilon[NUMVERTICES];
   Element::GetInputListOnVertices(&base[0],BaseEnum);
   Element::GetInputListOnVertices(&bed[0],BedEnum);
   Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
   Element::GetInputListOnVertices(&pressure[0],PressureEnum);
   Element::GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
   IssmDouble rho_ice   = FindParam(MaterialsRhoIceEnum);
   IssmDouble rho_water = FindParam(MaterialsRhoSeawaterEnum);
   IssmDouble gravity   = FindParam(ConstantsGEnum);
 
   /* Get node coordinates and dof list: */
   GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will be needing: */
   Input2* vx_input = this->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input = this->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vz_input = this->GetInput2(VzEnum); _assert_(vz_input);
 
   /*1. Recover stresses at the base*/
   GaussPenta* gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      /*Compute strain rate viscosity and pressure: */
      this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
      this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
      /*FIXME: this is for Hongju only*/
      //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
      //if (pressure[iv]/pressureice[iv]>1)  pressure[iv]=pressureice[iv];
 
      /*Compute Stress*/
      sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv]; // sigma = nu eps - pressure
      sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
      sigmazz[iv]=2*viscosity*epsilon[2]-pressure[iv];
      sigmaxy[iv]=2*viscosity*epsilon[3];
      sigmaxz[iv]=2*viscosity*epsilon[4];
      sigmayz[iv]=2*viscosity*epsilon[5];
   }
 
   /*2. compute contact condition*/
   for(int i=0;i<NUMVERTICES;i++){
 
      /*If was grounded*/
      if (phi[i]>=0.){
         NormalBase(&bed_normal[0],xyz_list);
         sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1] + sigmazz[i]*bed_normal[2]*bed_normal[2]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]+2*sigmaxz[i]*bed_normal[0]*bed_normal[2]+2*sigmayz[i]*bed_normal[1]*bed_normal[2]);
         water_pressure[i]=-gravity*rho_water*base[i];
         vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL);
         vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL);
      }
 
      /*If was floating*/
      else{
         /*Tricky part:
          * 1. if base is now touching, we put 1 for sigma_nn and leave water pressure at 0 so that the rest of the module will reground this vertex
          * 2. if base is still above bed, water pressure is set as 1, sigma_nn is left as 0, so the GL module will keep it afloat*/
         if(base[i]<bed[i]) vertex_sigmann->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
         else vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
      }
   }
 
   /*clean up*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
}
/*}}}*/
void       Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
   /*Computeportion of the element that is grounded*/
 
   int         i,j,k;
   IssmDouble  area_init,area_portion;
   IssmDouble  xyz_bis[3][3];
 
   area_init=fabs(xyz_list[1*3+0]*xyz_list[2*3+1] - xyz_list[1*3+1]*xyz_list[2*3+0] + xyz_list[0*3+0]*xyz_list[1*3+1] - xyz_list[0*3+1]*xyz_list[1*3+0] + xyz_list[2*3+0]*xyz_list[0*3+1] - xyz_list[2*3+1]*xyz_list[0*3+0])/2.;
 
   /*Initialize xyz_list with original xyz_list of triangle coordinates*/
   for(j=0;j<3;j++){
      for(k=0;k<3;k++){
         xyz_bis[j][k]=xyz_list[j*3+k];
      }
   }
   for(i=0;i<numpoints;i++){
      for(j=0;j<3;j++){
         for(k=0;k<3;k++){
            /*Change appropriate line*/
            xyz_bis[j][k]=xyz_zero[i*3+k];
         }
 
         /*Compute area fraction*/
         area_portion=fabs(xyz_bis[1][0]*xyz_bis[2][1] - xyz_bis[1][1]*xyz_bis[2][0] + xyz_bis[0][0]*xyz_bis[1][1] - xyz_bis[0][1]*xyz_bis[1][0] + xyz_bis[2][0]*xyz_bis[0][1] - xyz_bis[2][1]*xyz_bis[0][0])/2.;
         *(area_coordinates+3*i+j)=area_portion/area_init;
 
         /*Reinitialize xyz_list*/
         for(k=0;k<3;k++){
            /*Reinitialize xyz_list with original coordinates*/
            xyz_bis[j][k]=xyz_list[j*3+k];
         }
      }
   }
}
/*}}}*/
Element*   Penta::GetBasalElement(void){/*{{{*/
 
   /*Output*/
   Element* element=this->GetBasalPenta();
   return element;
}
/*}}}*/
Penta*     Penta::GetBasalPenta(void){/*{{{*/
 
   /*Output*/
   Penta* penta=NULL;
 
   /*Go through all pentas till the bed is reached*/
   penta=this;
   for(;;){
      /*Stop if we have reached the surface, else, take lower penta*/
      if (penta->IsOnBase()) break;
 
      /* get lower Penta*/
      penta=penta->GetLowerPenta();
      _assert_(penta->Id()!=this->id);
   }
 
   /*return output*/
   return penta;
}
/*}}}*/
int        Penta::GetElementType(){/*{{{*/
 
   /*return PentaRef field*/
   return this->element_type;
}
/*}}}*/
void       Penta::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
   /*Computeportion of the element that is grounded*/
 
   bool               floating=true;
   int                point;
   const IssmPDouble  epsilon= 1.e-15;
   IssmDouble         gl[NUMVERTICES];
   IssmDouble         f1,f2;
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&gl[0],MaskOceanLevelsetEnum);
 
   /*Be sure that values are not zero*/
   if(gl[0]==0.) gl[0]=gl[0]+epsilon;
   if(gl[1]==0.) gl[1]=gl[1]+epsilon;
   if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
   /*Check that not all nodes are grounded or floating*/
   if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
      point=0;
      f1=1.;
      f2=1.;
   }
   else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
      point=0;
      f1=0.;
      f2=0.;
   }
   else{
      if(gl[0]*gl[1]*gl[2]<0) floating=false;
 
      if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
         point=2;
         f1=gl[2]/(gl[2]-gl[0]);
         f2=gl[2]/(gl[2]-gl[1]);
      }
      else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
         point=0;
         f1=gl[0]/(gl[0]-gl[1]);
         f2=gl[0]/(gl[0]-gl[2]);
      }
      else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
         point=1;
         f1=gl[1]/(gl[1]-gl[2]);
         f2=gl[1]/(gl[1]-gl[0]);
      }
      else _error_("case not possible");
   }
   *point1=point;
   *fraction1=f1;
   *fraction2=f2;
   *mainlyfloating=floating;
}
/*}}}*/
IssmDouble Penta::GetGroundedPortion(IssmDouble* xyz_list){/*{{{*/
   /*Computeportion of the element that is grounded*/
 
   bool               mainlyfloating = true;
   const IssmPDouble  epsilon= 1.e-15;
   IssmDouble         phi,s1,s2,area_init,area_grounded;
   IssmDouble         gl[NUMVERTICES];
   IssmDouble         xyz_bis[NUMVERTICES2D][3];
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&gl[0],MaskOceanLevelsetEnum);
 
   /*Be sure that values are not zero*/
   if(gl[0]==0.) gl[0]=gl[0]+epsilon;
   if(gl[1]==0.) gl[1]=gl[1]+epsilon;
   if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
   /*Check that not all nodes are grounded or floating*/
   if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
      phi=1;
   }
   else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
      phi=0;
   }
   else{
      /*Figure out if two nodes are floating or grounded*/
      if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
 
      if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
         /*Coordinates of point 2: same as initial point 2*/
         xyz_bis[2][0]=xyz_list[3*2+0];
         xyz_bis[2][1]=xyz_list[3*2+1];
         xyz_bis[2][2]=xyz_list[3*2+2];
 
         /*Portion of the segments*/
         s1=gl[2]/(gl[2]-gl[1]);
         s2=gl[2]/(gl[2]-gl[0]);
 
         /*New point 1*/
         xyz_bis[1][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
         xyz_bis[1][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
         xyz_bis[1][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
 
         /*New point 0*/
         xyz_bis[0][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
         xyz_bis[0][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
         xyz_bis[0][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
      }
      else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
         /*Coordinates of point 0: same as initial point 2*/
         xyz_bis[0][0]=*(xyz_list+3*0+0);
         xyz_bis[0][1]=*(xyz_list+3*0+1);
         xyz_bis[0][2]=*(xyz_list+3*0+2);
 
         /*Portion of the segments*/
         s1=gl[0]/(gl[0]-gl[1]);
         s2=gl[0]/(gl[0]-gl[2]);
 
         /*New point 1*/
         xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
         xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
         xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
 
         /*New point 2*/
         xyz_bis[2][0]=*(xyz_list+3*0+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
         xyz_bis[2][1]=*(xyz_list+3*0+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
         xyz_bis[2][2]=*(xyz_list+3*0+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
      }
      else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
         /*Coordinates of point 1: same as initial point 2*/
         xyz_bis[1][0]=*(xyz_list+3*1+0);
         xyz_bis[1][1]=*(xyz_list+3*1+1);
         xyz_bis[1][2]=*(xyz_list+3*1+2);
 
         /*Portion of the segments*/
         s1=gl[1]/(gl[1]-gl[0]);
         s2=gl[1]/(gl[1]-gl[2]);
 
         /*New point 0*/
         xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
         xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
         xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
 
         /*New point 2*/
         xyz_bis[2][0]=*(xyz_list+3*1+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
         xyz_bis[2][1]=*(xyz_list+3*1+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
         xyz_bis[2][2]=*(xyz_list+3*1+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
      }
      else _error_("case not possible");
 
      /*Compute fraction of grounded element*/
      GetTriaJacobianDeterminant(&area_init, xyz_list,NULL);
      GetTriaJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
      if(mainlyfloating==true) area_grounded=area_init-area_grounded;
      phi=area_grounded/area_init;
   }
 
   if(phi>1. || phi<0.) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
 
   return phi;
}
/*}}}*/
IssmDouble Penta::GetIcefrontArea(){/*{{{*/
 
   IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
   IssmDouble  Haverage,frontarea;
   IssmDouble  x1,y1,x2,y2,distance;
   IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
   int* indices=NULL;
   IssmDouble* H=NULL;;
   int nrfrontbed,numiceverts;
 
   if(!this->IsOnBase()) return 0;
   if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
 
   /*Retrieve all inputs and parameters*/
   Element::GetInputListOnVertices(&bed[0],BedEnum);
   Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
   Element::GetInputListOnVertices(&bases[0],BaseEnum);
   Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
 
   nrfrontbed=0;
   for(int i=0;i<NUMVERTICES2D;i++){
      /*Find if bed<0*/
      if(bed[i]<0.) nrfrontbed++;
   }
 
   if(nrfrontbed==3){
      /*2. Find coordinates of where levelset crosses 0*/
      int         numiceverts;
      IssmDouble  s[2],x[2],y[2];
      this->GetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
      _assert_(numiceverts);
 
      /*3 Write coordinates*/
      IssmDouble  xyz_list[NUMVERTICES][3];
      ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
      int counter = 0;
      if((numiceverts>0) && (numiceverts<NUMVERTICES2D)){
         for(int i=0;i<numiceverts;i++){
            for(int n=numiceverts;n<NUMVERTICES2D;n++){ // iterate over no-ice vertices
               x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
               y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
               counter++;
            }
         }
      }
      else if(numiceverts==NUMVERTICES2D){ //NUMVERTICES ice vertices: calving front lies on element edge
 
         for(int i=0;i<NUMVERTICES2D;i++){
            if(lsf[indices[i]]==0.){
               x[counter]=xyz_list[indices[i]][0];
               y[counter]=xyz_list[indices[i]][1];
               counter++;
            }
            if(counter==2) break;
         }
         if(counter==1){
            /*We actually have only 1 vertex on levelset, write a single point as a segment*/
            x[counter]=x[0];
            y[counter]=y[0];
            counter++;
         }
      }
      else{
         _error_("not sure what's going on here...");
      }
      x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
      distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
 
      int numthk=numiceverts+2;
      H=xNew<IssmDouble>(numthk);
      for(int iv=0;iv<NUMVERTICES2D;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
 
      switch(numiceverts){
         case 1: // average over triangle
            H[0]=Haux[0];
            H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
            H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
            Haverage=(H[1]+H[2])/2;
            break;
         case 2: // average over quadrangle
            H[0]=Haux[0];
            H[1]=Haux[1];
            H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
            H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
            Haverage=(H[2]+H[3])/2;
            break;
         default:
            _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
            break;
      }
      frontarea=distance*Haverage;
   }
   else return 0;
 
   xDelete<int>(indices);
   xDelete<IssmDouble>(H);
 
   _assert_(frontarea>0);
   return frontarea;
}
/*}}}*/
void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
 
   /* Intermediaries */
   const int dim=3;
   int i, dir,nrfrontnodes;
   IssmDouble  levelset[NUMVERTICES];
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&levelset[0],levelsetenum);
 
   int* indicesfront = xNew<int>(NUMVERTICES);
   /* Get basal nodes where there is no ice */
   nrfrontnodes=0;
   for(i=0;i<NUMVERTICES2D;i++){
      if(levelset[i]>=0.){
         indicesfront[nrfrontnodes]=i;
         nrfrontnodes++;
      }
   }
   _assert_(nrfrontnodes==2);
 
   /* arrange order of basal frontnodes such that they are oriented counterclockwise */
   if((NUMVERTICES2D+indicesfront[0]-indicesfront[1])%NUMVERTICES2D!=NUMVERTICES2D-1){
      int index=indicesfront[0];
      indicesfront[0]=indicesfront[1];
      indicesfront[1]=index;
   }
 
   IssmDouble* xyz_front = xNew<IssmDouble>(2*dim*nrfrontnodes);
   /* Return basal and top front nodes */
   for(i=0;i<nrfrontnodes;i++){
      for(dir=0;dir<dim;dir++){
         int ind1=i*dim+dir, ind2=(2*nrfrontnodes-1-i)*dim+dir; // vertex structure front segment: base0, base1, top1, top0
         xyz_front[ind1]=xyz_list[dim*indicesfront[i]+dir];
         xyz_front[ind2]=xyz_list[dim*(indicesfront[i]+NUMVERTICES2D)+dir];
      }
   }
 
   *pxyz_front=xyz_front;
 
   xDelete<int>(indicesfront);
}/*}}}*/
Input2*    Penta::GetInput2(int inputenum){/*{{{*/
 
   /*Get Input from dataset*/
   PentaInput2* input = this->inputs2->GetPentaInput(inputenum);
   if(!input) return input;
 
   /*Intermediaries*/
   int numindices;
   int indices[30]; /*Max numnodes*/
 
   /*Check interpolation*/
   int interpolation = input->GetInterpolation();
   if(interpolation==P1Enum){
      numindices = 6;
      for(int i=0;i<6;i++) indices[i] = vertices[i]->lid;
      input->Serve(numindices,&indices[0]);
   }
   else{
      input->Serve(this->lid,this->GetNumberOfNodes(interpolation));
   }
 
   /*Tell input it is NOT collapsed*/
   //input->SetServeCollapsed(0); FIXME: not needed?
 
   /*Return*/
   return input;
}/*}}}*/
Input2*    Penta::GetInput2(int inputenum,IssmDouble time){/*{{{*/
 
   /*Get Input from dataset*/
   PentaInput2* input = this->inputs2->GetPentaInput(inputenum,time);
   if(!input) return input;
 
   /*Intermediaries*/
   int numindices;
   int indices[30]; /*Max numnodes*/
 
   /*Check interpolation*/
   int interpolation = input->GetInterpolation();
   if(interpolation==P1Enum){
      numindices = 6;
      for(int i=0;i<6;i++) indices[i] = vertices[i]->lid;
      input->Serve(numindices,&indices[0]);
   }
   else{
      input->Serve(this->lid,this->GetNumberOfNodes(interpolation));
   }
 
   /*Tell input it is NOT collapsed*/
   //input->SetServeCollapsed(0); FIXME: not needed?
 
   /*Return*/
   return input;
}/*}}}*/
void       Penta::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
 
   /*Checks in debugging mode*/
   _assert_(pvalue);
 
   /* Start looping on the number of vertices: */
   if(input){
      GaussPenta gauss;
      for(int iv=0;iv<NUMVERTICES;iv++){
         gauss.GaussVertex(iv);
         input->GetInputValue(&pvalue[iv],&gauss);
      }
   }
   else{
      for(int iv=0;iv<NUMVERTICES;iv++) pvalue[iv] = default_value;
   }
}
/*}}}*/
void       Penta::GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
 
   /*Checks in debugging mode*/
   _assert_(pvalue);
 
   /*What type of finite element are we dealing with?*/
   int fe       = this->FiniteElement();
   int numnodes = this->GetNumberOfNodes();
 
   /* Start looping on the number of vertices: */
   if(input){
      GaussPenta gauss;
      for(int iv=0;iv<numnodes;iv++){
         gauss.GaussNode(fe,iv);
         input->GetInputValue(&pvalue[iv],&gauss);
      }
   }
   else{
      for(int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
   }
}
/*}}}*/
DatasetInput2* Penta::GetDatasetInput2(int inputenum){/*{{{*/
 
   DatasetInput2* datasetinput = this->inputs2->GetDatasetInput2(inputenum);
   if(!datasetinput) return NULL;
 
   for(int i=0;i<datasetinput->GetNumIds();i++){
 
      PentaInput2* input = datasetinput->GetPentaInputByOffset(i); _assert_(input);
 
      /*Intermediaries*/
      int numindices;
      int indices[30]; /*Max numnodes*/
 
      /*Check interpolation*/
      int interpolation = input->GetInterpolation();
      if(interpolation==P1Enum){
         numindices = 6;
         for(int i=0;i<6;i++) indices[i] = vertices[i]->lid;
         input->Serve(numindices,&indices[0]);
      }
      else{
         input->Serve(this->lid,this->GetNumberOfNodes(interpolation));
      }
 
   }
 
   return datasetinput;
}/*}}}*/
void       Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
 
   Input2* input=this->GetInput2(enumtype);
   if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
 
   int index = this->GetNodeIndex(node);
 
   GaussPenta* gauss=new GaussPenta();
   gauss->GaussNode(this->element_type,index);
 
   input->GetInputValue(pvalue,gauss);
   delete gauss;
}
/*}}}*/
void       Penta::GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){/*{{{*/
 
   Input2* input=this->GetInput2(enumtype);
   if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
 
   int index = this->GetVertexIndex(vertex);
 
   GaussPenta* gauss=new GaussPenta();
   gauss->GaussVertex(index);
 
   input->GetInputValue(pvalue,gauss);
   delete gauss;
}
/*}}}*/
Penta*     Penta::GetLowerPenta(void){/*{{{*/
 
   Penta* lower_penta=NULL;
 
   lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
 
   return lower_penta;
}
/*}}}*/
void       Penta::GetLevelsetIntersectionBase(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
 
   /* GetLevelsetIntersection computes:
    * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
    * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
 
   /*Intermediaries*/
   int i, numiceverts, numnoiceverts;
   int ind0, ind1, lastindex;
   int indices_ice[NUMVERTICES2D],indices_noice[NUMVERTICES2D];
   IssmDouble lsf[NUMVERTICES];
   int* indices = xNew<int>(NUMVERTICES2D);
 
   /*Retrieve all inputs and parameters*/
   Element::GetInputListOnVertices(&lsf[0],levelset_enum);
 
   /* Determine distribution of ice over element.
    * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
   lastindex=0;
   for(i=0;i<NUMVERTICES2D;i++){ // go backwards along vertices, and check for sign change
      ind0=(NUMVERTICES2D-i)%NUMVERTICES2D;
      ind1=(ind0-1+NUMVERTICES2D)%NUMVERTICES2D;
      if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
         if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
            lastindex=ind1;
         else
            lastindex=ind0;
         break;
      }
   }
 
   numiceverts=0;
   numnoiceverts=0;
   for(i=0;i<NUMVERTICES2D;i++){
      ind0=(lastindex+i)%NUMVERTICES2D;
      if(lsf[i]<=level){
         indices_ice[numiceverts]=i;
         numiceverts++;
      }
      else{
         indices_noice[numnoiceverts]=i;
         numnoiceverts++;
      }
   }
   //merge indices
   for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
   for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
 
   switch (numiceverts){
      case 0: // no vertex has ice: element is ice free, no intersection
         for(i=0;i<2;i++)
            fraction[i]=0.;
         break;
      case 1: // one vertex has ice:
         for(i=0;i<2;i++){
            fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
         }
         break;
      case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
         for(i=0;i<2;i++){
            fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
         }
         break;
      case NUMVERTICES2D: // all vertices have ice: return triangle area
         for(i=0;i<2;i++)
            fraction[i]=1.;
         break;
      default:
         _error_("Wrong number of ice vertices in Penta::GetLevelsetIntersectionBase!");
         break;
   }
 
   *pindices=indices;
   *pnumiceverts=numiceverts;
}
/*}}}*/
int        Penta::GetVertexIndex(Vertex* vertex){/*{{{*/
   _assert_(vertices);
   for(int i=0;i<NUMVERTICES;i++){
      if(vertex==vertices[i])
       return i;
   }
   _error_("Vertex provided not found among element vertices");
}
/*}}}*/
int        Penta::GetNumberOfNodes(void){/*{{{*/
   return this->NumberofNodes(this->element_type);
}
/*}}}*/
int        Penta::GetNumberOfNodes(int enum_type){/*{{{*/
   return this->NumberofNodes(enum_type);
}
/*}}}*/
int        Penta::GetNumberOfVertices(void){/*{{{*/
   return NUMVERTICES;
}
/*}}}*/
Penta*     Penta::GetUpperPenta(void){/*{{{*/
 
   Penta* upper_penta=NULL;
 
   upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
 
   return upper_penta;
}
/*}}}*/
void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
 
   /*Get out if this is not an element input*/
   if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
 
   /*Prepare index list*/
   int idlist[NUMVERTICES];
   GradientIndexing(&idlist[0],control_index);
 
   /*Get input (either in element or material)*/
   if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
   ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);
 
   /*Intermediaries*/
   int numindices;
   int indices[NUMVERTICES];
 
   /*Check interpolation*/
   int interpolation = input->GetInterpolation();
   switch(interpolation){
      case P1Enum:
         numindices = NUMVERTICES;
         for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
         input->Serve(numindices,&indices[0]);
         break;
      default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
   }
 
   /* Start looping on the number of vertices: */
   IssmDouble values[NUMVERTICES];
   Gauss*gauss=this->NewGauss();
   for(int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
      input->GetInputValue(&values[iv],gauss);
   }
   delete gauss;
 
   vector->SetValues(NUMVERTICES,idlist,&values[0],INS_VAL);
}
/*}}}*/
void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset){/*{{{*/
 
   int* idlist = NULL;
   IssmDouble* values = NULL;
 
   /*Get input*/
   ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);
 
   /*Check what input we are dealing with*/
   switch(input->ObjectEnum()){
      case PentaInput2Enum:
              {
               PentaInput2* pentainput = xDynamicCast<PentaInput2*>(input);
               if(pentainput->GetInputInterpolationType()!=P1Enum) _error_("not supported yet");
 
               /*Create list of indices and values for global vector*/
               idlist = xNew<int>(NUMVERTICES);
               values = xNew<IssmDouble>(NUMVERTICES);
               GradientIndexing(&idlist[0],control_index);
               for(int i=0;i<NUMVERTICES;i++){
                  values[i] = pentainput->element_values[i];
               }
               vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
               break;
              }
 
            case TransientInputEnum:
              {
               _error_("not implemented (see Tria)");
               break;
              }
            default: _error_("input "<<input->ObjectEnum()<<" not supported yet");
 
      }
         xDelete<int>(idlist);
         xDelete<IssmDouble>(values);
   }
/*}}}*/
void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
 
   IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
   ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
 
   /*Assign output pointer*/
   *pxyz_list = xyz_list;
 
}/*}}}*/
void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
 
   IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
   ::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
 
   /*Assign output pointer*/
   *pxyz_list = xyz_list;
 
}/*}}}*/
IssmDouble Penta::GroundedArea(bool scaled){/*{{{*/
 
   /*Intermediaries*/
   int         domaintype;
   IssmDouble  phi,base_area,scalefactor,groundedarea;
   IssmDouble  xyz_list[NUMVERTICES][3];
 
   if(!IsIceInElement() || !IsOnBase())return 0.;
 
   /*Get problem dimension*/
   this->FindParam(&domaintype,DomainTypeEnum);
   if(domaintype!=Domain3DEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   phi=this->GetGroundedPortion(&xyz_list[0][0]);
   base_area= 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
 
   groundedarea=phi*base_area;
 
   if(scaled==true){
      Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
      scalefactor_input->GetInputAverage(&scalefactor);
      groundedarea=groundedarea*scalefactor;
   }
   /*Clean up and return*/
   return groundedarea;
}
/*}}}*/
IssmDouble Penta::GroundinglineMassFlux(bool scaled){/*{{{*/
 
   /*Make sure there is a grounding line here*/
   if(!IsOnBase()) return 0;
   if(!IsIceInElement()) return 0;
   if(!IsZeroLevelset(MaskOceanLevelsetEnum)) return 0;
 
   /*Scaled not implemented yet...*/
   _assert_(!scaled);
 
   int               domaintype,index1,index2;
   const IssmPDouble epsilon = 1.e-15;
   IssmDouble        s1,s2;
   IssmDouble        gl[NUMVERTICES];
   IssmDouble        xyz_front[2][3];
 
   IssmDouble *xyz_list = NULL;
   this->GetVerticesCoordinates(&xyz_list);
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&gl[0],MaskOceanLevelsetEnum);
 
   /*Be sure that values are not zero*/
   if(gl[0]==0.) gl[0]=gl[0]+epsilon;
   if(gl[1]==0.) gl[1]=gl[1]+epsilon;
   if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
   int pt1 = 0;
   int pt2 = 1;
   if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 
      /*Portion of the segments*/
      s1=gl[2]/(gl[2]-gl[1]);
      s2=gl[2]/(gl[2]-gl[0]);
      if(gl[2]<0.){
         pt1 = 1; pt2 = 0;
      }
      xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
      xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
      xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
      xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
      xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
      xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
   }
   else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
 
      /*Portion of the segments*/
      s1=gl[0]/(gl[0]-gl[1]);
      s2=gl[0]/(gl[0]-gl[2]);
      if(gl[0]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
      xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
      xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
      xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
      xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
      xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
   }
   else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
 
      /*Portion of the segments*/
      s1=gl[1]/(gl[1]-gl[0]);
      s2=gl[1]/(gl[1]-gl[2]);
      if(gl[1]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
      xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
      xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
      xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
      xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
      xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
   }
   else{
      _error_("case not possible");
   }
 
 
   /*Some checks in debugging mode*/
   _assert_(s1>=0 && s1<=1.);
   _assert_(s2>=0 && s2<=1.);
 
   /*Get normal vector*/
   IssmDouble normal[3];
   this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
   normal[0] = -normal[0];
   normal[1] = -normal[1];
 
   this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 
   /*Get inputs*/
   IssmDouble flux = 0.;
   IssmDouble vx,vy,thickness,Jdet;
   IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
   Input2 *thickness_input = this->GetInput2(ThicknessEnum); _assert_(thickness_input);
   Input2 *vx_input        = this->GetInput2(VxAverageEnum); _assert_(vx_input);
   Input2 *vy_input        = this->GetInput2(VyAverageEnum); _assert_(vy_input);
 
   /*Start looping on Gaussian points*/
   Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
      thickness_input->GetInputValue(&thickness,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
 
      flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
   }
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
   return flux;
 
}
/*}}}*/
IssmDouble Penta::IcefrontMassFluxLevelset(bool scaled){/*{{{*/
 
   /*Make sure there is an ice front here*/
   if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
 
   /*Scaled not implemented yet...*/
   _assert_(!scaled);
 
   int               domaintype,index1,index2;
   const IssmPDouble epsilon = 1.e-15;
   IssmDouble        s1,s2;
   IssmDouble        gl[NUMVERTICES];
   IssmDouble        xyz_front[2][3];
 
   IssmDouble *xyz_list = NULL;
   this->GetVerticesCoordinates(&xyz_list);
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
   /*Be sure that values are not zero*/
   if(gl[0]==0.) gl[0]=gl[0]+epsilon;
   if(gl[1]==0.) gl[1]=gl[1]+epsilon;
   if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
   int pt1 = 0;
   int pt2 = 1;
   if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 
      /*Portion of the segments*/
      s1=gl[2]/(gl[2]-gl[1]);
      s2=gl[2]/(gl[2]-gl[0]);
      if(gl[2]<0.){
         pt1 = 1; pt2 = 0;
      }
      xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
      xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
      xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
      xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
      xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
      xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
   }
   else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
 
      /*Portion of the segments*/
      s1=gl[0]/(gl[0]-gl[1]);
      s2=gl[0]/(gl[0]-gl[2]);
      if(gl[0]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
      xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
      xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
      xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
      xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
      xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
   }
   else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
 
      /*Portion of the segments*/
      s1=gl[1]/(gl[1]-gl[0]);
      s2=gl[1]/(gl[1]-gl[2]);
      if(gl[1]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
      xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
      xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
      xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
      xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
      xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
   }
   else{
      _error_("case not possible");
   }
 
 
   /*Some checks in debugging mode*/
   _assert_(s1>=0 && s1<=1.);
   _assert_(s2>=0 && s2<=1.);
 
   /*Get normal vector*/
   IssmDouble normal[3];
   this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
   normal[0] = -normal[0];
   normal[1] = -normal[1];
 
   this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 
   /*Get inputs*/
   IssmDouble flux = 0.;
   IssmDouble vx,vy,thickness,Jdet;
   IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
   Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
   Input2* vx_input=NULL;
   Input2* vy_input=NULL;
   vx_input=this->GetInput2(VxAverageEnum); _assert_(vx_input);
   vy_input=this->GetInput2(VyAverageEnum); _assert_(vy_input);
 
   /*Start looping on Gaussian points*/
   Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
      thickness_input->GetInputValue(&thickness,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
 
      flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   delete gauss;
   return flux;
 
}
/*}}}*/
IssmDouble Penta::IceVolume(bool scaled){/*{{{*/
 
   /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
   IssmDouble base,height,scalefactor;
   IssmDouble xyz_list[NUMVERTICES][3];
 
   if(!IsIceInElement())return 0;
 
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*First calculate the area of the base (cross section triangle)
    * http://en.wikipedia.org/wiki/Pentangle
    * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
   base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
 
   if(scaled==true){ //scale for area projection correction
      Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
      scalefactor_input->GetInputAverage(&scalefactor);
      base=base*scalefactor;
   }
 
   /*Now get the average height*/
   height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
 
   /*Return: */
   return base*height;
}
/*}}}*/
IssmDouble Penta::IceVolumeAboveFloatation(bool scaled){/*{{{*/
 
   /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
   IssmDouble rho_ice,rho_water;
   IssmDouble base,bed,surface,bathymetry,scalefactor;
   IssmDouble xyz_list[NUMVERTICES][3];
 
   if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
 
   rho_ice=FindParam(MaterialsRhoIceEnum);
   rho_water=FindParam(MaterialsRhoSeawaterEnum);
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*First calculate the area of the base (cross section triangle)
    * http://en.wikipedia.org/wiki/Pentangle
    * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
   base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
   if(scaled==true){
      Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
      scalefactor_input->GetInputAverage(&scalefactor);
      base=base*scalefactor;
   }
 
   /*Now get the average height above floatation*/
   Input2* surface_input    = this->GetInput2(SurfaceEnum);    _assert_(surface_input);
   Input2* base_input        = this->GetInput2(BaseEnum);        _assert_(base_input);
   Input2* bed_input = this->GetInput2(BedEnum); _assert_(bed_input);
   if(!bed_input) _error_("Could not find bed");
   surface_input->GetInputAverage(&surface);
   base_input->GetInputAverage(&bed);
   bed_input->GetInputAverage(&bathymetry);
 
   /*Return: */
   return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
}
/*}}}*/
void       Penta::InputDepthAverageAtBase(int original_enum,int average_enum){/*{{{*/
 
   IssmDouble  Jdet,value;
   IssmDouble  xyz_list[NUMVERTICES][3];
   IssmDouble  xyz_list_line[2][3];
   IssmDouble  total[NUMVERTICES]       = {0.};
   int         lidlist[NUMVERTICES];
   IssmDouble  intz[NUMVERTICES]        = {0.};
   Input2     *original_input           = NULL;
   Input2     *depth_averaged_input     = NULL;
 
   /*Are we on the base? If not, return*/
   if(!IsOnBase()) return;
 
   /*Now follow all the upper element from the base to the surface to integrate the input*/
   Penta* penta = this;
   int    step  = 0;
   Gauss* gauss[3];
   for(int iv=0;iv<3;iv++) gauss[iv] = penta->NewGaussLine(iv,iv+3,3);
 
   for(;;){
 
      /*Step1: Get original input (to be depth-avegaged): */
      original_input=penta->GetInput2(original_enum);
      if(!original_input) _error_("could not find input with enum " << EnumToStringx(original_enum));
 
      /*Step2: Create element thickness input*/
      ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
      for(int iv=0;iv<3;iv++){
         /*Get segment length*/
         for(int i=0;i<3;i++){
            xyz_list_line[0][i]=xyz_list[iv][i];
            xyz_list_line[1][i]=xyz_list[iv+3][i];
         }
         /*Integrate over edge*/
         for(int ig=gauss[iv]->begin();ig<gauss[iv]->end();ig++){
            gauss[iv]->GaussPoint(ig);
            penta->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss[iv]);
            original_input->GetInputValue(&value,gauss[iv]);
            total[iv] += value*Jdet*gauss[iv]->weight;
            intz[iv]  += Jdet*gauss[iv]->weight;
         }
      }
 
      /*Stop if we have reached the surface, else, take upper penta*/
      if(penta->IsOnSurface()) break;
 
      /* get upper Penta*/
      penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
      step++;
   }
   for(int iv=0;iv<3;iv++) delete gauss[iv];
 
   /*Now we only need to divide the depth integrated input by the total thickness!*/
   for(int iv=0;iv<3;iv++){
      total[iv  ] = total[iv]/intz[iv];
      total[iv+3] = total[iv];
   }
   GetVerticesLidList(&lidlist[0]);
   switch(original_input->ObjectEnum()){
      case PentaInputEnum:
      case PentaInput2Enum:
      case ControlInputEnum:
         this->inputs2->SetPentaInput(average_enum,P1Enum,NUMVERTICES,lidlist,&total[0]);
         break;
      default:
         _error_("Interpolation " << EnumToStringx(original_input->ObjectEnum()) << " not supported yet");
   }
}
/*}}}*/
void       Penta::DatasetInputExtrude(int enum_type,int start){/*{{{*/
 
   _assert_(start==-1 || start==+1);
   _assert_(this->inputs2);
 
   /*Are we on the the boundary we want to be?*/
   if(start==-1 && !IsOnBase())    return;
   if(start==+1 && !IsOnSurface()) return;
 
   /*Get original input*/
   DatasetInput2* dinput = this->inputs2->GetDatasetInput2(enum_type);
 
   int lidlist[NUMVERTICES];
   this->GetVerticesLidList(&lidlist[0]);
 
   for(int id=0;id<dinput->GetNumIds();id++){
 
      PentaInput2* pentainput = dinput->GetPentaInputByOffset(id);
      pentainput->Serve(NUMVERTICES,&lidlist[0]);
 
      if(pentainput->GetInterpolation()==P1Enum){
 
         /*Extrude values first*/
         IssmDouble extrudedvalues[NUMVERTICES];
         this->GetInputListOnVertices(&extrudedvalues[0],pentainput,0.);
 
         if(start==-1){
            for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i+NUMVERTICES2D]=extrudedvalues[i];
         }
         else{
            for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i]=extrudedvalues[i+NUMVERTICES2D];
         }
 
         /*Propagate to other Pentas*/
         Penta* penta=this;
         for(;;){
 
            /*Add input of the basal element to penta->inputs*/
            int vertexlids[NUMVERTICES];
            penta->GetVerticesLidList(&vertexlids[0]);
            pentainput->SetInput(P1Enum,NUMVERTICES,&vertexlids[0],&extrudedvalues[0]);
 
            /*Stop if we have reached the surface/base*/
            if(start==-1 && penta->IsOnSurface()) break;
            if(start==+1 && penta->IsOnBase())    break;
 
            /*get upper/lower Penta*/
            if(start==-1) penta=penta->GetUpperPenta();
            else          penta=penta->GetLowerPenta();
            _assert_(penta->Id()!=this->id);
         }
      }
      else{
         _error_("not implemented yet");
      }
   }
}
/*}}}*/
void       Penta::ControlInputExtrude(int enum_type,int start){/*{{{*/
 
   _assert_(start==-1 || start==+1);
   _assert_(this->inputs2);
 
   /*Are we on the the boundary we want to be?*/
   if(start==-1 && !IsOnBase())    return;
   if(start==+1 && !IsOnSurface()) return;
 
   /*Get original input*/
   ElementInput2* input  = this->inputs2->GetControlInput2Data(enum_type,"value");
   if(input->ObjectEnum()!=PentaInput2Enum) _error_("not supported yet");
   PentaInput2* pentainput = xDynamicCast<PentaInput2*>(input);
   ElementInput2* input2 = this->inputs2->GetControlInput2Data(enum_type,"savedvalues");
   if(input->ObjectEnum()!=PentaInput2Enum) _error_("not supported yet");
   PentaInput2* pentainput2= xDynamicCast<PentaInput2*>(input2);
   /*FIXME: this should not be necessary*/
   ElementInput2* input3 = this->inputs2->GetControlInput2Data(enum_type,"gradient");
   if(input->ObjectEnum()!=PentaInput2Enum) _error_("not supported yet");
   PentaInput2* pentainput3= xDynamicCast<PentaInput2*>(input3);
 
   int lidlist[NUMVERTICES];
   this->GetVerticesLidList(&lidlist[0]);
   pentainput->Serve(NUMVERTICES,&lidlist[0]);
   pentainput2->Serve(NUMVERTICES,&lidlist[0]);
   pentainput3->Serve(NUMVERTICES,&lidlist[0]);
 
   if(pentainput->GetInterpolation()==P1Enum){
 
      /*Extrude values first*/
      IssmDouble extrudedvalues[NUMVERTICES];
      IssmDouble extrudedvalues2[NUMVERTICES];
      IssmDouble extrudedvalues3[NUMVERTICES];
 
      this->GetInputListOnVertices(&extrudedvalues[0],pentainput,0.);
      this->GetInputListOnVertices(&extrudedvalues2[0],pentainput2,0.);
      this->GetInputListOnVertices(&extrudedvalues3[0],pentainput3,0.);
 
      if(start==-1){
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i+NUMVERTICES2D]=extrudedvalues[i];
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues2[i+NUMVERTICES2D]=extrudedvalues2[i];
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues3[i+NUMVERTICES2D]=extrudedvalues3[i]/2.; /*FIXME: this is just for NR*/
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues3[i]=extrudedvalues3[i]/2.; /*FIXME: this is just for NR*/
      }
      else{
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i]=extrudedvalues[i+NUMVERTICES2D];
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues2[i]=extrudedvalues2[i+NUMVERTICES2D];
      }
 
      /*Propagate to other Pentas*/
      Penta* penta=this;
      for(;;){
 
         if(penta->IsOnSurface() && start==-1){ /*FIXME: this is just for NR*/
            for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues3[i+NUMVERTICES2D]=0.;
         }
 
         /*Add input of the basal element to penta->inputs*/
         int vertexlids[NUMVERTICES];
         penta->GetVerticesLidList(&vertexlids[0]);
         pentainput->SetInput(P1Enum,NUMVERTICES,&vertexlids[0],&extrudedvalues[0]);
         pentainput2->SetInput(P1Enum,NUMVERTICES,&vertexlids[0],&extrudedvalues2[0]);
         if(start==-1 && !penta->IsOnBase()){
            pentainput3->SetInput(P1Enum,NUMVERTICES,&vertexlids[0],&extrudedvalues3[0]);
         }
 
         /*Stop if we have reached the surface/base*/
         if(start==-1 && penta->IsOnSurface()) break;
         if(start==+1 && penta->IsOnBase())    break;
 
         /*get upper/lower Penta*/
         if(start==-1) penta=penta->GetUpperPenta();
         else          penta=penta->GetLowerPenta();
         _assert_(penta->Id()!=this->id);
      }
   }
   else{
      _error_("not implemented yet");
   }
}
/*}}}*/
void       Penta::InputExtrude(int enum_type,int start){/*{{{*/
 
   _assert_(start==-1 || start==+1);
   _assert_(this->inputs2);
 
   /*Are we on the the boundary we want to be?*/
   if(start==-1 && !IsOnBase())    return;
   if(start==+1 && !IsOnSurface()) return;
 
 
   /*Get original input*/
   Input2* input = this->GetInput2(enum_type);
   if(input->ObjectEnum()!=PentaInput2Enum) _error_("not supported yet");
   PentaInput2* pentainput = xDynamicCast<PentaInput2*>(input);
 
   if(pentainput->GetInterpolation()==P1Enum || pentainput->GetInterpolation()==P1DGEnum){
      /*Extrude values first*/
      IssmDouble extrudedvalues[NUMVERTICES];
 
      Element::GetInputListOnVertices(&extrudedvalues[0],enum_type);
      if(start==-1){
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i+NUMVERTICES2D]=extrudedvalues[i];
      }
      else{
         for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i]=extrudedvalues[i+NUMVERTICES2D];
      }
 
      /*Propagate to other Pentas*/
      Penta* penta=this;
      for(;;){
 
         /*Add input of the basal element to penta->inputs*/
         penta->AddInput2(enum_type,&extrudedvalues[0],pentainput->GetInterpolation());
 
         /*Stop if we have reached the surface/base*/
         if(start==-1 && penta->IsOnSurface()) break;
         if(start==+1 && penta->IsOnBase())    break;
 
         /*get upper/lower Penta*/
         if(start==-1) penta=penta->GetUpperPenta();
         else          penta=penta->GetLowerPenta();
         _assert_(penta->Id()!=this->id);
      }
   }
   else{
      _error_("interpolation "<<EnumToStringx(pentainput->GetInterpolation())<<" not implemented yet (while trying to extrude "<<EnumToStringx(enum_type)<<")");
   }
}
/*}}}*/
void       Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ /*{{{*/
 
   /*Intermediaries*/
   int         i,j;
   int         penta_vertex_ids[NUMVERTICES];
   IssmDouble  nodeinputs[NUMVERTICES];
   IssmDouble  cmmininputs[NUMVERTICES];
   IssmDouble  cmmaxinputs[NUMVERTICES];
 
   IssmDouble  yts;
   bool    control_analysis;
   char**  controls = NULL;
   int     num_control_type,num_responses;
 
   /*Fetch parameters: */
   iomodel->FindConstant(&yts,"md.constants.yts");
   iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
   if(control_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters");
   if(control_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions");
 
   /*Recover vertices ids needed to initialize inputs*/
   _assert_(iomodel->elements);
   for(i=0;i<NUMVERTICES;i++){
      penta_vertex_ids[i]=iomodel->elements[NUMVERTICES*index+i]; //ids for vertices are in the elements array from Matlab
   }
}
/*}}}*/
void       Penta::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){/*{{{*/
 
   /*Intermediary*/
   int* doflist = NULL;
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = this->NumberofNodes(this->element_type);
 
   /*Fetch dof list and allocate solution vector*/
   GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
   IssmDouble* values = xNew<IssmDouble>(numnodes);
 
   /*Use the dof list to index into the solution vector: */
   for(int i=0;i<numnodes;i++){
      values[i]=solution[doflist[i]];
      if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
      if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
   }
 
   /*Add input to the element: */
   this->AddInput2(enum_type,values,this->element_type);
 
   /*Free ressources:*/
   xDelete<IssmDouble>(values);
   xDelete<int>(doflist);
}
/*}}}*/
void       Penta::InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){/*{{{*/
 
   const int  numdof   = NUMVERTICES;
   const int  numdof2d = NUMVERTICES2D;
 
   IssmDouble  values[numdof];
   int*    doflist = NULL;
   Penta  *penta   = NULL;
 
   /*If not on bed, return*/
   if (!IsOnBase()) return;
 
   /*Get dof list: */
   GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
 
   /*Use the dof list to index into the solution vector and extrude it */
   for(int i=0;i<numdof2d;i++){
      values[i]         =solution[doflist[i]];
      values[i+numdof2d]=values[i];
      if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
      if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
   }
 
   /*Start looping over all elements above current element and update all inputs*/
   penta=this;
   for(;;){
      /*Add input to the element: */
      penta->AddInput2(enum_type,values,P1Enum);
 
      /*Stop if we have reached the surface*/
      if (penta->IsOnSurface()) break;
 
      /* get upper Penta*/
      penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
   }
 
   /*Free ressources:*/
   xDelete<int>(doflist);
}
/*}}}*/
void       Penta::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
 
   const int   numdof         = NUMVERTICES;
   int        *doflist        = NULL;
   IssmDouble  values[numdof];
   int         lidlist[NUMVERTICES];
 
   GetVerticesLidList(&lidlist[0]);
 
   /*Check that name is an element input*/
   if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
 
   switch(type){
      case VertexLIdEnum:
         for (int i=0;i<NUMVERTICES;i++){
            values[i]=vector[this->vertices[i]->Lid()];
         }
         /*update input*/
         inputs2->SetPentaInput(name,P1Enum,NUMVERTICES,lidlist,values);
         return;
 
      case VertexPIdEnum:
         for (int i=0;i<NUMVERTICES;i++){
            values[i]=vector[this->vertices[i]->Pid()];
         }
         /*update input*/
         inputs2->SetPentaInput(name,P1Enum,NUMVERTICES,lidlist,values);
         return;
 
      case VertexSIdEnum:
         for (int i=0;i<NUMVERTICES;i++){
            values[i]=vector[this->vertices[i]->Sid()];
         }
         /*update input*/
         inputs2->SetPentaInput(name,P1Enum,NUMVERTICES,lidlist,values);
         return;
 
      case NodesEnum:
         /*Get dof list: */
         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
         /*Use the dof list to index into the vector: */
         for(int i=0;i<numdof;i++){
            values[i]=vector[doflist[i]];
            if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
            if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
         }
         /*Add input to the element: */
         inputs2->SetPentaInput(name,P1Enum,NUMVERTICES,lidlist,values);
 
         /*Free ressources:*/
         xDelete<int>(doflist);
         return;
 
     case NodeSIdEnum:
         for(int i=0;i<NUMVERTICES;i++){
            values[i]=vector[nodes[i]->Sid()];
            if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
            if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
         }
         /*Add input to the element: */
         inputs2->SetPentaInput(name,P1Enum,NUMVERTICES,lidlist,values);
 
         /*Free ressources:*/
         xDelete<int>(doflist);
         return;
 
     default:
         _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
   }
}
/*}}}*/
bool       Penta::IsIcefront(void){/*{{{*/
 
   bool isicefront;
   int i,nrice;
   IssmDouble ls[NUMVERTICES];
 
   /*Retrieve all inputs and parameters*/
   Element::GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
   /* If only one vertex has ice, there is an ice front here */
   isicefront=false;
   if(IsIceInElement()){
      nrice=0;
      for(i=0;i<NUMVERTICES2D;i++)
         if(ls[i]<0.) nrice++;
      if(nrice==1) isicefront= true;
   }
   return isicefront;
}/*}}}*/
bool       Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
 
   int  i;
   bool shelf=false;
 
   for(i=0;i<NUMVERTICES;i++){
      if (flags[vertices[i]->Pid()]<0.){
         shelf=true;
         break;
      }
   }
   return shelf;
}
/*}}}*/
bool       Penta::IsZeroLevelset(int levelset_enum){/*{{{*/
 
   bool        iszerols;
   IssmDouble  ls[NUMVERTICES];
 
   /*Retrieve all inputs and parameters*/
   Element::GetInputListOnVertices(&ls[0],levelset_enum);
 
   /*If the level set has always same sign, there is no ice front here*/
   iszerols = false;
   if(IsIceInElement()){
      if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
         iszerols = true;
      }
   }
   return iszerols;
}
/*}}}*/
void       Penta::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetJacobianDeterminant(pJdet,xyz_list,(GaussPenta*)gauss);
 
}
/*}}}*/
void       Penta::JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetTriaJacobianDeterminant(pJdet,xyz_list_base,(GaussPenta*)gauss);
 
}
/*}}}*/
void       Penta::JacobianDeterminantLine(IssmDouble* pJdet,IssmDouble* xyz_list_line,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetSegmentJacobianDeterminant(pJdet,xyz_list_line,(GaussPenta*)gauss);
 
}
/*}}}*/
void       Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
 
}
/*}}}*/
void       Penta::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetTriaJacobianDeterminant(pJdet,xyz_list_top,(GaussPenta*)gauss);
 
}
/*}}}*/
IssmDouble Penta::MassFlux(IssmDouble* segment){/*{{{*/
 
   IssmDouble mass_flux=0;
 
   if(!IsOnBase()) return mass_flux;
 
   /*Depth Averaging Vx and Vy*/
   this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 
   /*Spawn Tria element from the base of the Penta: */
   Tria* tria=(Tria*)SpawnTria(0,1,2);
   mass_flux=tria->MassFlux(segment);
   delete tria->material; delete tria;
 
   /*clean up and return*/
   return mass_flux;
}
/*}}}*/
IssmDouble Penta::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
 
   IssmDouble mass_flux=0;
 
   if(!IsOnBase()) return mass_flux;
 
   /*Depth Averaging Vx and Vy*/
   this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 
   /*Spawn Tria element from the base of the Penta: */
   Tria* tria=(Tria*)SpawnTria(0,1,2);
   mass_flux=tria->MassFlux(x1,y1,x2,y2,segment_id);
   delete tria->material; delete tria;
 
   /*clean up and return*/
   return mass_flux;
}
/*}}}*/
IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){/*{{{*/
   /*Return the minimum lenght of the nine egdes of the penta*/
 
   int    i,node0,node1;
   int    edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges
   IssmDouble length;
   IssmDouble minlength=-1;
 
   for(i=0;i<9;i++){
      /*Find the two nodes for this edge*/
      node0=edges[i][0];
      node1=edges[i][1];
 
      /*Compute the length of this edge and compare it to the minimal length*/
      length=sqrt(pow(xyz_list[node0*3+0]-xyz_list[node1*3+0],2)+pow(xyz_list[node0*3+1]-xyz_list[node1*3+1],2)+pow(xyz_list[node0*3+2]-xyz_list[node1*3+2],2));
      if(length<minlength || minlength<0) minlength=length;
   }
 
   return minlength;
}
/*}}}*/
Gauss*     Penta::NewGauss(void){/*{{{*/
   return new GaussPenta();
}
/*}}}*/
Gauss*     Penta::NewGauss(int order){/*{{{*/
   return new GaussPenta(order,order);
}
/*}}}*/
Gauss*     Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/
 
   IssmDouble  area_coordinates[4][3];
 
   GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
 
   return new GaussPenta(area_coordinates,order_horiz,order_vert);
}
/*}}}*/
Gauss*     Penta::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){/*{{{*/
   return new GaussPenta(point1,fraction1,fraction2,mainlyfloating,order);
}
/*}}}*/
Gauss*     Penta::NewGaussBase(int order){/*{{{*/
   return new GaussPenta(0,1,2,order);
}
/*}}}*/
Gauss*     Penta::NewGaussBase(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz){/*{{{*/
 
   IssmDouble  area_coordinates[2][3];
 
   GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
 
   return new GaussPenta(area_coordinates,order_horiz);
}
/*}}}*/
Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
   return new GaussPenta(vertex1,vertex2,order);
}
/*}}}*/
Gauss*     Penta::NewGaussTop(int order){/*{{{*/
   return new GaussPenta(3,4,5,order);
}
/*}}}*/
void       Penta::NodalFunctions(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
 
}
/*}}}*/
void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
 
}
/*}}}*/
void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
 
}
/*}}}*/
void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
 
}
/*}}}*/
void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
 
}
/*}}}*/
void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
 
}
/*}}}*/
void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
 
}
/*}}}*/
void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
 
}
/*}}}*/
void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
 
}
/*}}}*/
void       Penta::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->TensorInterpolation());
 
}
/*}}}*/
int        Penta::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){/*{{{*/
 
   int i;
   int found=0;
   IssmDouble value;
   GaussPenta* gauss=NULL;
 
   /*First, serarch the input: */
   Input2* data=this->GetInput2(natureofdataenum);
 
   /*figure out if we have the vertex id: */
   found=0;
   for(i=0;i<NUMVERTICES;i++){
      if(index==vertices[i]->Id()){
         /*Do we have natureofdataenum in our inputs? :*/
         if(data){
            /*ok, we are good. retrieve value of input at vertex :*/
            gauss=new GaussPenta(); gauss->GaussVertex(i);
            data->GetInputValue(&value,gauss);
            found=1;
            break;
         }
      }
   }
 
   delete gauss;
   if(found)*pvalue=value;
   return found;
}
/*}}}*/
void       Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){/*{{{*/
 
   IssmDouble v13[3],v23[3];
   IssmDouble normal[3];
   IssmDouble normal_norm;
 
   for(int i=0;i<3;i++){
      v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
      v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
   }
 
   normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
   normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
   normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
   normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
 
   /*Bed normal is opposite to surface normal*/
   bed_normal[0]=-normal[0]/normal_norm;
   bed_normal[1]=-normal[1]/normal_norm;
   bed_normal[2]=-normal[2]/normal_norm;
}
/*}}}*/
void       Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
 
   /*Build unit outward pointing vector*/
   IssmDouble AB[3];
   IssmDouble AC[3];
   IssmDouble norm;
 
   AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
   AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
   AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
   AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
   AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
   AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
 
   cross(normal,AB,AC);
   norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
 
   for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
}
/*}}}*/
void       Penta::NormalSectionBase(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
 
   /*Build unit outward pointing vector*/
   IssmDouble vector[2];
   IssmDouble norm;
 
   vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
   vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
 
   norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
 
   normal[0]= + vector[1]/norm;
   normal[1]= - vector[0]/norm;
}
/*}}}*/
void       Penta::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/
 
   int i;
   IssmDouble v13[3],v23[3];
   IssmDouble normal[3];
   IssmDouble normal_norm;
 
   for (i=0;i<3;i++){
      v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
      v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
   }
 
   normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
   normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
   normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
   normal_norm=sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
 
   top_normal[0]=normal[0]/normal_norm;
   top_normal[1]=normal[1]/normal_norm;
   top_normal[2]=normal[2]/normal_norm;
}
/*}}}*/
int        Penta::NumberofNodesPressure(void){/*{{{*/
   return PentaRef::NumberofNodes(this->PressureInterpolation());
}
/*}}}*/
int        Penta::NumberofNodesVelocity(void){/*{{{*/
   return PentaRef::NumberofNodes(this->VelocityInterpolation());
}
/*}}}*/
int        Penta::ObjectEnum(void){/*{{{*/
 
   return PentaEnum;
 
}
/*}}}*/
void       Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
 
   IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
   IssmDouble  bed_hydro;
   IssmDouble  rho_water,rho_ice,density;
 
   /*material parameters: */
   rho_water=FindParam(MaterialsRhoSeawaterEnum);
   rho_ice=FindParam(MaterialsRhoIceEnum);
   density=rho_ice/rho_water;
   Element::GetInputListOnVertices(&h[0],ThicknessEnum);
   Element::GetInputListOnVertices(&r[0],BedEnum);
   Element::GetInputListOnVertices(&gl[0],MaskOceanLevelsetEnum);
 
   /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
   for(int i=0;i<NUMVERTICES;i++){
      /*Find if grounded vertices want to start floating*/
      if (gl[i]>0.){
         bed_hydro=-density*h[i];
         if(bed_hydro>r[i]){
            /*Vertex that could potentially unground, flag it*/
            potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
         }
      }
   }
}
/*}}}*/
int        Penta::PressureInterpolation(void){/*{{{*/
   return PentaRef::PressureInterpolation(this->element_type);
}
/*}}}*/
void       Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/
 
   int analysis_type;
   parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
   if(pe){
      if(analysis_type==StressbalanceAnalysisEnum){
         if(this->element_type==MINIcondensedEnum){
            int approximation;
            this->GetInput2Value(&approximation,ApproximationEnum);
            if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
               //Do nothing, condensation already done in PVectorCoupling
            }
            else{
               int indices[3]={18,19,20};
               pe->StaticCondensation(Ke,3,&indices[0]);
            }
         }
         else if(this->element_type==P1bubblecondensedEnum){
            int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
            int offset = 0;
            for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
            int* indices=xNew<int>(size);
            for(int i=0;i<size;i++) indices[i] = offset+i;
            pe->StaticCondensation(Ke,size,indices);
            xDelete<int>(indices);
         }
      }
   }
 
   if(Ke){
      if(analysis_type==StressbalanceAnalysisEnum){
         int approximation;
         this->GetInput2Value(&approximation,ApproximationEnum);
         if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
            //Do nothing condensatino already done for Stokes part
         }
         else{
            if(this->element_type==MINIcondensedEnum){
               int indices[3]={18,19,20};
               Ke->StaticCondensation(3,&indices[0]);
            }
            else if(this->element_type==P1bubblecondensedEnum){
               int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
               int offset = 0;
               for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
               int* indices=xNew<int>(size);
               for(int i=0;i<size;i++) indices[i] = offset+i;
               Ke->StaticCondensation(size,indices);
               xDelete<int>(indices);
            }
         }
      }
   }
}
/*}}}*/
void       Penta::ResetFSBasalBoundaryCondition(void){/*{{{*/
 
   int          approximation;
   int          numindices;
   int         *indices = NULL;
   IssmDouble   slopex,slopey,groundedice;
   IssmDouble   xz_plane[6];
   IssmDouble*  vertexapproximation= NULL;
 
   /*For FS only: we want the CS to be tangential to the bedrock*/
   this->GetInput2Value(&approximation,ApproximationEnum);
   if(!IsOnBase() || (approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return;
 
   /*Get number of nodes for velocity only and base*/
   BasalNodeIndices(&numindices,&indices,this->VelocityInterpolation());
 
   /*Get inputs*/
   Input2* slopex_input=this->GetInput2(BedSlopeXEnum); _assert_(slopex_input);
   Input2* slopey_input=this->GetInput2(BedSlopeYEnum); _assert_(slopey_input);
   Input2* groundedicelevelset_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(groundedicelevelset_input);
 
   /*Loop over basal nodes and update their CS*/
   GaussPenta* gauss = new GaussPenta();
   for(int i=0;i<numindices;i++){//FIXME
      gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
 
      slopex_input->GetInputValue(&slopex,gauss);
      slopey_input->GetInputValue(&slopey,gauss);
      groundedicelevelset_input->GetInputValue(&groundedice,gauss);
 
      /*New X axis          New Z axis*/
      xz_plane[0]=1.;       xz_plane[3]=-slopex;
      xz_plane[1]=0.;       xz_plane[4]=-slopey;
      xz_plane[2]=slopex;   xz_plane[5]=1.;
 
      if(groundedice>=0){
         if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
            this->nodes[indices[i]]->DofInSSet(2); //vz
         }
         else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
            this->nodes[indices[i]]->DofInSSet(4); //vz
         }
         else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
      }
      else{
         if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
            this->nodes[indices[i]]->DofInFSet(2); //vz
         }
         else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
            this->nodes[indices[i]]->DofInFSet(4); //vz
         }
         else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
      }
 
      XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
   }
 
   /*cleanup*/
   xDelete<int>(indices);
   delete gauss;
}
/*}}}*/
void       Penta::ResetHooks(){/*{{{*/
 
   if(this->nodes) xDelete<Node*>(this->nodes);
   this->nodes=NULL;
   this->vertices=NULL;
   this->material=NULL;
   this->verticalneighbors=NULL;
   this->parameters=NULL;
 
   //deal with ElementHook mother class
   for(int i=0;i<this->numanalyses;i++) if(this->hnodes[i]) this->hnodes[i]->reset();
   this->hvertices->reset();
   this->hmaterial->reset();
   if(this->hneighbors) this->hneighbors->reset();
 
}
/*}}}*/
void       Penta::RignotMeltParameterization(){/*{{{*/
 
   if(!this->IsOnBase()) return;
 
   IssmDouble A, B, alpha, beta;
   IssmDouble bed,qsg,qsg_basin,TF,yts;
   int numbasins;
   IssmDouble basinid[NUMVERTICES];
   IssmDouble* basin_icefront_area=NULL;
 
   /* Coefficients */
   A    = 3e-4;
   B    = 0.15;
   alpha = 0.39;
   beta = 1.18;
 
   /*Get inputs*/
   Input2* bed_input = this->GetInput2(BedEnum);                     _assert_(bed_input);
   Input2* qsg_input = this->GetInput2(FrontalForcingsSubglacialDischargeEnum);      _assert_(qsg_input);
   Input2* TF_input  = this->GetInput2(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
   Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
 
   this->FindParam(&yts, ConstantsYtsEnum);
   this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
   this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
 
   IssmDouble meltrates[NUMVERTICES2D];  //frontal melt-rate
 
   /* Start looping on the number of vertices: */
   GaussPenta* gauss=new GaussPenta();
   for(int iv=0;iv<NUMVERTICES2D;iv++){
      gauss->GaussVertex(iv);
 
      /* Get variables */
      bed_input->GetInputValue(&bed,gauss);
      qsg_input->GetInputValue(&qsg,gauss);
      TF_input->GetInputValue(&TF,gauss);
 
      if(basinid[iv]==0 || basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
      else{
         /* change the unit of qsg (m^3/d -> m/d) with ice front area */
         qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
 
         /* calculate melt rates */
         meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
      }
 
      if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
      if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
   }
 
   /*Add input*/
   this->AddInput2(CalvingMeltingrateEnum,&meltrates[0],P1Enum);
 
   this->InputExtrude(CalvingMeltingrateEnum,-1);
 
   /*Cleanup and return*/
   xDelete<IssmDouble>(basin_icefront_area);
   delete gauss;
}
/*}}}*/
void       Penta::SetElementInput(int enum_in,IssmDouble value){/*{{{*/
 
   this->SetElementInput(this->inputs2,enum_in,value);
 
}
/*}}}*/
void       Penta::SetElementInput(Inputs2* inputs2,int enum_in,IssmDouble value){/*{{{*/
 
   _assert_(inputs2);
   inputs2->SetPentaInput(enum_in,P0Enum,this->lid,value);
 
}
/*}}}*/
void       Penta::SetElementInput(Inputs2* inputs2,int numindices,int* indices,IssmDouble* values,int enum_in){/*{{{*/
 
   _assert_(inputs2);
   inputs2->SetPentaInput(enum_in,P1Enum,numindices,indices,values);
 
}
/*}}}*/
void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N, int M){/*{{{*/
 
   IssmDouble  values[NUMVERTICES];
   int         vertexpidlist[NUMVERTICES],control_init;
 
   /*Specific case for depth averaged quantities*/
   control_init=control_enum;
   if(control_enum==MaterialsRheologyBbarEnum){
      control_enum=MaterialsRheologyBEnum;
      if(!IsOnBase()) return;
   }
   if(control_enum==DamageDbarEnum){
      control_enum=DamageDEnum;
      if(!IsOnBase()) return;
   }
 
   /*Get out if this is not an element input*/
   if(!IsInputEnum(control_enum)) return;
 
   /*Prepare index list*/
   GradientIndexing(&vertexpidlist[0],control_index);
 
   /*Get values on vertices*/
   for(int i=0;i<NUMVERTICES;i++){
      values[i]=vector[vertexpidlist[i]];
   }
   _error_("not implemented");
   //Input* new_input = new PentaInput(control_enum,values,P1Enum);
   Input2* input=(Input2*)this->GetInput2(control_enum);   _assert_(input);
   if(input->ObjectEnum()!=ControlInputEnum){
      _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
   }
 
   //((ControlInput*)input)->SetInput(new_input);
 
   if(control_init==MaterialsRheologyBbarEnum){
      this->InputExtrude(control_enum,-1);
   }
   if(control_init==DamageDbarEnum){
      this->InputExtrude(control_enum,-1);
   }
}
/*}}}*/
void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
 
   IssmDouble  values[NUMVERTICES];
   int         lidlist[NUMVERTICES];
   int         idlist[NUMVERTICES],control_init;
 
   /*Specific case for depth averaged quantities*/
   control_init=control_enum;
   if(control_enum==MaterialsRheologyBbarEnum){
      control_enum=MaterialsRheologyBEnum;
      if(!IsOnBase()) return;
   }
   if(control_enum==DamageDbarEnum){
      control_enum=DamageDEnum;
      if(!IsOnBase()) return;
   }
 
   /*Get Domain type*/
   int domaintype;
   parameters->FindParam(&domaintype,DomainTypeEnum);
 
   /*Specific case for depth averaged quantities*/
   if(domaintype==Domain2DverticalEnum){
      if(control_enum==MaterialsRheologyBbarEnum){
         control_enum=MaterialsRheologyBEnum;
         if(!IsOnBase()) return;
      }
      if(control_enum==DamageDbarEnum){
         control_enum=DamageDEnum;
         if(!IsOnBase()) return;
      }
   }
 
   /*Get out if this is not an element input*/
   if(!IsInputEnum(control_enum)) return;
 
   /*prepare index list*/
   this->GetVerticesLidList(&lidlist[0]);
   GradientIndexing(&idlist[0],control_index);
 
   /*Get values on vertices*/
   for(int i=0;i<NUMVERTICES;i++){
      values[i]=vector[idlist[i]];
   }
 
   /*Set Input*/
   ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"value");   _assert_(input);
   input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
   if(control_init==MaterialsRheologyBbarEnum){
      this->ControlInputExtrude(control_enum,-1);
   }
   if(control_init==DamageDbarEnum){
      this->ControlInputExtrude(control_enum,-1);
   }
}
/*}}}*/
void       Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
 
   /*go into parameters and get the analysis_counter: */
   int analysis_counter;
   parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
 
   /*Get Element type*/
   this->element_type=this->element_type_list[analysis_counter];
 
   /*Pick up nodes*/
   if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
   else this->nodes=NULL;
 
}
/*}}}*/
void       Penta::SetTemporaryElementType(int element_type_in){/*{{{*/
   this->element_type=element_type_in;
}
/*}}}*/
Element*   Penta::SpawnBasalElement(void){/*{{{*/
 
   _assert_(this->IsOnBase());
 
   switch(this->material->ObjectEnum()){
      case MaticeEnum:
         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
         if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
         if(this->material->IsEnhanced())this->InputDepthAverageAtBase(MaterialsRheologyEEnum,MaterialsRheologyEbarEnum);
         break;
      case MatestarEnum:
         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
         this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
         this->InputDepthAverageAtBase(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum);
         break;
      default:
         _error_("not supported yet");
   }
   if(this->GetInput2(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   if(this->GetInput2(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
   if(this->GetInput2(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
   if(this->GetInput2(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
 
   Tria* tria=(Tria*)SpawnTria(0,1,2);
 
   return tria;
}
/*}}}*/
Element*   Penta::SpawnTopElement(void){/*{{{*/
 
   _assert_(this->IsOnSurface());
 
   Tria* tria=(Tria*)SpawnTria(3,4,5);
 
   return tria;
}
/*}}}*/
Tria*      Penta::SpawnTria(int index1,int index2,int index3){/*{{{*/
 
   int analysis_counter;
 
   /*go into parameters and get the analysis_counter: */
   this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
 
   /*Create Tria*/
   Tria* tria=new Tria();
   tria->id=this->id;
   tria->sid=this->sid;
   tria->lid=this->lid;
   tria->parameters=this->parameters;
   tria->inputs2=this->inputs2;
   tria->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
   this->SpawnTriaHook(xDynamicCast<ElementHook*>(tria),index1,index2,index3);
 
   if(index1==0 && index2==1 && index3==2){
      tria->iscollapsed = 1;
   }
   else if(index1==3 && index2==4 && index3==5){
      tria->iscollapsed = 2;
   }
 
   /*Spawn material*/
   tria->material=(Material*)this->material->copy2(tria);
 
   /*recover nodes, material*/
   tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp();
   tria->vertices=(Vertex**)tria->hvertices->deliverp();
 
   /*Return new Tria*/
   return tria;
}
/*}}}*/
IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/
   /*Compute stabilization parameter*/
   /*kappa=thermalconductivity/(rho_ice*heatcapacity) for thermal model*/
   /*kappa=enthalpydiffusionparameter for enthalpy model*/
 
   IssmDouble normu;
   IssmDouble tau_parameter;
 
   normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
   if(normu*diameter/(3*2*kappa)<1){
      tau_parameter=pow(diameter,2)/(3*2*2*kappa);
   }
   else tau_parameter=diameter/(2*normu);
 
   return tau_parameter;
}/*}}}*/
void Penta::StabilizationParameterAnisotropic(IssmDouble* tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa){/*{{{*/
   /*Compute stabilization parameter*/
   /*kappa=thermalconductivity/(rho_ice*heatcapacity) for thermal model*/
   /*kappa=enthalpydiffusionparameter for enthalpy model*/
 
   IssmDouble normu,hk,C,area,p;
 
   /* compute tau for the horizontal direction */
   p=2.; C=3.;
   normu=pow(pow(u,p)+pow(v,p)+pow(w,p),1./p);
   hk=sqrt(pow(hx,2)+pow(hy,2));
 
   if(normu*hk/(C*2*kappa)<1){
      tau_parameter_anisotropic[0]=pow(hk,2)/(C*2*2*kappa);
   }
   else tau_parameter_anisotropic[0]=hk/(2*normu);
 
   /* compute tau for the vertical direction */
   hk=hz;
   if(normu*hk/(C*2*kappa)<1){
      tau_parameter_anisotropic[1]=pow(hk,2)/(C*2*2*kappa);
   }
   else{
      tau_parameter_anisotropic[1]=hk/(2*normu);
   }
}
/*}}}*/
void       Penta::StrainRateparallel(){/*{{{*/
 
   IssmDouble *xyz_list = NULL;
   IssmDouble  epsilon[6];
   GaussPenta* gauss=NULL;
   IssmDouble  vx,vy,vel;
   IssmDouble  strainxx;
   IssmDouble  strainxy;
   IssmDouble  strainyy;
   IssmDouble  strainparallel[NUMVERTICES];
 
   /* Get node coordinates and dof list: */
   this->GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will need*/
   Input2* vx_input=this->GetInput2(VxEnum);                                  _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum);                                  _assert_(vy_input);
   Input2* vz_input=this->GetInput2(VzEnum);                                  _assert_(vz_input);
 
   /* Start looping on the number of vertices: */
   gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      /* Get the value we need*/
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      vel=vx*vx+vy*vy;
 
      /*Compute strain rate and viscosity: */
      this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
      strainxx=epsilon[0];
      strainyy=epsilon[1];
      strainxy=epsilon[3];
 
      /*strainparallel= Strain rate along the ice flow direction */
      strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-14);
   }
 
   /*Add input*/
   this->AddInput2(StrainRateparallelEnum,&strainparallel[0],P1DGEnum);
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
}
/*}}}*/
void       Penta::StrainRateperpendicular(){/*{{{*/
 
   IssmDouble *xyz_list = NULL;
   IssmDouble  epsilon[6];
   GaussPenta* gauss=NULL;
   IssmDouble  vx,vy,vel;
   IssmDouble  strainxx;
   IssmDouble  strainxy;
   IssmDouble  strainyy;
   IssmDouble  strainperpendicular[NUMVERTICES];
 
   /* Get node coordinates and dof list: */
   this->GetVerticesCoordinates(&xyz_list);
 
   /*Retrieve all inputs we will need*/
   Input2* vx_input=this->GetInput2(VxEnum);                                  _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum);                                  _assert_(vy_input);
   Input2* vz_input=this->GetInput2(VzEnum);                                  _assert_(vz_input);
 
   /* Start looping on the number of vertices: */
   gauss=new GaussPenta();
   for (int iv=0;iv<NUMVERTICES;iv++){
      gauss->GaussVertex(iv);
 
      /* Get the value we need*/
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      vel=vx*vx+vy*vy;
 
      /*Compute strain rate and viscosity: */
      this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
      strainxx=epsilon[0];
      strainyy=epsilon[1];
      strainxy=epsilon[3];
 
      /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
      strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-14);
   }
 
   /*Add input*/
   this->AddInput2(StrainRateperpendicularEnum,&strainperpendicular[0],P1DGEnum);
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
}
/*}}}*/
void       Penta::StressIntensityFactor(){/*{{{*/
 
   /* Check if we are on the base */
   if(!IsOnBase()) return;
 
   IssmDouble  ki[6]={0.};
   IssmDouble  const_grav=9.81;
   IssmDouble  rho_ice=900;
   IssmDouble  rho_water=1000;
   IssmDouble  Jdet[3];
   IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
 
   Penta* penta=this;
   for(;;){
 
      IssmDouble  xyz_list[NUMVERTICES][3];
      /* Get node coordinates and dof list: */
      ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
 
      Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
      Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
      Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
 
      /*Retrieve all inputs we will need*/
      Input2* vx_input=this->GetInput2(VxEnum);                                  _assert_(vx_input);
      Input2* vy_input=this->GetInput2(VyEnum);                                  _assert_(vy_input);
      Input2* vel_input=this->GetInput2(VelEnum);                                _assert_(vel_input);
      Input2* pressure_input=this->GetInput2(PressureEnum);                      _assert_(pressure_input);
      Input2* deviaxx_input=this->GetInput2(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
      Input2* deviaxy_input=this->GetInput2(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
      Input2* deviayy_input=this->GetInput2(DeviatoricStressyyEnum);             _assert_(deviayy_input);
      Input2* surface_input=this->GetInput2(SurfaceEnum);                        _assert_(surface_input);
      Input2* thickness_input=this->GetInput2(ThicknessEnum);                    _assert_(thickness_input);
 
      /* Start looping on the number of 2D vertices: */
      for(int ig=0;ig<3;ig++){
         GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
         for (int iv=gauss->begin();iv<gauss->end();iv++){
            gauss->GaussPoint(iv);
 
            /* Get the value we need*/
            pressure_input->GetInputValue(&pressure,gauss);
            vx_input->GetInputValue(&vx,gauss);
            vy_input->GetInputValue(&vy,gauss);
            vel_input->GetInputValue(&vel,gauss);
            deviaxx_input->GetInputValue(&deviaxx,gauss);
            deviaxy_input->GetInputValue(&deviaxy,gauss);
            deviayy_input->GetInputValue(&deviayy,gauss);
            surface_input->GetInputValue(&water_depth,gauss);
            thickness_input->GetInputValue(&thickness,gauss);
            prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
 
            /*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
            stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
 
            if(prof<water_depth&prof<thickness){
               /* Compute the local stress intensity factor*/
               ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
            }
         }
         delete gauss;
      }
 
      /*Stop if we have reached the surface/base*/
      if(penta->IsOnSurface()) break;
 
      /*get upper Penta*/
      penta=penta->GetUpperPenta();
      _assert_(penta->Id()!=this->id);
   }
 
   /*Add input*/
   this->AddInput2(StressIntensityFactorEnum,&ki[0],P1Enum);
   this->InputExtrude(StressIntensityFactorEnum,-1);
}
/*}}}*/
IssmDouble Penta::SurfaceArea(void){/*{{{*/
 
   int    approximation;
   IssmDouble S;
   Tria*  tria=NULL;
 
   /*retrieve inputs :*/
   this->GetInput2Value(&approximation,ApproximationEnum);
 
   /*If on water, return 0: */
   if(!IsIceInElement())return 0;
 
   /*Bail out if this element if:
    * -> Non SSA not on the surface
    * -> SSA (2d model) and not on bed) */
   if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
      return 0;
   }
   else if (approximation==SSAApproximationEnum){
 
      /*This element should be collapsed into a tria element at its base. Create this tria element,
       * and compute SurfaceArea*/
      tria=(Tria*)SpawnTria(0,1,2);
      S=tria->SurfaceArea();
      delete tria->material; delete tria;
      return S;
   }
   else{
 
      tria=(Tria*)SpawnTria(3,4,5);
      S=tria->SurfaceArea();
      delete tria->material; delete tria;
      return S;
   }
}
/*}}}*/
IssmDouble Penta::TimeAdapt(void){/*{{{*/
 
   /*intermediary: */
   IssmDouble C;
   IssmDouble xyz_list[NUMVERTICES][3];
 
   /*get CFL coefficient:*/
   this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum);
 
   /*Get for Vx and Vy, the max of abs value: */
   Input2* vx_input = this->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input = this->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vz_input = this->GetInput2(VzEnum); _assert_(vz_input);
   IssmDouble maxabsvx = vx_input->GetInputMaxAbs();
   IssmDouble maxabsvy = vy_input->GetInputMaxAbs();
   IssmDouble maxabsvz = vz_input->GetInputMaxAbs();
 
   /* Get node coordinates and dof list: */
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   IssmDouble minx=xyz_list[0][0];
   IssmDouble maxx=xyz_list[0][0];
   IssmDouble miny=xyz_list[0][1];
   IssmDouble maxy=xyz_list[0][1];
   IssmDouble minz=xyz_list[0][2];
   IssmDouble maxz=xyz_list[0][2];
 
   for(int i=1;i<NUMVERTICES;i++){
      if(xyz_list[i][0]<minx) minx=xyz_list[i][0];
      if(xyz_list[i][0]>maxx) maxx=xyz_list[i][0];
      if(xyz_list[i][1]<miny) miny=xyz_list[i][1];
      if(xyz_list[i][1]>maxy) maxy=xyz_list[i][1];
      if(xyz_list[i][2]<minz) minz=xyz_list[i][2];
      if(xyz_list[i][2]>maxz) maxz=xyz_list[i][2];
   }
   IssmDouble dx=maxx-minx;
   IssmDouble dy=maxy-miny;
   IssmDouble dz=maxz-minz;
 
   /*CFL criterion: */
   IssmDouble dt = C/(maxabsvx/dx+maxabsvy/dy+maxabsvz/dz);
 
   return dt;
}/*}}}*/
IssmDouble Penta::TotalCalvingFluxLevelset(bool scaled){/*{{{*/
 
   /*Make sure there is an ice front here*/
   if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
 
   /*Scaled not implemented yet...*/
   _assert_(!scaled);
 
   int               domaintype,index1,index2;
   const IssmPDouble epsilon = 1.e-15;
   IssmDouble        s1,s2;
   IssmDouble        gl[NUMVERTICES];
   IssmDouble        xyz_front[2][3];
 
   IssmDouble *xyz_list = NULL;
   this->GetVerticesCoordinates(&xyz_list);
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
   /*Be sure that values are not zero*/
   if(gl[0]==0.) gl[0]=gl[0]+epsilon;
   if(gl[1]==0.) gl[1]=gl[1]+epsilon;
   if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
   int pt1 = 0;
   int pt2 = 1;
   if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 
      /*Portion of the segments*/
      s1=gl[2]/(gl[2]-gl[1]);
      s2=gl[2]/(gl[2]-gl[0]);
      if(gl[2]<0.){
         pt1 = 1; pt2 = 0;
      }
      xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
      xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
      xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
      xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
      xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
      xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
   }
   else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
 
      /*Portion of the segments*/
      s1=gl[0]/(gl[0]-gl[1]);
      s2=gl[0]/(gl[0]-gl[2]);
      if(gl[0]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
      xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
      xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
      xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
      xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
      xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
   }
   else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
 
      /*Portion of the segments*/
      s1=gl[1]/(gl[1]-gl[0]);
      s2=gl[1]/(gl[1]-gl[2]);
      if(gl[1]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
      xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
      xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
      xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
      xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
      xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
   }
   else{
      _error_("case not possible");
   }
 
 
   /*Some checks in debugging mode*/
   _assert_(s1>=0 && s1<=1.);
   _assert_(s2>=0 && s2<=1.);
 
   /*Get normal vector*/
   IssmDouble normal[3];
   this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
   normal[0] = -normal[0];
   normal[1] = -normal[1];
 
   /*Get inputs*/
   IssmDouble flux = 0.;
   IssmDouble calvingratex,calvingratey,thickness,Jdet;
   IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
   Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
   Input2* calvingratex_input=NULL;
   Input2* calvingratey_input=NULL;
   calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
   calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
 
   /*Start looping on Gaussian points*/
   Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
      thickness_input->GetInputValue(&thickness,gauss);
      calvingratex_input->GetInputValue(&calvingratex,gauss);
      calvingratey_input->GetInputValue(&calvingratey,gauss);
      this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
 
      flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
   }
 
   /*Clean up and return*/
   delete gauss;
   return flux;
 
}
/*}}}*/
IssmDouble Penta::TotalCalvingMeltingFluxLevelset(bool scaled){/*{{{*/
 
   /*Make sure there is an ice front here*/
   if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
 
   /*Scaled not implemented yet...*/
   _assert_(!scaled);
 
   int               domaintype,index1,index2;
   const IssmPDouble epsilon = 1.e-15;
   IssmDouble        s1,s2;
   IssmDouble        gl[NUMVERTICES];
   IssmDouble        xyz_front[2][3];
 
   IssmDouble *xyz_list = NULL;
   this->GetVerticesCoordinates(&xyz_list);
 
   /*Recover parameters and values*/
   Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
   /*Be sure that values are not zero*/
   if(gl[0]==0.) gl[0]=gl[0]+epsilon;
   if(gl[1]==0.) gl[1]=gl[1]+epsilon;
   if(gl[2]==0.) gl[2]=gl[2]+epsilon;
 
   int pt1 = 0;
   int pt2 = 1;
   if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 
      /*Portion of the segments*/
      s1=gl[2]/(gl[2]-gl[1]);
      s2=gl[2]/(gl[2]-gl[0]);
      if(gl[2]<0.){
         pt1 = 1; pt2 = 0;
      }
      xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
      xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
      xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
      xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
      xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
      xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
   }
   else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
 
      /*Portion of the segments*/
      s1=gl[0]/(gl[0]-gl[1]);
      s2=gl[0]/(gl[0]-gl[2]);
      if(gl[0]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
      xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
      xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
      xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
      xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
      xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
   }
   else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
 
      /*Portion of the segments*/
      s1=gl[1]/(gl[1]-gl[0]);
      s2=gl[1]/(gl[1]-gl[2]);
      if(gl[1]<0.){
         pt1 = 1; pt2 = 0;
      }
 
      xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
      xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
      xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
      xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
      xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
      xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
   }
   else{
      _error_("case not possible");
   }
 
 
   /*Some checks in debugging mode*/
   _assert_(s1>=0 && s1<=1.);
   _assert_(s2>=0 && s2<=1.);
 
   /*Get normal vector*/
   IssmDouble normal[3];
   this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
   normal[0] = -normal[0];
   normal[1] = -normal[1];
 
   this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
   this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
 
   /*Get inputs*/
   IssmDouble flux = 0.;
   IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
   IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
   Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
   Input2* calvingratex_input=NULL;
   Input2* calvingratey_input=NULL;
   Input2* vx_input=NULL;
   Input2* vy_input=NULL;
   Input2* meltingrate_input=NULL;
   calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
   calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
   vx_input=this->GetInput2(VxAverageEnum); _assert_(vx_input);
   vy_input=this->GetInput2(VyAverageEnum); _assert_(vy_input);
   meltingrate_input=this->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
 
   /*Start looping on Gaussian points*/
   Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
      thickness_input->GetInputValue(&thickness,gauss);
      calvingratex_input->GetInputValue(&calvingratex,gauss);
      calvingratey_input->GetInputValue(&calvingratey,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      vel=vx*vx+vy*vy;
      meltingrate_input->GetInputValue(&meltingrate,gauss);
      meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
      meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
      this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
 
      flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
   }
 
   /*Clean up and return*/
   delete gauss;
   return flux;
 
}
/*}}}*/
IssmDouble Penta::TotalFloatingBmb(bool scaled){/*{{{*/
 
   /*The fbmb[kg yr-1] of one element is area[m2] * melting_rate [kg m^-2 yr^-1]*/
   int        point1;
   bool       mainlyfloating;
   IssmDouble fbmb=0;
   IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet,scalefactor;
   IssmDouble Total_Fbmb=0;
   IssmDouble xyz_list[NUMVERTICES][3];
   Gauss*     gauss     = NULL;
 
   if(!IsIceInElement() || !IsOnBase())return 0;
 
   /*Get material parameters :*/
   rho_ice=FindParam(MaterialsRhoIceEnum);
   Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
   Input2* gllevelset_input = this->GetInput2(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
   Input2* scalefactor_input = NULL;
   if(scaled==true){
      scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
   }
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
   /* Start  looping on the number of gaussian points: */
   gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
      this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
      floatingmelt_input->GetInputValue(&floatingmelt,gauss);
      if(scaled==true){
         scalefactor_input->GetInputValue(&scalefactor,gauss);
      }
      else scalefactor=1;
      fbmb+=floatingmelt*Jdet*gauss->weight*scalefactor;
   }
 
   Total_Fbmb=rho_ice*fbmb;           // from volume to mass
 
   /*Return: */
   delete gauss;
   return Total_Fbmb;
}
/*}}}*/
IssmDouble Penta::TotalGroundedBmb(bool scaled){/*{{{*/
 
   /*The gbmb[kg yr-1] of one element is area[m2] * gounded melting rate [kg m^-2 yr^-1]*/
   int        point1;
   bool       mainlyfloating;
   IssmDouble gbmb=0;
   IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet,scalefactor;
   IssmDouble Total_Gbmb=0;
   IssmDouble xyz_list[NUMVERTICES][3];
   Gauss*     gauss     = NULL;
 
   if(!IsIceInElement() || !IsOnBase())return 0;
 
   /*Get material parameters :*/
   rho_ice=FindParam(MaterialsRhoIceEnum);
   Input2* groundedmelt_input = this->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input);
   Input2* gllevelset_input   = this->GetInput2(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
   Input2* scalefactor_input  = NULL;
   if(scaled==true){
      scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
   }
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
   /* Start  looping on the number of gaussian points: */
   gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
      this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
      groundedmelt_input->GetInputValue(&groundedmelt,gauss);
      if(scaled==true){
         scalefactor_input->GetInputValue(&scalefactor,gauss);
      }
      else scalefactor=1;
      gbmb+=groundedmelt*Jdet*gauss->weight*scalefactor;
   }
 
   Total_Gbmb=rho_ice*gbmb;           // from volume to mass
 
   /*Return: */
   delete gauss;
   return Total_Gbmb;
}
/*}}}*/
IssmDouble Penta::TotalSmb(bool scaled){/*{{{*/
 
   /*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
   IssmDouble base,smb,rho_ice,scalefactor;
   IssmDouble Total_Smb=0;
   IssmDouble xyz_list[NUMVERTICES][3];
 
   /*Get material parameters :*/
   rho_ice=FindParam(MaterialsRhoIceEnum);
 
   if(!IsIceInElement() || !IsOnSurface()) return 0.;
 
   ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
   /*First calculate the area of the base (cross section triangle)
    * http://en.wikipedia.org/wiki/Triangle
    * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
   base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
 
   /*Now get the average SMB over the element*/
   Input2* smb_input = this->GetInput2(SmbMassBalanceEnum); _assert_(smb_input);
 
   smb_input->GetInputAverage(&smb);
   if(scaled==true){
      Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
      scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
   }
   else{
      scalefactor=1.;
   }
   Total_Smb=rho_ice*base*smb*scalefactor;// smb on element in kg s-1
 
   /*Return: */
   return Total_Smb;
}
/*}}}*/
void       Penta::Update(Inputs2* inputs2,int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ /*{{{*/
 
   /*Intermediaries*/
   int        i;
   int        penta_vertex_ids[6];
   IssmDouble nodeinputs[6];
   IssmDouble yts;
   bool       dakota_analysis;
   int        numnodes;
   int*       penta_node_ids = NULL;
 
   /*Fetch parameters: */
   iomodel->FindConstant(&yts,"md.constants.yts");
   iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
 
   /*Checks if debuging*/
   _assert_(iomodel->elements);
   _assert_(index==this->sid);
 
   /*Recover element type*/
   this->element_type_list[analysis_counter]=finiteelement_type;
 
   /*Recover vertices ids needed to initialize inputs*/
   for(i=0;i<6;i++) penta_vertex_ids[i]=iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
 
   /*Recover nodes ids needed to initialize the node hook.*/
   switch(finiteelement_type){
      case P1Enum:
         numnodes         = 6;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[0]=iomodel->elements[6*index+0];
         penta_node_ids[1]=iomodel->elements[6*index+1];
         penta_node_ids[2]=iomodel->elements[6*index+2];
         penta_node_ids[3]=iomodel->elements[6*index+3];
         penta_node_ids[4]=iomodel->elements[6*index+4];
         penta_node_ids[5]=iomodel->elements[6*index+5];
         break;
      case P1bubbleEnum: case P1bubblecondensedEnum:
         numnodes         = 7;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[0]=iomodel->elements[6*index+0];
         penta_node_ids[1]=iomodel->elements[6*index+1];
         penta_node_ids[2]=iomodel->elements[6*index+2];
         penta_node_ids[3]=iomodel->elements[6*index+3];
         penta_node_ids[4]=iomodel->elements[6*index+4];
         penta_node_ids[5]=iomodel->elements[6*index+5];
         penta_node_ids[6]=iomodel->numberofvertices+index+1;
         break;
      case P1xP2Enum:
         numnodes         = 9;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+2]+1;
         break;
      case P1xP3Enum:
         numnodes         = 12;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2;
         penta_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2;
         penta_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2;
         break;
      case P2xP1Enum:
         numnodes         = 12;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+5]+1;
         break;
      case P1xP4Enum:
         numnodes         = 15;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0]; /*Vertex 1*/
         penta_node_ids[ 1]=iomodel->elements[6*index+1]; /*Vertex 2*/
         penta_node_ids[ 2]=iomodel->elements[6*index+2]; /*Vertex 3*/
         penta_node_ids[ 3]=iomodel->elements[6*index+3]; /*Vertex 4*/
         penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/
         penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/
         penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /*mid vertical edge 1*/
         penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /*mid vertical edge 2*/
         penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /*mid vertical edge 3*/
         penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 1/4 vertical edge 1*/
         penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 1/4 vertical edge 2*/
         penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 1/4 vertical edge 3*/
         penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+3; /* 3/4 vertical edge 1*/
         penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+3; /* 3/4 vertical edge 2*/
         penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+3; /* 3/4 vertical edge 3*/
         break;
      case P2xP4Enum:
         numnodes         = 30;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0]; /*Vertex 1*/
         penta_node_ids[ 1]=iomodel->elements[6*index+1]; /*Vertex 2*/
         penta_node_ids[ 2]=iomodel->elements[6*index+2]; /*Vertex 3*/
         penta_node_ids[ 3]=iomodel->elements[6*index+3]; /*Vertex 4*/
         penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/
         penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /* 1/4 vertical edge 1*/
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /* 1/4 vertical edge 2*/
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /* 1/4 vertical edge 3*/
         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 3/4 vertical edge 1*/
         penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 3/4 vertical edge 2*/
         penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 3/4 vertical edge 3*/
         penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; /* 1/4 vertical face 1*/
         penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; /* 1/4 vertical face 2*/
         penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; /* 1/4 vertical face 3*/
         penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+2; /* 2/4 vertical face 1*/
         penta_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+2; /* 2/4 vertical face 2*/
         penta_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+2; /* 2/4 vertical face 3*/
         penta_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+3; /* 3/4 vertical face 1*/
         penta_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+3; /* 3/4 vertical face 2*/
         penta_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+3; /* 3/4 vertical face 3*/
         break;
      case P2Enum:
         numnodes         = 18;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+0]+1;
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+1]+1;
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+2]+1;
         break;
      case P2bubbleEnum: case P2bubblecondensedEnum:
         numnodes         = 19;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
         break;
      case P1P1Enum: case P1P1GLSEnum:
         numnodes         = 12;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
 
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elements[6*index+0];
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elements[6*index+1];
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elements[6*index+2];
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elements[6*index+3];
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elements[6*index+4];
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elements[6*index+5];
         break;
      case MINIEnum: case MINIcondensedEnum:
         numnodes         = 13;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+index+1;
 
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+0];
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+1];
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+2];
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+3];
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+4];
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+5];
         break;
      case TaylorHoodEnum:
         numnodes         = 24;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+0]+1;
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+1]+1;
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+2]+1;
 
         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+0];
         penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+1];
         penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+2];
         penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+3];
         penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+4];
         penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+5];
         break;
      case LATaylorHoodEnum:
         numnodes         = 18;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
         break;
      case OneLayerP4zEnum:
         numnodes         = 30+6;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0]; /*Vertex 1*/
         penta_node_ids[ 1]=iomodel->elements[6*index+1]; /*Vertex 2*/
         penta_node_ids[ 2]=iomodel->elements[6*index+2]; /*Vertex 3*/
         penta_node_ids[ 3]=iomodel->elements[6*index+3]; /*Vertex 4*/
         penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/
         penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /* 1/4 vertical edge 1*/
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /* 1/4 vertical edge 2*/
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /* 1/4 vertical edge 3*/
         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 3/4 vertical edge 1*/
         penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 3/4 vertical edge 2*/
         penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 3/4 vertical edge 3*/
         penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; /* 1/4 vertical face 1*/
         penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; /* 1/4 vertical face 2*/
         penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; /* 1/4 vertical face 3*/
         penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+2; /* 2/4 vertical face 1*/
         penta_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+2; /* 2/4 vertical face 2*/
         penta_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+2; /* 2/4 vertical face 3*/
         penta_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+3; /* 3/4 vertical face 1*/
         penta_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+3; /* 3/4 vertical face 2*/
         penta_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+3; /* 3/4 vertical face 3*/
 
         penta_node_ids[30]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+0];
         penta_node_ids[31]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+1];
         penta_node_ids[32]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+2];
         penta_node_ids[33]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+3];
         penta_node_ids[34]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+4];
         penta_node_ids[35]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+5];
         break;
      case CrouzeixRaviartEnum:
         numnodes         = 25;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
 
         penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+1;
         penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+2;
         penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+3;
         penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+4;
         penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+5;
         penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+6;
         break;
      case LACrouzeixRaviartEnum:
         numnodes         = 19;
         penta_node_ids   = xNew<int>(numnodes);
         penta_node_ids[ 0]=iomodel->elements[6*index+0];
         penta_node_ids[ 1]=iomodel->elements[6*index+1];
         penta_node_ids[ 2]=iomodel->elements[6*index+2];
         penta_node_ids[ 3]=iomodel->elements[6*index+3];
         penta_node_ids[ 4]=iomodel->elements[6*index+4];
         penta_node_ids[ 5]=iomodel->elements[6*index+5];
         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
         penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
         penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
         penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
         break;
      default:
         _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
   }
 
   /*hooks: */
   this->SetHookNodes(penta_node_ids,numnodes,analysis_counter); this->nodes=NULL;
   xDelete<int>(penta_node_ids);
 
   /*Defaults if not provided in iomodel*/
   switch(analysis_type){
 
      case StressbalanceAnalysisEnum:
         _assert_(iomodel->Data("md.flowequation.element_equation"));
 
         if((IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[index])))==HOFSApproximationEnum){
            int vertexlids[NUMVERTICES];
            for(i=0;i<NUMVERTICES;i++) vertexlids[i]=iomodel->my_vertices_lids[penta_vertex_ids[i]-1];
            /*Create VzHO and VzFS Enums*/
            if(iomodel->Data("md.initialization.vz") && iomodel->Data("md.flowequation.borderFS")){
               for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1];
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzFSEnum);
               for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*(1-iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1]);
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzHOEnum);
            }
            else{
               for(i=0;i<6;i++)nodeinputs[i]=0;
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzFSEnum);
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzHOEnum);
            }
         }
         if((IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[index])))==SSAFSApproximationEnum){
            int vertexlids[NUMVERTICES];
            for(i=0;i<NUMVERTICES;i++) vertexlids[i]=iomodel->my_vertices_lids[penta_vertex_ids[i]-1];
            /*Create VzSSA and VzFS Enums*/
            if(iomodel->Data("md.initialization.vz") && iomodel->Data("md.flowequation.borderFS")){
               for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1];
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzFSEnum);
               for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*(1-iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1]);
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzSSAEnum);
            }
            else{
               for(i=0;i<6;i++)nodeinputs[i]=0;
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzFSEnum);
               this->SetElementInput(inputs2,NUMVERTICES,vertexlids,nodeinputs,VzSSAEnum);
            }
         }
         break;
      default:
         /*No update for other solution types*/
         break;
   }
}
/*}}}*/
void       Penta::UpdateConstraintsExtrudeFromBase(void){/*{{{*/
 
   if(!IsOnBase()) return;
 
   int        extrusioninput;
   IssmDouble value,isonbase;
 
   this->parameters->FindParam(&extrusioninput,InputToExtrudeEnum);
   Input2* input = this->GetInput2(extrusioninput);      _assert_(extrusioninput);
   Input2* onbase = this->GetInput2(MeshVertexonbaseEnum); _assert_(onbase);
 
   GaussPenta* gauss=new GaussPenta();
   for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
      gauss->GaussNode(this->element_type,iv);
      onbase->GetInputValue(&isonbase,gauss);
      if(isonbase==1.){
         input->GetInputValue(&value,gauss);
         this->nodes[iv]->ApplyConstraint(0,value);
      }
   }
   delete gauss;
 
}
/*}}}*/
void       Penta::UpdateConstraintsExtrudeFromTop(void){/*{{{*/
 
   if(!IsOnSurface()) return;
 
   int extrusioninput;
   int indices[3]={3,4,5};
   IssmDouble value;
 
   this->parameters->FindParam(&extrusioninput,InputToExtrudeEnum);
   Input2* input = this->GetInput2(extrusioninput); _assert_(extrusioninput);
 
   GaussPenta* gauss=new GaussPenta();
   for(int i=0;i<3;i++){
      gauss->GaussNode(P1Enum,indices[i]);
      input->GetInputValue(&value,gauss);
      this->nodes[indices[i]]->ApplyConstraint(0,value);
   }
   delete gauss;
 
}
/*}}}*/
int        Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
 
   int i;
   int nflipped=0;
 
   /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
   for(i=0;i<NUMVERTICES;i++){
      if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
         vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
 
         /*If node was not on ice shelf, we flipped*/
         if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
            nflipped++;
         }
      }
   }
   return nflipped;
}
/*}}}*/
void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
   PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
}
/*}}}*/
void       Penta::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
   PentaRef::GetInputValue(pvalue,values,gauss,P1Enum);
}
/*}}}*/
int        Penta::VelocityInterpolation(void){/*{{{*/
   return PentaRef::VelocityInterpolation(this->element_type);
}
/*}}}*/
int        Penta::VertexConnectivity(int vertexindex){/*{{{*/
   _assert_(this->vertices);
   return this->vertices[vertexindex]->Connectivity();
}
/*}}}*/
void       Penta::VerticalSegmentIndices(int** pindices,int* pnumseg){/*{{{*/
 
   /*output*/
   int *indices = xNew<int>(3*2);
   indices[0*2 + 0] = 0; indices[0*2 + 1] = 3;
   indices[1*2 + 0] = 1; indices[1*2 + 1] = 4;
   indices[2*2 + 0] = 2; indices[2*2 + 1] = 5;
 
   /*Assign output pointers*/
   *pindices = indices;
   *pnumseg  = 3;
}
/*}}}*/
void       Penta::VerticalSegmentIndicesBase(int** pindices,int* pnumseg){/*{{{*/
 
   PentaRef::VerticalSegmentIndicesBase(pindices,pnumseg,this->GetElementType());
 
   return;
}
/*}}}*/
void       Penta::ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input,Input2* vz_input){/*{{{*/
 
   /*Intermediaries*/
   IssmDouble phi;
   IssmDouble viscosity;
   IssmDouble epsilon[6];
 
   _assert_(gauss->Enum()==GaussPentaEnum);
   this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
   this->material->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
   GetPhi(&phi,&epsilon[0],viscosity);
 
   /*Assign output pointer*/
   *pphi = phi;
}
/*}}}*/
 
#ifdef _HAVE_GIA_
void       Penta::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y){/*{{{*/
   _error_("GIA deflection not implemented yet!");
}
/*}}}*/
#endif
 
#ifdef _HAVE_DAKOTA_
void       Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
 
   /*Check that name is an element input*/
   if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
   TransientInput2* transientinput = inputs2->GetTransientInput(name);
 
   switch(type){
 
      case VertexEnum:
 
         /*Get LID lists once for all*/
         IssmDouble  values[NUMVERTICES];
         int         lidlist[NUMVERTICES];
         this->GetVerticesLidList(&lidlist[0]);
 
         /*Create transient input: */
         for(int t=0;t<ncols;t++){ //ncols is the number of times
 
            /*create input values: */
            for(int i=0;i<6;i++){
               int row=this->vertices[i]->Sid();
               values[i]=matrix[ncols*row+t];
            }
 
            /*time:*/
            IssmDouble time=matrix[(nrows-1)*ncols+t];
 
            transientinput->AddPentaTimeInput(t,NUMVERTICES,&lidlist[0],&values[0],P1Enum);
         }
         break;
 
      case ElementEnum:
         /*Get value for the element: */
         for(int t=0;t<ncols;t++){ //ncols is the number of times
            IssmDouble value=matrix[ncols*(this->Sid())+t];
            IssmDouble time=matrix[(nrows-1)*ncols+t];
            transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum);
         }
         break;
 
      default:
         _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
   }
 
}
/*}}}*/
void       Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
 
   int i,j;
 
   /*Check that name is an element input*/
   if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
 
   switch(type){
 
      case VertexEnum:
 
         /*New PentaInput*/
         IssmDouble values[6];
 
         /*Get values on the 6 vertices*/
         for (i=0;i<6;i++){
            values[i]=vector[this->vertices[i]->Sid()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
         }
 
         /*Branch on the specified type of update: */
         switch(name){
            case ThicknessEnum:
               /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium*/
               IssmDouble  thickness[6];
               IssmDouble  thickness_init[6];
               IssmDouble  hydrostatic_ratio[6];
               IssmDouble  surface[6];
               IssmDouble  bed[6];
 
               /*retrieve inputs: */
               Element::GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
               Element::GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
               Element::GetInputListOnVertices(&bed[0],BaseEnum);
               Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
 
               /*build new thickness: */
//             for(j=0;j<6;j++)thickness[j]=values[j];
 
               /*build new bed and surface: */
               if (this->IsFloating()){
                  /*hydrostatic equilibrium: */
                  IssmDouble rho_ice,rho_water,di;
                  rho_ice=this->FindParam(MaterialsRhoIceEnum);
                  rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
 
                  di=rho_ice/rho_water;
 
                  /*build new thickness: */
                  for (j=0; j<6; j++) {
                  /*  for observed/interpolated/hydrostatic thickness, remove scaling from any hydrostatic thickness  */
                     if     (hydrostatic_ratio[j] >= 0.)
                        thickness[j]=values[j]-(values[j]/thickness_init[j]-1.)*hydrostatic_ratio[j]*surface[j]/(1.-di);
                  /*  for minimum thickness, don't scale  */
                     else
                        thickness[j]=thickness_init[j];
 
                  /*  check the computed thickness and update bed  */
                     if (thickness[j] < 0.)
                        thickness[j]=1./(1.-di);
                     bed[j]=surface[j]-thickness[j];
                  }
 
//                for(j=0;j<6;j++){
//                   surface[j]=(1-di)*thickness[j];
//                   bed[j]=-di*thickness[j];
//                }
               }
               else{
                  /*build new thickness: */
                  for (j=0; j<6; j++) {
                  /*  for observed thickness, use scaled value  */
                     if(hydrostatic_ratio[j] >= 0.)
                        thickness[j]=values[j];
                  /*  for minimum thickness, don't scale  */
                     else
                        thickness[j]=thickness_init[j];
                  }
 
                  /*update bed on grounded ice: */
//                for(j=0;j<6;j++)surface[j]=bed[j]+thickness[j];
                  for(j=0;j<6;j++)bed[j]=surface[j]-thickness[j];
               }
 
               /*Add new inputs: */
               this->AddInput2(ThicknessEnum,thickness,P1Enum);
               this->AddInput2(BaseEnum,bed,P1Enum);
               this->AddInput2(SurfaceEnum,surface,P1Enum);
               break;
 
            default:
               this->AddInput2(name,values,P1Enum);
         }
         break;
 
      case ElementEnum:
         IssmDouble value;
         /*Get value for the element: */
         value=vector[this->Sid()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
         this->AddInput2(name,&value,P0Enum);
         break;
 
      default:
         _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
   }
 
}
/*}}}*/
#endif