Changeset 4250

Timestamp:

06/26/10 23:10:55 (15 years ago)

Author:

Eric.Larour

Message:

reorganized element routines

Location:

issm/trunk/src/c/objects/Elements

Files:

• Beam.cpp (modified) (12 diffs)
• Beam.h (modified) (1 diff)
• Penta.cpp (modified) (19 diffs)
• Penta.h (modified) (3 diffs)

issm/trunk/src/c/objects/Elements/Beam.cpp ¶

@@ -31 +31 @@
 }
 /*}}}*/
-/*FUNCTION void Beam::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/
-void Beam::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-
-/*Object virtual functions definitions: */
+
+/*Object virtual functions definitions:*/
+/*FUNCTION Beam::copy{{{1*/
+Object* Beam::copy() {
+
+        int i;
+        Beam* beam=NULL;
+
+        beam=new Beam();
+
+        /*copy fields: */
+        beam->id=this->id;
+        if(this->inputs){
+                beam->inputs=(Inputs*)this->inputs->Copy();
+        }
+        else{
+                beam->inputs=new Inputs();
+        }
+        /*point parameters: */
+        beam->parameters=this->parameters;
+
+        /*pointers: */
+        beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); 
+        for(i=0;i<2;i++)beam->nodes[i]=this->nodes[i];
+        beam->matice=this->matice;
+        beam->matpar=this->matpar;
+
+        return beam;
+}
+/*}}}*/
 /*FUNCTION Beam::DeepEcho{{{1*/
 void Beam::DeepEcho(void){
@@ -58 +81 @@
 }
 /*}}}*/
-/*FUNCTION Beam::Id{{{1*/
-int    Beam::Id(void){ return id; }
-/*}}}*/
-/*FUNCTION Beam::MyRank{{{1*/
-int    Beam::MyRank(void){ 
-        extern int my_rank;
-        return my_rank; 
-}
-/*}}}*/
-/*FUNCTION Beam::Marshall{{{1*/
-void  Beam::Marshall(char** pmarshalled_dataset){
-        ISSMERROR("not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::MarshallSize{{{1*/
-int   Beam::MarshallSize(){
-        ISSMERROR("not supported yet!");
-}
-/*}}}*/
 /*FUNCTION Beam::Demarshall{{{1*/
 void  Beam::Demarshall(char** pmarshalled_dataset){
         ISSMERROR("not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::Enum{{{1*/
-int Beam::Enum(void){
-
-        return BeamEnum;
-
-}
-/*}}}*/
-/*FUNCTION Beam::copy{{{1*/
-Object* Beam::copy() {
-
-        int i;
-        Beam* beam=NULL;
-
-        beam=new Beam();
-
-        /*copy fields: */
-        beam->id=this->id;
-        if(this->inputs){
-                beam->inputs=(Inputs*)this->inputs->Copy();
-        }
-        else{
-                beam->inputs=new Inputs();
-        }
-        /*point parameters: */
-        beam->parameters=this->parameters;
-
-        /*pointers: */
-        beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); 
-        for(i=0;i<2;i++)beam->nodes[i]=this->nodes[i];
-        beam->matice=this->matice;
-        beam->matpar=this->matpar;
-
-        return beam;
-}
-/*}}}*/
-
-/*Beam management*/
-/*FUNCTION Beam::Configure {{{1*/
-void  Beam::Configure(Elements* elementsin,Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
-
-        ISSMERROR(" not supported yet!");
-
 }
 /*}}}*/
@@ -146 +106 @@
 }
 /*}}}*/
-/*FUNCTION Beam::IsInput{{{1*/
-bool Beam::IsInput(int name){
-        if (name==SurfaceSlopeXEnum ||
-                                name==SurfaceSlopeYEnum){
-                return true;
-        }
-        else return false;
+/*FUNCTION Beam::Enum{{{1*/
+int Beam::Enum(void){
+
+        return BeamEnum;
+
+}
+/*}}}*/
+/*FUNCTION Beam::Id{{{1*/
+int    Beam::Id(void){ return id; }
+/*}}}*/
+/*FUNCTION Beam::Marshall{{{1*/
+void  Beam::Marshall(char** pmarshalled_dataset){
+        ISSMERROR("not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::MarshallSize{{{1*/
+int   Beam::MarshallSize(){
+        ISSMERROR("not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::MyRank{{{1*/
+int    Beam::MyRank(void){ 
+        extern int my_rank;
+        return my_rank; 
+}
+/*}}}*/
+
+/*Update virtual functions definitions:*/
+/*FUNCTION Beam::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
+void  Beam::InputUpdateFromVector(double* vector, int name, int type){
+
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+        switch(type){
+
+                case VertexEnum:
+
+                        /*New PentaVertexInpu*/
+                        double values[2];
+
+                        /*Get values on the 6 vertices*/
+                        for (int i=0;i<2;i++){
+                                values[i]=vector[nodes[i]->GetVertexDof()];
+                        }
+
+                        /*update input*/
+                        this->inputs->AddInput(new BeamVertexInput(name,values));
+                        return;
+
+                default:
+
+                        ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
+        }
+}
+/*}}}*/
+/*FUNCTION Beam::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
+void  Beam::InputUpdateFromVector(int* vector, int name, int type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
+void  Beam::InputUpdateFromVector(bool* vector, int name, int type){
+        ISSMERROR(" not supported yet!");
 }
 /*}}}*/
@@ -160 +176 @@
 }
 /*}}}*/
-/*FUNCTION Beam::GetSolutionFromInputs(Vec solution);{{{1*/
-void  Beam::GetSolutionFromInputs(Vec solution){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::InputToResult(int enum_type,int step,double time){{{1*/
-void  Beam::InputToResult(int enum_type,int step,double time){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::ProcessResultsUnits(void){{{1*/
-void  Beam::ProcessResultsUnits(void){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-
-/*Object functions*/
+
+/*Element virtual functions definitions:*/
 /*FUNCTION Beam::ComputeBasalStress{{{1*/
 void  Beam::ComputeBasalStress(Vec eps){
@@ -232 +233 @@
 }
 /*}}}*/
+/*FUNCTION Beam::Configure {{{1*/
+void  Beam::Configure(Elements* elementsin,Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
+
+        ISSMERROR(" not supported yet!");
+
+}
+/*}}}*/
 /*FUNCTION Beam::CostFunction{{{1*/
 double Beam::CostFunction(void){
@@ -255 +263 @@
 }
 /*}}}*/
+/*FUNCTION Beam::CreatePVector{{{1*/
+void  Beam::CreatePVector(Vec pg){
+
+        int analysis_type,sub_analysis_type;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
+        
+        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
+        if (analysis_type==DiagnosticHutterAnalysisEnum) {
+                CreatePVectorDiagnosticHutter( pg);
+        }
+        else{
+                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
+        }
+
+}
+/*}}}*/
+/*FUNCTION Beam::Du{{{1*/
+void  Beam::Du(Vec){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GetBedList{{{1*/
+void  Beam::GetBedList(double*){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GetMatPar{{{1*/
+void* Beam::GetMatPar(){
+
+        return matpar;
+}
+/*}}}*/
+/*FUNCTION Beam::GetNodes{{{1*/
+void  Beam::GetNodes(void** vpnodes){
+        int i;
+        Node** pnodes=NULL;
+
+        /*recover nodes: */
+        pnodes=(Node**)vpnodes;
+
+        for(i=0;i<3;i++){
+                pnodes[i]=nodes[i];
+        }
+}
+/*}}}*/
+/*FUNCTION Beam::GetOnBed{{{1*/
+bool   Beam::GetOnBed(){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GetShelf{{{1*/
+bool   Beam::GetShelf(){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GetSolutionFromInputs(Vec solution);{{{1*/
+void  Beam::GetSolutionFromInputs(Vec solution){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GetThicknessList{{{1*/
+void  Beam::GetThicknessList(double* thickness_list){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
+void  Beam::GetVectorFromInputs(Vec vector,int NameEnum){
+
+        int i;
+        const int numvertices=2;
+        int doflist1[numvertices];
+
+        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
+        for(i=0;i<this->inputs->Size();i++){
+                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+                if(input->EnumType()==NameEnum){
+                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+                        this->GetDofList1(&doflist1[0]);
+                        input->GetVectorFromInputs(vector,&doflist1[0]);
+                        break;
+                }
+        }
+}
+/*}}}*/
+/*FUNCTION Beam::Gradj{{{1*/
+void  Beam::Gradj(Vec gradient,int control_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GradjB{{{1*/
+void  Beam::GradjB(Vec gradient){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::GradjDrag{{{1*/
+void  Beam::GradjDrag(Vec gradient){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/
+void  Beam::InputAXPY(int YEnum, double scalar, int XEnum){
+
+        Input* xinput=NULL;
+        Input* yinput=NULL;
+
+        /*Find x and y inputs: */
+        xinput=(Input*)this->inputs->GetInput(XEnum);
+        yinput=(Input*)this->inputs->GetInput(YEnum);
+
+        /*some checks: */
+        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
+        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
+
+        /*Scale: */
+        yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
+/*FUNCTION Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
+void  Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){
+
+        Input* input=NULL;
+
+        /*Find input: */
+        input=(Input*)this->inputs->GetInput(control_type);
+        
+        /*Do nothing if we  don't find it: */
+        if(!input)return;
+
+        /*Constrain input using cm_min and cm_max: */
+        input->Constrain(cm_min,cm_max);
+
+}
+/*}}}*/
+/*FUNCTION Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
+void  Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
+
+        int     i;
+        Input** new_inputs=NULL;
+        Input** old_inputs=NULL;
+        int     converged=1;
+
+        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
+        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
+        
+        for(i=0;i<num_enums/2;i++){
+                new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
+                old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
+                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
+                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
+        }
+
+        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
+        for(i=0;i<num_criterionenums;i++){
+                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
+                if(eps[i]>criterionvalues[i]) converged=0; 
+        }
+
+        /*Assign output pointers:*/
+        *pconverged=converged;
+
+}
+/*}}}*/
+/*FUNCTION Beam::InputDuplicate(int original_enum,int new_enum){{{1*/
+void  Beam::InputDuplicate(int original_enum,int new_enum){
+
+        Input* original=NULL;
+        Input* copy=NULL;
+
+        /*Make a copy of the original input: */
+        original=(Input*)this->inputs->GetInput(original_enum);
+        copy=(Input*)original->copy();
+
+        /*Change copy enum to reinitialized_enum: */
+        copy->ChangeEnum(new_enum);
+
+        /*Add copy into inputs, it will wipe off the one already there: */
+        inputs->AddObject((Input*)copy);
+}
+/*}}}*/
+/*FUNCTION Beam::InputScale(int enum_type,double scale_factor){{{1*/
+void  Beam::InputScale(int enum_type,double scale_factor){
+
+        Input* input=NULL;
+
+        /*Make a copy of the original input: */
+        input=(Input*)this->inputs->GetInput(enum_type);
+
+        /*Scale: */
+        input->Scale(scale_factor);
+}
+/*}}}*/
+/*FUNCTION Beam::InputToResult(int enum_type,int step,double time){{{1*/
+void  Beam::InputToResult(int enum_type,int step,double time){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::MassFlux{{{1*/
+double Beam::MassFlux( double* segment){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/
+void  Beam::MaxAbsVx(double* pmaxabsvx, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vx_values[numgrids];
+        double  maxabsvx;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+
+        /*now, compute maximum:*/
+        maxabsvx=fabs(vx_values[0]);
+        for(i=1;i<numgrids;i++){
+                if (fabs(vx_values[i])>maxabsvx)maxabsvx=fabs(vx_values[i]);
+        }
+
+        /*Assign output pointers:*/
+        *pmaxabsvx=maxabsvx;
+}
+/*}}}*/
+/*FUNCTION Beam::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/
+void  Beam::MaxAbsVy(double* pmaxabsvy, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vy_values[numgrids];
+        double  maxabsvy;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+
+        /*now, compute maximum:*/
+        maxabsvy=fabs(vy_values[0]);
+        for(i=1;i<numgrids;i++){
+                if (fabs(vy_values[i])>maxabsvy)maxabsvy=fabs(vy_values[i]);
+        }
+
+        /*Assign output pointers:*/
+        *pmaxabsvy=maxabsvy;
+}
+/*}}}*/
+/*FUNCTION Beam::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/
+void  Beam::MaxAbsVz(double* pmaxabsvz, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vz_values[numgrids];
+        double  maxabsvz;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute maximum:*/
+        maxabsvz=fabs(vz_values[0]);
+        for(i=1;i<numgrids;i++){
+                if (fabs(vz_values[i])>maxabsvz)maxabsvz=fabs(vz_values[i]);
+        }
+
+        /*Assign output pointers:*/
+        *pmaxabsvz=maxabsvz;
+}
+/*}}}*/
+/*FUNCTION Beam::MaxVel(double* pmaxvel, bool process_units);{{{1*/
+void  Beam::MaxVel(double* pmaxvel, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vx_values[numgrids];
+        double  vy_values[numgrids];
+        double  vz_values[numgrids];
+        double  vel_values[numgrids];
+        double  maxvel;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute maximum of velocity :*/
+        if(dim==2){
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
+        }
+        else{
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
+        }
+
+        /*now, compute maximum:*/
+        maxvel=vel_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vel_values[i]>maxvel)maxvel=vel_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvel=maxvel;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MaxVx(double* pmaxvx, bool process_units);{{{1*/
+void  Beam::MaxVx(double* pmaxvx, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vx_values[numgrids];
+        double  maxvx;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+
+        /*now, compute maximum:*/
+        maxvx=vx_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vx_values[i]>maxvx)maxvx=vx_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvx=maxvx;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MaxVy(double* pmaxvy, bool process_units);{{{1*/
+void  Beam::MaxVy(double* pmaxvy, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vy_values[numgrids];
+        double  maxvy;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+
+        /*now, compute maximum:*/
+        maxvy=vy_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vy_values[i]>maxvy)maxvy=vy_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvy=maxvy;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MaxVz(double* pmaxvz, bool process_units);{{{1*/
+void  Beam::MaxVz(double* pmaxvz, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vz_values[numgrids];
+        double  maxvz;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute maximum:*/
+        maxvz=vz_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vz_values[i]>maxvz)maxvz=vz_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvz=maxvz;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MinVel(double* pminvel, bool process_units);{{{1*/
+void  Beam::MinVel(double* pminvel, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vx_values[numgrids];
+        double  vy_values[numgrids];
+        double  vz_values[numgrids];
+        double  vel_values[numgrids];
+        double  minvel;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute minimum of velocity :*/
+        if(dim==2){
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
+        }
+        else{
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
+        }
+
+        /*now, compute minimum:*/
+        minvel=vel_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vel_values[i]<minvel)minvel=vel_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvel=minvel;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MinVx(double* pminvx, bool process_units);{{{1*/
+void  Beam::MinVx(double* pminvx, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vx_values[numgrids];
+        double  minvx;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+
+        /*now, compute minimum:*/
+        minvx=vx_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vx_values[i]<minvx)minvx=vx_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvx=minvx;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MinVy(double* pminvy, bool process_units);{{{1*/
+void  Beam::MinVy(double* pminvy, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vy_values[numgrids];
+        double  minvy;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+
+        /*now, compute minimum:*/
+        minvy=vy_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vy_values[i]<minvy)minvy=vy_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvy=minvy;
+
+}
+/*}}}*/
+/*FUNCTION Beam::MinVz(double* pminvz, bool process_units);{{{1*/
+void  Beam::MinVz(double* pminvz, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=2;
+        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
+        double  vz_values[numgrids];
+        double  minvz;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute minimum:*/
+        minvz=vz_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vz_values[i]<minvz)minvz=vz_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvz=minvz;
+
+}
+/*}}}*/
+/*FUNCTION Beam::Misfit{{{1*/
+double Beam::Misfit(void){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::PatchFill(int* pcount, Patch* patch){{{1*/
+void  Beam::PatchFill(int* pcount, Patch* patch){
+
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/
+void  Beam::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
+
+        ISSMERROR(" not supported yet!");
+        
+}
+/*}}}*/
+/*FUNCTION Beam::ProcessResultsUnits(void){{{1*/
+void  Beam::ProcessResultsUnits(void){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::SurfaceArea{{{1*/
+double Beam::SurfaceArea(void){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Beam::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/
+void Beam::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+
+/*Beam specific routines: */
 /*FUNCTION Beam::CreateKMatrixDiagnosticHutter{{{1*/
 
@@ -304 +871 @@
         /*Add Ke_gg to global matrix Kgg: */
         MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-}
-/*}}}*/
-/*FUNCTION Beam::CreatePVector{{{1*/
-void  Beam::CreatePVector(Vec pg){
-
-        int analysis_type,sub_analysis_type;
-
-        /*retrive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
-        
-        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
-        if (analysis_type==DiagnosticHutterAnalysisEnum) {
-                CreatePVectorDiagnosticHutter( pg);
-        }
-        else{
-                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
-        }
 
 }
@@ -455 +1003 @@
 }
 /*}}}*/
-/*FUNCTION Beam::Du{{{1*/
-void  Beam::Du(Vec){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::GetBedList{{{1*/
-void  Beam::GetBedList(double*){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
 /*FUNCTION Beam::GetDofList{{{1*/
 void  Beam::GetDofList(int* doflist,int* pnumberofdofspernode){
@@ -509 +1047 @@
 }
 /*}}}*/
-/*FUNCTION Beam::GetMatPar{{{1*/
-void* Beam::GetMatPar(){
-
-        return matpar;
-}
-/*}}}*/
 /*FUNCTION Beam::GetNodalFunctions{{{1*/
 void Beam::GetNodalFunctions(double* l1l2, double gauss_coord){
@@ -527 +1059 @@
 }
 /*}}}*/
-/*FUNCTION Beam::GetNodes{{{1*/
-void  Beam::GetNodes(void** vpnodes){
-        int i;
-        Node** pnodes=NULL;
-
-        /*recover nodes: */
-        pnodes=(Node**)vpnodes;
-
-        for(i=0;i<3;i++){
-                pnodes[i]=nodes[i];
-        }
-}
-/*}}}*/
-/*FUNCTION Beam::GetOnBed{{{1*/
-bool   Beam::GetOnBed(){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
 /*FUNCTION Beam::GetParameterValue{{{1*/
 void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){
@@ -555 +1069 @@
 }
 /*}}}*/
-/*FUNCTION Beam::GetShelf{{{1*/
-bool   Beam::GetShelf(){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::GetThicknessList{{{1*/
-void  Beam::GetThicknessList(double* thickness_list){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::Gradj{{{1*/
-void  Beam::Gradj(Vec gradient,int control_type){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::GradjB{{{1*/
-void  Beam::GradjB(Vec gradient){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::GradjDrag{{{1*/
-void  Beam::GradjDrag(Vec gradient){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::MassFlux{{{1*/
-double Beam::MassFlux( double* segment){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::Misfit{{{1*/
-double Beam::Misfit(void){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::SurfaceArea{{{1*/
-double Beam::SurfaceArea(void){
-        ISSMERROR(" not supported yet!");
+/*FUNCTION Beam::IsInput{{{1*/
+bool Beam::IsInput(int name){
+        if (name==SurfaceSlopeXEnum ||
+                                name==SurfaceSlopeYEnum){
+                return true;
+        }
+        else return false;
 }
 /*}}}*/
@@ -601 +1084 @@
 }
 /*}}}1*/
-/*FUNCTION Beam::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
-void  Beam::InputUpdateFromVector(double* vector, int name, int type){
-
-        /*Check that name is an element input*/
-        if (!IsInput(name)) return;
-        switch(type){
-
-                case VertexEnum:
-
-                        /*New PentaVertexInpu*/
-                        double values[2];
-
-                        /*Get values on the 6 vertices*/
-                        for (int i=0;i<2;i++){
-                                values[i]=vector[nodes[i]->GetVertexDof()];
-                        }
-
-                        /*update input*/
-                        this->inputs->AddInput(new BeamVertexInput(name,values));
-                        return;
-
-                default:
-
-                        ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
-        }
-}
-/*}}}*/
-/*FUNCTION Beam::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
-void  Beam::InputUpdateFromVector(int* vector, int name, int type){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
-void  Beam::InputUpdateFromVector(bool* vector, int name, int type){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/
-void  Beam::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
-
-        ISSMERROR(" not supported yet!");
-        
-}
-/*}}}*/
-/*FUNCTION Beam::PatchFill(int* pcount, Patch* patch){{{1*/
-void  Beam::PatchFill(int* pcount, Patch* patch){
-
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::MinVel(double* pminvel, bool process_units);{{{1*/
-void  Beam::MinVel(double* pminvel, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vx_values[numgrids];
-        double  vy_values[numgrids];
-        double  vz_values[numgrids];
-        double  vel_values[numgrids];
-        double  minvel;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute minimum of velocity :*/
-        if(dim==2){
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
-        }
-        else{
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
-        }
-
-        /*now, compute minimum:*/
-        minvel=vel_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vel_values[i]<minvel)minvel=vel_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvel=minvel;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxVel(double* pmaxvel, bool process_units);{{{1*/
-void  Beam::MaxVel(double* pmaxvel, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vx_values[numgrids];
-        double  vy_values[numgrids];
-        double  vz_values[numgrids];
-        double  vel_values[numgrids];
-        double  maxvel;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute maximum of velocity :*/
-        if(dim==2){
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
-        }
-        else{
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
-        }
-
-        /*now, compute maximum:*/
-        maxvel=vel_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vel_values[i]>maxvel)maxvel=vel_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvel=maxvel;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MinVx(double* pminvx, bool process_units);{{{1*/
-void  Beam::MinVx(double* pminvx, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vx_values[numgrids];
-        double  minvx;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-
-        /*now, compute minimum:*/
-        minvx=vx_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vx_values[i]<minvx)minvx=vx_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvx=minvx;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxVx(double* pmaxvx, bool process_units);{{{1*/
-void  Beam::MaxVx(double* pmaxvx, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vx_values[numgrids];
-        double  maxvx;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-
-        /*now, compute maximum:*/
-        maxvx=vx_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vx_values[i]>maxvx)maxvx=vx_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvx=maxvx;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/
-void  Beam::MaxAbsVx(double* pmaxabsvx, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vx_values[numgrids];
-        double  maxabsvx;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-
-        /*now, compute maximum:*/
-        maxabsvx=fabs(vx_values[0]);
-        for(i=1;i<numgrids;i++){
-                if (fabs(vx_values[i])>maxabsvx)maxabsvx=fabs(vx_values[i]);
-        }
-
-        /*Assign output pointers:*/
-        *pmaxabsvx=maxabsvx;
-}
-/*}}}*/
-/*FUNCTION Beam::MinVy(double* pminvy, bool process_units);{{{1*/
-void  Beam::MinVy(double* pminvy, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vy_values[numgrids];
-        double  minvy;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-
-        /*now, compute minimum:*/
-        minvy=vy_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vy_values[i]<minvy)minvy=vy_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvy=minvy;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxVy(double* pmaxvy, bool process_units);{{{1*/
-void  Beam::MaxVy(double* pmaxvy, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vy_values[numgrids];
-        double  maxvy;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-
-        /*now, compute maximum:*/
-        maxvy=vy_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vy_values[i]>maxvy)maxvy=vy_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvy=maxvy;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/
-void  Beam::MaxAbsVy(double* pmaxabsvy, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vy_values[numgrids];
-        double  maxabsvy;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-
-        /*now, compute maximum:*/
-        maxabsvy=fabs(vy_values[0]);
-        for(i=1;i<numgrids;i++){
-                if (fabs(vy_values[i])>maxabsvy)maxabsvy=fabs(vy_values[i]);
-        }
-
-        /*Assign output pointers:*/
-        *pmaxabsvy=maxabsvy;
-}
-/*}}}*/
-/*FUNCTION Beam::MinVz(double* pminvz, bool process_units);{{{1*/
-void  Beam::MinVz(double* pminvz, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vz_values[numgrids];
-        double  minvz;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute minimum:*/
-        minvz=vz_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vz_values[i]<minvz)minvz=vz_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvz=minvz;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxVz(double* pmaxvz, bool process_units);{{{1*/
-void  Beam::MaxVz(double* pmaxvz, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vz_values[numgrids];
-        double  maxvz;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute maximum:*/
-        maxvz=vz_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vz_values[i]>maxvz)maxvz=vz_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvz=maxvz;
-
-}
-/*}}}*/
-/*FUNCTION Beam::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/
-void  Beam::MaxAbsVz(double* pmaxabsvz, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=2;
-        double  gaussgrids[numgrids][2]={{0,1},{1,0}};
-        double  vz_values[numgrids];
-        double  maxabsvz;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute maximum:*/
-        maxabsvz=fabs(vz_values[0]);
-        for(i=1;i<numgrids;i++){
-                if (fabs(vz_values[i])>maxabsvz)maxabsvz=fabs(vz_values[i]);
-        }
-
-        /*Assign output pointers:*/
-        *pmaxabsvz=maxabsvz;
-}
-/*}}}*/
-/*FUNCTION Beam::InputDuplicate(int original_enum,int new_enum){{{1*/
-void  Beam::InputDuplicate(int original_enum,int new_enum){
-
-        Input* original=NULL;
-        Input* copy=NULL;
-
-        /*Make a copy of the original input: */
-        original=(Input*)this->inputs->GetInput(original_enum);
-        copy=(Input*)original->copy();
-
-        /*Change copy enum to reinitialized_enum: */
-        copy->ChangeEnum(new_enum);
-
-        /*Add copy into inputs, it will wipe off the one already there: */
-        inputs->AddObject((Input*)copy);
-}
-/*}}}*/
-/*FUNCTION Beam::InputScale(int enum_type,double scale_factor){{{1*/
-void  Beam::InputScale(int enum_type,double scale_factor){
-
-        Input* input=NULL;
-
-        /*Make a copy of the original input: */
-        input=(Input*)this->inputs->GetInput(enum_type);
-
-        /*Scale: */
-        input->Scale(scale_factor);
-}
-/*}}}*/
-/*FUNCTION Beam::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/
-void  Beam::InputAXPY(int YEnum, double scalar, int XEnum){
-
-        Input* xinput=NULL;
-        Input* yinput=NULL;
-
-        /*Find x and y inputs: */
-        xinput=(Input*)this->inputs->GetInput(XEnum);
-        yinput=(Input*)this->inputs->GetInput(YEnum);
-
-        /*some checks: */
-        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
-        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
-
-        /*Scale: */
-        yinput->AXPY(xinput,scalar);
-}
-/*}}}*/
-/*FUNCTION Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
-void  Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){
-
-        Input* input=NULL;
-
-        /*Find input: */
-        input=(Input*)this->inputs->GetInput(control_type);
-        
-        /*Do nothing if we  don't find it: */
-        if(!input)return;
-
-        /*Constrain input using cm_min and cm_max: */
-        input->Constrain(cm_min,cm_max);
-
-}
-/*}}}*/
-/*FUNCTION Beam::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
-void  Beam::GetVectorFromInputs(Vec vector,int NameEnum){
-
-        int i;
-        const int numvertices=2;
-        int doflist1[numvertices];
-
-        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
-        for(i=0;i<this->inputs->Size();i++){
-                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
-                if(input->EnumType()==NameEnum){
-                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
-                        this->GetDofList1(&doflist1[0]);
-                        input->GetVectorFromInputs(vector,&doflist1[0]);
-                        break;
-                }
-        }
-}
-/*}}}*/
-/*FUNCTION Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
-void  Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
-
-        int     i;
-        Input** new_inputs=NULL;
-        Input** old_inputs=NULL;
-        int     converged=1;
-
-        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
-        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
-        
-        for(i=0;i<num_enums/2;i++){
-                new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
-                old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
-                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
-                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
-        }
-
-        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
-        for(i=0;i<num_criterionenums;i++){
-                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
-                if(eps[i]>criterionvalues[i]) converged=0; 
-        }
-
-        /*Assign output pointers:*/
-        *pconverged=converged;
-
-}
-/*}}}*/

issm/trunk/src/c/objects/Elements/Beam.h ¶

@@ -42 +42 @@
                 /*}}}*/
                 /*Object virtual functions definitions: {{{1*/
-                void  Echo();
-                void  DeepEcho();
-                int   Id(); 
-                int   MyRank();
-                void  Marshall(char** pmarshalled_dataset);
-                int   MarshallSize();
-                void  Demarshall(char** pmarshalled_dataset);
-                int   Enum();
                 Object* copy();
+                void    DeepEcho();
+                void    Demarshall(char** pmarshalled_dataset);
+                void    Echo();
+                int     Enum();
+                int     Id(); 
+                void    Marshall(char** pmarshalled_dataset);
+                int     MarshallSize();
+                int     MyRank();
                 /*}}}*/
                 /*Update virtual functions resolution: {{{1*/

issm/trunk/src/c/objects/Elements/Penta.cpp ¶

@@ -88 +88 @@
 }
 /*}}}*/
-/*FUNCTION Penta::Update(IoModel* iomodel,int analysis_counter,int analysis_type) {{{1*/
-void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type){ 
-
-        /*Intermediaries*/
-        IssmInt i;
-        int penta_node_ids[6];
-        int penta_vertex_ids[6];
-        double nodeinputs[6];
-
-        /*Checks if debuging*/
-        /*{{{2*/
-        ISSMASSERT(iomodel->elements);
-        /*}}}*/
-
-        /*Recover vertices ids needed to initialize inputs*/
-        for(i=0;i<6;i++){ 
-                penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
-        }
-
-        /*Recover nodes ids needed to initialize the node hook.*/
-        for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
-                penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
-        }
-
-        /*hooks: */
-        this->SetHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
-
-        //add as many inputs per element as requested: 
-        if (iomodel->thickness) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));
-        }
-        if (iomodel->surface) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));
-        }
-        if (iomodel->bed) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs));
-        }
-        if (iomodel->drag_coefficient) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs));
-
-                if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
-                if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
-                this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
-
-        }
-        if (iomodel->melting_rate) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));
-        }
-        if (iomodel->accumulation_rate) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));
-        }
-        if (iomodel->geothermalflux) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs));
-        }       
-        if (iomodel->pressure) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));
-        }
-        if (iomodel->temperature) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs));
-        }
-        if (iomodel->dhdt) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs));
-        }
-        /*vx,vy and vz: */
-        if (iomodel->vx) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
-                this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));
-        }
-        if (iomodel->vy) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
-                this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));
-        }
-        if (iomodel->vz) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
-                this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
-        }
-        if (iomodel->vx_obs) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs));
-        }
-        if (iomodel->vy_obs) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs));
-        }
-        if (iomodel->vz_obs) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
-                this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs));
-        }
-        if (iomodel->weights) {
-                for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];
-                this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs));
-        }
-
-        /*default vx,vy and vz: either observation or 0 */
-        if (!iomodel->vx){
-                if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
-                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
-                this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
-                this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));
-        }
-        if (!iomodel->vy){
-                if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
-                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
-                this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
-                this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));
-        }
-        if (!iomodel->vz){
-                if (iomodel->vz_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
-                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
-                this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
-                this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
-        }
-        if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
-        if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
-        if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
-        if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
-
-        //elements of type 3 are macayeal type penta. we collapse the formulation on their base.
-        if(iomodel->elements_type){
-                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ 
-                        collapse[analysis_counter]=true;
-                }
-                else{
-                        collapse[analysis_counter]=false;
-                }
-        }
-}
-/*}}}*/
 
 /*Object virtual functions definitions: */
-/*FUNCTION Penta::Echo{{{1*/
-
-void Penta::Echo(void){
-        this->DeepEcho();
+/*FUNCTION Penta::copy {{{1*/
+Object* Penta::copy() {
+
+        int i;
+
+        Penta* penta=NULL;
+
+        penta=new Penta();
+
+        /*copy fields: */
+        penta->id=this->id;
+        if(this->inputs){
+                penta->inputs=(Inputs*)this->inputs->Copy();
+        }
+        else{
+                penta->inputs=new Inputs();
+        }
+        if(this->results){
+                penta->results=(Results*)this->results->Copy();
+        }
+        else{
+                penta->results=new Results();
+        }
+        /*point parameters: */
+        penta->parameters=this->parameters;
+
+        /*now deal with hooks and objects: */
+        penta->InitHookNodes(this->numanalyses);
+        for(i=0;i<this->numanalyses;i++)penta->hnodes[i].copy(&this->hnodes[i]);
+        penta->hmatice.copy(&this->hmatice);
+        penta->hmatpar.copy(&this->hmatpar);
+        penta->hneighbors.copy(&this->hneighbors);
+
+        /*recover objects: */
+        penta->nodes=(Node**)xmalloc(6*sizeof(Node*)); //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<6;i++)penta->nodes[i]=this->nodes[i];
+        penta->matice=(Matice*)penta->hmatice.delivers();
+        penta->matpar=(Matpar*)penta->hmatpar.delivers();
+        penta->neighbors=(Penta**)penta->hneighbors.deliverp();
+
+        return penta;
 }
 /*}}}*/
@@ -264 +160 @@
 }
 /*}}}*/
+/*FUNCTION Penta::Demarshall {{{1*/
+void  Penta::Demarshall(char** pmarshalled_dataset){
+
+        char* marshalled_dataset=NULL;
+        int   i;
+
+        /*recover marshalled_dataset: */
+        marshalled_dataset=*pmarshalled_dataset;
+
+        /*this time, no need to get enum type, the pointer directly points to the beginning of the 
+         *object data (thanks to DataSet::Demarshall):*/
+
+        memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
+        memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
+
+        /*allocate dynamic memory: */
+        collapse=(bool*)xmalloc(numanalyses*sizeof(bool));
+        InitHookNodes(numanalyses);
+
+        /*demarshall hooks: */
+        for(i=0;i<numanalyses;i++)hnodes[i].Demarshall(&marshalled_dataset);
+        hmatice.Demarshall(&marshalled_dataset);
+        hmatpar.Demarshall(&marshalled_dataset);
+        hneighbors.Demarshall(&marshalled_dataset);
+
+        /*pointers are garbage, until configuration is carried out: */
+        nodes=NULL;
+        matice=NULL;
+        matpar=NULL;
+        neighbors=NULL;
+        
+        /*demarshall inputs and results: */
+        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
+        results=(Results*)DataSetDemarshallRaw(&marshalled_dataset); 
+
+        /*demarshall internal parameters: */
+        memcpy(collapse,marshalled_dataset,numanalyses*sizeof(bool));marshalled_dataset+=numanalyses*sizeof(bool);
+
+        /*parameters: may not exist even yet, so let Configure handle it: */
+        this->parameters=NULL;
+
+        /*return: */
+        *pmarshalled_dataset=marshalled_dataset;
+        return;
+}
+/*}}}*/
+/*FUNCTION Penta::Echo{{{1*/
+
+void Penta::Echo(void){
+        this->DeepEcho();
+}
+/*}}}*/
+/*FUNCTION Penta::Enum {{{1*/
+int Penta::Enum(void){
+
+        return PentaEnum;
+
+}
+/*}}}*/
 /*FUNCTION Penta::Id {{{1*/
 int    Penta::Id(void){
         return id; 
-}
-/*}}}*/
-/*FUNCTION Penta::MyRank {{{1*/
-int    Penta::MyRank(void){ 
-        extern int my_rank;
-        return my_rank; 
 }
 /*}}}*/
@@ -347 +296 @@
 }
 /*}}}*/
-/*FUNCTION Penta::Demarshall {{{1*/
-void  Penta::Demarshall(char** pmarshalled_dataset){
-
-        char* marshalled_dataset=NULL;
-        int   i;
-
-        /*recover marshalled_dataset: */
-        marshalled_dataset=*pmarshalled_dataset;
-
-        /*this time, no need to get enum type, the pointer directly points to the beginning of the 
-         *object data (thanks to DataSet::Demarshall):*/
-
-        memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
-        memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);
-
-        /*allocate dynamic memory: */
-        collapse=(bool*)xmalloc(numanalyses*sizeof(bool));
-        InitHookNodes(numanalyses);
-
-        /*demarshall hooks: */
-        for(i=0;i<numanalyses;i++)hnodes[i].Demarshall(&marshalled_dataset);
-        hmatice.Demarshall(&marshalled_dataset);
-        hmatpar.Demarshall(&marshalled_dataset);
-        hneighbors.Demarshall(&marshalled_dataset);
-
-        /*pointers are garbage, until configuration is carried out: */
-        nodes=NULL;
-        matice=NULL;
-        matpar=NULL;
-        neighbors=NULL;
+/*FUNCTION Penta::MyRank {{{1*/
+int    Penta::MyRank(void){ 
+        extern int my_rank;
+        return my_rank; 
+}
+/*}}}*/
+
+/*Update virtual functions definitions: */
+/*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{1*/
+void  Penta::InputUpdateFromConstant(bool constant, int name){
+        /*Nothing updated for now*/
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromConstant(double value, int name);{{{1*/
+void  Penta::InputUpdateFromConstant(double constant, int name){
+        /*Nothing updated for now*/
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromConstant(int value, int name);{{{1*/
+void  Penta::InputUpdateFromConstant(int constant, int name){
+        /*Nothing updated for now*/
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolution {{{1*/
+void  Penta::InputUpdateFromSolution(double* solution){
+
+        int analysis_type;
+
+        /*retreive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
+        if (analysis_type==ControlAnalysisEnum){
+                InputUpdateFromSolutionDiagnosticHoriz( solution);
+        }
+        else if (analysis_type==DiagnosticHorizAnalysisEnum){
+                InputUpdateFromSolutionDiagnosticHoriz( solution);
+        }
+        else if (analysis_type==DiagnosticStokesAnalysisEnum){
+                InputUpdateFromSolutionDiagnosticStokes( solution);
+        }
+        else if (analysis_type==SlopeAnalysisEnum){
+                InputUpdateFromSolutionSlopeCompute( solution);
+        }
+        else if (analysis_type==PrognosticAnalysisEnum){
+                InputUpdateFromSolutionPrognostic( solution);
+        }
+        else if (analysis_type==Prognostic2AnalysisEnum){
+                InputUpdateFromSolutionPrognostic2(solution);
+        }
+        else if (analysis_type==BalancedthicknessAnalysisEnum){
+                InputUpdateFromSolutionBalancedthickness( solution);
+        }
+        else if (analysis_type==Balancedthickness2AnalysisEnum){
+                InputUpdateFromSolutionBalancedthickness2( solution);
+        }
+        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
+                InputUpdateFromSolutionBalancedvelocities( solution);
+        }
+        else{
+                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
+void  Penta::InputUpdateFromVector(double* vector, int name, int type){
+
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        switch(type){
+
+                case VertexEnum:
+
+                        /*New PentaVertexInpu*/
+                        double values[6];
+
+                        /*Get values on the 6 vertices*/
+                        for (int i=0;i<6;i++){
+                                values[i]=vector[this->nodes[i]->GetVertexDof()];
+                        }
+
+                        /*update input*/
+                        this->inputs->AddInput(new PentaVertexInput(name,values));
+                        return;
+
+                default:
+
+                        ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
+void  Penta::InputUpdateFromVector(int* vector, int name, int type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
+void  Penta::InputUpdateFromVector(bool* vector, int name, int type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+
+/*Element virtual functions definitions: */
+/*FUNCTION Penta::ComputeBasalStress {{{1*/
+void  Penta::ComputeBasalStress(Vec sigma_b){
+
+        int i,j;
+        const int numgrids=6;
+        int    doflist[numgrids];
+        double xyz_list[numgrids][3];
+        double xyz_list_tria[3][3];
+
+        /*Parameters*/
+        double rho_ice,gravity;
+        double surface_normal[3];
+        double bed_normal[3];
+        double bed;
+        double basalforce[3];
+        double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+        double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+        double stresstensor[6]={0.0};
+        double viscosity;
+        int analysis_type,sub_analysis_type;
+
+        int  dofv[3]={0,1,2};
+        int  dofp[1]={3};
+        double Jdet2d;
+        Tria* tria=NULL;
+
+        /*Gauss*/
+        int     num_gauss,ig;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_coord[4];
+
+        /*Output*/
+        double pressure;
+        double sigma_xx,sigma_yy,sigma_zz;
+        double sigma_xy,sigma_xz,sigma_yz;
+        double surface=0;
+        double value=0;
         
-        /*demarshall inputs and results: */
-        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
-        results=(Results*)DataSetDemarshallRaw(&marshalled_dataset); 
-
-        /*demarshall internal parameters: */
-        memcpy(collapse,marshalled_dataset,numanalyses*sizeof(bool));marshalled_dataset+=numanalyses*sizeof(bool);
-
-        /*parameters: may not exist even yet, so let Configure handle it: */
-        this->parameters=NULL;
-
-        /*return: */
-        *pmarshalled_dataset=marshalled_dataset;
-        return;
-}
-/*}}}*/
-/*FUNCTION Penta::Enum {{{1*/
-int Penta::Enum(void){
-
-        return PentaEnum;
-
-}
-/*}}}*/
-/*FUNCTION Penta::copy {{{1*/
-Object* Penta::copy() {
+        /*flags: */
+        bool onbed;
+
+        /*parameters: */
+        double stokesreconditioning;
+
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
+
+        /*Check analysis_types*/
+        if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!");
+
+        /*recover some inputs: */
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+        
+        if(!onbed){
+                //put zero
+                VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES);
+                return;
+        }
+
+        /*recovre material parameters: */
+        rho_ice=matpar->GetRhoIce();
+        gravity=matpar->GetG();
+
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+        for(i=0;i<3;i++){
+                for(j=0;j<3;j++){
+                        xyz_list_tria[i][j]=xyz_list[i][j];
+                }
+        }
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+
+                        /*Pick up the gaussian point: */
+                        gauss_weight=*(gauss_weights+ig);
+                        gauss_coord[0]=*(first_gauss_area_coord+ig); 
+                        gauss_coord[1]=*(second_gauss_area_coord+ig);
+                        gauss_coord[2]=*(third_gauss_area_coord+ig);
+                        gauss_coord[3]=-1.0; //take the base
+
+                        /*Compute strain rate viscosity and pressure: */
+                        inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
+                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+                        inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum);
+
+                        /*Compute Stress*/
+                        sigma_xx=viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure
+                        sigma_yy=viscosity*epsilon[1]-pressure*stokesreconditioning;
+                        sigma_zz=viscosity*epsilon[2]-pressure*stokesreconditioning;
+                        sigma_xy=viscosity*epsilon[3];
+                        sigma_xz=viscosity*epsilon[4];
+                        sigma_yz=viscosity*epsilon[5];
+
+                        /*Get normal vector to the bed */
+                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
+                        bed_normal[0] = - surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
+                        bed_normal[1] = - surface_normal[1];
+                        bed_normal[2] = - surface_normal[2];
+
+                        /*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];
+
+                        /*Get the Jacobian determinant */
+                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord);
+                        value+=sigma_zz*Jdet2d*gauss_weight;
+                        surface+=Jdet2d*gauss_weight;
+        }
+        value=value/surface;
+
+        /*Add value to output*/
+        VecSetValue(sigma_b,id-1,(const double)value,INSERT_VALUES);
+}
+/*}}}*/
+/*FUNCTION Penta::ComputePressure {{{1*/
+void  Penta::ComputePressure(Vec pg){
 
         int i;
-
-        Penta* penta=NULL;
-
-        penta=new Penta();
-
-        /*copy fields: */
-        penta->id=this->id;
-        if(this->inputs){
-                penta->inputs=(Inputs*)this->inputs->Copy();
-        }
-        else{
-                penta->inputs=new Inputs();
-        }
-        if(this->results){
-                penta->results=(Results*)this->results->Copy();
-        }
-        else{
-                penta->results=new Results();
-        }
-        /*point parameters: */
-        penta->parameters=this->parameters;
-
-        /*now deal with hooks and objects: */
-        penta->InitHookNodes(this->numanalyses);
-        for(i=0;i<this->numanalyses;i++)penta->hnodes[i].copy(&this->hnodes[i]);
-        penta->hmatice.copy(&this->hmatice);
-        penta->hmatpar.copy(&this->hmatpar);
-        penta->hneighbors.copy(&this->hneighbors);
-
-        /*recover objects: */
-        penta->nodes=(Node**)xmalloc(6*sizeof(Node*)); //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<6;i++)penta->nodes[i]=this->nodes[i];
-        penta->matice=(Matice*)penta->hmatice.delivers();
-        penta->matpar=(Matpar*)penta->hmatpar.delivers();
-        penta->neighbors=(Penta**)penta->hneighbors.deliverp();
-
-        return penta;
-}
-/*}}}*/
-
-
-/*Penta functionality: */
+        const int numgrids=6;
+        int    doflist[numgrids];
+        double pressure[numgrids];
+        double rho_ice,g;
+        double surface[numgrids];
+        double xyz_list[numgrids][3];
+        double gauss[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+        /*inputs: */
+        bool onwater;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+
+        /*If on water, skip: */
+        if(onwater)return;
+
+        /*Get node data: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+
+        /*pressure is lithostatic: */
+        //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab
+
+        /*Get dof list on which we will plug the pressure values: */
+        GetDofList1(&doflist[0]);
+
+        /*recover value of surface at grids: */
+        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
+
+        /*pressure is lithostatic: */
+        rho_ice=matpar->GetRhoIce();
+        g=matpar->GetG();
+        for(i=0;i<numgrids;i++){
+                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+        }
+
+        /*plug local pressure values into global pressure vector: */
+        VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION Penta::ComputeStrainRate {{{1*/
+void  Penta::ComputeStrainRate(Vec eps){
+
+        ISSMERROR("Not implemented yet");
+
+}
+/*}}}*/
                 /*FUNCTION Penta::Configure {{{1*/
 void  Penta::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
@@ -471 +605 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CostFunction {{{1*/
+double Penta::CostFunction(void){
+
+        double J;
+        Tria* tria=NULL;
+
+        /*flags: */
+        bool onbed;
+        bool onwater;
+        bool collapse;
+        bool onsurface;
+
+        /*recover some inputs: */
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+        /*If on water, return 0: */
+        if(onwater)return 0;
+
+        /*Bail out if this element if:
+         * -> Non collapsed and not on the surface
+         * -> collapsed (2d model) and not on bed) */
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
+                return 0;
+        }
+        else if (collapse){
+
+                /*This element should be collapsed into a tria element at its base. Create this tria element, 
+                 * and compute CostFunction*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
+                J=tria->CostFunction();
+                delete tria;
+                return J;
+        }
+        else{
+
+                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+                J=tria->CostFunction();
+                delete tria;
+                return J;
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrix {{{1*/
+void  Penta::CreateKMatrix(Mat Kgg){
+
+        int analysis_type;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        /*if debugging mode, check that all pointers exist*/
+        ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
+
+        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
+        if (analysis_type==ControlAnalysisEnum){
+                CreateKMatrixDiagnosticHoriz( Kgg);
+        }
+        else if (analysis_type==DiagnosticHorizAnalysisEnum){
+                CreateKMatrixDiagnosticHoriz( Kgg);
+        }
+        else if (analysis_type==DiagnosticHutterAnalysisEnum){
+                CreateKMatrixDiagnosticHutter( Kgg);
+        }
+        else if (analysis_type==DiagnosticVertAnalysisEnum){
+                CreateKMatrixDiagnosticVert( Kgg);
+        }
+        else if (analysis_type==DiagnosticStokesAnalysisEnum){
+                CreateKMatrixDiagnosticStokes( Kgg);
+        }
+        else if (analysis_type==SlopeAnalysisEnum){
+                CreateKMatrixSlopeCompute( Kgg);
+        }
+        else if (analysis_type==PrognosticAnalysisEnum){
+                CreateKMatrixPrognostic( Kgg);
+        }
+        else if (analysis_type==BalancedthicknessAnalysisEnum){
+                CreateKMatrixBalancedthickness( Kgg);
+        }
+        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
+                CreateKMatrixBalancedvelocities( Kgg);
+        }
+        else if (analysis_type==ThermalAnalysisEnum){
+                CreateKMatrixThermal( Kgg);
+        }
+        else if (analysis_type==MeltingAnalysisEnum){
+                CreateKMatrixMelting( Kgg);
+        }
+        else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
+
+}
+/*}}}*/
+/*FUNCTION Penta::CreatePVector {{{1*/
+void  Penta::CreatePVector(Vec pg){
+
+        int analysis_type,sub_analysis_type;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
+
+        /*if debugging mode, check that all pointers exist*/
+        ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
+
+        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
+        if (analysis_type==ControlAnalysisEnum){
+                CreatePVectorDiagnosticHoriz( pg);
+        }
+        else if (analysis_type==DiagnosticHorizAnalysisEnum){
+                CreatePVectorDiagnosticHoriz( pg);
+        }
+        else if (analysis_type==DiagnosticHutterAnalysisEnum){
+                CreatePVectorDiagnosticHutter( pg);
+        }
+        else if (analysis_type==DiagnosticVertAnalysisEnum){
+                CreatePVectorDiagnosticVert( pg);
+        }
+        else if (analysis_type==DiagnosticStokesAnalysisEnum){
+                CreatePVectorDiagnosticStokes( pg);
+        }
+        else if (analysis_type==SlopeAnalysisEnum){
+                CreatePVectorSlopeCompute( pg);
+        }
+        else if (analysis_type==PrognosticAnalysisEnum){
+                CreatePVectorPrognostic( pg);
+        }
+        else if (analysis_type==BalancedthicknessAnalysisEnum){
+                CreatePVectorBalancedthickness( pg);
+        }
+        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
+                CreatePVectorBalancedvelocities( pg);
+        }
+        else if (analysis_type==ThermalAnalysisEnum){
+                CreatePVectorThermal( pg);
+        }
+        else if (analysis_type==MeltingAnalysisEnum){
+                CreatePVectorMelting( pg);
+        }
+        else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
+
+}
+/*}}}*/
+/*FUNCTION Penta::Du {{{1*/
+void  Penta::Du(Vec du_g){
+
+        int i;
+        Tria* tria=NULL;
+
+        /*inputs: */
+        bool onwater;
+        bool collapse;
+        bool onsurface;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+        /*If on water, skip: */
+        if(onwater)return;
+
+        /*Bail out if this element if:
+         * -> Non collapsed and not on the surface
+         * -> collapsed (2d model) and not on bed) */
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
+                return;
+        }
+        else if (collapse){
+
+                /*This element should be collapsed into a tria element at its base. Create this tria element, 
+                 * and compute Du*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
+                tria->Du(du_g);
+                delete tria;
+                return;
+        }
+        else{
+
+                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+                tria->Du(du_g);
+                delete tria;
+                return;
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::GetBedList {{{1*/
+void Penta::GetBedList(double* bedlist){
+
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        
+        inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum);
+
+}
+/*}}}*/
+/*FUNCTION Penta::GetMatPar {{{1*/
+void* Penta::GetMatPar(){
+
+        return matpar;
+}
+/*}}}*/
+/*FUNCTION Penta::GetNodes {{{1*/
+void  Penta::GetNodes(void** vpnodes){
+
+        int i;
+        Node** pnodes=NULL;
+        
+        /*virtual object: */
+        pnodes=(Node**)vpnodes;
+
+        for(i=0;i<6;i++){
+                pnodes[i]=nodes[i];
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::GetOnBed {{{1*/
+bool Penta::GetOnBed(){
+        
+        bool onbed;
+
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        return onbed;
+}
+/*}}}*/
+/*FUNCTION Penta::GetShelf {{{1*/
+bool   Penta::GetShelf(){
+        bool onshelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
+
+        return onshelf;
+}
+/*}}}*/
+/*FUNCTION Penta::GetSolutionFromInputs(Vec solution){{{1*/
+void  Penta::GetSolutionFromInputs(Vec solution){
+
+        int analysis_type,sub_analysis_type;
+
+        /*retrive parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
+
+
+        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
+        if (analysis_type==DiagnosticAnalysisEnum){
+                if (sub_analysis_type==HorizAnalysisEnum){
+                        GetSolutionFromInputsDiagnosticHoriz(solution);
+                }
+                else if(sub_analysis_type==VertAnalysisEnum){
+                        GetSolutionFromInputsDiagnosticVert(solution);
+                }
+                else if(sub_analysis_type==StokesAnalysisEnum){
+                        GetSolutionFromInputsDiagnosticStokes(solution);
+                }
+                else ISSMERROR("sub_analysis: %i (%s) not supported yet",sub_analysis_type,EnumAsString(sub_analysis_type));
+        }
+        else{
+                ISSMERROR("analysis: %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::GetThicknessList {{{1*/
+void Penta::GetThicknessList(double* thickness_list){
+
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum);
+
+}
+/*}}}*/
+/*FUNCTION Penta::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
+void  Penta::GetVectorFromInputs(Vec vector,int NameEnum){
+
+        int i;
+        const int numvertices=6;
+        int doflist1[numvertices];
+
+        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
+        for(i=0;i<this->inputs->Size();i++){
+                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+                if(input->EnumType()==NameEnum){
+                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+                        this->GetDofList1(&doflist1[0]);
+                        input->GetVectorFromInputs(vector,&doflist1[0]);
+                        break;
+                }
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::Gradj {{{1*/
+void  Penta::Gradj(Vec gradient,int control_type){
+
+        /*inputs: */
+        bool onwater;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+
+        /*If on water, skip grad (=0): */
+        if(onwater)return;
+
+        if (control_type==DragCoefficientEnum){
+                GradjDrag(gradient);
+        }
+        else if (control_type=RheologyBEnum){
+                GradjB(gradient);
+        }
+        else ISSMERROR("%s%i","control type not supported yet: ",control_type);
+}
+/*}}}*/
+/*FUNCTION Penta::GradjDrag {{{1*/
+void  Penta::GradjDrag(Vec gradient){
+
+        int i;
+        Tria* tria=NULL;
+        TriaVertexInput* triavertexinput=NULL;
+        double temp_gradient[6]={0,0,0,0,0,0};
+
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+        int analysis_type,sub_analysis_type;
+
+        /*retrieve parameters: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
+        /*If on water, skip: */
+        if(onwater)return;
+
+        /*If on shelf, skip: */
+        if(shelf)return;
+
+        /*Bail out if this element does not touch the bedrock: */
+        if (!onbed) return;
+
+        if (sub_analysis_type==HorizAnalysisEnum){
+
+                /*MacAyeal or Pattyn*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->GradjDrag(gradient);
+                delete tria;
+        }
+        else if (sub_analysis_type==StokesAnalysisEnum){
+
+                /*Stokes*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+                tria->GradjDragStokes(gradient);
+                delete tria;
+        }
+        else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
+
+
+}
+/*}}}*/
+/*FUNCTION Penta::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/
+void  Penta::InputAXPY(int YEnum, double scalar, int XEnum){
+
+        Input* xinput=NULL;
+        Input* yinput=NULL;
+
+        /*Find x and y inputs: */
+        xinput=(Input*)this->inputs->GetInput(XEnum);
+        yinput=(Input*)this->inputs->GetInput(YEnum);
+
+        /*some checks: */
+        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
+        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
+
+        /*Scale: */
+        yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
+/*FUNCTION Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
+void  Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){
+
+        Input* input=NULL;
+
+        /*Find input: */
+        input=(Input*)this->inputs->GetInput(control_type);
+        
+        /*Do nothing if we  don't find it: */
+        if(!input)return;
+
+        /*Constrain input using cm_min and cm_max: */
+        input->Constrain(cm_min,cm_max);
+
+}
+/*}}}*/
+/*FUNCTION Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
+void  Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
+
+        int     i;
+        Input** new_inputs=NULL;
+        Input** old_inputs=NULL;
+        int     converged=1;
+
+        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
+        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
+        
+        for(i=0;i<num_enums/2;i++){
+                new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
+                old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
+                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
+                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
+        }
+
+        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
+        for(i=0;i<num_criterionenums;i++){
+                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
+                if(eps[i]>criterionvalues[i]) converged=0; 
+        }
+
+        /*Assign output pointers:*/
+        *pconverged=converged;
+
+}
+/*}}}*/
+/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
+void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
+        ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)");
+}
+/*}}}*/
+/*FUNCTION Penta::InputDuplicate(int original_enum,int new_enum){{{1*/
+void  Penta::InputDuplicate(int original_enum,int new_enum){
+
+        Input* original=NULL;
+        Input* copy=NULL;
+
+        /*Make a copy of the original input: */
+        original=(Input*)this->inputs->GetInput(original_enum);
+        copy=(Input*)original->copy();
+
+        /*Change copy enum to reinitialized_enum: */
+        copy->ChangeEnum(new_enum);
+
+        /*Add copy into inputs, it will wipe off the one already there: */
+        inputs->AddObject((Input*)copy);
+}
+/*}}}*/
+/*FUNCTION Penta::InputScale(int enum_type,double scale_factor){{{1*/
+void  Penta::InputScale(int enum_type,double scale_factor){
+
+        Input* input=NULL;
+
+        /*Make a copy of the original input: */
+        input=(Input*)this->inputs->GetInput(enum_type);
+
+        /*Scale: */
+        input->Scale(scale_factor);
+}
+/*}}}*/
+/*FUNCTION Penta::InputToResult(int enum_type,int step,double time){{{1*/
+void  Penta::InputToResult(int enum_type,int step,double time){
+
+        int    i;
+        bool   found = false;
+        Input *input = NULL;
+
+        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
+        for (i=0;i<this->inputs->Size();i++){
+                input=(Input*)this->inputs->GetObjectByOffset(i);
+                if (input->EnumType()==enum_type){
+                        found=true;
+                        break;
+                }
+        }
+
+        /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 
+         * object out of the input, with the additional step and time information: */
+        this->results->AddObject((Object*)input->SpawnResult(step,time));
+
+}
+/*}}}*/
+/*FUNCTION Penta::MassFlux {{{1*/
+double Penta::MassFlux( double* segment){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/
+void  Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vx_values[numgrids];
+        double  maxabsvx;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+
+        /*now, compute maximum:*/
+        maxabsvx=fabs(vx_values[0]);
+        for(i=1;i<numgrids;i++){
+                if (fabs(vx_values[i])>maxabsvx)maxabsvx=fabs(vx_values[i]);
+        }
+
+        /*Assign output pointers:*/
+        *pmaxabsvx=maxabsvx;
+}
+/*}}}*/
+/*FUNCTION Penta::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/
+void  Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vy_values[numgrids];
+        double  maxabsvy;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+
+        /*now, compute maximum:*/
+        maxabsvy=fabs(vy_values[0]);
+        for(i=1;i<numgrids;i++){
+                if (fabs(vy_values[i])>maxabsvy)maxabsvy=fabs(vy_values[i]);
+        }
+
+        /*Assign output pointers:*/
+        *pmaxabsvy=maxabsvy;
+}
+/*}}}*/
+/*FUNCTION Penta::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/
+void  Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vz_values[numgrids];
+        double  maxabsvz;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute maximum:*/
+        maxabsvz=fabs(vz_values[0]);
+        for(i=1;i<numgrids;i++){
+                if (fabs(vz_values[i])>maxabsvz)maxabsvz=fabs(vz_values[i]);
+        }
+
+        /*Assign output pointers:*/
+        *pmaxabsvz=maxabsvz;
+}
+/*}}}*/
+/*FUNCTION Penta::MaxVel(double* pmaxvel, bool process_units);{{{1*/
+void  Penta::MaxVel(double* pmaxvel, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vx_values[numgrids];
+        double  vy_values[numgrids];
+        double  vz_values[numgrids];
+        double  vel_values[numgrids];
+        double  maxvel;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute maximum of velocity :*/
+        if(dim==2){
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
+        }
+        else{
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
+        }
+
+        /*now, compute maximum:*/
+        maxvel=vel_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vel_values[i]>maxvel)maxvel=vel_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvel=maxvel;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MaxVx(double* pmaxvx, bool process_units);{{{1*/
+void  Penta::MaxVx(double* pmaxvx, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vx_values[numgrids];
+        double  maxvx;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+
+        /*now, compute maximum:*/
+        maxvx=vx_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vx_values[i]>maxvx)maxvx=vx_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvx=maxvx;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MaxVy(double* pmaxvy, bool process_units);{{{1*/
+void  Penta::MaxVy(double* pmaxvy, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vy_values[numgrids];
+        double  maxvy;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+
+        /*now, compute maximum:*/
+        maxvy=vy_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vy_values[i]>maxvy)maxvy=vy_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvy=maxvy;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MaxVz(double* pmaxvz, bool process_units);{{{1*/
+void  Penta::MaxVz(double* pmaxvz, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vz_values[numgrids];
+        double  maxvz;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute maximum:*/
+        maxvz=vz_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vz_values[i]>maxvz)maxvz=vz_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pmaxvz=maxvz;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MinVel(double* pminvel, bool process_units);{{{1*/
+void  Penta::MinVel(double* pminvel, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vx_values[numgrids];
+        double  vy_values[numgrids];
+        double  vz_values[numgrids];
+        double  vel_values[numgrids];
+        double  minvel;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute minimum of velocity :*/
+        if(dim==2){
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
+        }
+        else{
+                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
+        }
+
+        /*now, compute minimum:*/
+        minvel=vel_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vel_values[i]<minvel)minvel=vel_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvel=minvel;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MinVx(double* pminvx, bool process_units);{{{1*/
+void  Penta::MinVx(double* pminvx, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vx_values[numgrids];
+        double  minvx;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
+
+        /*now, compute minimum:*/
+        minvx=vx_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vx_values[i]<minvx)minvx=vx_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvx=minvx;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MinVy(double* pminvy, bool process_units);{{{1*/
+void  Penta::MinVy(double* pminvy, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vy_values[numgrids];
+        double  minvy;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
+
+        /*now, compute minimum:*/
+        minvy=vy_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vy_values[i]<minvy)minvy=vy_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvy=minvy;
+
+}
+/*}}}*/
+/*FUNCTION Penta::MinVz(double* pminvz, bool process_units);{{{1*/
+void  Penta::MinVz(double* pminvz, bool process_units){
+
+        int i;
+        int dim;
+        const int numgrids=6;
+        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+        double  vz_values[numgrids];
+        double  minvz;
+
+        /*retrieve dim parameter: */
+        parameters->FindParam(&dim,DimEnum);
+
+        /*retrive velocity values at nodes */
+        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
+
+        /*now, compute minimum:*/
+        minvz=vz_values[0];
+        for(i=1;i<numgrids;i++){
+                if (vz_values[i]<minvz)minvz=vz_values[i];
+        }
+
+        /*Assign output pointers:*/
+        *pminvz=minvz;
+
+}
+/*}}}*/
+/*FUNCTION Penta::Misfit {{{1*/
+double Penta::Misfit(void){
+
+        double J;
+        Tria* tria=NULL;
+
+        /*inputs: */
+        bool onwater;
+        bool collapse;
+        bool onsurface;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+        /*If on water, return 0: */
+        if(onwater)return 0;
+
+        /*Bail out if this element if:
+         * -> Non collapsed and not on the surface
+         * -> collapsed (2d model) and not on bed) */
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
+                return 0;
+        }
+        else if (collapse){
+
+                /*This element should be collapsed into a tria element at its base. Create this tria element, 
+                 * and compute Misfit*/
+                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
+                J=tria->Misfit();
+                delete tria;
+                return J;
+        }
+        else{
+
+                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+                J=tria->Misfit();
+                delete tria;
+                return J;
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::PatchFill(int* pcount, Patch* patch){{{1*/
+void  Penta::PatchFill(int* pcount, Patch* patch){
+
+        int i;
+        int count;
+        int vertices_ids[6];
+
+
+        /*recover pointer: */
+        count=*pcount;
+                
+        /*will be needed later: */
+        for(i=0;i<6;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
+
+        for(i=0;i<this->results->Size();i++){
+                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
+
+                /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand 
+                 *it to the result object, to fill the rest: */
+                patch->fillelementinfo(count,this->id,vertices_ids,6);
+                elementresult->PatchFill(count,patch);
+
+                /*increment counter: */
+                count++;
+        }
+
+        /*Assign output pointers:*/
+        *pcount=count;
+}
+/*FUNCTION Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/
+void  Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
+
+        int     i;
+        
+        /*output: */
+        int     numrows     = 0;
+        int     numvertices = 0;
+        int     numnodes    = 0;
+
+        /*Go through all the results objects, and update the counters: */
+        for (i=0;i<this->results->Size();i++){
+                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
+                /*first, we have one more result: */
+                numrows++;
+                /*now, how many vertices and how many nodal values for this result? :*/
+                numvertices=6; //this is a penta element, with 6 vertices
+                numnodes=elementresult->NumberOfNodalValues(); //ask result object.
+        }
+
+        /*Assign output pointers:*/
+        *pnumrows=numrows;
+        *pnumvertices=numvertices;
+        *pnumnodes=numnodes;
+        
+}
+/*}}}*/
+/*FUNCTION Penta::ProcessResultsUnits(void){{{1*/
+void  Penta::ProcessResultsUnits(void){
+
+        int i;
+
+        for(i=0;i<this->results->Size();i++){
+                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
+                elementresult->ProcessUnits(this->parameters);
+        }
+
+}
+/*}}}*/
+/*FUNCTION Penta::SurfaceArea {{{1*/
+double Penta::SurfaceArea(void){
+
+        double S;
+        Tria* tria=NULL;
+
+        /*inputs: */
+        bool onwater;
+        bool collapse;
+        bool onsurface;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+        /*If on water, return 0: */
+        if(onwater)return 0;
+
+        /*Bail out if this element if:
+         * -> Non collapsed and not on the surface
+         * -> collapsed (2d model) and not on bed) */
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
+                return 0;
+        }
+        else if (collapse){
+
+                /*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); //grids 0, 1 and 2 make the new tria (lower face).
+                S=tria->SurfaceArea();
+                delete tria;
+                return S;
+        }
+        else{
+
+                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
+                S=tria->SurfaceArea();
+                delete tria;
+                return S;
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::Update(IoModel* iomodel,int analysis_counter,int analysis_type) {{{1*/
+void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type){ 
+
+        /*Intermediaries*/
+        IssmInt i;
+        int penta_node_ids[6];
+        int penta_vertex_ids[6];
+        double nodeinputs[6];
+
+        /*Checks if debuging*/
+        /*{{{2*/
+        ISSMASSERT(iomodel->elements);
+        /*}}}*/
+
+        /*Recover vertices ids needed to initialize inputs*/
+        for(i=0;i<6;i++){ 
+                penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+        }
+
+        /*Recover nodes ids needed to initialize the node hook.*/
+        for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
+                penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+        }
+
+        /*hooks: */
+        this->SetHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
+
+        //add as many inputs per element as requested: 
+        if (iomodel->thickness) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));
+        }
+        if (iomodel->surface) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));
+        }
+        if (iomodel->bed) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs));
+        }
+        if (iomodel->drag_coefficient) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs));
+
+                if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
+                if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
+                this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
+
+        }
+        if (iomodel->melting_rate) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));
+        }
+        if (iomodel->accumulation_rate) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));
+        }
+        if (iomodel->geothermalflux) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs));
+        }       
+        if (iomodel->pressure) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));
+        }
+        if (iomodel->temperature) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs));
+        }
+        if (iomodel->dhdt) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs));
+        }
+        /*vx,vy and vz: */
+        if (iomodel->vx) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
+                this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));
+        }
+        if (iomodel->vy) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
+                this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));
+        }
+        if (iomodel->vz) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
+                this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
+        }
+        if (iomodel->vx_obs) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs));
+        }
+        if (iomodel->vy_obs) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs));
+        }
+        if (iomodel->vz_obs) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
+                this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs));
+        }
+        if (iomodel->weights) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs));
+        }
+
+        /*default vx,vy and vz: either observation or 0 */
+        if (!iomodel->vx){
+                if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
+                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
+                this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
+                this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs));
+        }
+        if (!iomodel->vy){
+                if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
+                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
+                this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
+                this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs));
+        }
+        if (!iomodel->vz){
+                if (iomodel->vz_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts;
+                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
+                this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs));
+                this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs));
+        }
+        if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
+        if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
+        if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
+        if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
+
+        //elements of type 3 are macayeal type penta. we collapse the formulation on their base.
+        if(iomodel->elements_type){
+                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ 
+                        collapse[analysis_counter]=true;
+                }
+                else{
+                        collapse[analysis_counter]=false;
+                }
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateGeometry{{{1*/
+void  Penta::UpdateGeometry(void){
+
+        /*Intermediaries*/
+        int i;
+        double rho_ice,rho_water;
+
+        /*If shelf: hydrostatic equilibrium*/
+        if (this->GetShelf()){
+
+                /*recover material parameters: */
+                rho_ice=matpar->GetRhoIce();
+                rho_water=matpar->GetRhoWater();
+
+                /*Create New Surface: s = (1-rho_ice/rho_water) h*/
+                InputDuplicate(ThicknessEnum,SurfaceEnum);     //1: copy thickness into surface
+                InputScale(SurfaceEnum,(1-rho_ice/rho_water)); //2: surface = surface * (1-di)
+
+                /*Create New Bed b = -rho_ice/rho_water h*/
+                InputDuplicate(ThicknessEnum,BedEnum);         //1: copy thickness into bed
+                InputScale(BedEnum, -rho_ice/rho_water);       //2: bed = bed * (-di)
+        }
+
+        /*If sheet: surface = bed + thickness*/
+        else{
+
+                /*The bed does not change, update surface only s = b + h*/
+                InputDuplicate(BedEnum,SurfaceEnum);          //1: copy bed into surface
+                InputAXPY(SurfaceEnum,1.0,ThicknessEnum);     //2: surface = surface + 1 * thickness
+        }
+
+}
+/*}}}*/
+
+/*Penta specific routines: */
 /*FUNCTION Penta::SpawnTria {{{1*/
 void*  Penta::SpawnTria(int g0, int g1, int g2){
@@ -602 +1881 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputToResult(int enum_type,int step,double time){{{1*/
-void  Penta::InputToResult(int enum_type,int step,double time){
-
-        int    i;
-        bool   found = false;
-        Input *input = NULL;
-
-        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
-        for (i=0;i<this->inputs->Size();i++){
-                input=(Input*)this->inputs->GetObjectByOffset(i);
-                if (input->EnumType()==enum_type){
-                        found=true;
-                        break;
-                }
-        }
-
-        /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 
-         * object out of the input, with the additional step and time information: */
-        this->results->AddObject((Object*)input->SpawnResult(step,time));
-
-}
-/*}}}*/
-/*FUNCTION Penta::ProcessResultsUnits(void){{{1*/
-void  Penta::ProcessResultsUnits(void){
-
-        int i;
-
-        for(i=0;i<this->results->Size();i++){
-                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
-                elementresult->ProcessUnits(this->parameters);
-        }
-
-}
-/*}}}*/
-
-/*updates: */
 /*FUNCTION Penta::UpdateFromDakota {{{1*/
 void  Penta::UpdateFromDakota(void* vinputs){
@@ -645 +1888 @@
 }
 /*}}}*/
-/*FUNCTION Penta::UpdateGeometry{{{1*/
-void  Penta::UpdateGeometry(void){
-
-        /*Intermediaries*/
-        int i;
-        double rho_ice,rho_water;
-
-        /*If shelf: hydrostatic equilibrium*/
-        if (this->GetShelf()){
-
-                /*recover material parameters: */
-                rho_ice=matpar->GetRhoIce();
-                rho_water=matpar->GetRhoWater();
-
-                /*Create New Surface: s = (1-rho_ice/rho_water) h*/
-                InputDuplicate(ThicknessEnum,SurfaceEnum);     //1: copy thickness into surface
-                InputScale(SurfaceEnum,(1-rho_ice/rho_water)); //2: surface = surface * (1-di)
-
-                /*Create New Bed b = -rho_ice/rho_water h*/
-                InputDuplicate(ThicknessEnum,BedEnum);         //1: copy thickness into bed
-                InputScale(BedEnum, -rho_ice/rho_water);       //2: bed = bed * (-di)
-        }
-
-        /*If sheet: surface = bed + thickness*/
-        else{
-
-                /*The bed does not change, update surface only s = b + h*/
-                InputDuplicate(BedEnum,SurfaceEnum);          //1: copy bed into surface
-                InputAXPY(SurfaceEnum,1.0,ThicknessEnum);     //2: surface = surface + 1 * thickness
-        }
-
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolution {{{1*/
-void  Penta::InputUpdateFromSolution(double* solution){
-
-        int analysis_type;
-
-        /*retreive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-        if (analysis_type==ControlAnalysisEnum){
-                InputUpdateFromSolutionDiagnosticHoriz( solution);
-        }
-        else if (analysis_type==DiagnosticHorizAnalysisEnum){
-                InputUpdateFromSolutionDiagnosticHoriz( solution);
-        }
-        else if (analysis_type==DiagnosticStokesAnalysisEnum){
-                InputUpdateFromSolutionDiagnosticStokes( solution);
-        }
-        else if (analysis_type==SlopeAnalysisEnum){
-                InputUpdateFromSolutionSlopeCompute( solution);
-        }
-        else if (analysis_type==PrognosticAnalysisEnum){
-                InputUpdateFromSolutionPrognostic( solution);
-        }
-        else if (analysis_type==Prognostic2AnalysisEnum){
-                InputUpdateFromSolutionPrognostic2(solution);
-        }
-        else if (analysis_type==BalancedthicknessAnalysisEnum){
-                InputUpdateFromSolutionBalancedthickness( solution);
-        }
-        else if (analysis_type==Balancedthickness2AnalysisEnum){
-                InputUpdateFromSolutionBalancedthickness2( solution);
-        }
-        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
-                InputUpdateFromSolutionBalancedvelocities( solution);
-        }
-        else{
-                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
-        }
-}
-/*Object functions*/
 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
 void  Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
@@ -867 +2036 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetSolutionFromInputs(Vec solution){{{1*/
-void  Penta::GetSolutionFromInputs(Vec solution){
-
-        int analysis_type,sub_analysis_type;
-
-        /*retrive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
-
-
-        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-        if (analysis_type==DiagnosticAnalysisEnum){
-                if (sub_analysis_type==HorizAnalysisEnum){
-                        GetSolutionFromInputsDiagnosticHoriz(solution);
-                }
-                else if(sub_analysis_type==VertAnalysisEnum){
-                        GetSolutionFromInputsDiagnosticVert(solution);
-                }
-                else if(sub_analysis_type==StokesAnalysisEnum){
-                        GetSolutionFromInputsDiagnosticStokes(solution);
-                }
-                else ISSMERROR("sub_analysis: %i (%s) not supported yet",sub_analysis_type,EnumAsString(sub_analysis_type));
-        }
-        else{
-                ISSMERROR("analysis: %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
-        }
-}
-/*}}}*/
 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){{{1*/
 void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
@@ -1002 +2143 @@
 }
 /*}}}*/
-
-/*Object functions*/
 /*FUNCTION Penta::IsInput{{{1*/
 bool Penta::IsInput(int name){
@@ -1042 +2181 @@
         upper_penta=(Penta*)neighbors[1]; //first one under, second one above
         return upper_penta;
-
-}
-/*}}}*/
-/*FUNCTION Penta::ComputeBasalStress {{{1*/
-void  Penta::ComputeBasalStress(Vec sigma_b){
-
-        int i,j;
-        const int numgrids=6;
-        int    doflist[numgrids];
-        double xyz_list[numgrids][3];
-        double xyz_list_tria[3][3];
-
-        /*Parameters*/
-        double rho_ice,gravity;
-        double surface_normal[3];
-        double bed_normal[3];
-        double bed;
-        double basalforce[3];
-        double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-        double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-        double stresstensor[6]={0.0};
-        double viscosity;
-        int analysis_type,sub_analysis_type;
-
-        int  dofv[3]={0,1,2};
-        int  dofp[1]={3};
-        double Jdet2d;
-        Tria* tria=NULL;
-
-        /*Gauss*/
-        int     num_gauss,ig;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* gauss_weights           =  NULL;
-        double  gauss_weight;
-        double  gauss_coord[4];
-
-        /*Output*/
-        double pressure;
-        double sigma_xx,sigma_yy,sigma_zz;
-        double sigma_xy,sigma_xz,sigma_yz;
-        double surface=0;
-        double value=0;
-        
-        /*flags: */
-        bool onbed;
-
-        /*parameters: */
-        double stokesreconditioning;
-
-
-        /*retrive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
-
-        /*Check analysis_types*/
-        if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!");
-
-        /*recover some inputs: */
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-
-        /*retrieve some parameters: */
-        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
-        
-        if(!onbed){
-                //put zero
-                VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES);
-                return;
-        }
-
-        /*recovre material parameters: */
-        rho_ice=matpar->GetRhoIce();
-        gravity=matpar->GetG();
-
-        /* Get node coordinates and dof list: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-        for(i=0;i<3;i++){
-                for(j=0;j<3;j++){
-                        xyz_list_tria[i][j]=xyz_list[i][j];
-                }
-        }
-
-        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-        GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
-
-        /* Start  looping on the number of gaussian points: */
-        for (ig=0; ig<num_gauss; ig++){
-
-                        /*Pick up the gaussian point: */
-                        gauss_weight=*(gauss_weights+ig);
-                        gauss_coord[0]=*(first_gauss_area_coord+ig); 
-                        gauss_coord[1]=*(second_gauss_area_coord+ig);
-                        gauss_coord[2]=*(third_gauss_area_coord+ig);
-                        gauss_coord[3]=-1.0; //take the base
-
-                        /*Compute strain rate viscosity and pressure: */
-                        inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
-                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-                        inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum);
-
-                        /*Compute Stress*/
-                        sigma_xx=viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure
-                        sigma_yy=viscosity*epsilon[1]-pressure*stokesreconditioning;
-                        sigma_zz=viscosity*epsilon[2]-pressure*stokesreconditioning;
-                        sigma_xy=viscosity*epsilon[3];
-                        sigma_xz=viscosity*epsilon[4];
-                        sigma_yz=viscosity*epsilon[5];
-
-                        /*Get normal vector to the bed */
-                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
-                        bed_normal[0] = - surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
-                        bed_normal[1] = - surface_normal[1];
-                        bed_normal[2] = - surface_normal[2];
-
-                        /*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];
-
-                        /*Get the Jacobian determinant */
-                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord);
-                        value+=sigma_zz*Jdet2d*gauss_weight;
-                        surface+=Jdet2d*gauss_weight;
-        }
-        value=value/surface;
-
-        /*Add value to output*/
-        VecSetValue(sigma_b,id-1,(const double)value,INSERT_VALUES);
-}
-/*}}}*/
-/*FUNCTION Penta::ComputePressure {{{1*/
-void  Penta::ComputePressure(Vec pg){
-
-        int i;
-        const int numgrids=6;
-        int    doflist[numgrids];
-        double pressure[numgrids];
-        double rho_ice,g;
-        double surface[numgrids];
-        double xyz_list[numgrids][3];
-        double gauss[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-
-        /*inputs: */
-        bool onwater;
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-
-        /*If on water, skip: */
-        if(onwater)return;
-
-        /*Get node data: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-
-        /*pressure is lithostatic: */
-        //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab
-
-        /*Get dof list on which we will plug the pressure values: */
-        GetDofList1(&doflist[0]);
-
-        /*recover value of surface at grids: */
-        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
-
-        /*pressure is lithostatic: */
-        rho_ice=matpar->GetRhoIce();
-        g=matpar->GetG();
-        for(i=0;i<numgrids;i++){
-                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
-        }
-
-        /*plug local pressure values into global pressure vector: */
-        VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
-
-}
-/*}}}*/
-/*FUNCTION Penta::ComputeStrainRate {{{1*/
-void  Penta::ComputeStrainRate(Vec eps){
-
-        ISSMERROR("Not implemented yet");
-
-}
-/*}}}*/
-/*FUNCTION Penta::CostFunction {{{1*/
-double Penta::CostFunction(void){
-
-        double J;
-        Tria* tria=NULL;
-
-        /*flags: */
-        bool onbed;
-        bool onwater;
-        bool collapse;
-        bool onsurface;
-
-        /*recover some inputs: */
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-        inputs->GetParameterValue(&collapse,CollapseEnum);
-        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
-
-        /*If on water, return 0: */
-        if(onwater)return 0;
-
-        /*Bail out if this element if:
-         * -> Non collapsed and not on the surface
-         * -> collapsed (2d model) and not on bed) */
-        if ((!collapse && !onsurface) || (collapse && !onbed)){
-                return 0;
-        }
-        else if (collapse){
-
-                /*This element should be collapsed into a tria element at its base. Create this tria element, 
-                 * and compute CostFunction*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-                J=tria->CostFunction();
-                delete tria;
-                return J;
-        }
-        else{
-
-                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                J=tria->CostFunction();
-                delete tria;
-                return J;
-        }
-}
-/*}}}*/
-/*FUNCTION Penta::CreateKMatrix {{{1*/
-void  Penta::CreateKMatrix(Mat Kgg){
-
-        int analysis_type;
-
-        /*retrive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-        /*if debugging mode, check that all pointers exist*/
-        ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
-
-        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-        if (analysis_type==ControlAnalysisEnum){
-                CreateKMatrixDiagnosticHoriz( Kgg);
-        }
-        else if (analysis_type==DiagnosticHorizAnalysisEnum){
-                CreateKMatrixDiagnosticHoriz( Kgg);
-        }
-        else if (analysis_type==DiagnosticHutterAnalysisEnum){
-                CreateKMatrixDiagnosticHutter( Kgg);
-        }
-        else if (analysis_type==DiagnosticVertAnalysisEnum){
-                CreateKMatrixDiagnosticVert( Kgg);
-        }
-        else if (analysis_type==DiagnosticStokesAnalysisEnum){
-                CreateKMatrixDiagnosticStokes( Kgg);
-        }
-        else if (analysis_type==SlopeAnalysisEnum){
-                CreateKMatrixSlopeCompute( Kgg);
-        }
-        else if (analysis_type==PrognosticAnalysisEnum){
-                CreateKMatrixPrognostic( Kgg);
-        }
-        else if (analysis_type==BalancedthicknessAnalysisEnum){
-                CreateKMatrixBalancedthickness( Kgg);
-        }
-        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
-                CreateKMatrixBalancedvelocities( Kgg);
-        }
-        else if (analysis_type==ThermalAnalysisEnum){
-                CreateKMatrixThermal( Kgg);
-        }
-        else if (analysis_type==MeltingAnalysisEnum){
-                CreateKMatrixMelting( Kgg);
-        }
-        else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
 
 }
@@ -2353 +3218 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVector {{{1*/
-void  Penta::CreatePVector(Vec pg){
-
-        int analysis_type,sub_analysis_type;
-
-        /*retrive parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
-
-        /*if debugging mode, check that all pointers exist*/
-        ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs);
-
-        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-        if (analysis_type==ControlAnalysisEnum){
-                CreatePVectorDiagnosticHoriz( pg);
-        }
-        else if (analysis_type==DiagnosticHorizAnalysisEnum){
-                CreatePVectorDiagnosticHoriz( pg);
-        }
-        else if (analysis_type==DiagnosticHutterAnalysisEnum){
-                CreatePVectorDiagnosticHutter( pg);
-        }
-        else if (analysis_type==DiagnosticVertAnalysisEnum){
-                CreatePVectorDiagnosticVert( pg);
-        }
-        else if (analysis_type==DiagnosticStokesAnalysisEnum){
-                CreatePVectorDiagnosticStokes( pg);
-        }
-        else if (analysis_type==SlopeAnalysisEnum){
-                CreatePVectorSlopeCompute( pg);
-        }
-        else if (analysis_type==PrognosticAnalysisEnum){
-                CreatePVectorPrognostic( pg);
-        }
-        else if (analysis_type==BalancedthicknessAnalysisEnum){
-                CreatePVectorBalancedthickness( pg);
-        }
-        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
-                CreatePVectorBalancedvelocities( pg);
-        }
-        else if (analysis_type==ThermalAnalysisEnum){
-                CreatePVectorThermal( pg);
-        }
-        else if (analysis_type==MeltingAnalysisEnum){
-                CreatePVectorMelting( pg);
-        }
-        else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
-
-}
-/*}}}*/
 /*FUNCTION Penta::CreatePVectorBalancedthickness {{{1*/
 void Penta::CreatePVectorBalancedthickness( Vec pg){
@@ -3242 +4057 @@
 }
 /*}}}*/
-/*FUNCTION Penta::Du {{{1*/
-void  Penta::Du(Vec du_g){
-
-        int i;
-        Tria* tria=NULL;
-
-        /*inputs: */
-        bool onwater;
-        bool collapse;
-        bool onsurface;
-        bool onbed;
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-        inputs->GetParameterValue(&collapse,CollapseEnum);
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
-
-        /*If on water, skip: */
-        if(onwater)return;
-
-        /*Bail out if this element if:
-         * -> Non collapsed and not on the surface
-         * -> collapsed (2d model) and not on bed) */
-        if ((!collapse && !onsurface) || (collapse && !onbed)){
-                return;
-        }
-        else if (collapse){
-
-                /*This element should be collapsed into a tria element at its base. Create this tria element, 
-                 * and compute Du*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-                tria->Du(du_g);
-                delete tria;
-                return;
-        }
-        else{
-
-                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                tria->Du(du_g);
-                delete tria;
-                return;
-        }
-}
-/*}}}*/
 /*FUNCTION Penta::VecExtrude {{{1*/
 void  Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){
@@ -3389 +4159 @@
 
         return;
-}
-/*}}}*/
-/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
-void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
-        ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)");
 }
 /*}}}*/
@@ -3558 +4323 @@
                 B[i]=dh1dh6[2][i];  
         }
-
-}
-/*}}}*/
-/*FUNCTION Penta::GetBedList {{{1*/
-void Penta::GetBedList(double* bedlist){
-
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        
-        inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum);
 
 }
@@ -4127 +4882 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetMatPar {{{1*/
-void* Penta::GetMatPar(){
-
-        return matpar;
-}
-/*}}}*/
 /*FUNCTION Penta::GetMatrixInvert {{{1*/
 void Penta::GetMatrixInvert(double*  Ke_invert, double* Ke){
@@ -4376 +5125 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetNodes {{{1*/
-void  Penta::GetNodes(void** vpnodes){
-
-        int i;
-        Node** pnodes=NULL;
-        
-        /*virtual object: */
-        pnodes=(Node**)vpnodes;
-
-        for(i=0;i<6;i++){
-                pnodes[i]=nodes[i];
-        }
-}
-/*}}}*/
-/*FUNCTION Penta::GetOnBed {{{1*/
-bool Penta::GetOnBed(){
-        
-        bool onbed;
-
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-
-        return onbed;
-}
-/*}}}*/
 /*FUNCTION Penta::GetParameterDerivativeValue {{{1*/
 void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord){
@@ -4471 +5196 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetShelf {{{1*/
-bool   Penta::GetShelf(){
-        bool onshelf;
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
-
-        return onshelf;
-}
-/*}}}*/
-/*FUNCTION Penta::GetThicknessList {{{1*/
-void Penta::GetThicknessList(double* thickness_list){
-
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum);
-
-}
-/*}}}*/
-/*FUNCTION Penta::Gradj {{{1*/
-void  Penta::Gradj(Vec gradient,int control_type){
-
-        /*inputs: */
-        bool onwater;
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-
-        /*If on water, skip grad (=0): */
-        if(onwater)return;
-
-        if (control_type==DragCoefficientEnum){
-                GradjDrag(gradient);
-        }
-        else if (control_type=RheologyBEnum){
-                GradjB(gradient);
-        }
-        else ISSMERROR("%s%i","control type not supported yet: ",control_type);
-}
-/*}}}*/
-/*FUNCTION Penta::GradjDrag {{{1*/
-void  Penta::GradjDrag(Vec gradient){
-
-        int i;
-        Tria* tria=NULL;
-        TriaVertexInput* triavertexinput=NULL;
-        double temp_gradient[6]={0,0,0,0,0,0};
-
-        /*inputs: */
-        bool onwater;
-        bool onbed;
-        bool shelf;
-        int analysis_type,sub_analysis_type;
-
-        /*retrieve parameters: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum);
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
-
-        /*If on water, skip: */
-        if(onwater)return;
-
-        /*If on shelf, skip: */
-        if(shelf)return;
-
-        /*Bail out if this element does not touch the bedrock: */
-        if (!onbed) return;
-
-        if (sub_analysis_type==HorizAnalysisEnum){
-
-                /*MacAyeal or Pattyn*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->GradjDrag(gradient);
-                delete tria;
-        }
-        else if (sub_analysis_type==StokesAnalysisEnum){
-
-                /*Stokes*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->GradjDragStokes(gradient);
-                delete tria;
-        }
-        else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
-
-
-}
-/*}}}*/
 /*FUNCTION Penta::GradjB {{{1*/
 void  Penta::GradjB(Vec gradient){
@@ -4600 +5234 @@
         
 
-}
-/*}}}*/
-/*FUNCTION Penta::MassFlux {{{1*/
-double Penta::MassFlux( double* segment){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Penta::Misfit {{{1*/
-double Penta::Misfit(void){
-
-        double J;
-        Tria* tria=NULL;
-
-        /*inputs: */
-        bool onwater;
-        bool collapse;
-        bool onsurface;
-        bool onbed;
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-        inputs->GetParameterValue(&collapse,CollapseEnum);
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
-
-        /*If on water, return 0: */
-        if(onwater)return 0;
-
-        /*Bail out if this element if:
-         * -> Non collapsed and not on the surface
-         * -> collapsed (2d model) and not on bed) */
-        if ((!collapse && !onsurface) || (collapse && !onbed)){
-                return 0;
-        }
-        else if (collapse){
-
-                /*This element should be collapsed into a tria element at its base. Create this tria element, 
-                 * and compute Misfit*/
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-                J=tria->Misfit();
-                delete tria;
-                return J;
-        }
-        else{
-
-                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                J=tria->Misfit();
-                delete tria;
-                return J;
-        }
 }
 /*}}}*/
@@ -4746 +5330 @@
 }
 /*}}}1*/
-/*FUNCTION Penta::SurfaceArea {{{1*/
-double Penta::SurfaceArea(void){
-
-        double S;
-        Tria* tria=NULL;
-
-        /*inputs: */
-        bool onwater;
-        bool collapse;
-        bool onsurface;
-        bool onbed;
-
-        /*retrieve inputs :*/
-        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-        inputs->GetParameterValue(&collapse,CollapseEnum);
-        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
-
-        /*If on water, return 0: */
-        if(onwater)return 0;
-
-        /*Bail out if this element if:
-         * -> Non collapsed and not on the surface
-         * -> collapsed (2d model) and not on bed) */
-        if ((!collapse && !onsurface) || (collapse && !onbed)){
-                return 0;
-        }
-        else if (collapse){
-
-                /*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); //grids 0, 1 and 2 make the new tria (lower face).
-                S=tria->SurfaceArea();
-                delete tria;
-                return S;
-        }
-        else{
-
-                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                S=tria->SurfaceArea();
-                delete tria;
-                return S;
-        }
-}
-/*}}}*/
 /*FUNCTION Penta::SurfaceNormal {{{1*/
 void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
@@ -4817 +5356 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
-void  Penta::InputUpdateFromVector(double* vector, int name, int type){
-
-        /*Check that name is an element input*/
-        if (!IsInput(name)) return;
-
-        switch(type){
-
-                case VertexEnum:
-
-                        /*New PentaVertexInpu*/
-                        double values[6];
-
-                        /*Get values on the 6 vertices*/
-                        for (int i=0;i<6;i++){
-                                values[i]=vector[this->nodes[i]->GetVertexDof()];
-                        }
-
-                        /*update input*/
-                        this->inputs->AddInput(new PentaVertexInput(name,values));
-                        return;
-
-                default:
-
-                        ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
-        }
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
-void  Penta::InputUpdateFromVector(int* vector, int name, int type){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
-void  Penta::InputUpdateFromVector(bool* vector, int name, int type){
-        ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromConstant(int value, int name);{{{1*/
-void  Penta::InputUpdateFromConstant(int constant, int name){
-        /*Nothing updated for now*/
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{1*/
-void  Penta::InputUpdateFromConstant(bool constant, int name){
-        /*Nothing updated for now*/
-}
-/*}}}*/
-/*FUNCTION Penta::InputUpdateFromConstant(double value, int name);{{{1*/
-void  Penta::InputUpdateFromConstant(double constant, int name){
-        /*Nothing updated for now*/
-}
-/*}}}*/
-/*FUNCTION Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/
-void  Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
-
-        int     i;
-        
-        /*output: */
-        int     numrows     = 0;
-        int     numvertices = 0;
-        int     numnodes    = 0;
-
-        /*Go through all the results objects, and update the counters: */
-        for (i=0;i<this->results->Size();i++){
-                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
-                /*first, we have one more result: */
-                numrows++;
-                /*now, how many vertices and how many nodal values for this result? :*/
-                numvertices=6; //this is a penta element, with 6 vertices
-                numnodes=elementresult->NumberOfNodalValues(); //ask result object.
-        }
-
-        /*Assign output pointers:*/
-        *pnumrows=numrows;
-        *pnumvertices=numvertices;
-        *pnumnodes=numnodes;
-        
-}
-/*}}}*/
-/*FUNCTION Penta::PatchFill(int* pcount, Patch* patch){{{1*/
-void  Penta::PatchFill(int* pcount, Patch* patch){
-
-        int i;
-        int count;
-        int vertices_ids[6];
-
-
-        /*recover pointer: */
-        count=*pcount;
-                
-        /*will be needed later: */
-        for(i=0;i<6;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
-
-        for(i=0;i<this->results->Size();i++){
-                ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
-
-                /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand 
-                 *it to the result object, to fill the rest: */
-                patch->fillelementinfo(count,this->id,vertices_ids,6);
-                elementresult->PatchFill(count,patch);
-
-                /*increment counter: */
-                count++;
-        }
-
-        /*Assign output pointers:*/
-        *pcount=count;
-}
-/*FUNCTION Penta::MinVel(double* pminvel, bool process_units);{{{1*/
-void  Penta::MinVel(double* pminvel, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vx_values[numgrids];
-        double  vy_values[numgrids];
-        double  vz_values[numgrids];
-        double  vel_values[numgrids];
-        double  minvel;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute minimum of velocity :*/
-        if(dim==2){
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
-        }
-        else{
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
-        }
-
-        /*now, compute minimum:*/
-        minvel=vel_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vel_values[i]<minvel)minvel=vel_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvel=minvel;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxVel(double* pmaxvel, bool process_units);{{{1*/
-void  Penta::MaxVel(double* pmaxvel, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vx_values[numgrids];
-        double  vy_values[numgrids];
-        double  vz_values[numgrids];
-        double  vel_values[numgrids];
-        double  maxvel;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-        if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute maximum of velocity :*/
-        if(dim==2){
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2));
-        }
-        else{
-                for(i=0;i<numgrids;i++)vel_values[i]=sqrt(pow(vx_values[i],2)+pow(vy_values[i],2)+pow(vz_values[i],2));
-        }
-
-        /*now, compute maximum:*/
-        maxvel=vel_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vel_values[i]>maxvel)maxvel=vel_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvel=maxvel;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MinVx(double* pminvx, bool process_units);{{{1*/
-void  Penta::MinVx(double* pminvx, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vx_values[numgrids];
-        double  minvx;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-
-        /*now, compute minimum:*/
-        minvx=vx_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vx_values[i]<minvx)minvx=vx_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvx=minvx;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxVx(double* pmaxvx, bool process_units);{{{1*/
-void  Penta::MaxVx(double* pmaxvx, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vx_values[numgrids];
-        double  maxvx;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-
-        /*now, compute maximum:*/
-        maxvx=vx_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vx_values[i]>maxvx)maxvx=vx_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvx=maxvx;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/
-void  Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vx_values[numgrids];
-        double  maxabsvx;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
-
-        /*now, compute maximum:*/
-        maxabsvx=fabs(vx_values[0]);
-        for(i=1;i<numgrids;i++){
-                if (fabs(vx_values[i])>maxabsvx)maxabsvx=fabs(vx_values[i]);
-        }
-
-        /*Assign output pointers:*/
-        *pmaxabsvx=maxabsvx;
-}
-/*}}}*/
-/*FUNCTION Penta::MinVy(double* pminvy, bool process_units);{{{1*/
-void  Penta::MinVy(double* pminvy, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vy_values[numgrids];
-        double  minvy;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-
-        /*now, compute minimum:*/
-        minvy=vy_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vy_values[i]<minvy)minvy=vy_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvy=minvy;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxVy(double* pmaxvy, bool process_units);{{{1*/
-void  Penta::MaxVy(double* pmaxvy, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vy_values[numgrids];
-        double  maxvy;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-
-        /*now, compute maximum:*/
-        maxvy=vy_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vy_values[i]>maxvy)maxvy=vy_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvy=maxvy;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/
-void  Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vy_values[numgrids];
-        double  maxabsvy;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
-
-        /*now, compute maximum:*/
-        maxabsvy=fabs(vy_values[0]);
-        for(i=1;i<numgrids;i++){
-                if (fabs(vy_values[i])>maxabsvy)maxabsvy=fabs(vy_values[i]);
-        }
-
-        /*Assign output pointers:*/
-        *pmaxabsvy=maxabsvy;
-}
-/*}}}*/
-/*FUNCTION Penta::MinVz(double* pminvz, bool process_units);{{{1*/
-void  Penta::MinVz(double* pminvz, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vz_values[numgrids];
-        double  minvz;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute minimum:*/
-        minvz=vz_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vz_values[i]<minvz)minvz=vz_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pminvz=minvz;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxVz(double* pmaxvz, bool process_units);{{{1*/
-void  Penta::MaxVz(double* pmaxvz, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vz_values[numgrids];
-        double  maxvz;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute maximum:*/
-        maxvz=vz_values[0];
-        for(i=1;i<numgrids;i++){
-                if (vz_values[i]>maxvz)maxvz=vz_values[i];
-        }
-
-        /*Assign output pointers:*/
-        *pmaxvz=maxvz;
-
-}
-/*}}}*/
-/*FUNCTION Penta::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/
-void  Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){
-
-        int i;
-        int dim;
-        const int numgrids=6;
-        double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-        double  vz_values[numgrids];
-        double  maxabsvz;
-
-        /*retrieve dim parameter: */
-        parameters->FindParam(&dim,DimEnum);
-
-        /*retrive velocity values at nodes */
-        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
-
-        /*now, compute maximum:*/
-        maxabsvz=fabs(vz_values[0]);
-        for(i=1;i<numgrids;i++){
-                if (fabs(vz_values[i])>maxabsvz)maxabsvz=fabs(vz_values[i]);
-        }
-
-        /*Assign output pointers:*/
-        *pmaxabsvz=maxabsvz;
-}
-/*}}}*/
-/*FUNCTION Penta::InputDuplicate(int original_enum,int new_enum){{{1*/
-void  Penta::InputDuplicate(int original_enum,int new_enum){
-
-        Input* original=NULL;
-        Input* copy=NULL;
-
-        /*Make a copy of the original input: */
-        original=(Input*)this->inputs->GetInput(original_enum);
-        copy=(Input*)original->copy();
-
-        /*Change copy enum to reinitialized_enum: */
-        copy->ChangeEnum(new_enum);
-
-        /*Add copy into inputs, it will wipe off the one already there: */
-        inputs->AddObject((Input*)copy);
-}
-/*}}}*/
-/*FUNCTION Penta::InputScale(int enum_type,double scale_factor){{{1*/
-void  Penta::InputScale(int enum_type,double scale_factor){
-
-        Input* input=NULL;
-
-        /*Make a copy of the original input: */
-        input=(Input*)this->inputs->GetInput(enum_type);
-
-        /*Scale: */
-        input->Scale(scale_factor);
-}
-/*}}}*/
-/*FUNCTION Penta::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/
-void  Penta::InputAXPY(int YEnum, double scalar, int XEnum){
-
-        Input* xinput=NULL;
-        Input* yinput=NULL;
-
-        /*Find x and y inputs: */
-        xinput=(Input*)this->inputs->GetInput(XEnum);
-        yinput=(Input*)this->inputs->GetInput(YEnum);
-
-        /*some checks: */
-        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
-        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
-
-        /*Scale: */
-        yinput->AXPY(xinput,scalar);
-}
-/*}}}*/
-/*FUNCTION Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
-void  Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){
-
-        Input* input=NULL;
-
-        /*Find input: */
-        input=(Input*)this->inputs->GetInput(control_type);
-        
-        /*Do nothing if we  don't find it: */
-        if(!input)return;
-
-        /*Constrain input using cm_min and cm_max: */
-        input->Constrain(cm_min,cm_max);
-
-}
-/*}}}*/
-/*FUNCTION Penta::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
-void  Penta::GetVectorFromInputs(Vec vector,int NameEnum){
-
-        int i;
-        const int numvertices=6;
-        int doflist1[numvertices];
-
-        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
-        for(i=0;i<this->inputs->Size();i++){
-                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
-                if(input->EnumType()==NameEnum){
-                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
-                        this->GetDofList1(&doflist1[0]);
-                        input->GetVectorFromInputs(vector,&doflist1[0]);
-                        break;
-                }
-        }
-}
-/*}}}*/
-/*FUNCTION Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
-void  Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
-
-        int     i;
-        Input** new_inputs=NULL;
-        Input** old_inputs=NULL;
-        int     converged=1;
-
-        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
-        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
-        
-        for(i=0;i<num_enums/2;i++){
-                new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
-                old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
-                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
-                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
-        }
-
-        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
-        for(i=0;i<num_criterionenums;i++){
-                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
-                if(eps[i]>criterionvalues[i]) converged=0; 
-        }
-
-        /*Assign output pointers:*/
-        *pconverged=converged;
-
-}
-/*}}}*/

issm/trunk/src/c/objects/Elements/Penta.h ¶

@@ -62 +62 @@
                 void  InputUpdateFromConstant(int constant, int name);
                 void  InputUpdateFromSolution(double* solutiong);
-                void  InputUpdateFromSolutionBalancedthickness( double* solutiong);
-                void  InputUpdateFromSolutionBalancedthickness2( double* solutiong);
-                void  InputUpdateFromSolutionBalancedvelocities( double* solutiong);
-                void  InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
-                void  InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
-                void  InputUpdateFromSolutionPrognostic( double* solutiong);
-                void  InputUpdateFromSolutionPrognostic2(double* solutiong);
-                void  InputUpdateFromSolutionSlopeCompute( double* solutiong);
                 void  InputUpdateFromVector(bool* vector, int name, int type);
                 void  InputUpdateFromVector(double* vector, int name, int type);
                 void  InputUpdateFromVector(int* vector, int name, int type);
-                void  UpdateFromDakota(void* inputs);
                 /*}}}*/
                 /*Element virtual functions definitions: {{{1*/
@@ -178 +169 @@
                 Penta*    GetUpperElement(void);
                 void      InputExtrude(int enum_type);
+                void      InputUpdateFromSolutionBalancedthickness( double* solutiong);
+                void      InputUpdateFromSolutionBalancedthickness2( double* solutiong);
+                void      InputUpdateFromSolutionBalancedvelocities( double* solutiong);
+                void      InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
+                void      InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
+                void      InputUpdateFromSolutionPrognostic( double* solutiong);
+                void      InputUpdateFromSolutionPrognostic2(double* solutiong);
+                void      InputUpdateFromSolutionSlopeCompute( double* solutiong);
                 bool      IsInput(int name);
                 bool      IsOnSurface(void);
@@ -187 +186 @@
                 void*     SpawnTria(int g0, int g1, int g2);
                 void      SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
+                void      UpdateFromDakota(void* inputs);
                 void      VecExtrude(Vec vector,double* vector_serial,int iscollapsed);
                 /*}}}*/