Changeset 3620

Timestamp:

04/26/10 16:46:08 (15 years ago)

Author:

Eric.Larour

Message:

Big rewrite of elements

Location:

issm/trunk/src/c/objects

Files:

• Beam.cpp (modified) (28 diffs)
• Beam.h (modified) (4 diffs)
• BeamVertexInput.cpp (added)
• BeamVertexInput.h (added)
• Element.h (modified) (2 diffs)
• FemModel.cpp (modified) (2 diffs)
• FemModel.h (modified) (4 diffs)
• Model.cpp (modified) (2 diffs)
• Penta.cpp (modified) (119 diffs)
• Penta.h (modified) (5 diffs)
• PentaVertexInput.cpp (added)
• PentaVertexInput.h (added)
• Sing.cpp (modified) (24 diffs)
• Sing.h (modified) (4 diffs)
• SingVertexInput.cpp (added)
• SingVertexInput.h (added)
• Tria.cpp (modified) (128 diffs)
• Tria.h (modified) (2 diffs)

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

@@ -11 +11 @@
 #include "stdio.h"
 #include "./Beam.h"
+#include "./BeamVertexInput.h"
 #include <string.h>
 #include "../EnumDefinitions/EnumDefinitions.h"
@@ -21 +22 @@
 /*FUNCTION Beam::Beam(){{{1*/
 Beam::Beam(){
+        this->inputs=NULL;
         return;
 }
 /*}}}*/
 /*FUNCTION Beam::Beam(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id, ElementProperties* properties){{{1*/
-Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id, ElementProperties* beamproperties): 
+Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id):
         hnodes(beam_node_ids,2),
         hmatice(&beam_matice_id,1),
         hmatpar(&beam_matpar_id,1),
-        hnumpar(&beam_numpar_id,1),
-        properties(beamproperties)
+        hnumpar(&beam_numpar_id,1)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=beam_id;
-
-        return;
+        this->inputs=new Inputs();
 }
 /*}}}*/
 /*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, ElementProperties* properties){{{1*/
-Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, ElementProperties* beam_properties):
+Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs):
         hnodes(beam_hnodes),
         hmatice(beam_hmatice),
         hmatpar(beam_hmatpar),
-        hnumpar(beam_hnumpar),
-        properties(beam_properties)
+        hnumpar(beam_hnumpar)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=beam_id;
-
+        if(beam_inputs){
+                this->inputs=(Inputs*)beam_inputs->Copy();
+        }
+        else{
+                this->inputs=new Inputs();
+        }
         return;
 }
 /*}}}*/
 /*FUNCTION Beam::Beam(int i, IoModel* iomodel){{{1*/
-Beam::Beam(int i,IoModel* iomodel){
-
+Beam::Beam(int index,IoModel* iomodel){
+
+
+        int i;
 
         /*beam constructor input: */
@@ -64 +70 @@
         int   beam_numpar_id;
         int   beam_node_ids[2];
-        double beam_h[2];
-        double beam_s[2];
-        double beam_b[2];
-        double beam_k[2];
-        int    beam_onbed;
-        
-        
+        double nodeinputs[2];
+
         /*id: */
-        beam_id=i+1; 
+        beam_id=index+1; 
         this->id=beam_id;
 
         /*hooks: */
-        beam_matice_id=i+1; //refers to the corresponding material property card
+        beam_matice_id=index+1; //refers to the corresponding material property card
         beam_matpar_id=iomodel->numberofvertices2d*(iomodel->numlayers-1)+1;//refers to the corresponding matpar property card
         beam_numpar_id=1;
-        beam_node_ids[0]=i+1;
-        beam_node_ids[1]=(int)iomodel->uppernodes[i]; //grid that lays right on top
+        beam_node_ids[0]=index+1;
+        beam_node_ids[1]=(int)iomodel->uppernodes[index]; //grid that lays right on top
         
         this->hnodes.Init(beam_node_ids,2);
@@ -87 +88 @@
         this->hnumpar.Init(&beam_numpar_id,1);
 
-        /*properties: */
-        beam_h[0]=iomodel->thickness[i];
-        beam_h[1]=iomodel->thickness[(int)(iomodel->uppernodes[i]-1)];
-        beam_s[0]=iomodel->surface[i];
-        beam_s[1]=iomodel->surface[(int)(iomodel->uppernodes[i]-1)];
-        beam_b[0]=iomodel->bed[i];
-        beam_b[1]=iomodel->bed[(int)(iomodel->uppernodes[i]-1)];
-        beam_k[0]=iomodel->drag[i];
-        beam_k[1]=iomodel->drag[(int)(iomodel->uppernodes[i]-1)];
-        beam_onbed=(bool)iomodel->gridonbed[i];
-
-        this->properties.Init(2,beam_h, beam_s, beam_b, beam_k,NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, beam_onbed, UNDEF, UNDEF, UNDEF, UNDEF);
+        //intialize inputs, and add as many inputs per element as requested: 
+        this->inputs=new Inputs();
+
+        if (iomodel->thickness) {
+                nodeinputs[0]=iomodel->thickness[index];
+                nodeinputs[1]=iomodel->thickness[(int)(iomodel->uppernodes[index]-1)];
+                this->inputs->AddInput(new BeamVertexInput(ThicknessEnum,nodeinputs));
+        }
+        if (iomodel->surface) {
+                nodeinputs[0]=iomodel->surface[index];
+                nodeinputs[1]=iomodel->surface[(int)(iomodel->uppernodes[index]-1)];
+                this->inputs->AddInput(new BeamVertexInput(SurfaceEnum,nodeinputs));
+        }       
+        if (iomodel->bed) {
+                nodeinputs[0]=iomodel->bed[index];
+                nodeinputs[1]=iomodel->bed[(int)(iomodel->uppernodes[index]-1)];
+                this->inputs->AddInput(new BeamVertexInput(BedEnum,nodeinputs));
+        }       
+        if (iomodel->drag_coefficient) {
+                nodeinputs[0]=iomodel->drag_coefficient[index];
+                nodeinputs[1]=iomodel->drag_coefficient[(int)(iomodel->uppernodes[index]-1)];
+                this->inputs->AddInput(new BeamVertexInput(DragCoefficientEnum,nodeinputs));
+        }       
+        if (iomodel->gridonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->gridonbed[index]));
+
 }
 /*}}}*/
 /*FUNCTION Beam::~Beam(){{{1*/
 Beam::~Beam(){
+        delete inputs;
         return;
 }
@@ -136 +151 @@
 Object* Beam::copy() {
         
-        return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties);
+        return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs);
 
 }
@@ -149 +164 @@
         hmatpar.DeepEcho();
         hnumpar.DeepEcho();
-        properties.DeepEcho();
+        printf("   inputs\n");
+        inputs->DeepEcho();
 
         return;
@@ -173 +189 @@
         hnumpar.Demarshall(&marshalled_dataset);
 
-        /*demarshall properties: */
-        properties.Demarshall(&marshalled_dataset);
+        /*demarshall inputs: */
+        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
 
         /*return: */
@@ -190 +206 @@
         hmatpar.Echo();
         hnumpar.Echo();
-        properties.Echo();
+        printf("   inputs\n");
+        inputs->Echo();
 
         return;
@@ -200 +217 @@
         char* marshalled_dataset=NULL;
         int   enum_type=0;
+        char* marshalled_inputs=NULL;
+        int   marshalled_inputs_size;
 
         /*recover marshalled_dataset: */
@@ -219 +238 @@
         hnumpar.Marshall(&marshalled_dataset);
 
-        /*Marshall properties: */
-        properties.Marshall(&marshalled_dataset);
+        /*Marshall inputs: */
+        marshalled_inputs_size=inputs->MarshallSize();
+        marshalled_inputs=inputs->Marshall();
+        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
+        marshalled_dataset+=marshalled_inputs_size;
+
+        xfree((void**)&marshalled_inputs);
 
         *pmarshalled_dataset=marshalled_dataset;
@@ -234 +258 @@
                 +hmatpar.MarshallSize()
                 +hnumpar.MarshallSize()
-                +properties.MarshallSize()
+                +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
 }
@@ -259 +283 @@
         int doflist[numgrids];
         double pressure[numgrids];
+        double surface[numgrids];
         double rho_ice,g;
         double xyz_list[numgrids][3];
+        double gauss[numgrids][numgrids]={{1,0},{0,1}};
 
         /*dynamic objects pointed to by hooks: */
@@ -286 +312 @@
         rho_ice=matpar->GetRhoIce();
         g=matpar->GetG();
+
+        /*recover value of surface at gauss points (0,1) and (1,0): */
+        inputs->GetParameterValues(&surface[0],&gauss[0][0],2,SurfaceEnum);
         for(i=0;i<numgrids;i++){
-                pressure[i]=rho_ice*g*(this->properties.s[i]-xyz_list[i][2]);
-        }
-        
+                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+        }
+
         /*plug local pressure values into global pressure vector: */
         VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
@@ -303 +332 @@
 /*}}}*/
 /*FUNCTION Beam::CostFunction{{{1*/
-double Beam::CostFunction(void*,int,int){
+double Beam::CostFunction(int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -316 +345 @@
                 if (sub_analysis_type==HutterAnalysisEnum) {
 
-                        CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
+                        CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type);
                 }
                 else 
@@ -338 +367 @@
         double    Ke_gg[numdofs][numdofs]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
         int       numberofdofspernode;
-        
+        bool onbed;
+        
+        
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         GetDofList(&doflist[0],&numberofdofspernode);
 
-        if (this->properties.onbed){
+        if (onbed){
                 Ke_gg[0][0]=1;
                 Ke_gg[1][1]=1;
@@ -367 +400 @@
         if (analysis_type==DiagnosticAnalysisEnum) {
                 if (sub_analysis_type==HutterAnalysisEnum) {
-                        CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
+                        CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type);
                 }
                 else
@@ -403 +436 @@
         double*  gauss_weights=NULL;
         double   gauss_weight;
+        double   gauss1[2]={1,0};
         int      ig;
         double   l1l2[2];
@@ -412 +446 @@
         double   Jdet;
         double   ub,vb;
+        double   surface,thickness;
+        
+        bool onbed;
 
         /*dynamic objects pointed to by hooks: */
@@ -419 +456 @@
         Numpar* numpar=NULL;
 
-        ParameterInputs* inputs=NULL;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-
         /*recover objects from hooks: */
         nodes=(Node**)hnodes.deliverp();
@@ -433 +465 @@
         GetDofList(&doflist[0],&numberofdofspernode);
 
+        /*recover some inputs: */
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*recover parameters: */
         rho_ice=matpar->GetRhoIce();
@@ -440 +475 @@
 
         //recover extra inputs
-        found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,1,(void**)&nodes[0]); //only recover from fist node, second node will have the same slope
-        if(!found)ISSMERROR(" surfaceslopex missing from inputs!");
-        
-        found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,1,(void**)&nodes[0]);//only recover from fist node, second node will have the same slope
-        if(!found)ISSMERROR(" surfaceslopey missing from inputs!");
+        inputs->GetParameterValue(&slope[0],&gauss1[0],SurfaceSlopexEnum);
+        inputs->GetParameterValue(&slope[1],&gauss1[0],SurfaceSlopeyEnum);
+        inputs->GetParameterValue(&surface,&gauss1[0],SurfaceEnum);
+        inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum);
 
         //Get all element grid data:
@@ -474 +508 @@
                 
                 for(j=0;j<NDOF2;j++){
-                        pe_g_gaussian[NDOF2+j]=constant_part*pow((this->properties.s[0]-z_g)/B,n)*slope[j]*Jdet*gauss_weight;
+                        pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight;
                 }
                 
@@ -484 +518 @@
 
         //Deal with lower surface
-        if (this->properties.onbed){
+        if (onbed){
 
                 //compute ub
-                constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*this->properties.h[0];
+                constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
                 ub=constant_part*slope[0];
                 vb=constant_part*slope[1];
@@ -505 +539 @@
 /*}}}*/
 /*FUNCTION Beam::Du{{{1*/
-void  Beam::Du(_p_Vec*,void*,int,int){
+void  Beam::Du(Vec,int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -654 +688 @@
 /*}}}*/
 /*FUNCTION Beam::Gradj{{{1*/
-void  Beam::Gradj(_p_Vec*, void*, int, int,char*){
+void  Beam::Gradj(Vec, int, int,char*){
         ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Beam::GradjB{{{1*/
-void  Beam::GradjB(_p_Vec*, void*, int, int){
+void  Beam::GradjB(Vec, int, int){
         ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Beam::GradjDrag{{{1*/
-void  Beam::GradjDrag(_p_Vec*, void*, int,int ){
+void  Beam::GradjDrag(Vec, int,int ){
         ISSMERROR(" not supported yet!");
 }
@@ -674 +708 @@
 /*}}}*/
 /*FUNCTION Beam::Misfit{{{1*/
-double Beam::Misfit(void*,int,int){
+double Beam::Misfit(int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -685 +719 @@
 /*}}}*/
 /*FUNCTION Beam::SurfaceArea{{{1*/
-double Beam::SurfaceArea(void*,int,int){
+double Beam::SurfaceArea(int,int){
         ISSMERROR(" not supported yet!");
 }

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

@@ -12 +12 @@
 #include "./Matpar.h"
 #include "../ModelProcessorx/IoModel.h"
+#include "../DataSet/Inputs.h"
 #include "./Hook.h"
 
@@ -30 +31 @@
                 Hook hmatpar; //hook to 1 matpar
                 Hook hnumpar; //hook to 1 numpar
-
                 Inputs* inputs;
         
@@ -52 +52 @@
                 int   GetId(); 
                 int   MyRank();
-                void  Configure(void* loads,void* nodes,void* materials,void* parameters);
+                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
                 Object* copy();
                 void  SetClone(int* minranks);
@@ -78 +78 @@
                 void  GetBedList(double*);
                 void  GetThicknessList(double* thickness_list);
-                void  Du(_p_Vec*,void*, int,int);
-                void  Gradj(_p_Vec*, void*, int, int,char*);
-                void  GradjDrag(_p_Vec*, void*, int,int );
-                void  GradjB(_p_Vec*, void*, int,int );
-                double Misfit(void*,int,int);
-                double SurfaceArea(void*,int,int);
-                double CostFunction(void*,int,int);
+                void  Du(Vec, int,int);
+                void  Gradj(Vec,  int, int,char*);
+                void  GradjDrag(Vec,  int,int );
+                void  GradjB(Vec,  int,int );
+                double Misfit(int,int);
+                double SurfaceArea(int,int);
+                double CostFunction(int,int);
                 void  GetNodalFunctions(double* l1l2, double gauss_coord);
                 void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);

issm/trunk/src/c/objects/Element.h

@@ -10 +10 @@
 
 #include "./Object.h"
+#include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 #include "../toolkits/toolkits.h"
 
@@ -18 +20 @@
                 virtual        ~Element(){};
                 virtual int    Enum()=0;
-                virtual void   Configure(void* loads,void* nodes,void* materials,void* parameters)=0;
+                virtual void   Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters)=0;
                 
                 virtual void   CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type)=0;

issm/trunk/src/c/objects/FemModel.cpp

@@ -39 +39 @@
 /*FUNCTION FemModel::FemModel {{{1*/
 FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_vertices, DataSet* femmodel_constraints,DataSet* femmodel_loads,
-                DataSet* femmodel_materials,DataSet* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg,
+                DataSet* femmodel_materials,Parameters* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg,
                 Mat femmodel_Rmg,Mat femmodel_Gmn,NodeSets* femmodel_nodesets,Vec femmodel_ys,Vec femmodel_ys0){
 
@@ -266 +266 @@
 /*}}}*/
 /*FUNCTION FemModel::get_parameters {{{1*/
-DataSet* FemModel::get_parameters(void){
+Parameters* FemModel::get_parameters(void){
         return parameters;
 }

issm/trunk/src/c/objects/FemModel.h

@@ -11 +11 @@
 #include "./Object.h"
 #include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 #include "./DofVec.h"
 #include "../toolkits/toolkits.h"
@@ -27 +28 @@
                 DataSet*            loads;
                 DataSet*            materials;
-                DataSet*            parameters;
+                Parameters*            parameters;
 
                 DofVec*             partition;
@@ -41 +42 @@
                 FemModel();
                 ~FemModel();
-                FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters,
+                FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,Parameters* parameters,
                                       DofVec* partition,DofVec* tpartition,DofVec* yg,Mat Rmg,Mat Gmn,NodeSets* nodesets,Vec ys,Vec ys0);
      
@@ -69 +70 @@
                 DataSet* get_loads(void);
                 DataSet* get_materials(void);
-                DataSet* get_parameters(void);
+                Parameters* get_parameters(void);
                 DofVec*      get_partition(void);
                 DofVec*      get_tpartition(void);

issm/trunk/src/c/objects/Model.cpp

@@ -47 +47 @@
         DataSet*            loads=NULL;
         DataSet*            materials=NULL;
-        DataSet*            parameters=NULL;
+        Parameters*            parameters=NULL;
         DofVec*                 partition=NULL;
         DofVec*                 tpartition=NULL;
@@ -121 +121 @@
         DataSet*            loads=NULL;
         DataSet*            materials=NULL;
-        DataSet*            parameters=NULL;
+        Parameters*            parameters=NULL;
         DofVec*                 partition=NULL;
         DofVec*                 tpartition=NULL;

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

@@ -10 +10 @@
 
 #include "stdio.h"
+#include "./PentaVertexInput.h"
 #include "./Penta.h"
 #include <string.h>
@@ -20 +21 @@
 /*FUNCTION Penta default constructor {{{1*/
 Penta::Penta(){
+        this->inputs=NULL;
         return;
 }
 /*}}}*/
 /*FUNCTION Penta constructor {{{1*/
-Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id, ElementProperties* pentaproperties): 
+Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id):
         hnodes(penta_node_ids,6),
         hmatice(&penta_matice_id,1),
         hmatpar(&penta_matpar_id,1),
-        hnumpar(&penta_numpar_id,1),
-        properties(pentaproperties)
+        hnumpar(&penta_numpar_id,1)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=penta_id;
-
-        return;
+        this->inputs=new Inputs();
+
 }
 /*}}}*/
 /*FUNCTION Penta other constructor {{{1*/
-Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties):
+Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* penta_inputs):
         hnodes(penta_hnodes),
         hmatice(penta_hmatice),
         hmatpar(penta_hmatpar),
-        hnumpar(penta_hnumpar),
-        properties(penta_properties)
+        hnumpar(penta_hnumpar)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=penta_id;
-
-        return;
+        if(penta_inputs){
+                this->inputs=(Inputs*)penta_inputs->Copy();
+        }
+        else{
+                this->inputs=new Inputs();
+        }
 }
 /*}}}*/
 /*FUNCTION Penta other constructor {{{1*/
-Penta::Penta(int i, IoModel* iomodel){ //i is the element index
-
-        int j;
+Penta::Penta(int index, IoModel* iomodel){ //i is the element index
+
+        IssmInt i;
         int penta_node_ids[6];
         int penta_matice_id;
         int penta_matpar_id;
         int penta_numpar_id;
+        double nodeinputs[6];
+        bool collapse;
 
         /*id: */
-        this->id=i+1;
+        this->id=index+1;
 
         /*hooks: */
-        for(j=0;j<6;j++){ //go recover node ids, needed to initialize the node hook.
-                penta_node_ids[j]=(int)*(iomodel->elements+6*i+j); //ids for vertices are in the elements array from Matlab
-        }
-        penta_matice_id=i+1; //refers to the corresponding ice material object
+        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
+        }
+        penta_matice_id=index+1; //refers to the corresponding ice material object
         penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
         penta_numpar_id=1; //refers to numerical parameters object
 
-        this->hnodes.Init(penta_node_ids,6);
+        this->hnodes.Init(&penta_node_ids[0],6);
         this->hmatice.Init(&penta_matice_id,1);
         this->hmatpar.Init(&penta_matpar_id,1);
         this->hnumpar.Init(&penta_numpar_id,1);
-        this->properties.Init(6,penta_node_ids,this->id,iomodel);
+
+        //intialize inputs, and add as many inputs per element as requested: 
+        this->inputs=new Inputs();
+
+        if (iomodel->thickness) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));
+        }
+        if (iomodel->surface) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_node_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));
+        }
+        if (iomodel->bed) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_node_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_node_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_node_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));
+        }
+        if (iomodel->accumulation_rate) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_node_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));
+        }
+        if (iomodel->geothermalflux) {
+                for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_node_ids[i]-1];
+                this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,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){ 
-                        this->properties.collapse=1;
+                        collapse=1;
                 }
                 else{
-                        this->properties.collapse=0;
+                        collapse=0;
                 }
         }
+        this->inputs->AddInput(new BoolInput(CollapseEnum,(IssmBool)collapse));
 
 }
@@ -93 +141 @@
 /*FUNCTION Penta destructor {{{1*/
 Penta::~Penta(){
+        delete inputs;
         return;
 }
@@ -125 +174 @@
 /*FUNCTION copy {{{1*/
 Object* Penta::copy() {
-        return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties); 
+        return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); 
 }
 /*}}}*/
@@ -148 +197 @@
         hnumpar.Demarshall(&marshalled_dataset);
 
-        /*demarshall properties: */
-        properties.Demarshall(&marshalled_dataset);
+        /*demarshall inputs: */
+        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
 
         /*return: */
@@ -166 +215 @@
         hmatpar.DeepEcho();
         hnumpar.DeepEcho();
-        properties.DeepEcho();
+        printf("   inputs\n");
+        inputs->DeepEcho();
 
         return;
@@ -181 +231 @@
         hmatpar.Echo();
         hnumpar.Echo();
-        properties.Echo();
-
-        return;
+        printf("   inputs\n");
+        inputs->Echo();
+
+}
+/*}}}*/
+/*FUNCTION Enum {{{1*/
+int Penta::Enum(void){
+
+        return PentaEnum;
+
 }
 /*}}}*/
@@ -191 +248 @@
         char* marshalled_dataset=NULL;
         int   enum_type=0;
+        char* marshalled_inputs=NULL;
+        int   marshalled_inputs_size;
 
         /*recover marshalled_dataset: */
@@ -210 +269 @@
         hnumpar.Marshall(&marshalled_dataset);
 
-        /*Marshall properties: */
-        properties.Marshall(&marshalled_dataset);
+        /*Marshall inputs: */
+        marshalled_inputs_size=inputs->MarshallSize();
+        marshalled_inputs=inputs->Marshall();
+        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
+        marshalled_dataset+=marshalled_inputs_size;
+
+        xfree((void**)&marshalled_inputs);
+
 
         *pmarshalled_dataset=marshalled_dataset;
@@ -225 +290 @@
                 +hmatpar.MarshallSize()
                 +hnumpar.MarshallSize()
-                +properties.MarshallSize()
+                +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
 }
@@ -240 +305 @@
         Hook* tria_hmatpar=NULL;
         Hook* tria_hnumpar=NULL;
-        ElementProperties* tria_properties=NULL;
-        DataSet* tria_inputs=NULL;
+        Inputs* tria_inputs=NULL;
 
         indices[0]=g0;
@@ -251 +315 @@
         tria_hmatpar=this->hmatpar.Spawn(&zero,1);
         tria_hnumpar=this->hnumpar.Spawn(&zero,1);
-        tria_properties=(ElementProperties*)this->properties.Spawn(indices,3);
-        tria_inputs=(DataSet*)this->inputs->Spawn(indices,3);
-
-        tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_properties,tria_inputs);
+        tria_inputs=(Inputs*)this->inputs->Spawn(indices,3);
+
+        tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_inputs);
 
         delete tria_hnodes;
@@ -260 +323 @@
         delete tria_hmatpar;
         delete tria_hnumpar;
-        delete tria_properties;
         delete tria_inputs;
 
@@ -266 +328 @@
 }
 /*}}}*/
+
+/*updates: */
 /*FUNCTION UpdateFromDakota {{{1*/
 void  Penta::UpdateFromDakota(void* vinputs){
@@ -275 +339 @@
 /*FUNCTION Penta::UpdateInputs {{{1*/
 void  Penta::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){
+
+        /*Just branch to the correct UpdateInputs generator, according to the type of analysis we are carrying out: */
+        if (analysis_type==ControlAnalysisEnum){
+                
+                UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==DiagnosticAnalysisEnum){
+        
+                if (sub_analysis_type==HorizAnalysisEnum){
+
+                        UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type);
+                }
+                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
+
+        }
+        else if (analysis_type==SlopecomputeAnalysisEnum){
+
+                UpdateInputsSlopeCompute( solution,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==PrognosticAnalysisEnum){
+
+                UpdateInputsPrognostic( solution,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==Prognostic2AnalysisEnum){
+
+                UpdateInputsPrognostic2(solution,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==BalancedthicknessAnalysisEnum){
+
+                UpdateInputsBalancedthickness( solution,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==Balancedthickness2AnalysisEnum){
+
+                UpdateInputsBalancedthickness2( solution,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
+
+                UpdateInputsBalancedvelocities( solution,analysis_type,sub_analysis_type);
+        }
+        else{
+
+                ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet");
+        }
+}
+/*Object functions*/
+/*FUNCTION Penta::UpdateInputsDiagnosticHoriz {{{1*/
+void  Penta::UpdateInputsDiagnosticHoriz(double* solution, int analysis_type, int sub_analysis_type){
+        
+        
+        int i;
+
+        const int    numvertices=6;
+        const int    numdofpervertex=2;
+        const int    numdof=numdofpervertex*numvertices;
+        
+        int          doflist[numdof];
+        double       values[numdof];
+        double       vx[numvertices];
+        double       vy[numvertices];
+
+        int          dummy;
+        
+        /*Get dof list: */
+        GetDofList(&doflist[0],&dummy);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof;i++){
+                values[i]=solution[doflist[i]];
+        }
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<numvertices;i++){
+                vx[i]=values[i*numdofpervertex+0];
+                vy[i]=values[i*numdofpervertex+1];
+        }
+
+        /*Add vx and vy as inputs to the tria element: */
+        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
+        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
+}
+
+/*}}}*/
+/*FUNCTION Penta::UpdateInputsSlopeCompute {{{1*/
+void  Penta::UpdateInputsSlopeCompute(double* solution, int analysis_type, int sub_analysis_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateInputsPrognostic {{{1*/
+void  Penta::UpdateInputsPrognostic(double* solution, int analysis_type, int sub_analysis_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateInputsPrognostic2 {{{1*/
+void  Penta::UpdateInputsPrognostic2(double* solution, int analysis_type, int sub_analysis_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateInputsBalancedthickness {{{1*/
+void  Penta::UpdateInputsBalancedthickness(double* solution, int analysis_type, int sub_analysis_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateInputsBalancedthickness2 {{{1*/
+void  Penta::UpdateInputsBalancedthickness2(double* solution, int analysis_type, int sub_analysis_type){
+        ISSMERROR(" not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateInputsBalancedvelocities {{{1*/
+void  Penta::UpdateInputsBalancedvelocities(double* solution, int analysis_type, int sub_analysis_type){
         ISSMERROR(" not supported yet!");
 }
@@ -295 +468 @@
         double bed;
         double basalforce[3];
-        double vxvyvz_list[numgrids][3];
-        double pressure_list[numgrids];
         double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -304 +475 @@
         int  dofv[3]={0,1,2};
         int  dofp[1]={3};
-        ParameterInputs* inputs=NULL;
         double Jdet2d;
         Tria* tria=NULL;
@@ -323 +493 @@
         double surface=0;
         double value=0;
+        
+        /*flags: */
+        bool onbed;
 
         /*dynamic objects pointed to by hooks: */
@@ -333 +506 @@
         if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!");
 
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-
         /*recover objects from hooks: */
         nodes=(Node**)hnodes.deliverp();
@@ -342 +512 @@
         numpar=(Numpar*)hnumpar.delivers();
 
-        if(!this->properties.onbed){
+        /*recover some inputs: */
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        if(!onbed){
                 //put zero
                 VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES);
@@ -351 +524 @@
         rho_ice=matpar->GetRhoIce();
         gravity=matpar->GetG();
-
-        /*recover extra inputs from users, at current convergence iteration: */
-        inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofv,numgrids,(void**)nodes);
-        inputs->Recover("velocity",&pressure_list[0] ,1,dofp,numgrids,(void**)nodes);
 
         /* Get node coordinates and dof list: */
@@ -378 +547 @@
 
                         /*Compute strain rate viscosity and pressure: */
-                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-                        GetParameterValue(&pressure,&pressure_list[0],&gauss_coord[0]);
+                        inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum);
 
                         /*Compute Stress*/
@@ -417 +586 @@
         int i;
         const int numgrids=6;
-        int doflist[numgrids];
+        int    doflist[numgrids];
         double pressure[numgrids];
         double rho_ice,g;
-        double       xyz_list[numgrids][3];
+        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;
 
         /*dynamic objects pointed to by hooks: */
@@ -430 +604 @@
         matpar=(Matpar*)hmatpar.delivers();
 
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*Get node data: */
@@ -441 +618 @@
         /*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: */
@@ -446 +626 @@
         g=matpar->GetG();
         for(i=0;i<numgrids;i++){
-                pressure[i]=rho_ice*g*(this->properties.s[i]-xyz_list[i][2]);
+                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
         }
 
@@ -467 +647 @@
         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(this->properties.onwater)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 ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
                 return 0;
         }
-        else if (this->properties.collapse){
+        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(inputs,analysis_type,sub_analysis_type);
+                J=tria->CostFunction(analysis_type,sub_analysis_type);
                 delete tria;
                 return J;
@@ -488 +680 @@
 
                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                J=tria->CostFunction(inputs,analysis_type,sub_analysis_type);
+                J=tria->CostFunction(analysis_type,sub_analysis_type);
                 delete tria;
                 return J;
@@ -501 +693 @@
         if (analysis_type==ControlAnalysisEnum){
 
-                CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
+                CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type);
 
         }
@@ -508 +700 @@
                 if (sub_analysis_type==HorizAnalysisEnum){
 
-                        CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
+                        CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==VertAnalysisEnum){
 
-                        CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
+                        CreateKMatrixDiagnosticVert( Kgg,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==StokesAnalysisEnum){
 
-                        CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
+                        CreateKMatrixDiagnosticStokes( Kgg,analysis_type,sub_analysis_type);
 
                 }
@@ -523 +715 @@
         else if (analysis_type==SlopecomputeAnalysisEnum){
 
-                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
+                CreateKMatrixSlopeCompute( Kgg,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==PrognosticAnalysisEnum){
 
-                CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
+                CreateKMatrixPrognostic( Kgg,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==BalancedthicknessAnalysisEnum){
@@ -539 +731 @@
         else if (analysis_type==ThermalAnalysisEnum){
 
-                CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
+                CreateKMatrixThermal( Kgg,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==MeltingAnalysisEnum){
 
-                CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
+                CreateKMatrixMelting( Kgg,analysis_type,sub_analysis_type);
         }
         else{
@@ -645 +837 @@
         double Bprime[5][numdof];
         double L[2][numdof];
-        double D[5][5]={{ 0,0,0,0,0 },{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}};              // material matrix, simple scalar matrix.
+        double D[5][5]={0.0};            // material matrix, simple scalar matrix.
         double D_scalar;
-        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+        double DL[2][2]={0.0}; //for basal drag
         double DL_scalar;
 
         /* local element matrices: */
-        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
+
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
         double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
@@ -657 +850 @@
 
         /*slope: */
-        double  slope[2]={0.0,0.0};
+        double  slope[2]={0.0};
         double  slope_magnitude;
-
-        /*input parameters for structural analysis (diagnostic): */
-        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
-        double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
-        int     dofs[2]={0,1};
-        double  thickness;
 
         /*friction: */
@@ -672 +859 @@
         double MAXSLOPE=.06; // 6 %
         double MOUNTAINKEXPONENT=10;
-
-        ParameterInputs* inputs=NULL;
 
         /*Collapsed formulation: */
@@ -688 +873 @@
         numpar=(Numpar*)hnumpar.delivers();
 
+        /*inputs: */
+        bool onwater;
+        bool collapse;
+        bool onbed;
+        bool shelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
         /*If on water, skip stiffness: */
-        if(this->properties.onwater)return;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
+        if(onwater)return;
 
         /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 
@@ -699 +893 @@
 
 
-        if ((this->properties.collapse==1) && (this->properties.onbed==0)){
+        if ((collapse==1) && (onbed==0)){
                 /*This element should be collapsed, but this element is not on the bedrock, therefore all its 
                  * dofs have already been frozen! Do nothing: */
                 return;
         }
-        else if ((this->properties.collapse==1) && (this->properties.onbed==1)){
+        else if ((collapse==1) && (onbed==1)){
 
                 /*This element should be collapsed into a tria element at its base. Create this tria element, 
                  *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
+                tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -716 +910 @@
 
                 /*Implement standard penta element: */
-
-                /*recover extra inputs from users, at current convergence iteration: */
-                inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
-                inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
 
                 /* Get node coordinates and dof list: */
@@ -725 +915 @@
                 GetDofList(&doflist[0],&numberofdofspernode);
 
-                /* Set Ke_gg to 0: */
-                for(i=0;i<numdof;i++){
-                        for(j=0;j<numdof;j++){
-                                Ke_gg[i][j]=0.0;
-                        }
-                }
 
                 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -760 +944 @@
 
                                 /*Get strain rate from velocity: */
-                                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_coord);
-                                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_coord);
+                                inputs->GetStrainRate(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum);
+                                inputs->GetStrainRate(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum);
 
                                 /*Get viscosity: */
@@ -774 +958 @@
                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
 
-                                /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
+                                /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 
                                   onto this scalar matrix, so that we win some computational time: */
 
@@ -803 +987 @@
 
                 //Deal with 2d friction at the bedrock interface
-                if((this->properties.onbed && !this->properties.shelf)){
+                if((onbed && !shelf)){
 
                         /*Build a tria element using the 3 grids of the base of the penta. Then use 
@@ -810 +994 @@
 
                         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                        tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
+                        tria->CreateKMatrixDiagnosticHorizFriction(Kgg,analysis_type,sub_analysis_type);
                         delete tria;
                 }
@@ -845 +1029 @@
         int   dofs[3]={0,1,2};
 
-        double K_terms[numdof][numdof];
+        double K_terms[numdof][numdof]={0.0};
 
         /*Material properties: */
@@ -860 +1044 @@
         double             surface_normal[3];
         double             bed_normal[3];
-        double         vxvyvz_list[numgrids][3];
         double         thickness;
 
@@ -889 +1072 @@
         double* area_gauss_weights  =  NULL;
         double* vert_gauss_weights  =  NULL;
+        double  gaussgrids[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}};
 
         /* specific gaussian point: */
@@ -902 +1086 @@
         double  alpha2_gauss;
 
-        ParameterInputs* inputs=NULL;
+        double  vx_list[numgrids];
+        double  vy_list[numgrids];
+        double  vz_list[numgrids];
+        double  thickness_list[numgrids];
+        double  bed_list[numgrids];
+        double  dragcoefficient_list[numgrids];
+        double  drag_p,drag_q;
 
         /*dynamic objects pointed to by hooks: */
@@ -910 +1100 @@
         Numpar* numpar=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
         /*If on water, skip stiffness: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*recover objects from hooks: */
@@ -919 +1119 @@
         numpar=(Numpar*)hnumpar.delivers();
 
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-
-        /* Set K_terms  to 0: */
-        for(i=0;i<numdof;i++){
-                for(j=0;j<numdof;j++){
-                        K_terms[i][j]=0.0;
-                }
-        }
 
         /*recovre material parameters: */
@@ -937 +1128 @@
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /*recover extra inputs from users, at current convergence iteration: */
-        inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -965 +1153 @@
                         gauss_coord[3]=*(vert_gauss_coord+igvert);
 
-                        /*Compute thickness: */
-
                         /*Compute strain rate: */
-                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
 
                         /*Get viscosity: */
@@ -1003 +1189 @@
         }
 
-        if((this->properties.onbed==1) && (this->properties.shelf==0)){
+        if((onbed==1) && (shelf==0)){
 
                 /*Build alpha2_list used by drag stiffness matrix*/
@@ -1012 +1198 @@
                 strcpy(friction->element_type,"2d");
 
+                inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],6,VxEnum);
+                inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],6,VyEnum);
+                inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],6,VzEnum);
+                inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],6,DragCoefficientEnum);
+                inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],6,BedEnum);
+                inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],6,ThicknessEnum);
+                inputs->GetParameterValue(&drag_p,DragPEnum);
+                inputs->GetParameterValue(&drag_q,DragQEnum);
+
                 friction->gravity=matpar->GetG();
                 friction->rho_ice=matpar->GetRhoIce();
                 friction->rho_water=matpar->GetRhoWater();
-                friction->K=&this->properties.k[0];
-                friction->bed=&this->properties.b[0];
-                friction->thickness=&this->properties.h[0];
-                friction->velocities=&vxvyvz_list[0][0];
-                friction->p=this->properties.p;
-                friction->q=this->properties.q;
-                friction->analysis_type=analysis_type;
+                friction->K=&dragcoefficient_list[0];
+                friction->bed=&bed_list[0];
+                friction->thickness=&thickness_list[0];
+                friction->vx=&vx_list[0];
+                friction->vy=&vy_list[0];
+                friction->vz=&vz_list[0];
+                friction->p=drag_p;
+                friction->q=drag_q;
 
                 /*Compute alpha2_list: */
@@ -1058 +1254 @@
 
                         /*Compute strain rate: */
-                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
 
                         /*Get viscosity at last iteration: */
@@ -1098 +1294 @@
                         }
                 }
-        } //if ( (this->properties.onbed==1) && (shelf==0))
+        } //if ( (onbed==1) && (shelf==0))
 
         /*Reduce the matrix */
@@ -1155 +1351 @@
 
         /* matrices: */
-        double  Ke_gg[numdof][numdof];
+        double  Ke_gg[numdof][numdof]={0.0};
         double  Ke_gg_gaussian[numdof][numdof];
         double  Jdet;
@@ -1168 +1364 @@
         nodes=(Node**)hnodes.deliverp();
 
-        ParameterInputs* inputs=NULL;
-
         /*Collapsed formulation: */
         Tria*  tria=NULL;
 
+        /*inputs: */
+        bool onwater; 
+        bool onsurface;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
         /*If on water, skip stiffness: */
-        if(this->properties.onwater)return;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-
+        if(onwater)return;
 
         /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness 
          * matrix: */
-        if(this->properties.onsurface){
+        if(onsurface){
                 tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
-                tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type);
+                tria->CreateKMatrixDiagnosticSurfaceVert(Kgg, analysis_type,sub_analysis_type);
                 delete tria;
         }
@@ -1194 +1392 @@
         GetDofList(&doflist[0],&numberofdofspernode);
 
-        /* Set Ke_gg to 0: */
-        for(i=0;i<numdof;i++){
-                for(j=0;j<numdof;j++){
-                        Ke_gg[i][j]=0.0;
-                }
-        }
 
         /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -1265 +1457 @@
         Tria* tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
-
-        if (!this->properties.onbed){
+        if(onwater)return;
+
+        if (!onbed){
                 return;
         }
@@ -1274 +1474 @@
 
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
+                tria->CreateKMatrixMelting(Kgg, analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -1287 +1487 @@
         Tria*  tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*Is this element on the bed? :*/
-        if(!this->properties.onbed)return;
+        if(!onbed)return;
 
         /*Spawn Tria element from the base of the Penta: */
         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-        tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
+        tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type);
         delete tria;
         return;
@@ -1308 +1516 @@
         Tria*  tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*Is this element on the bed? :*/
-        if(!this->properties.onbed)return;
+        if(!onbed)return;
 
         /*Spawn Tria element from the base of the Penta: */
         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-        tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
+        tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type);
         delete tria;
         return;
@@ -1356 +1573 @@
         double  K[2][2]={0.0};
 
-        double  vxvyvz_list[numgrids][3];
-        double  vx_list[numgrids];
-        int     vx_list_exists;
-        double  u;
-        double  vy_list[numgrids];
-        int     vy_list_exists;
-        double  v;
-        double  vz_list[numgrids];
-        int     vz_list_exists;
-        double  w;
+        double  u,v,w;
 
         /*matrices: */
@@ -1404 +1612 @@
         /*Collapsed formulation: */
         Tria*  tria=NULL;
-        ParameterInputs* inputs=NULL;
+
+
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
 
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*recover objects from hooks: */
@@ -1413 +1631 @@
         matpar=(Matpar*)hmatpar.delivers();
         numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
 
         /* Get node coordinates and dof list: */
@@ -1427 +1642 @@
         heatcapacity=matpar->GetHeatCapacity();
         thermalconductivity=matpar->GetThermalConductivity();
-
-        /*recover extra inputs from users, dt and velocity: */
-        found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
-        if(!found)ISSMERROR(" could not find velocity in inputs!");
-        found=inputs->Recover("dt",&dt);
-        if(!found)ISSMERROR(" could not find dt in inputs!");
-
-        for(i=0;i<numgrids;i++){
-                vx_list[i]=vxvyvz_list[i][0];
-                vy_list[i]=vxvyvz_list[i][1];
-                vz_list[i]=vxvyvz_list[i][2];
-        }
-
+        dt=numpar->dt;
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -1495 +1698 @@
 
                         //Build the D matrix
-                        GetParameterValue(&u, &vx_list[0],gauss_coord);
-                        GetParameterValue(&v, &vy_list[0],gauss_coord);
-                        GetParameterValue(&w, &vz_list[0],gauss_coord);
+                        inputs->GetParameterValue(&u, gauss_coord,VxEnum);
+                        inputs->GetParameterValue(&v, gauss_coord,VyEnum);
+                        inputs->GetParameterValue(&w, gauss_coord,VzEnum);
 
                         D_scalar=gauss_weight*Jdet;
@@ -1577 +1780 @@
 
         //Ice/ocean heat exchange flux on ice shelf base 
-        if(this->properties.onbed && this->properties.shelf){
+        if(onbed && shelf){
 
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
+                tria->CreateKMatrixThermal(Kgg, analysis_type,sub_analysis_type);
                 delete tria;
         }
@@ -1591 +1794 @@
         if (analysis_type==ControlAnalysisEnum){
 
-                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+                CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==DiagnosticAnalysisEnum){
@@ -1597 +1800 @@
                 if (sub_analysis_type==HorizAnalysisEnum){
 
-                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+                        CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==VertAnalysisEnum){
 
-                        CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
+                        CreatePVectorDiagnosticVert( pg,analysis_type,sub_analysis_type);
                 }
                 else if (sub_analysis_type==StokesAnalysisEnum){
 
-                        CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
+                        CreatePVectorDiagnosticStokes( pg,analysis_type,sub_analysis_type);
                 }
                 else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
@@ -1611 +1814 @@
         else if (analysis_type==SlopecomputeAnalysisEnum){
 
-                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
+                CreatePVectorSlopeCompute( pg,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==PrognosticAnalysisEnum){
 
+                CreatePVectorPrognostic( pg,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==BalancedthicknessAnalysisEnum){
+
                 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
         }
-        else if (analysis_type==BalancedthicknessAnalysisEnum){
+        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
 
                 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
         }
-        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
-
-                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
-        }
         else if (analysis_type==ThermalAnalysisEnum){
 
-                CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
+                CreatePVectorThermal( pg,analysis_type,sub_analysis_type);
         }
         else if (analysis_type==MeltingAnalysisEnum){
 
-                CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
+                CreatePVectorMelting( pg,analysis_type,sub_analysis_type);
         }
         else{
@@ -1720 +1923 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdof];
+        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
-
-        ParameterInputs* inputs=NULL;
 
         /*dynamic objects pointed to by hooks: */
@@ -1732 +1933 @@
         Tria* tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool collapse;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*If on water, skip load: */
-        if(this->properties.onwater)return;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
+        if(onwater)return;
 
         /*recover objects from hooks: */
@@ -1747 +1955 @@
           the load vector. */
 
-        if ((this->properties.collapse==1) && (this->properties.onbed==0)){
+        if ((collapse==1) && (onbed==0)){
                 /*This element should be collapsed, but this element is not on the bedrock, therefore all its 
                  * dofs have already been frozen! Do nothing: */
                 return;
         }
-        else if ((this->properties.collapse==1) && (this->properties.onbed==1)){
+        else if ((collapse==1) && (onbed==1)){
 
                 /*This element should be collapsed into a tria element at its base. Create this tria element, 
                  *and use its CreatePVector functionality to return an elementary load vector: */
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
+                tria->CreatePVector(pg, analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -1768 +1976 @@
                 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
                 GetDofList(&doflist[0],&numberofdofspernode);
-
-                /* Set pe_g to 0: */
-                for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
                 /*Get gaussian points and weights :*/
@@ -1793 +1998 @@
 
                                 /*Compute thickness at gaussian point: */
-                                GetParameterValue(&thickness, &this->properties.h[0],gauss_coord);
+                                inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum);
 
                                 /*Compute slope at gaussian point: */
-                                GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_coord);
+                                inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum);
 
                                 /* Get Jacobian determinant: */
@@ -1820 +2025 @@
                 } //for (ig1=0; ig1<num_area_gauss; ig1++)
 
-        } //else if ((collapse==1) && (this->properties.onbed==1))
+        } //else if ((collapse==1) && (onbed==1))
 
         /*Add pe_g to global vector pg: */
@@ -1857 +2062 @@
         double             bed_normal[3];
         double         bed;
-        double         vxvyvz_list[numgrids][3];
 
         /* gaussian points: */
@@ -1898 +2102 @@
         double     D_scalar;
         double     tBD[27][8];
-        double     P_terms[numdof];
-
-        ParameterInputs* inputs=NULL;
+        double     P_terms[numdof]={0.0};
+
         Tria*            tria=NULL;
 
@@ -1909 +2112 @@
         Numpar* numpar=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
         /*If on water, skip load: */
-        if(this->properties.onwater)return;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
+        if(onwater)return;
 
         /*recover objects from hooks: */
@@ -1921 +2131 @@
         numpar=(Numpar*)hnumpar.delivers();
 
-        /* Set P_terms to 0: */
-        for(i=0;i<numdof;i++){
-                P_terms[i]=0;
-        }       
 
         /*recovre material parameters: */
@@ -1930 +2136 @@
         rho_ice=matpar->GetRhoIce();
         gravity=matpar->GetG();
-
-        /*recover extra inputs from users, at current convergence iteration: */
-        inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
 
         /* Get node coordinates and dof list: */
@@ -1960 +2163 @@
 
                         /*Compute strain rate and viscosity: */
-                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
@@ -2014 +2217 @@
 
         /*Deal with 2d friction at the bedrock interface: */
-        if ( (this->properties.onbed==1) && (this->properties.shelf==1)){
+        if ( (onbed==1) && (shelf==1)){
 
                 for(i=0;i<numgrids2d;i++){
@@ -2041 +2244 @@
 
                         /* Get bed at gaussian point */
-                        GetParameterValue(&bed,&this->properties.b[0],gauss_coord);
+                        inputs->GetParameterValue(&bed, gauss_coord,BedEnum);
 
                         /*Get L matrix : */
@@ -2062 +2265 @@
                         }
                 }
-        } //if ( (this->properties.onbed==1) && (this->properties.shelf==1))
+        } //if ( (onbed==1) && (shelf==1))
 
         /*Reduce the matrix */
@@ -2118 +2321 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdof];
+        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
         double l1l6[6];
@@ -2126 +2329 @@
 
         /*input parameters for structural analysis (diagnostic): */
-        double vx_list[numgrids]={0,0,0,0,0,0};
-        double vy_list[numgrids]={0,0,0,0,0,0};
         double du[3];
         double dv[3];
@@ -2137 +2338 @@
         Node**  nodes=NULL;
 
-        ParameterInputs* inputs=NULL;
+        /*inputs: */
+        bool onwater;
+        bool onbed;
 
         /*recover objects from hooks: */
         nodes=(Node**)hnodes.deliverp();
 
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the 
          *diagnostic base vertical stifness: */
-        if(this->properties.onbed){
+        if(onbed){
                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
-                tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type,sub_analysis_type);
+                tria->CreatePVectorDiagnosticBaseVert(pg, analysis_type,sub_analysis_type);
                 delete tria;
         }
@@ -2161 +2366 @@
         GetDofList(&doflist[0],&numberofdofspernode);
 
-        /* Set pe_g to 0: */
-        for(i=0;i<numdof;i++) pe_g[i]=0.0;
-
-        /* recover input parameters: */
-        if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))ISSMERROR(" cannot compute vertical velocity without horizontal velocity!");
-        inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
-
         /*Get gaussian points and weights :*/
         order_area_gauss=2;
@@ -2189 +2387 @@
 
                         /*Get velocity derivative, with respect to x and y: */
-                        GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_coord);
-                        GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_coord);
+
+                        inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum);
+                        inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VyEnum);
                         dudx=du[0];
                         dvdy=dv[1];
@@ -2235 +2434 @@
         Tria*  tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*Is this element on the bed? :*/
-        if(!this->properties.onbed)return;
+        if(!onbed)return;
 
         /*Spawn Tria element from the base of the Penta: */
         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-        tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
+        tria->CreatePVector(pg, analysis_type,sub_analysis_type);
         delete tria;
         return;
@@ -2255 +2462 @@
         Tria*  tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*Is this element on the bed? :*/
-        if(!this->properties.onbed)return;
+        if(!onbed)return;
 
         /*Spawn Tria element from the base of the Penta: */
         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-        tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
+        tria->CreatePVector(pg, analysis_type,sub_analysis_type);
         delete tria;
         return;
@@ -2299 +2514 @@
 
         double dt;
-        double vx_list[numgrids];
-        double vy_list[numgrids];
-        double vz_list[numgrids];
-        double vxvyvz_list[numgrids][3];
         double temperature_list[numgrids];
         double temperature;
@@ -2340 +2551 @@
         /*Collapsed formulation: */
         Tria*  tria=NULL;
-        ParameterInputs* inputs=NULL;
 
         /*dynamic objects pointed to by hooks: */
@@ -2346 +2556 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
+        Numpar* numpar=NULL;
+
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
 
         /*If on water, skip: */
-        if(this->properties.onwater)return;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
+        if(onwater)return;
 
         /*recover objects from hooks: */
@@ -2357 +2575 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
+        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -2369 +2588 @@
         beta=matpar->GetBeta();
         meltingpoint=matpar->GetMeltingPoint();
-
-        /*recover extra inputs from users, dt and velocity: */
-        found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
-        if(!found)ISSMERROR(" could not find velocity in inputs!");
-        found=inputs->Recover("dt",&dt);
-        if(!found)ISSMERROR(" could not find dt in inputs!");
-
-        if(dt){
-                found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes);
-                if(!found)ISSMERROR(" could not find temperature in inputs!");
-        }
-
-        for(i=0;i<numgrids;i++){
-                vx_list[i]=vxvyvz_list[i][0];
-                vy_list[i]=vxvyvz_list[i][1];
-                vz_list[i]=vxvyvz_list[i][2];
-        }       
+        dt=numpar->dt;
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -2409 +2612 @@
 
                         /*Compute strain rate and viscosity: */
-                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
+                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
@@ -2433 +2636 @@
                         /* Build transient now */
                         if(dt){
-                                GetParameterValue(&temperature, &temperature_list[0],gauss_coord);
+                                inputs->GetParameterValue(&temperature, gauss_coord,TemperatureEnum);
                                 scalar_transient=temperature*Jdet*gauss_weight;
                                 for(i=0;i<numgrids;i++){
@@ -2446 +2649 @@
 
         /* Ice/ocean heat exchange flux on ice shelf base */
-        if(this->properties.onbed && this->properties.shelf){
+        if(onbed && shelf){
 
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type);
+                tria->CreatePVectorThermalShelf(pg, analysis_type,sub_analysis_type);
                 delete tria;
         }
 
         /* Geothermal flux on ice sheet base and basal friction */
-        if(this->properties.onbed && !this->properties.shelf){
+        if(onbed && !shelf){
 
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type);
+                tria->CreatePVectorThermalSheet(pg, analysis_type,sub_analysis_type);
                 delete tria;
         }
@@ -2477 +2680 @@
         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(this->properties.onwater)return;
+        if(onwater)return;
 
         /*Bail out if this element if:
          * -> Non collapsed and not on the surface
          * -> collapsed (2d model) and not on bed) */
-        if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
                 return;
         }
-        else if (this->properties.collapse){
+        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,inputs,analysis_type,sub_analysis_type);
+                tria->Du(du_g,analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -2498 +2713 @@
 
                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                tria->Du(du_g,inputs,analysis_type,sub_analysis_type);
+                tria->Du(du_g,analysis_type,sub_analysis_type);
                 delete tria;
                 return;
         }
-}
-/*}}}*/
-/*FUNCTION Enum {{{1*/
-int Penta::Enum(void){
-
-        return PentaEnum;
-
 }
 /*}}}*/
@@ -2527 +2735 @@
         nodes=(Node**)hnodes.deliverp();
 
+        /*inputs: */
+        bool collapse;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*Figure out if we should extrude for this element: */
         if (iscollapsed){
                 /*From higher level, we are told to extrude only elements that have the collapse flag on: */
-                if (this->properties.collapse)extrude=1;
+                if (collapse)extrude=1;
                 else extrude=0;
         }
@@ -2540 +2756 @@
         /*Now, extrusion starts from the bed on, so double check this element is on 
          * the bedrock: */
-        if(this->properties.onbed==0)extrude=0;
+        if(onbed==0)extrude=0;
 
         /*Go on and extrude field: */
@@ -2822 +3038 @@
 /*}}}*/
 /*FUNCTION GetBedList {{{1*/
-void Penta::GetBedList(double* bed_list){
-
-        int i;
-        for(i=0;i<6;i++)bed_list[i]=this->properties.b[i];
+void Penta::GetBedList(double* bedlist){
+
+        const int numgrids=6;
+        double  gaussgrids[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->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum);
 
 }
@@ -3682 +3900 @@
 /*FUNCTION GetOnBed {{{1*/
 bool Penta::GetOnBed(){
-        return this->properties.onbed;
+        
+        bool onbed;
+
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        return onbed;
 }
 /*}}}*/
@@ -3757 +3980 @@
 /*}}}*/
 /*FUNCTION GetShelf {{{1*/
-int   Penta::GetShelf(){
-        return this->properties.shelf;
+bool   Penta::GetShelf(){
+        bool onshelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
+
+        return onshelf;
 }
 /*}}}*/
@@ -3821 +4049 @@
 void Penta::GetThicknessList(double* thickness_list){
 
-        int i;
-        for(i=0;i<6;i++)thickness_list[i]=this->properties.h[i];
+        const int numgrids=6;
+        double  gaussgrids[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->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum);
+
 }
 /*}}}*/
@@ -3828 +4058 @@
 void  Penta::Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type){
 
+        /*inputs: */
+        bool onwater;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+
         /*If on water, skip grad (=0): */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         if (strcmp(control_type,"drag")==0){
-                GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type);
+                GradjDrag( grad_g,analysis_type,sub_analysis_type);
         }
         else if (strcmp(control_type,"B")==0){
-                GradjB( grad_g, inputs,analysis_type,sub_analysis_type);
+                GradjB( grad_g, analysis_type,sub_analysis_type);
         }
         else ISSMERROR("%s%s","control type not supported yet: ",control_type);
@@ -3845 +4081 @@
         Tria* tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool onbed;
+        bool shelf;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
+        if(onwater)return;
 
         /*If on shelf, skip: */
-        if(this->properties.shelf)return;
+        if(shelf)return;
 
         /*Bail out if this element does not touch the bedrock: */
-        if (!this->properties.onbed) return;
+        if (!onbed) return;
 
         if (sub_analysis_type==HorizAnalysisEnum){
@@ -3858 +4104 @@
                 /*MacAyeal or Pattyn*/
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type);
+                tria->GradjDrag( grad_g,analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -3866 +4112 @@
                 /*Stokes*/
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type);
+                tria->GradjDragStokes( grad_g,analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -3878 +4124 @@
         Tria* tria=NULL;
 
+        /*inputs: */
+        bool onwater;
+        bool collapse;
+        bool onbed;
+
+        /*retrieve inputs :*/
+        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
+        inputs->GetParameterValue(&collapse,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
         /*If on water, skip: */
-        if(this->properties.onwater)return;
-
-        if (this->properties.collapse){
+        if(onwater)return;
+
+        if (collapse){
                 /*Bail out element if collapsed (2d) and not on bed*/
-                if (!this->properties.onbed) return;
+                if (!onbed) return;
 
                 /*This element should be collapsed into a tria element at its base. Create this tria element, 
                  * and compute gardj*/
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-                tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type);
+                tria->GradjB(grad_g,analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -3895 +4151 @@
                 /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-                tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type);
+                tria->GradjB(grad_g,analysis_type,sub_analysis_type);
                 delete tria;
                 return;
@@ -3912 +4168 @@
         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(this->properties.onwater)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 ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
                 return 0;
         }
-        else if (this->properties.collapse){
+        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(inputs,analysis_type,sub_analysis_type);
+                J=tria->Misfit(analysis_type,sub_analysis_type);
                 delete tria;
                 return J;
@@ -3933 +4201 @@
 
                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                J=tria->Misfit(inputs,analysis_type,sub_analysis_type);
+                J=tria->Misfit(analysis_type,sub_analysis_type);
                 delete tria;
                 return J;
@@ -4045 +4313 @@
         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(this->properties.onwater)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 ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
+        if ((!collapse && !onsurface) || (collapse && !onbed)){
                 return 0;
         }
-        else if (this->properties.collapse){
+        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(inputs,analysis_type,sub_analysis_type);
+                S=tria->SurfaceArea(analysis_type,sub_analysis_type);
                 delete tria;
                 return S;
@@ -4066 +4346 @@
 
                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
-                S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type);
+                S=tria->SurfaceArea(analysis_type,sub_analysis_type);
                 delete tria;
                 return S;

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

@@ -16 +16 @@
 #include "./Element.h"
 #include "./Matpar.h"
-#include "./Numpar.h"
 #include "./Matice.h"
 #include "./Tria.h"
@@ -45 +44 @@
                 /*}}}*/
                 /*FUNCTION object management {{{1*/
-                void  Configure(void* loads,void* nodes,void* materials,void* parameters);
+                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
                 Object* copy();
                 void  DeepEcho();
@@ -57 +56 @@
                 int   MyRank();
                 void*  SpawnTria(int g0, int g1, int g2);
-                void  UpdateFromDakota(void* inputs);
-                void  UpdateInputs(double* solution, int analysis_type, int sub_analysis_type);
                 void  SetClone(int* minranks);
 
@@ -107 +104 @@
                 void  CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type);
                 void  CreatePVectorPrognostic( Vec pg,  int analysis_type,int sub_analysis_type);
-                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
-                void  CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
-                void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
-                void  CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
 
                 void  CreateKMatrixDiagnosticStokes( Mat Kgg,  int analysis_type,int sub_analysis_type);
@@ -136 +129 @@
                 void  GetPhi(double* phi, double*  epsilon, double viscosity);
                 double MassFlux(double* segment,double* ug);
+
+                /*updates: */
+                void  UpdateFromDakota(void* inputs);
+                void  UpdateInputs(double* solution, int analysis_type, int sub_analysis_type);
+                void  UpdateInputsDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type);
+                void  UpdateInputsSlopeCompute( double* solution,int analysis_type,int sub_analysis_type);
+                void  UpdateInputsPrognostic( double* solution,int analysis_type,int sub_analysis_type);
+                void  UpdateInputsPrognostic2(double* solution,int analysis_type,int sub_analysis_type);
+                void  UpdateInputsBalancedthickness( double* solution,int analysis_type,int sub_analysis_type);
+                void  UpdateInputsBalancedthickness2( double* solution,int analysis_type,int sub_analysis_type);
+                void  UpdateInputsBalancedvelocities( double* solution,int analysis_type,int sub_analysis_type);
+
                 /*}}}*/
 

issm/trunk/src/c/objects/Sing.cpp

@@ -11 +11 @@
 #include "stdio.h"
 #include "./Sing.h"
+#include "./SingVertexInput.h"
 #include <string.h>
 #include "../EnumDefinitions/EnumDefinitions.h"
 #include "../shared/shared.h"
 #include "../DataSet/DataSet.h"
+#include "../DataSet/Inputs.h"
 #include "../include/typedefs.h"
 #include "../include/macros.h"
@@ -21 +23 @@
 /*FUNCTION Sing::constructor {{{1*/
 Sing::Sing(){
+        this->inputs=NULL;
         return;
 }
 /*}}}*/
 /*FUNCTION Sing constructor {{{1*/
-Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id, ElementProperties* singproperties): 
+Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id):
         hnodes(sing_node_ids,1),
         hmatice(&sing_matice_id,1),
         hmatpar(&sing_matpar_id,1),
-        hnumpar(&sing_numpar_id,1),
-        properties(singproperties)
+        hnumpar(&sing_numpar_id,1)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=sing_id;
+        this->inputs=new Inputs();
 
         return;
@@ -40 +43 @@
 /*}}}*/
 /*FUNCTION Sing other constructor {{{1*/
-Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, ElementProperties* sing_properties):
+Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs):
         hnodes(sing_hnodes),
         hmatice(sing_hmatice),
         hmatpar(sing_hmatpar),
-        hnumpar(sing_hnumpar),
-        properties(sing_properties)
+        hnumpar(sing_hnumpar)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=sing_id;
-
+        if(sing_inputs){
+                this->inputs=(Inputs*)sing_inputs->Copy();
+        }
+        else{
+                this->inputs=new Inputs();
+        }
         return;
 }
@@ -80 +87 @@
         this->hnumpar.Init(&sing_numpar_id,1);
 
-        /*properties: */
-        sing_h=iomodel->thickness[i];
-        sing_k=iomodel->drag[i];
-
-        this->properties.Init(1,&sing_h, NULL, NULL, &sing_k, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, UNDEF,UNDEF, UNDEF,UNDEF,UNDEF);
+        //intialize inputs, and add as many inputs per element as requested: 
+        this->inputs=new Inputs();
+
+        if (iomodel->thickness) this->inputs->AddInput(new SingVertexInput(ThicknessEnum,iomodel->thickness[i]));
+        if (iomodel->drag_coefficient) this->inputs->AddInput(new SingVertexInput(DragCoefficientEnum,iomodel->drag_coefficient[i]));
 
 }
@@ -90 +97 @@
 /*FUNCTION Sing::destructor {{{1*/
 Sing::~Sing(){
+        delete inputs;
         return;
 }
@@ -136 +144 @@
         hmatpar.DeepEcho();
         hnumpar.DeepEcho();
-        properties.DeepEcho();
+        printf("   inputs\n");
+        inputs->DeepEcho();
 
         return;
@@ -160 +169 @@
         hnumpar.Demarshall(&marshalled_dataset);
 
-        /*demarshall properties: */
-        properties.Demarshall(&marshalled_dataset);
+        /*demarshall inputs: */
+        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
 
         /*return: */
@@ -178 +187 @@
         hmatpar.Echo();
         hnumpar.Echo();
-        properties.Echo();
-
-        return;
+        printf("   inputs\n");
+        inputs->Echo();
 }
 /*}}}*/
@@ -188 +196 @@
         char* marshalled_dataset=NULL;
         int   enum_type=0;
+        char* marshalled_inputs=NULL;
+        int   marshalled_inputs_size;
 
         /*recover marshalled_dataset: */
@@ -207 +217 @@
         hnumpar.Marshall(&marshalled_dataset);
 
-        /*Marshall properties: */
-        properties.Marshall(&marshalled_dataset);
+        /*Marshall inputs: */
+        marshalled_inputs_size=inputs->MarshallSize();
+        marshalled_inputs=inputs->Marshall();
+        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
+        marshalled_dataset+=marshalled_inputs_size;
+        
+        xfree((void**)&marshalled_inputs);
 
         *pmarshalled_dataset=marshalled_dataset;
@@ -222 +237 @@
                 +hmatpar.MarshallSize()
                 +hnumpar.MarshallSize()
-                +properties.MarshallSize()
+                +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
 }
@@ -243 +258 @@
 void  Sing::ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type){
 
-        int i;
-        const int numgrids=1;
-        int doflist[numgrids];
-        double pressure[numgrids];
+        int    dof;
+        double pressure;
+        double thickness;
         double rho_ice,g;
 
         /*dynamic objects pointed to by hooks: */
-        Node**  nodes=NULL;
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*Get dof list on which we will plug the pressure values: */
-        GetDofList1(&doflist[0]);
+        GetDofList1(&dof);
 
         /*recover objects from hooks: */
-        nodes=(Node**)hnodes.deliverp();
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*pressure is lithostatic: */
         rho_ice=matpar->GetRhoIce();
         g=matpar->GetG();
-        pressure[0]=rho_ice*g*this->properties.h[0];
+        inputs->GetParameterValue(&thickness,ThicknessEnum);
+        pressure=rho_ice*g*thickness;
         
         /*plug local pressure values into global pressure vector: */
-        VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
+        VecSetValue(p_g,dof,pressure,INSERT_VALUES);
 
 }
@@ -282 +293 @@
 /*}}}*/
 /*FUNCTION Sing::CostFunction {{{1*/
-double Sing::CostFunction(void*, int,int){
+double Sing::CostFunction( int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -293 +304 @@
         if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){
 
-                CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
+                CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type);
 
         }
@@ -313 +324 @@
         int    numberofdofspernode;
 
-        ParameterInputs* inputs=NULL;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-        
         GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -330 +336 @@
         if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){
         
-                        CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
+                        CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type);
 
         }
@@ -365 +371 @@
         Numpar* numpar=NULL;
 
-        ParameterInputs* inputs=NULL;
-
-        /*recover pointers: */
-        inputs=(ParameterInputs*)vinputs;
-
         /*recover objects from hooks: */
         node=(Node**)hnodes.deliverp();
@@ -376 +377 @@
         numpar=(Numpar*)hnumpar.delivers();
 
-        found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,numgrids,(void**)&node[0]);
-        if(!found)ISSMERROR(" surfaceslopex missing from inputs!");
-        found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,numgrids,(void**)&node[0]);
-        if(!found)ISSMERROR(" surfaceslopey missing from inputs!");
+        inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum);
+        inputs->GetParameterValue(&slope[1],SurfaceSlopeyEnum);
 
         GetDofList(&doflist[0],&numberofdofspernode);
@@ -391 +390 @@
         n=matice->GetN();
         B=matice->GetB();
-        thickness=this->properties.h[0];
+        inputs->GetParameterValue(&thickness,ThicknessEnum);
 
         ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
@@ -408 +407 @@
 /*}}}*/
 /*FUNCTION Sing::Du {{{1*/
-void  Sing::Du(_p_Vec*,void*,int,int){
+void  Sing::Du(Vec,int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -504 +503 @@
 /*}}}*/
 /*FUNCTION Sing::GetShelf {{{1*/
-int   Sing::GetShelf(){
+bool   Sing::GetShelf(){
         ISSMERROR(" not supported yet!");
 }
@@ -514 +513 @@
 /*}}}*/
 /*FUNCTION Sing::Gradj {{{1*/
-void  Sing::Gradj(_p_Vec*, void*, int, int ,char*){
+void  Sing::Gradj(Vec,  int, int ,char*){
         ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Sing::GradB {{{1*/
-void  Sing::GradjB(_p_Vec*, void*, int,int){
+void  Sing::GradjB(Vec,  int,int){
         ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Sing::GradjDrag {{{1*/
-void  Sing::GradjDrag(_p_Vec*, void*, int,int){
+void  Sing::GradjDrag(Vec,  int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -534 +533 @@
 /*}}}*/
 /*FUNCTION Sing::Misfit {{{1*/
-double Sing::Misfit(void*, int,int){
+double Sing::Misfit( int,int){
         ISSMERROR(" not supported yet!");
 }
@@ -545 +544 @@
 /*}}}*/
 /*FUNCTION Sing::SurfaceArea {{{1*/
-double Sing::SurfaceArea(void*, int,int){
+double Sing::SurfaceArea( int,int){
         ISSMERROR(" not supported yet!");
 }

issm/trunk/src/c/objects/Sing.h

@@ -12 +12 @@
 #include "../ModelProcessorx/IoModel.h"
 #include "./Hook.h"
-
-class Hook;
+#include "../DataSet/Inputs.h"
 
 class Sing: public Element{
@@ -26 +25 @@
                 Hook hmatpar; //hook to 1 matpar
                 Hook hnumpar; //hook to 1 numpar
-
                 Inputs* inputs;
 
@@ -39 +37 @@
                 /*}}}*/
                 /*object management: {{{1*/
-                void  Configure(void* loads,void* nodes,void* materials,void* parameters);
+                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
                 Object* copy();
                 void  DeepEcho();
@@ -72 +70 @@
                 void  GetBedList(double*);
                 void  GetThicknessList(double* thickness_list);
-                void  Du(_p_Vec*,void*,int,int);
-                void  Gradj(_p_Vec*, void*, int, int,char*);
-                void  GradjDrag(_p_Vec*, void*, int,int);
-                void  GradjB(_p_Vec*, void*, int,int);
-                double Misfit(void*,int,int);
-                double SurfaceArea(void*,int,int);
-                double CostFunction(void*,int,int);
+                void  Du(Vec,int,int);
+                void  Gradj(Vec, int, int,char*);
+                void  GradjDrag(Vec , int,int);
+                void  GradjB(Vec,  int,int);
+                double Misfit(int,int);
+                double SurfaceArea(int,int);
+                double CostFunction(int,int);
                 double MassFlux(double* segment,double* ug);
                 /*}}}*/

issm/trunk/src/c/objects/Tria.cpp

@@ -30 +30 @@
 }
 /*}}}*/
-/*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id){{{1*/
-Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id): 
+/*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id){{{1*/
+Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id): 
         hnodes(tria_node_ids,3),
         hmatice(&tria_matice_id,1),
-        hmatpar(&tria_matpar_id,1),
-        hnumpar(&tria_numpar_id,1)
+        hmatpar(&tria_matpar_id,1)
 {
 
         /*all the initialization has been done by the initializer, just fill in the id: */
         this->id=tria_id;
+        this->parameters=NULL;
         this->inputs=new Inputs();
 
-        return;
-}
-/*}}}*/
-/*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, Inputs* tria_inputs) {{{1*/
-Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs):
+}
+/*}}}*/
+/*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* tria_inputs) {{{1*/
+Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* tria_parameters, Inputs* tria_inputs):
         hnodes(tria_hnodes),
         hmatice(tria_hmatice),
-        hmatpar(tria_hmatpar),
-        hnumpar(tria_hnumpar)
+        hmatpar(tria_hmatpar)
 {
 
@@ -61 +59 @@
                 this->inputs=new Inputs();
         }
-
-        return;
+        /*point parameters: */
+        this->parameters=tria_parameters;
 }
 /*}}}*/
@@ -72 +70 @@
         int tria_matice_id;
         int tria_matpar_id;
-        int tria_numpar_id;
         double nodeinputs[3];
 
@@ -94 +91 @@
         tria_matice_id=index+1; //refers to the corresponding ice material object
         tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
-        tria_numpar_id=1; //refers to numerical parameters object
 
         this->hnodes.Init(tria_node_ids,3);
         this->hmatice.Init(&tria_matice_id,1);
         this->hmatpar.Init(&tria_matpar_id,1);
-        this->hnumpar.Init(&tria_numpar_id,1);
 
         //intialize inputs, and add as many inputs per element as requested: 
@@ -143 +138 @@
         if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
 
+        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
+        this->parameters=NULL;
+
 
 }
@@ -149 +147 @@
 Tria::~Tria(){
         delete inputs;
+        this->parameters=NULL;
         return;
 }
@@ -155 +154 @@
 /*Object management: */
 /*FUNCTION Tria::Configure {{{1*/
-void  Tria::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){
+void  Tria::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){
 
         int i;
-
-        DataSet* loadsin=NULL;
-        DataSet* nodesin=NULL;
-        DataSet* materialsin=NULL;
-        DataSet* parametersin=NULL;
-
-        /*Recover pointers :*/
-        loadsin=(DataSet*)ploadsin;
-        nodesin=(DataSet*)pnodesin;
-        materialsin=(DataSet*)pmaterialsin;
-        parametersin=(DataSet*)pparametersin;
 
         /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
@@ -175 +163 @@
         hmatice.configure(materialsin);
         hmatpar.configure(materialsin);
-        hnumpar.configure(parametersin);
+
+        /*point parameters to real dataset: */
+        this->parameters=parametersin;
 
 }
@@ -182 +172 @@
 Object* Tria::copy() {
 
-        return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs);
+        return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs);
 
 }
@@ -205 +195 @@
         hmatice.Demarshall(&marshalled_dataset);
         hmatpar.Demarshall(&marshalled_dataset);
-        hnumpar.Demarshall(&marshalled_dataset);
         
         /*demarshall inputs: */
         inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
+
+        /*parameters: may not exist even yet, so let Configure handle it: */
+        this->parameters=NULL;
 
         /*return: */
@@ -224 +216 @@
         hmatice.DeepEcho();
         hmatpar.DeepEcho();
-        hnumpar.DeepEcho();
+        printf("   parameters\n");
+        parameters->DeepEcho();
         printf("   inputs\n");
         inputs->DeepEcho();
-
+        
         return;
 }
@@ -240 +233 @@
         hmatice.Echo();
         hmatpar.Echo();
-        hnumpar.Echo();
+        printf("   parameters\n");
+        parameters->Echo();
         printf("   inputs\n");
         inputs->Echo();
-
-        return;
 }
 /*}}}*/
@@ -271 +263 @@
         hmatice.Marshall(&marshalled_dataset);
         hmatpar.Marshall(&marshalled_dataset);
-        hnumpar.Marshall(&marshalled_dataset);
 
         /*Marshall inputs: */
@@ -279 +270 @@
         marshalled_dataset+=marshalled_inputs_size;
 
+        /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/
+
+        xfree((void**)&marshalled_inputs);
+
         *pmarshalled_dataset=marshalled_dataset;
         return;
@@ -290 +285 @@
                 +hmatice.MarshallSize()
                 +hmatpar.MarshallSize()
-                +hnumpar.MarshallSize()
                 +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
@@ -312 +306 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -318 +311 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*Update internal data if inputs holds new values: */
@@ -468 +460 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -474 +465 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*plug local pressure values into global pressure vector: */
@@ -497 +487 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -503 +492 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*Get dof list on which we will plug the pressure values: */
@@ -536 +524 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -542 +529 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*plug local pressure values into global pressure vector: */
@@ -583 +569 @@
         double  dk[NDOF2]; 
         double  dB[NDOF2]; 
+        char*   control_type=NULL;
+        double  cm_noisedmp;
 
         /* Jacobian: */
@@ -591 +579 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*inputs: */
@@ -608 +595 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&control_type,"control_type");
+        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
 
         /* Get node coordinates and dof list: */
@@ -631 +621 @@
 
                 /*Add Tikhonov regularization term to misfit*/
-                if (strcmp(numpar->control_type,"drag")==0){
+                if (strcmp(control_type,"drag")==0){
                         if (shelf){
                                 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
-                                Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
+                                Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
 
                         }
                 }
-                else if (strcmp(numpar->control_type,"B")==0){
+                else if (strcmp(control_type,"B")==0){
                         inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBEnum);
-                        Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
+                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
                 }
                 else{
-                        ISSMERROR("%s%s","unsupported control type: ",numpar->control_type);
+                        ISSMERROR("%s%s","unsupported control type: ",control_type);
                 }
 
@@ -652 +642 @@
         xfree((void**)&third_gauss_area_coord);
         xfree((void**)&gauss_weights);
+        xfree((void**)&control_type);
 
         /*Return: */
@@ -763 +754 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
+
+        /*parameters: */
+        double artdiff;
 
 
@@ -770 +763 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -780 +772 @@
         inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
 
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&artdiff,"artdiff");
+
         //Create Artificial diffusivity once for all if requested
-        if(numpar->artdiff){
+        if(artdiff){
                 //Get the Jacobian determinant
                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
@@ -848 +843 @@
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
 
-                if(numpar->artdiff){
+                if(artdiff){
 
                         /* Compute artificial diffusivity */
@@ -918 +913 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
 
@@ -925 +919 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -1034 +1027 @@
         int     found=0;
 
+        /*parameters: */
+        double artdiff;
+
         /*dynamic objects pointed to by hooks: */
         Node**  nodes=NULL;
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -1044 +1039 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*Recover velocity: */
         inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
         inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&artdiff,"artdiff");
 
         /* Get node coordinates and dof list: */
@@ -1075 +1072 @@
 
         //Create Artificial diffusivity once for all if requested
-        if(numpar->artdiff){
+        if(artdiff){
                 //Get the Jacobian determinant
                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
@@ -1131 +1128 @@
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
 
-                if(numpar->artdiff){
+                if(artdiff){
 
                         /* Compute artificial diffusivity */
@@ -1195 +1192 @@
         double B[3][numdof];
         double Bprime[3][numdof];
-        double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
+        double D[3][3]={0.0};  // material matrix, simple scalar matrix.
         double D_scalar;
 
+        /*parameters: */
+        double viscosity_overshoot;
+
         /* local element matrices: */
-        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg[numdof][numdof]={0.0};
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
 
@@ -1212 +1212 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*inputs: */
@@ -1220 +1219 @@
         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&viscosity_overshoot,"viscosity_overshoot");
 
         /*First, if we are on water, return empty matrix: */
@@ -1228 +1230 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set Ke_gg to 0: */
-        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -1265 +1263 @@
                 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
                         onto this scalar matrix, so that we win some computational time: */
-                newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity);
+                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
                 D_scalar=newviscosity*thickness*gauss_weight*Jdet;
 
@@ -1333 +1331 @@
 
         /* local element matrices: */
-        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg[numdof][numdof]={0.0};
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
         
@@ -1365 +1363 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
 
@@ -1372 +1369 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*retrieve inputs :*/
@@ -1381 +1377 @@
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set Ke_gg to 0: */
-        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         if (shelf){
@@ -1522 +1515 @@
 
         /* local element matrices: */
-        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
 
@@ -1529 +1522 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -1535 +1527 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set Ke_gg to 0: */
-        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -1655 +1643 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -1661 +1648 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*Recover constants of ice */
@@ -1759 +1745 @@
         double  K[2][2]={0.0};
         double  KDL[2][2]={0.0};
-        double  dt;
         int     dofs[2]={0,1};
         int     found;
+
+        /*parameters: */
+        double dt;
+        double artdiff;
 
         /*dynamic objects pointed to by hooks: */
@@ -1767 +1756 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -1773 +1761 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover dt: */
-        dt=numpar->dt;
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
+        this->parameters->FindParam(&artdiff,"artdiff");
 
         /* Get node coordinates and dof list: */
@@ -1787 +1775 @@
 
         //Create Artificial diffusivity once for all if requested
-        if(numpar->artdiff==1){
+        if(artdiff==1){
                 //Get the Jacobian determinant
                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
@@ -1867 +1855 @@
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
 
-                if(numpar->artdiff==1){
+                if(artdiff==1){
 
                         /* Compute artificial diffusivity */
@@ -1934 +1922 @@
         /*input parameters for structural analysis (diagnostic): */
         double  vx,vy;
-        double  dt;
         int     dofs[1]={0};
         int     found;
+
+        /*parameters: */
+        double dt;
 
         /*dynamic objects pointed to by hooks: */
@@ -1942 +1932 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -1948 +1937 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover dt: */
-        dt=numpar->dt;
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
 
         /* Get node coordinates and dof list: */
@@ -2050 +2038 @@
 
         /* local element matrices: */
-        double Ke_gg[numdof][numdof]; //local element stiffness matrix 
+        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
         
@@ -2059 +2047 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2065 +2052 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set Ke_gg to 0: */
-        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -2133 +2116 @@
         double rho_ice;
         double heatcapacity;
-        double dt;
 
         int     num_gauss,ig;
@@ -2151 +2133 @@
         double  D_scalar;
 
+        /*parameters: */
+        double dt;
+
         /*dynamic objects pointed to by hooks: */
         Node**  nodes=NULL;
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2161 +2145 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover extra inputs from users, dt: */
-        dt=numpar->dt;
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
 
         /* Get node coordinates and dof list: */
@@ -2304 +2287 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2310 +2292 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -2391 +2372 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2397 +2377 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -2478 +2457 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2484 +2462 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover objects from hooks: */
-        nodes=(Node**)hnodes.deliverp();
-        matpar=(Matpar*)hmatpar.delivers();
-        matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -2562 +2533 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdof];
+        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
 
@@ -2578 +2549 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2584 +2554 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set pe_g to 0: */
-        for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -2680 +2646 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdof];
+        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
 
@@ -2690 +2656 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*inputs: */
@@ -2707 +2672 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set pe_g to 0: */
-        for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
 
@@ -2812 +2773 @@
         double  melting_g;
         double  thickness_g;
+
+        /*parameters: */
         double  dt;
 
@@ -2818 +2781 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2824 +2786 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover dt: */
-        dt=numpar->dt;
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
 
         /* Get node coordinates and dof list: */
@@ -2909 +2870 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -2915 +2875 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
-
-        /*recover dt: */
-        dt=numpar->dt;
-
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
 
         /* Get node coordinates and dof list: */
@@ -2993 +2951 @@
 
         /*element vector at the gaussian points: */
-        double  pe_g[numdof];
+        double  pe_g[numdof]={0.0};
         double  pe_g_gaussian[numdof];
         double  slope[2];
@@ -3001 +2959 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -3007 +2964 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList(&doflist[0],&numberofdofspernode);
-
-        /* Set pe_g to 0: */
-        for(i=0;i<numdof;i++) pe_g[i]=0.0;
 
 
@@ -3112 +3065 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -3118 +3070 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
         
         /* Get node coordinates and dof list: */
@@ -3132 +3083 @@
         beta=matpar->GetBeta();
         meltingpoint=matpar->GetMeltingPoint();
-        dt=numpar->dt;
-
+        
+        /*retrieve some solution parameters: */
+        this->parameters->FindParam(&dt,"dt");
 
         /* Ice/ocean heat exchange flux on ice shelf base */
@@ -3231 +3183 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -3237 +3188 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*retrieve inputs :*/
@@ -3249 +3199 @@
         rho_ice=matpar->GetRhoIce();
         heatcapacity=matpar->GetHeatCapacity();
-        dt=numpar->dt;
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
 
 
@@ -3369 +3321 @@
         double  obs_velocity_mag,velocity_mag;
         double  dux,duy;
+        double  meanvel, epsvel;
 
         /*element vector : */
-        double  due_g[numdof]={0,0,0,0,0,0};
+        double  due_g[numdof]={0.0};
         double  due_g_gaussian[numdof];
 
@@ -3391 +3344 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -3397 +3349 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&meanvel,"meanvel");
+        this->parameters->FindParam(&epsvel,"epsvel");
 
         /* Get node coordinates and dof list: */
@@ -3460 +3415 @@
                  */
                 for (i=0;i<numgrids;i++){
-                        scalex=pow(numpar->meanvel/(obs_vx_list[i]+numpar->epsvel),2);
-                        scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),2);
+                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
+                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
                         if(obs_vx_list[i]==0)scalex=0;
                         if(obs_vy_list[i]==0)scaley=0;
@@ -3482 +3437 @@
                  */
                 for (i=0;i<numgrids;i++){
-                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+numpar->epsvel; //epsvel to avoid velocity being nil.
-                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+numpar->epsvel; //epsvel to avoid observed velocity being nil.
-                        scale=-8*pow(numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
+                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
+                        scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
                         dux_list[i]=scale*vx_list[i];
                         duy_list[i]=scale*vy_list[i];
@@ -3501 +3456 @@
                  */
                 for (i=0;i<numgrids;i++){
-                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+numpar->epsvel);
+                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
                         dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
                         duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
@@ -3518 +3473 @@
                  */
                 for (i=0;i<numgrids;i++){
-                        dux_list[i] = - pow(numpar->meanvel,(double)2)*(
-                                                log((fabs(vx_list[i])+numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)) * 1/(vx_list[i]+numpar->epsvel));
-                        duy_list[i] = - pow(numpar->meanvel,(double)2)*(
-                                                log((fabs(vy_list[i])+numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)) * 1/(vy_list[i]+numpar->epsvel));
+                        dux_list[i] = - pow(meanvel,(double)2)*(
+                                                log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
+                        duy_list[i] = - pow(meanvel,(double)2)*(
+                                                log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
                 }
         }
@@ -4228 +4183 @@
 
         /*element vector at the gaussian points: */
-        double  grade_g[numgrids];
+        double  grade_g[numgrids]={0.0};
         double  grade_g_gaussian[numgrids];
 
@@ -4252 +4207 @@
         double  dB[NDOF2]; 
         double  B_gauss;
+        
+        /*parameters: */
+        double  cm_noisedmp;
+        double  cm_mindmp_slope;
+        double  cm_mindmp_value;
+        double  cm_maxdmp_value;
+        double  cm_maxdmp_slope;
 
         /*dynamic objects pointed to by hooks: */
@@ -4257 +4219 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -4263 +4224 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
+        this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value");
+        this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope");
+        this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value");
+        this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope");
 
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
         GetDofList1(&doflist1[0]);
-
-        /* Set grade_g to 0: */
-        for(i=0;i<numgrids;i++) grade_g[i]=0.0;
 
         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -4317 +4281 @@
 
                         //Add regularization term
-                        grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);
+                        grade_g_gaussian[i]-=cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);
 
                         //min dampening
-                        if(B_gauss<numpar->cm_mindmp_value){ 
-                                grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        if(B_gauss<cm_mindmp_value){ 
+                                grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
                         }
 
                         //max dampening
-                        if(B_gauss>numpar->cm_maxdmp_value){ 
-                                grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        if(B_gauss>cm_maxdmp_value){ 
+                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
                         }
 
@@ -4394 +4358 @@
 
         /*element vector at the gaussian points: */
-        double  grade_g[numgrids]={0,0,0};
+        double  grade_g[numgrids]={0.0};
         double  grade_g_gaussian[numgrids];
 
@@ -4408 +4372 @@
         /*inputs: */
         bool shelf;
+
+        /*parameters: */
+        double  cm_noisedmp;
+        double  cm_mindmp_slope;
+        double  cm_mindmp_value;
+        double  cm_maxdmp_value;
+        double  cm_maxdmp_slope;
 
         /*dynamic objects pointed to by hooks: */
@@ -4413 +4384 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -4419 +4389 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*retrieve inputs :*/
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
+        this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value");
+        this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope");
+        this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value");
+        this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope");
 
         /*Get out if shelf*/
@@ -4520 +4496 @@
 
                         //noise dampening d/dki(1/2*(dk/dx)^2)
-                        grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
+                        grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
                         
                         //min dampening
-                        if(drag<numpar->cm_mindmp_value){ 
-                                grade_g_gaussian[i]+=numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        if(drag<cm_mindmp_value){ 
+                                grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
                         }
 
                         //max dampening
-                        if(drag>numpar->cm_maxdmp_value){ 
-                                grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        if(drag>cm_maxdmp_value){ 
+                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
                         }
                 }
@@ -4564 +4540 @@
         double vy_list[numgrids];
         double vz_list[numgrids];
-        double vxvyvz_list[numgrids][3];
-        double adjx_list[numgrids];
-        double adjy_list[numgrids];
-        double adjz_list[numgrids];
-        double adjxyz_list[numgrids][3];
-
         double drag;
-        int    dofs1[1]={0};
-        int    dofs3[3]={0,1,2};
+        double  thickness_list[numgrids];
+        double  bed_list[numgrids];
+        double  dragcoefficient_list[numgrids];
+        double  drag_p,drag_q;
+        double alpha_complement_list[numgrids];
+        double alpha_complement;
 
         /* gaussian points: */
@@ -4582 +4556 @@
         double  gauss_weight;
         double  gauss_l1l2l3[3];
+        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
 
         /* parameters: */
@@ -4591 +4566 @@
         double  dk[NDOF2]; 
 
-        /*drag: */
-        double  pcoeff,qcoeff;
-        double alpha_complement_list[numgrids];
-        double alpha_complement;
-
         /*element vector at the gaussian points: */
-        double  grade_g[numgrids];
+        double  grade_g[numgrids]={0.0};
         double  grade_g_gaussian[numgrids];
 
@@ -4612 +4582 @@
         bool shelf;
         int  drag_type;
+
+        /*parameters: */
+        double  cm_noisedmp;
+        double  cm_mindmp_slope;
+        double  cm_mindmp_value;
+        double  cm_maxdmp_value;
+        double  cm_maxdmp_slope;
 
         /*dynamic objects pointed to by hooks: */
@@ -4617 +4594 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*retrieve inputs :*/
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
         inputs->GetParameterValue(&drag_type,DragTypeEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
+        this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value");
+        this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope");
+        this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value");
+        this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope");
 
         /*Get out if shelf*/
@@ -4630 +4613 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -4636 +4618 @@
         GetDofList1(&doflist1[0]);
 
-        /* Set grade_g to 0: */
-        for(i=0;i<numgrids;i++) grade_g[i]=0.0;
-
-        /* recover input parameters: */
-        inputs->Recover("drag",&this->properties.k[0],1,dofs1,numgrids,(void**)nodes);
-        inputs->Recover("bed",&this->properties.b[0],1,dofs1,numgrids,(void**)nodes);
-        inputs->Recover("thickness",&this->properties.h[0],1,dofs1,numgrids,(void**)nodes);
-        if(!inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
-                ISSMERROR("missing velocity input parameter");
-        }
-        if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
-                ISSMERROR("missing adjoint input parameter");
-        }
-
-        /*Initialize parameter lists: */
-        for(i=0;i<numgrids;i++){
-                vx_list[i]=vxvyvz_list[i][0];
-                vy_list[i]=vxvyvz_list[i][1];
-                vz_list[i]=vxvyvz_list[i][2];
-                adjx_list[i]=adjxyz_list[i][0];
-                adjy_list[i]=adjxyz_list[i][1];
-                adjz_list[i]=adjxyz_list[i][2];
+        /*Build alpha_complement_list: */
+        if (drag_type==2){
+
+                /*Allocate friction object: */
+                Friction* friction=NewFriction();
+
+                inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
+                inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
+                inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],3,VzAverageEnum);
+                inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum);
+                inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum);
+                inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum);
+                inputs->GetParameterValue(&drag_p,DragPEnum);
+                inputs->GetParameterValue(&drag_q,DragQEnum);
+
+
+                /*Initialize all fields: */
+                friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
+                strcpy(friction->element_type,"2d");
+
+                friction->gravity=matpar->GetG();
+                friction->rho_ice=matpar->GetRhoIce();
+                friction->rho_water=matpar->GetRhoWater();
+                friction->K=&dragcoefficient_list[0];
+                friction->bed=&bed_list[0];
+                friction->thickness=&thickness_list[0];
+                friction->vx=&vx_list[0];
+                friction->vy=&vy_list[0];
+                friction->p=drag_p;
+                friction->q=drag_q;
+
+
+                if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");
+
+                /*Compute alpha complement: */
+                FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
+
+                /*Erase friction object: */
+                DeleteFriction(&friction);
+        }
+        else{
+                alpha_complement_list[0]=0;
+                alpha_complement_list[1]=0;
+                alpha_complement_list[2]=0;
         }
 
@@ -4671 +4675 @@
                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
 
-                /*Build alpha_complement_list: */
-                if (drag_type==2){
-
-                        /*Allocate friction object: */
-                        Friction* friction=NewFriction();
-
-                        /*Initialize all fields: */
-                        friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
-                        strcpy(friction->element_type,"2d");
-
-                        friction->gravity=matpar->GetG();
-                        friction->rho_ice=matpar->GetRhoIce();
-                        friction->rho_water=matpar->GetRhoWater();
-                        friction->K=&this->properties.k[0];
-                        friction->bed=&this->properties.b[0];
-                        friction->thickness=&this->properties.h[0];
-                        friction->velocities=&vxvyvz_list[0][0];
-                        friction->p=this->properties.p;
-                        friction->q=this->properties.q;
-
-                        if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");
-
-                        /*Compute alpha complement: */
-                        FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
-
-                        /*Erase friction object: */
-                        DeleteFriction(&friction);
-                }
-                else{
-                        alpha_complement_list[0]=0;
-                        alpha_complement_list[1]=0;
-                        alpha_complement_list[2]=0;
-                }
-
                 /*Recover alpha_complement and k: */
                 GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
-                GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);
+                inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
 
                 /*recover lambda mu and xi: */
-                GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
-                GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
-                GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);
+                inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum);
+                inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum);
+                inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum);
 
                 /*recover vx vy and vz: */
-                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
-                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
-                GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);
+                inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum);
+                inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum);
+                inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum);
 
                 /*Get normal vecyor to the bed */
@@ -4736 +4706 @@
 
                 /*Get k derivative: dk/dx */
-                GetParameterDerivativeValue(&dk[0], &this->properties.k[0],&xyz_list[0][0], gauss_l1l2l3);
+                inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
 
                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
@@ -4748 +4718 @@
 
                         //Add regularization term
-                        grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
+                        grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
 
                         //min dampening
-                        if(drag<numpar->cm_mindmp_value){ 
-                                grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        if(drag<cm_mindmp_value){ 
+                                grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
                         }
 
                         //max dampening
-                        if(drag>numpar->cm_maxdmp_value){ 
-                                grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
+                        if(drag>cm_maxdmp_value){ 
+                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
                         }
                 }
@@ -4839 +4809 @@
 
         /*get thickness and velocity at two segment extremities: */
-        GetParameterValue(&h1, &this->properties.h[0],gauss_1);
-        GetParameterValue(&h2, &this->properties.h[0],gauss_2);
-        GetParameterValue(&vx1, &vx_list[0],gauss_1);
-        GetParameterValue(&vy1, &vy_list[0],gauss_1);
-        GetParameterValue(&vx2, &vx_list[0],gauss_2);
-        GetParameterValue(&vy2, &vy_list[0],gauss_2);
-
+        inputs->GetParameterValue(&h1, &gauss_1[0],ThicknessEnum);
+        inputs->GetParameterValue(&h2, &gauss_2[0],ThicknessEnum);
+        inputs->GetParameterValue(&vx1, &gauss_1[0],VxEnum);
+        inputs->GetParameterValue(&vx2, &gauss_2[0],VxEnum);
+        inputs->GetParameterValue(&vy1, &gauss_1[0],VyEnum);
+        inputs->GetParameterValue(&vy2, &gauss_2[0],VyEnum);
 
         mass_flux= rho_ice*length*(  
@@ -4871 +4840 @@
 
         /* grid data: */
-        double vxvy_list[numgrids][2];
         double vx_list[numgrids];
         double vy_list[numgrids];
-        double obs_vxvy_list[numgrids][2];
         double obs_vx_list[numgrids];
         double obs_vy_list[numgrids];
@@ -4888 +4855 @@
         double  gauss_weight;
         double  gauss_l1l2l3[3];
+        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
 
         /* parameters: */
@@ -4906 +4874 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*inputs: */
         bool onwater;
+        
+        double  meanvel, epsvel;
 
         /*retrieve inputs :*/
@@ -4921 +4890 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -4927 +4895 @@
 
         /* Recover input data: */
-        if(!inputs->Recover("fit",&fit)) ISSMERROR(" missing fit input parameter");
-        if(fit==3 && !inputs->Recover("surfacearea",&S)){
-                ISSMERROR("missing surface area input parameter");
-        }
-        if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
-                ISSMERROR("missing velocity_obs input parameter");
-        }
-        if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
-                ISSMERROR("missing velocity input parameter");
-        }
-        if(!inputs->Recover("weights",&weights_list[0],1,dofs1,numgrids,(void**)nodes)){
-                ISSMERROR("missing weights input parameter");
-        }
-
-        /*Initialize velocities: */
-        for(i=0;i<numgrids;i++){
-                obs_vx_list[i]=obs_vxvy_list[i][0];
-                obs_vy_list[i]=obs_vxvy_list[i][1];
-                vx_list[i]=vxvy_list[i][0];
-                vy_list[i]=vxvy_list[i][1];
-        }
-
+        inputs->GetParameterValue(&fit,FitEnum);
+        if(fit==3){
+                inputs->GetParameterValue(&S,SurfaceAreaEnum);
+        }
+        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
+        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
+        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
+        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
+        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&meanvel,"meanvel");
+        this->parameters->FindParam(&epsvel,"epsvel");
+        
         /* Compute Misfit at the 3 nodes
          * Here we integrate linearized functions:
@@ -4979 +4939 @@
                  */
                 for (i=0;i<numgrids;i++){
-                        scalex=pow(numpar->meanvel/(obs_vx_list[i]+numpar->epsvel),(double)2);
-                        scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),(double)2);
+                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
+                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
                         if(obs_vx_list[i]==0)scalex=0;
                         if(obs_vy_list[i]==0)scaley=0;
@@ -4995 +4955 @@
                 */
                 for (i=0;i<numgrids;i++){
-                        velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid velocity being nil.
-                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid observed velocity being nil.
-                        misfit_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
+                        velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
+                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
+                        misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
                 }
         }
@@ -5020 +4980 @@
                  */
                 for (i=0;i<numgrids;i++){
-                        misfit_list[i]=0.5*pow(numpar->meanvel,(double)2)*(
-                          pow(log((fabs(vx_list[i])+numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)),(double)2) +
-                          pow(log((fabs(vy_list[i])+numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)),(double)2) );
+                        misfit_list[i]=0.5*pow(meanvel,(double)2)*(
+                          pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
+                          pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) );
                 }
         }

issm/trunk/src/c/objects/Tria.h

@@ -31 +31 @@
                 Hook hmatice; //hook to 1 matice
                 Hook hmatpar; //hook to 1 matpar
-                Hook hnumpar; //hook to 1 numpar
 
-                Inputs* inputs;
+                Parameters* parameters; //pointer to solution parameters
+                Inputs*  inputs;
         
         public:
@@ -39 +39 @@
                 /*FUNCTION constructors, destructors {{{1*/
                 Tria();
-                Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id);
-                Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs);
+                Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id);
+                Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* parameters, Inputs* tria_inputs);
                 Tria(int i, IoModel* iomodel);
                 ~Tria();
                 /*}}}*/
                 /*FUNCTION object management {{{1*/
-                void  Configure(void* loads,void* nodes,void* materials,void* parameters);
+                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
                 Object* copy();
                 void  DeepEcho();