Changeset 3620


Ignore:
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:
6 added
12 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Beam.cpp

    r3612 r3620  
    1111#include "stdio.h"
    1212#include "./Beam.h"
     13#include "./BeamVertexInput.h"
    1314#include <string.h>
    1415#include "../EnumDefinitions/EnumDefinitions.h"
     
    2122/*FUNCTION Beam::Beam(){{{1*/
    2223Beam::Beam(){
     24        this->inputs=NULL;
    2325        return;
    2426}
    2527/*}}}*/
    2628/*FUNCTION Beam::Beam(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id, ElementProperties* properties){{{1*/
    27 Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id, ElementProperties* beamproperties):
     29Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id):
    2830        hnodes(beam_node_ids,2),
    2931        hmatice(&beam_matice_id,1),
    3032        hmatpar(&beam_matpar_id,1),
    31         hnumpar(&beam_numpar_id,1),
    32         properties(beamproperties)
     33        hnumpar(&beam_numpar_id,1)
    3334{
    3435
    3536        /*all the initialization has been done by the initializer, just fill in the id: */
    3637        this->id=beam_id;
    37 
    38         return;
     38        this->inputs=new Inputs();
    3939}
    4040/*}}}*/
    4141/*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, ElementProperties* properties){{{1*/
    42 Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, ElementProperties* beam_properties):
     42Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs):
    4343        hnodes(beam_hnodes),
    4444        hmatice(beam_hmatice),
    4545        hmatpar(beam_hmatpar),
    46         hnumpar(beam_hnumpar),
    47         properties(beam_properties)
     46        hnumpar(beam_hnumpar)
    4847{
    4948
    5049        /*all the initialization has been done by the initializer, just fill in the id: */
    5150        this->id=beam_id;
    52 
     51        if(beam_inputs){
     52                this->inputs=(Inputs*)beam_inputs->Copy();
     53        }
     54        else{
     55                this->inputs=new Inputs();
     56        }
    5357        return;
    5458}
    5559/*}}}*/
    5660/*FUNCTION Beam::Beam(int i, IoModel* iomodel){{{1*/
    57 Beam::Beam(int i,IoModel* iomodel){
    58 
     61Beam::Beam(int index,IoModel* iomodel){
     62
     63
     64        int i;
    5965
    6066        /*beam constructor input: */
     
    6470        int   beam_numpar_id;
    6571        int   beam_node_ids[2];
    66         double beam_h[2];
    67         double beam_s[2];
    68         double beam_b[2];
    69         double beam_k[2];
    70         int    beam_onbed;
    71        
    72        
     72        double nodeinputs[2];
     73
    7374        /*id: */
    74         beam_id=i+1;
     75        beam_id=index+1;
    7576        this->id=beam_id;
    7677
    7778        /*hooks: */
    78         beam_matice_id=i+1; //refers to the corresponding material property card
     79        beam_matice_id=index+1; //refers to the corresponding material property card
    7980        beam_matpar_id=iomodel->numberofvertices2d*(iomodel->numlayers-1)+1;//refers to the corresponding matpar property card
    8081        beam_numpar_id=1;
    81         beam_node_ids[0]=i+1;
    82         beam_node_ids[1]=(int)iomodel->uppernodes[i]; //grid that lays right on top
     82        beam_node_ids[0]=index+1;
     83        beam_node_ids[1]=(int)iomodel->uppernodes[index]; //grid that lays right on top
    8384       
    8485        this->hnodes.Init(beam_node_ids,2);
     
    8788        this->hnumpar.Init(&beam_numpar_id,1);
    8889
    89         /*properties: */
    90         beam_h[0]=iomodel->thickness[i];
    91         beam_h[1]=iomodel->thickness[(int)(iomodel->uppernodes[i]-1)];
    92         beam_s[0]=iomodel->surface[i];
    93         beam_s[1]=iomodel->surface[(int)(iomodel->uppernodes[i]-1)];
    94         beam_b[0]=iomodel->bed[i];
    95         beam_b[1]=iomodel->bed[(int)(iomodel->uppernodes[i]-1)];
    96         beam_k[0]=iomodel->drag[i];
    97         beam_k[1]=iomodel->drag[(int)(iomodel->uppernodes[i]-1)];
    98         beam_onbed=(bool)iomodel->gridonbed[i];
    99 
    100         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);
     90        //intialize inputs, and add as many inputs per element as requested:
     91        this->inputs=new Inputs();
     92
     93        if (iomodel->thickness) {
     94                nodeinputs[0]=iomodel->thickness[index];
     95                nodeinputs[1]=iomodel->thickness[(int)(iomodel->uppernodes[index]-1)];
     96                this->inputs->AddInput(new BeamVertexInput(ThicknessEnum,nodeinputs));
     97        }
     98        if (iomodel->surface) {
     99                nodeinputs[0]=iomodel->surface[index];
     100                nodeinputs[1]=iomodel->surface[(int)(iomodel->uppernodes[index]-1)];
     101                this->inputs->AddInput(new BeamVertexInput(SurfaceEnum,nodeinputs));
     102        }       
     103        if (iomodel->bed) {
     104                nodeinputs[0]=iomodel->bed[index];
     105                nodeinputs[1]=iomodel->bed[(int)(iomodel->uppernodes[index]-1)];
     106                this->inputs->AddInput(new BeamVertexInput(BedEnum,nodeinputs));
     107        }       
     108        if (iomodel->drag_coefficient) {
     109                nodeinputs[0]=iomodel->drag_coefficient[index];
     110                nodeinputs[1]=iomodel->drag_coefficient[(int)(iomodel->uppernodes[index]-1)];
     111                this->inputs->AddInput(new BeamVertexInput(DragCoefficientEnum,nodeinputs));
     112        }       
     113        if (iomodel->gridonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->gridonbed[index]));
     114
    101115}
    102116/*}}}*/
    103117/*FUNCTION Beam::~Beam(){{{1*/
    104118Beam::~Beam(){
     119        delete inputs;
    105120        return;
    106121}
     
    136151Object* Beam::copy() {
    137152       
    138         return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties);
     153        return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs);
    139154
    140155}
     
    149164        hmatpar.DeepEcho();
    150165        hnumpar.DeepEcho();
    151         properties.DeepEcho();
     166        printf("   inputs\n");
     167        inputs->DeepEcho();
    152168
    153169        return;
     
    173189        hnumpar.Demarshall(&marshalled_dataset);
    174190
    175         /*demarshall properties: */
    176         properties.Demarshall(&marshalled_dataset);
     191        /*demarshall inputs: */
     192        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);
    177193
    178194        /*return: */
     
    190206        hmatpar.Echo();
    191207        hnumpar.Echo();
    192         properties.Echo();
     208        printf("   inputs\n");
     209        inputs->Echo();
    193210
    194211        return;
     
    200217        char* marshalled_dataset=NULL;
    201218        int   enum_type=0;
     219        char* marshalled_inputs=NULL;
     220        int   marshalled_inputs_size;
    202221
    203222        /*recover marshalled_dataset: */
     
    219238        hnumpar.Marshall(&marshalled_dataset);
    220239
    221         /*Marshall properties: */
    222         properties.Marshall(&marshalled_dataset);
     240        /*Marshall inputs: */
     241        marshalled_inputs_size=inputs->MarshallSize();
     242        marshalled_inputs=inputs->Marshall();
     243        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
     244        marshalled_dataset+=marshalled_inputs_size;
     245
     246        xfree((void**)&marshalled_inputs);
    223247
    224248        *pmarshalled_dataset=marshalled_dataset;
     
    234258                +hmatpar.MarshallSize()
    235259                +hnumpar.MarshallSize()
    236                 +properties.MarshallSize()
     260                +inputs->MarshallSize()
    237261                +sizeof(int); //sizeof(int) for enum type
    238262}
     
    259283        int doflist[numgrids];
    260284        double pressure[numgrids];
     285        double surface[numgrids];
    261286        double rho_ice,g;
    262287        double xyz_list[numgrids][3];
     288        double gauss[numgrids][numgrids]={{1,0},{0,1}};
    263289
    264290        /*dynamic objects pointed to by hooks: */
     
    286312        rho_ice=matpar->GetRhoIce();
    287313        g=matpar->GetG();
     314
     315        /*recover value of surface at gauss points (0,1) and (1,0): */
     316        inputs->GetParameterValues(&surface[0],&gauss[0][0],2,SurfaceEnum);
    288317        for(i=0;i<numgrids;i++){
    289                 pressure[i]=rho_ice*g*(this->properties.s[i]-xyz_list[i][2]);
    290         }
    291        
     318                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     319        }
     320
    292321        /*plug local pressure values into global pressure vector: */
    293322        VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
     
    303332/*}}}*/
    304333/*FUNCTION Beam::CostFunction{{{1*/
    305 double Beam::CostFunction(void*,int,int){
     334double Beam::CostFunction(int,int){
    306335        ISSMERROR(" not supported yet!");
    307336}
     
    316345                if (sub_analysis_type==HutterAnalysisEnum) {
    317346
    318                         CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
     347                        CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type);
    319348                }
    320349                else
     
    338367        double    Ke_gg[numdofs][numdofs]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
    339368        int       numberofdofspernode;
    340        
     369        bool onbed;
     370       
     371       
     372        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     373
    341374        GetDofList(&doflist[0],&numberofdofspernode);
    342375
    343         if (this->properties.onbed){
     376        if (onbed){
    344377                Ke_gg[0][0]=1;
    345378                Ke_gg[1][1]=1;
     
    367400        if (analysis_type==DiagnosticAnalysisEnum) {
    368401                if (sub_analysis_type==HutterAnalysisEnum) {
    369                         CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
     402                        CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type);
    370403                }
    371404                else
     
    403436        double*  gauss_weights=NULL;
    404437        double   gauss_weight;
     438        double   gauss1[2]={1,0};
    405439        int      ig;
    406440        double   l1l2[2];
     
    412446        double   Jdet;
    413447        double   ub,vb;
     448        double   surface,thickness;
     449       
     450        bool onbed;
    414451
    415452        /*dynamic objects pointed to by hooks: */
     
    419456        Numpar* numpar=NULL;
    420457
    421         ParameterInputs* inputs=NULL;
    422 
    423         /*recover pointers: */
    424         inputs=(ParameterInputs*)vinputs;
    425 
    426458        /*recover objects from hooks: */
    427459        nodes=(Node**)hnodes.deliverp();
     
    433465        GetDofList(&doflist[0],&numberofdofspernode);
    434466
     467        /*recover some inputs: */
     468        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     469
    435470        /*recover parameters: */
    436471        rho_ice=matpar->GetRhoIce();
     
    440475
    441476        //recover extra inputs
    442         found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,1,(void**)&nodes[0]); //only recover from fist node, second node will have the same slope
    443         if(!found)ISSMERROR(" surfaceslopex missing from inputs!");
    444        
    445         found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,1,(void**)&nodes[0]);//only recover from fist node, second node will have the same slope
    446         if(!found)ISSMERROR(" surfaceslopey missing from inputs!");
     477        inputs->GetParameterValue(&slope[0],&gauss1[0],SurfaceSlopexEnum);
     478        inputs->GetParameterValue(&slope[1],&gauss1[0],SurfaceSlopeyEnum);
     479        inputs->GetParameterValue(&surface,&gauss1[0],SurfaceEnum);
     480        inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum);
    447481
    448482        //Get all element grid data:
     
    474508               
    475509                for(j=0;j<NDOF2;j++){
    476                         pe_g_gaussian[NDOF2+j]=constant_part*pow((this->properties.s[0]-z_g)/B,n)*slope[j]*Jdet*gauss_weight;
     510                        pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight;
    477511                }
    478512               
     
    484518
    485519        //Deal with lower surface
    486         if (this->properties.onbed){
     520        if (onbed){
    487521
    488522                //compute ub
    489                 constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*this->properties.h[0];
     523                constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
    490524                ub=constant_part*slope[0];
    491525                vb=constant_part*slope[1];
     
    505539/*}}}*/
    506540/*FUNCTION Beam::Du{{{1*/
    507 void  Beam::Du(_p_Vec*,void*,int,int){
     541void  Beam::Du(Vec,int,int){
    508542        ISSMERROR(" not supported yet!");
    509543}
     
    654688/*}}}*/
    655689/*FUNCTION Beam::Gradj{{{1*/
    656 void  Beam::Gradj(_p_Vec*, void*, int, int,char*){
     690void  Beam::Gradj(Vec, int, int,char*){
    657691        ISSMERROR(" not supported yet!");
    658692}
    659693/*}}}*/
    660694/*FUNCTION Beam::GradjB{{{1*/
    661 void  Beam::GradjB(_p_Vec*, void*, int, int){
     695void  Beam::GradjB(Vec, int, int){
    662696        ISSMERROR(" not supported yet!");
    663697}
    664698/*}}}*/
    665699/*FUNCTION Beam::GradjDrag{{{1*/
    666 void  Beam::GradjDrag(_p_Vec*, void*, int,int ){
     700void  Beam::GradjDrag(Vec, int,int ){
    667701        ISSMERROR(" not supported yet!");
    668702}
     
    674708/*}}}*/
    675709/*FUNCTION Beam::Misfit{{{1*/
    676 double Beam::Misfit(void*,int,int){
     710double Beam::Misfit(int,int){
    677711        ISSMERROR(" not supported yet!");
    678712}
     
    685719/*}}}*/
    686720/*FUNCTION Beam::SurfaceArea{{{1*/
    687 double Beam::SurfaceArea(void*,int,int){
     721double Beam::SurfaceArea(int,int){
    688722        ISSMERROR(" not supported yet!");
    689723}
  • issm/trunk/src/c/objects/Beam.h

    r3612 r3620  
    1212#include "./Matpar.h"
    1313#include "../ModelProcessorx/IoModel.h"
     14#include "../DataSet/Inputs.h"
    1415#include "./Hook.h"
    1516
     
    3031                Hook hmatpar; //hook to 1 matpar
    3132                Hook hnumpar; //hook to 1 numpar
    32 
    3333                Inputs* inputs;
    3434       
     
    5252                int   GetId();
    5353                int   MyRank();
    54                 void  Configure(void* loads,void* nodes,void* materials,void* parameters);
     54                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
    5555                Object* copy();
    5656                void  SetClone(int* minranks);
     
    7878                void  GetBedList(double*);
    7979                void  GetThicknessList(double* thickness_list);
    80                 void  Du(_p_Vec*,void*, int,int);
    81                 void  Gradj(_p_Vec*, void*, int, int,char*);
    82                 void  GradjDrag(_p_Vec*, void*, int,int );
    83                 void  GradjB(_p_Vec*, void*, int,int );
    84                 double Misfit(void*,int,int);
    85                 double SurfaceArea(void*,int,int);
    86                 double CostFunction(void*,int,int);
     80                void  Du(Vec, int,int);
     81                void  Gradj(Vec, int, int,char*);
     82                void  GradjDrag(Vec, int,int );
     83                void  GradjB(Vec, int,int );
     84                double Misfit(int,int);
     85                double SurfaceArea(int,int);
     86                double CostFunction(int,int);
    8787                void  GetNodalFunctions(double* l1l2, double gauss_coord);
    8888                void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);
  • issm/trunk/src/c/objects/Element.h

    r3612 r3620  
    1010
    1111#include "./Object.h"
     12#include "../DataSet/DataSet.h"
     13#include "../DataSet/Parameters.h"
    1214#include "../toolkits/toolkits.h"
    1315
     
    1820                virtual        ~Element(){};
    1921                virtual int    Enum()=0;
    20                 virtual void   Configure(void* loads,void* nodes,void* materials,void* parameters)=0;
     22                virtual void   Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters)=0;
    2123               
    2224                virtual void   CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type)=0;
  • issm/trunk/src/c/objects/FemModel.cpp

    r3496 r3620  
    3939/*FUNCTION FemModel::FemModel {{{1*/
    4040FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_vertices, DataSet* femmodel_constraints,DataSet* femmodel_loads,
    41                 DataSet* femmodel_materials,DataSet* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg,
     41                DataSet* femmodel_materials,Parameters* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg,
    4242                Mat femmodel_Rmg,Mat femmodel_Gmn,NodeSets* femmodel_nodesets,Vec femmodel_ys,Vec femmodel_ys0){
    4343
     
    266266/*}}}*/
    267267/*FUNCTION FemModel::get_parameters {{{1*/
    268 DataSet* FemModel::get_parameters(void){
     268Parameters* FemModel::get_parameters(void){
    269269        return parameters;
    270270}
  • issm/trunk/src/c/objects/FemModel.h

    r3463 r3620  
    1111#include "./Object.h"
    1212#include "../DataSet/DataSet.h"
     13#include "../DataSet/Parameters.h"
    1314#include "./DofVec.h"
    1415#include "../toolkits/toolkits.h"
     
    2728                DataSet*            loads;
    2829                DataSet*            materials;
    29                 DataSet*            parameters;
     30                Parameters*            parameters;
    3031
    3132                DofVec*             partition;
     
    4142                FemModel();
    4243                ~FemModel();
    43                 FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters,
     44                FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,Parameters* parameters,
    4445                                      DofVec* partition,DofVec* tpartition,DofVec* yg,Mat Rmg,Mat Gmn,NodeSets* nodesets,Vec ys,Vec ys0);
    4546     
     
    6970                DataSet* get_loads(void);
    7071                DataSet* get_materials(void);
    71                 DataSet* get_parameters(void);
     72                Parameters* get_parameters(void);
    7273                DofVec*      get_partition(void);
    7374                DofVec*      get_tpartition(void);
  • issm/trunk/src/c/objects/Model.cpp

    r3567 r3620  
    4747        DataSet*            loads=NULL;
    4848        DataSet*            materials=NULL;
    49         DataSet*            parameters=NULL;
     49        Parameters*            parameters=NULL;
    5050        DofVec*                 partition=NULL;
    5151        DofVec*                 tpartition=NULL;
     
    121121        DataSet*            loads=NULL;
    122122        DataSet*            materials=NULL;
    123         DataSet*            parameters=NULL;
     123        Parameters*            parameters=NULL;
    124124        DofVec*                 partition=NULL;
    125125        DofVec*                 tpartition=NULL;
  • issm/trunk/src/c/objects/Penta.cpp

    r3613 r3620  
    1010
    1111#include "stdio.h"
     12#include "./PentaVertexInput.h"
    1213#include "./Penta.h"
    1314#include <string.h>
     
    2021/*FUNCTION Penta default constructor {{{1*/
    2122Penta::Penta(){
     23        this->inputs=NULL;
    2224        return;
    2325}
    2426/*}}}*/
    2527/*FUNCTION Penta constructor {{{1*/
    26 Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id, ElementProperties* pentaproperties):
     28Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id):
    2729        hnodes(penta_node_ids,6),
    2830        hmatice(&penta_matice_id,1),
    2931        hmatpar(&penta_matpar_id,1),
    30         hnumpar(&penta_numpar_id,1),
    31         properties(pentaproperties)
     32        hnumpar(&penta_numpar_id,1)
    3233{
    3334
    3435        /*all the initialization has been done by the initializer, just fill in the id: */
    3536        this->id=penta_id;
    36 
    37         return;
     37        this->inputs=new Inputs();
     38
    3839}
    3940/*}}}*/
    4041/*FUNCTION Penta other constructor {{{1*/
    41 Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties):
     42Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* penta_inputs):
    4243        hnodes(penta_hnodes),
    4344        hmatice(penta_hmatice),
    4445        hmatpar(penta_hmatpar),
    45         hnumpar(penta_hnumpar),
    46         properties(penta_properties)
     46        hnumpar(penta_hnumpar)
    4747{
    4848
    4949        /*all the initialization has been done by the initializer, just fill in the id: */
    5050        this->id=penta_id;
    51 
    52         return;
     51        if(penta_inputs){
     52                this->inputs=(Inputs*)penta_inputs->Copy();
     53        }
     54        else{
     55                this->inputs=new Inputs();
     56        }
    5357}
    5458/*}}}*/
    5559/*FUNCTION Penta other constructor {{{1*/
    56 Penta::Penta(int i, IoModel* iomodel){ //i is the element index
    57 
    58         int j;
     60Penta::Penta(int index, IoModel* iomodel){ //i is the element index
     61
     62        IssmInt i;
    5963        int penta_node_ids[6];
    6064        int penta_matice_id;
    6165        int penta_matpar_id;
    6266        int penta_numpar_id;
     67        double nodeinputs[6];
     68        bool collapse;
    6369
    6470        /*id: */
    65         this->id=i+1;
     71        this->id=index+1;
    6672
    6773        /*hooks: */
    68         for(j=0;j<6;j++){ //go recover node ids, needed to initialize the node hook.
    69                 penta_node_ids[j]=(int)*(iomodel->elements+6*i+j); //ids for vertices are in the elements array from Matlab
    70         }
    71         penta_matice_id=i+1; //refers to the corresponding ice material object
     74        for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
     75                penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab
     76        }
     77        penta_matice_id=index+1; //refers to the corresponding ice material object
    7278        penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
    7379        penta_numpar_id=1; //refers to numerical parameters object
    7480
    75         this->hnodes.Init(penta_node_ids,6);
     81        this->hnodes.Init(&penta_node_ids[0],6);
    7682        this->hmatice.Init(&penta_matice_id,1);
    7783        this->hmatpar.Init(&penta_matpar_id,1);
    7884        this->hnumpar.Init(&penta_numpar_id,1);
    79         this->properties.Init(6,penta_node_ids,this->id,iomodel);
     85
     86        //intialize inputs, and add as many inputs per element as requested:
     87        this->inputs=new Inputs();
     88
     89        if (iomodel->thickness) {
     90                for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1];
     91                this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs));
     92        }
     93        if (iomodel->surface) {
     94                for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_node_ids[i]-1];
     95                this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs));
     96        }
     97        if (iomodel->bed) {
     98                for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_node_ids[i]-1];
     99                this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs));
     100        }
     101        if (iomodel->drag_coefficient) {
     102                for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_node_ids[i]-1];
     103                this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs));
     104
     105                if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
     106                if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
     107                this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
     108
     109        }
     110        if (iomodel->melting_rate) {
     111                for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_node_ids[i]-1];
     112                this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs));
     113        }
     114        if (iomodel->accumulation_rate) {
     115                for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_node_ids[i]-1];
     116                this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs));
     117        }
     118        if (iomodel->geothermalflux) {
     119                for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_node_ids[i]-1];
     120                this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs));
     121        }       
     122
     123        if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
     124        if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
     125        if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
     126        if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
    80127
    81128        //elements of type 3 are macayeal type penta. we collapse the formulation on their base.
    82129        if(iomodel->elements_type){
    83130                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){
    84                         this->properties.collapse=1;
     131                        collapse=1;
    85132                }
    86133                else{
    87                         this->properties.collapse=0;
     134                        collapse=0;
    88135                }
    89136        }
     137        this->inputs->AddInput(new BoolInput(CollapseEnum,(IssmBool)collapse));
    90138
    91139}
     
    93141/*FUNCTION Penta destructor {{{1*/
    94142Penta::~Penta(){
     143        delete inputs;
    95144        return;
    96145}
     
    125174/*FUNCTION copy {{{1*/
    126175Object* Penta::copy() {
    127         return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties);
     176        return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs);
    128177}
    129178/*}}}*/
     
    148197        hnumpar.Demarshall(&marshalled_dataset);
    149198
    150         /*demarshall properties: */
    151         properties.Demarshall(&marshalled_dataset);
     199        /*demarshall inputs: */
     200        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);
    152201
    153202        /*return: */
     
    166215        hmatpar.DeepEcho();
    167216        hnumpar.DeepEcho();
    168         properties.DeepEcho();
     217        printf("   inputs\n");
     218        inputs->DeepEcho();
    169219
    170220        return;
     
    181231        hmatpar.Echo();
    182232        hnumpar.Echo();
    183         properties.Echo();
    184 
    185         return;
     233        printf("   inputs\n");
     234        inputs->Echo();
     235
     236}
     237/*}}}*/
     238/*FUNCTION Enum {{{1*/
     239int Penta::Enum(void){
     240
     241        return PentaEnum;
     242
    186243}
    187244/*}}}*/
     
    191248        char* marshalled_dataset=NULL;
    192249        int   enum_type=0;
     250        char* marshalled_inputs=NULL;
     251        int   marshalled_inputs_size;
    193252
    194253        /*recover marshalled_dataset: */
     
    210269        hnumpar.Marshall(&marshalled_dataset);
    211270
    212         /*Marshall properties: */
    213         properties.Marshall(&marshalled_dataset);
     271        /*Marshall inputs: */
     272        marshalled_inputs_size=inputs->MarshallSize();
     273        marshalled_inputs=inputs->Marshall();
     274        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
     275        marshalled_dataset+=marshalled_inputs_size;
     276
     277        xfree((void**)&marshalled_inputs);
     278
    214279
    215280        *pmarshalled_dataset=marshalled_dataset;
     
    225290                +hmatpar.MarshallSize()
    226291                +hnumpar.MarshallSize()
    227                 +properties.MarshallSize()
     292                +inputs->MarshallSize()
    228293                +sizeof(int); //sizeof(int) for enum type
    229294}
     
    240305        Hook* tria_hmatpar=NULL;
    241306        Hook* tria_hnumpar=NULL;
    242         ElementProperties* tria_properties=NULL;
    243         DataSet* tria_inputs=NULL;
     307        Inputs* tria_inputs=NULL;
    244308
    245309        indices[0]=g0;
     
    251315        tria_hmatpar=this->hmatpar.Spawn(&zero,1);
    252316        tria_hnumpar=this->hnumpar.Spawn(&zero,1);
    253         tria_properties=(ElementProperties*)this->properties.Spawn(indices,3);
    254         tria_inputs=(DataSet*)this->inputs->Spawn(indices,3);
    255 
    256         tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_properties,tria_inputs);
     317        tria_inputs=(Inputs*)this->inputs->Spawn(indices,3);
     318
     319        tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_inputs);
    257320
    258321        delete tria_hnodes;
     
    260323        delete tria_hmatpar;
    261324        delete tria_hnumpar;
    262         delete tria_properties;
    263325        delete tria_inputs;
    264326
     
    266328}
    267329/*}}}*/
     330
     331/*updates: */
    268332/*FUNCTION UpdateFromDakota {{{1*/
    269333void  Penta::UpdateFromDakota(void* vinputs){
     
    275339/*FUNCTION Penta::UpdateInputs {{{1*/
    276340void  Penta::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){
     341
     342        /*Just branch to the correct UpdateInputs generator, according to the type of analysis we are carrying out: */
     343        if (analysis_type==ControlAnalysisEnum){
     344               
     345                UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type);
     346        }
     347        else if (analysis_type==DiagnosticAnalysisEnum){
     348       
     349                if (sub_analysis_type==HorizAnalysisEnum){
     350
     351                        UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type);
     352                }
     353                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
     354
     355        }
     356        else if (analysis_type==SlopecomputeAnalysisEnum){
     357
     358                UpdateInputsSlopeCompute( solution,analysis_type,sub_analysis_type);
     359        }
     360        else if (analysis_type==PrognosticAnalysisEnum){
     361
     362                UpdateInputsPrognostic( solution,analysis_type,sub_analysis_type);
     363        }
     364        else if (analysis_type==Prognostic2AnalysisEnum){
     365
     366                UpdateInputsPrognostic2(solution,analysis_type,sub_analysis_type);
     367        }
     368        else if (analysis_type==BalancedthicknessAnalysisEnum){
     369
     370                UpdateInputsBalancedthickness( solution,analysis_type,sub_analysis_type);
     371        }
     372        else if (analysis_type==Balancedthickness2AnalysisEnum){
     373
     374                UpdateInputsBalancedthickness2( solution,analysis_type,sub_analysis_type);
     375        }
     376        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
     377
     378                UpdateInputsBalancedvelocities( solution,analysis_type,sub_analysis_type);
     379        }
     380        else{
     381
     382                ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet");
     383        }
     384}
     385/*Object functions*/
     386/*FUNCTION Penta::UpdateInputsDiagnosticHoriz {{{1*/
     387void  Penta::UpdateInputsDiagnosticHoriz(double* solution, int analysis_type, int sub_analysis_type){
     388       
     389       
     390        int i;
     391
     392        const int    numvertices=6;
     393        const int    numdofpervertex=2;
     394        const int    numdof=numdofpervertex*numvertices;
     395       
     396        int          doflist[numdof];
     397        double       values[numdof];
     398        double       vx[numvertices];
     399        double       vy[numvertices];
     400
     401        int          dummy;
     402       
     403        /*Get dof list: */
     404        GetDofList(&doflist[0],&dummy);
     405
     406        /*Use the dof list to index into the solution vector: */
     407        for(i=0;i<numdof;i++){
     408                values[i]=solution[doflist[i]];
     409        }
     410
     411        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     412        for(i=0;i<numvertices;i++){
     413                vx[i]=values[i*numdofpervertex+0];
     414                vy[i]=values[i*numdofpervertex+1];
     415        }
     416
     417        /*Add vx and vy as inputs to the tria element: */
     418        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
     419        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     420}
     421
     422/*}}}*/
     423/*FUNCTION Penta::UpdateInputsSlopeCompute {{{1*/
     424void  Penta::UpdateInputsSlopeCompute(double* solution, int analysis_type, int sub_analysis_type){
     425        ISSMERROR(" not supported yet!");
     426}
     427/*}}}*/
     428/*FUNCTION Penta::UpdateInputsPrognostic {{{1*/
     429void  Penta::UpdateInputsPrognostic(double* solution, int analysis_type, int sub_analysis_type){
     430        ISSMERROR(" not supported yet!");
     431}
     432/*}}}*/
     433/*FUNCTION Penta::UpdateInputsPrognostic2 {{{1*/
     434void  Penta::UpdateInputsPrognostic2(double* solution, int analysis_type, int sub_analysis_type){
     435        ISSMERROR(" not supported yet!");
     436}
     437/*}}}*/
     438/*FUNCTION Penta::UpdateInputsBalancedthickness {{{1*/
     439void  Penta::UpdateInputsBalancedthickness(double* solution, int analysis_type, int sub_analysis_type){
     440        ISSMERROR(" not supported yet!");
     441}
     442/*}}}*/
     443/*FUNCTION Penta::UpdateInputsBalancedthickness2 {{{1*/
     444void  Penta::UpdateInputsBalancedthickness2(double* solution, int analysis_type, int sub_analysis_type){
     445        ISSMERROR(" not supported yet!");
     446}
     447/*}}}*/
     448/*FUNCTION Penta::UpdateInputsBalancedvelocities {{{1*/
     449void  Penta::UpdateInputsBalancedvelocities(double* solution, int analysis_type, int sub_analysis_type){
    277450        ISSMERROR(" not supported yet!");
    278451}
     
    295468        double bed;
    296469        double basalforce[3];
    297         double vxvyvz_list[numgrids][3];
    298         double pressure_list[numgrids];
    299470        double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    300471        double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     
    304475        int  dofv[3]={0,1,2};
    305476        int  dofp[1]={3};
    306         ParameterInputs* inputs=NULL;
    307477        double Jdet2d;
    308478        Tria* tria=NULL;
     
    323493        double surface=0;
    324494        double value=0;
     495       
     496        /*flags: */
     497        bool onbed;
    325498
    326499        /*dynamic objects pointed to by hooks: */
     
    333506        if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!");
    334507
    335         /*recover pointers: */
    336         inputs=(ParameterInputs*)vinputs;
    337 
    338508        /*recover objects from hooks: */
    339509        nodes=(Node**)hnodes.deliverp();
     
    342512        numpar=(Numpar*)hnumpar.delivers();
    343513
    344         if(!this->properties.onbed){
     514        /*recover some inputs: */
     515        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     516
     517        if(!onbed){
    345518                //put zero
    346519                VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES);
     
    351524        rho_ice=matpar->GetRhoIce();
    352525        gravity=matpar->GetG();
    353 
    354         /*recover extra inputs from users, at current convergence iteration: */
    355         inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofv,numgrids,(void**)nodes);
    356         inputs->Recover("velocity",&pressure_list[0] ,1,dofp,numgrids,(void**)nodes);
    357526
    358527        /* Get node coordinates and dof list: */
     
    378547
    379548                        /*Compute strain rate viscosity and pressure: */
    380                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     549                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
    381550                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    382                         GetParameterValue(&pressure,&pressure_list[0],&gauss_coord[0]);
     551                        inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum);
    383552
    384553                        /*Compute Stress*/
     
    417586        int i;
    418587        const int numgrids=6;
    419         int doflist[numgrids];
     588        int    doflist[numgrids];
    420589        double pressure[numgrids];
    421590        double rho_ice,g;
    422         double       xyz_list[numgrids][3];
     591        double surface[numgrids];
     592        double xyz_list[numgrids][3];
     593        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}};
     594
     595        /*inputs: */
     596        bool onwater;
    423597
    424598        /*dynamic objects pointed to by hooks: */
     
    430604        matpar=(Matpar*)hmatpar.delivers();
    431605
     606        /*retrieve inputs :*/
     607        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     608
    432609        /*If on water, skip: */
    433         if(this->properties.onwater)return;
     610        if(onwater)return;
    434611
    435612        /*Get node data: */
     
    441618        /*Get dof list on which we will plug the pressure values: */
    442619        GetDofList1(&doflist[0]);
     620
     621        /*recover value of surface at grids: */
     622        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
    443623
    444624        /*pressure is lithostatic: */
     
    446626        g=matpar->GetG();
    447627        for(i=0;i<numgrids;i++){
    448                 pressure[i]=rho_ice*g*(this->properties.s[i]-xyz_list[i][2]);
     628                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    449629        }
    450630
     
    467647        Tria* tria=NULL;
    468648
     649        /*flags: */
     650        bool onbed;
     651        bool onwater;
     652        bool collapse;
     653        bool onsurface;
     654
     655        /*recover some inputs: */
     656        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     657        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     658        inputs->GetParameterValue(&collapse,CollapseEnum);
     659        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     660
    469661        /*If on water, return 0: */
    470         if(this->properties.onwater)return 0;
     662        if(onwater)return 0;
    471663
    472664        /*Bail out if this element if:
    473665         * -> Non collapsed and not on the surface
    474666         * -> collapsed (2d model) and not on bed) */
    475         if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
     667        if ((!collapse && !onsurface) || (collapse && !onbed)){
    476668                return 0;
    477669        }
    478         else if (this->properties.collapse){
     670        else if (collapse){
    479671
    480672                /*This element should be collapsed into a tria element at its base. Create this tria element,
    481673                 * and compute CostFunction*/
    482674                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    483                 J=tria->CostFunction(inputs,analysis_type,sub_analysis_type);
     675                J=tria->CostFunction(analysis_type,sub_analysis_type);
    484676                delete tria;
    485677                return J;
     
    488680
    489681                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    490                 J=tria->CostFunction(inputs,analysis_type,sub_analysis_type);
     682                J=tria->CostFunction(analysis_type,sub_analysis_type);
    491683                delete tria;
    492684                return J;
     
    501693        if (analysis_type==ControlAnalysisEnum){
    502694
    503                 CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
     695                CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type);
    504696
    505697        }
     
    508700                if (sub_analysis_type==HorizAnalysisEnum){
    509701
    510                         CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
     702                        CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type);
    511703                }
    512704                else if (sub_analysis_type==VertAnalysisEnum){
    513705
    514                         CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
     706                        CreateKMatrixDiagnosticVert( Kgg,analysis_type,sub_analysis_type);
    515707                }
    516708                else if (sub_analysis_type==StokesAnalysisEnum){
    517709
    518                         CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
     710                        CreateKMatrixDiagnosticStokes( Kgg,analysis_type,sub_analysis_type);
    519711
    520712                }
     
    523715        else if (analysis_type==SlopecomputeAnalysisEnum){
    524716
    525                 CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
     717                CreateKMatrixSlopeCompute( Kgg,analysis_type,sub_analysis_type);
    526718        }
    527719        else if (analysis_type==PrognosticAnalysisEnum){
    528720
    529                 CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
     721                CreateKMatrixPrognostic( Kgg,analysis_type,sub_analysis_type);
    530722        }
    531723        else if (analysis_type==BalancedthicknessAnalysisEnum){
     
    539731        else if (analysis_type==ThermalAnalysisEnum){
    540732
    541                 CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
     733                CreateKMatrixThermal( Kgg,analysis_type,sub_analysis_type);
    542734        }
    543735        else if (analysis_type==MeltingAnalysisEnum){
    544736
    545                 CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
     737                CreateKMatrixMelting( Kgg,analysis_type,sub_analysis_type);
    546738        }
    547739        else{
     
    645837        double Bprime[5][numdof];
    646838        double L[2][numdof];
    647         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.
     839        double D[5][5]={0.0};            // material matrix, simple scalar matrix.
    648840        double D_scalar;
    649         double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
     841        double DL[2][2]={0.0}; //for basal drag
    650842        double DL_scalar;
    651843
    652844        /* local element matrices: */
    653         double Ke_gg[numdof][numdof]; //local element stiffness matrix
     845        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
     846
    654847        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    655848        double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
     
    657850
    658851        /*slope: */
    659         double  slope[2]={0.0,0.0};
     852        double  slope[2]={0.0};
    660853        double  slope_magnitude;
    661 
    662         /*input parameters for structural analysis (diagnostic): */
    663         double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    664         double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    665         int     dofs[2]={0,1};
    666         double  thickness;
    667854
    668855        /*friction: */
     
    672859        double MAXSLOPE=.06; // 6 %
    673860        double MOUNTAINKEXPONENT=10;
    674 
    675         ParameterInputs* inputs=NULL;
    676861
    677862        /*Collapsed formulation: */
     
    688873        numpar=(Numpar*)hnumpar.delivers();
    689874
     875        /*inputs: */
     876        bool onwater;
     877        bool collapse;
     878        bool onbed;
     879        bool shelf;
     880
     881        /*retrieve inputs :*/
     882        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     883        inputs->GetParameterValue(&collapse,CollapseEnum);
     884        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     885        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     886
    690887        /*If on water, skip stiffness: */
    691         if(this->properties.onwater)return;
    692 
    693         /*recover pointers: */
    694         inputs=(ParameterInputs*)vinputs;
     888        if(onwater)return;
    695889
    696890        /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
     
    699893
    700894
    701         if ((this->properties.collapse==1) && (this->properties.onbed==0)){
     895        if ((collapse==1) && (onbed==0)){
    702896                /*This element should be collapsed, but this element is not on the bedrock, therefore all its
    703897                 * dofs have already been frozen! Do nothing: */
    704898                return;
    705899        }
    706         else if ((this->properties.collapse==1) && (this->properties.onbed==1)){
     900        else if ((collapse==1) && (onbed==1)){
    707901
    708902                /*This element should be collapsed into a tria element at its base. Create this tria element,
    709903                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
    710904                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    711                 tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
     905                tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type);
    712906                delete tria;
    713907                return;
     
    716910
    717911                /*Implement standard penta element: */
    718 
    719                 /*recover extra inputs from users, at current convergence iteration: */
    720                 inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    721                 inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    722912
    723913                /* Get node coordinates and dof list: */
     
    725915                GetDofList(&doflist[0],&numberofdofspernode);
    726916
    727                 /* Set Ke_gg to 0: */
    728                 for(i=0;i<numdof;i++){
    729                         for(j=0;j<numdof;j++){
    730                                 Ke_gg[i][j]=0.0;
    731                         }
    732                 }
    733917
    734918                /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     
    760944
    761945                                /*Get strain rate from velocity: */
    762                                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_coord);
    763                                 GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_coord);
     946                                inputs->GetStrainRate(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum);
     947                                inputs->GetStrainRate(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum);
    764948
    765949                                /*Get viscosity: */
     
    774958                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    775959
    776                                 /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     960                                /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
    777961                                  onto this scalar matrix, so that we win some computational time: */
    778962
     
    803987
    804988                //Deal with 2d friction at the bedrock interface
    805                 if((this->properties.onbed && !this->properties.shelf)){
     989                if((onbed && !shelf)){
    806990
    807991                        /*Build a tria element using the 3 grids of the base of the penta. Then use
     
    810994
    811995                        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    812                         tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
     996                        tria->CreateKMatrixDiagnosticHorizFriction(Kgg,analysis_type,sub_analysis_type);
    813997                        delete tria;
    814998                }
     
    8451029        int   dofs[3]={0,1,2};
    8461030
    847         double K_terms[numdof][numdof];
     1031        double K_terms[numdof][numdof]={0.0};
    8481032
    8491033        /*Material properties: */
     
    8601044        double             surface_normal[3];
    8611045        double             bed_normal[3];
    862         double         vxvyvz_list[numgrids][3];
    8631046        double         thickness;
    8641047
     
    8891072        double* area_gauss_weights  =  NULL;
    8901073        double* vert_gauss_weights  =  NULL;
     1074        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}};
    8911075
    8921076        /* specific gaussian point: */
     
    9021086        double  alpha2_gauss;
    9031087
    904         ParameterInputs* inputs=NULL;
     1088        double  vx_list[numgrids];
     1089        double  vy_list[numgrids];
     1090        double  vz_list[numgrids];
     1091        double  thickness_list[numgrids];
     1092        double  bed_list[numgrids];
     1093        double  dragcoefficient_list[numgrids];
     1094        double  drag_p,drag_q;
    9051095
    9061096        /*dynamic objects pointed to by hooks: */
     
    9101100        Numpar* numpar=NULL;
    9111101
     1102        /*inputs: */
     1103        bool onwater;
     1104        bool onbed;
     1105        bool shelf;
     1106
     1107        /*retrieve inputs :*/
     1108        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1109        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1110        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     1111
    9121112        /*If on water, skip stiffness: */
    913         if(this->properties.onwater)return;
     1113        if(onwater)return;
    9141114
    9151115        /*recover objects from hooks: */
     
    9191119        numpar=(Numpar*)hnumpar.delivers();
    9201120
    921         /*recover pointers: */
    922         inputs=(ParameterInputs*)vinputs;
    923 
    924         /* Set K_terms  to 0: */
    925         for(i=0;i<numdof;i++){
    926                 for(j=0;j<numdof;j++){
    927                         K_terms[i][j]=0.0;
    928                 }
    929         }
    9301121
    9311122        /*recovre material parameters: */
     
    9371128        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    9381129        GetDofList(&doflist[0],&numberofdofspernode);
    939 
    940         /*recover extra inputs from users, at current convergence iteration: */
    941         inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
    9421130
    9431131        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     
    9651153                        gauss_coord[3]=*(vert_gauss_coord+igvert);
    9661154
    967                         /*Compute thickness: */
    968 
    9691155                        /*Compute strain rate: */
    970                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     1156                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
    9711157
    9721158                        /*Get viscosity: */
     
    10031189        }
    10041190
    1005         if((this->properties.onbed==1) && (this->properties.shelf==0)){
     1191        if((onbed==1) && (shelf==0)){
    10061192
    10071193                /*Build alpha2_list used by drag stiffness matrix*/
     
    10121198                strcpy(friction->element_type,"2d");
    10131199
     1200                inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],6,VxEnum);
     1201                inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],6,VyEnum);
     1202                inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],6,VzEnum);
     1203                inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],6,DragCoefficientEnum);
     1204                inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],6,BedEnum);
     1205                inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],6,ThicknessEnum);
     1206                inputs->GetParameterValue(&drag_p,DragPEnum);
     1207                inputs->GetParameterValue(&drag_q,DragQEnum);
     1208
    10141209                friction->gravity=matpar->GetG();
    10151210                friction->rho_ice=matpar->GetRhoIce();
    10161211                friction->rho_water=matpar->GetRhoWater();
    1017                 friction->K=&this->properties.k[0];
    1018                 friction->bed=&this->properties.b[0];
    1019                 friction->thickness=&this->properties.h[0];
    1020                 friction->velocities=&vxvyvz_list[0][0];
    1021                 friction->p=this->properties.p;
    1022                 friction->q=this->properties.q;
    1023                 friction->analysis_type=analysis_type;
     1212                friction->K=&dragcoefficient_list[0];
     1213                friction->bed=&bed_list[0];
     1214                friction->thickness=&thickness_list[0];
     1215                friction->vx=&vx_list[0];
     1216                friction->vy=&vy_list[0];
     1217                friction->vz=&vz_list[0];
     1218                friction->p=drag_p;
     1219                friction->q=drag_q;
    10241220
    10251221                /*Compute alpha2_list: */
     
    10581254
    10591255                        /*Compute strain rate: */
    1060                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     1256                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
    10611257
    10621258                        /*Get viscosity at last iteration: */
     
    10981294                        }
    10991295                }
    1100         } //if ( (this->properties.onbed==1) && (shelf==0))
     1296        } //if ( (onbed==1) && (shelf==0))
    11011297
    11021298        /*Reduce the matrix */
     
    11551351
    11561352        /* matrices: */
    1157         double  Ke_gg[numdof][numdof];
     1353        double  Ke_gg[numdof][numdof]={0.0};
    11581354        double  Ke_gg_gaussian[numdof][numdof];
    11591355        double  Jdet;
     
    11681364        nodes=(Node**)hnodes.deliverp();
    11691365
    1170         ParameterInputs* inputs=NULL;
    1171 
    11721366        /*Collapsed formulation: */
    11731367        Tria*  tria=NULL;
    11741368
     1369        /*inputs: */
     1370        bool onwater;
     1371        bool onsurface;
     1372
     1373        /*retrieve inputs :*/
     1374        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1375        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     1376
    11751377        /*If on water, skip stiffness: */
    1176         if(this->properties.onwater)return;
    1177 
    1178         /*recover pointers: */
    1179         inputs=(ParameterInputs*)vinputs;
    1180 
     1378        if(onwater)return;
    11811379
    11821380        /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness
    11831381         * matrix: */
    1184         if(this->properties.onsurface){
     1382        if(onsurface){
    11851383                tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
    1186                 tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type);
     1384                tria->CreateKMatrixDiagnosticSurfaceVert(Kgg, analysis_type,sub_analysis_type);
    11871385                delete tria;
    11881386        }
     
    11941392        GetDofList(&doflist[0],&numberofdofspernode);
    11951393
    1196         /* Set Ke_gg to 0: */
    1197         for(i=0;i<numdof;i++){
    1198                 for(j=0;j<numdof;j++){
    1199                         Ke_gg[i][j]=0.0;
    1200                 }
    1201         }
    12021394
    12031395        /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     
    12651457        Tria* tria=NULL;
    12661458
     1459        /*inputs: */
     1460        bool onwater;
     1461        bool onbed;
     1462
     1463        /*retrieve inputs :*/
     1464        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1465        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1466
    12671467        /*If on water, skip: */
    1268         if(this->properties.onwater)return;
    1269 
    1270         if (!this->properties.onbed){
     1468        if(onwater)return;
     1469
     1470        if (!onbed){
    12711471                return;
    12721472        }
     
    12741474
    12751475                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1276                 tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
     1476                tria->CreateKMatrixMelting(Kgg, analysis_type,sub_analysis_type);
    12771477                delete tria;
    12781478                return;
     
    12871487        Tria*  tria=NULL;
    12881488
     1489        /*inputs: */
     1490        bool onwater;
     1491        bool onbed;
     1492
     1493        /*retrieve inputs :*/
     1494        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1495        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1496
    12891497        /*If on water, skip: */
    1290         if(this->properties.onwater)return;
     1498        if(onwater)return;
    12911499
    12921500        /*Is this element on the bed? :*/
    1293         if(!this->properties.onbed)return;
     1501        if(!onbed)return;
    12941502
    12951503        /*Spawn Tria element from the base of the Penta: */
    12961504        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1297         tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
     1505        tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type);
    12981506        delete tria;
    12991507        return;
     
    13081516        Tria*  tria=NULL;
    13091517
     1518        /*inputs: */
     1519        bool onwater;
     1520        bool onbed;
     1521
     1522        /*retrieve inputs :*/
     1523        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1524        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1525
     1526
    13101527        /*If on water, skip: */
    1311         if(this->properties.onwater)return;
     1528        if(onwater)return;
    13121529
    13131530        /*Is this element on the bed? :*/
    1314         if(!this->properties.onbed)return;
     1531        if(!onbed)return;
    13151532
    13161533        /*Spawn Tria element from the base of the Penta: */
    13171534        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1318         tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
     1535        tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type);
    13191536        delete tria;
    13201537        return;
     
    13561573        double  K[2][2]={0.0};
    13571574
    1358         double  vxvyvz_list[numgrids][3];
    1359         double  vx_list[numgrids];
    1360         int     vx_list_exists;
    1361         double  u;
    1362         double  vy_list[numgrids];
    1363         int     vy_list_exists;
    1364         double  v;
    1365         double  vz_list[numgrids];
    1366         int     vz_list_exists;
    1367         double  w;
     1575        double  u,v,w;
    13681576
    13691577        /*matrices: */
     
    14041612        /*Collapsed formulation: */
    14051613        Tria*  tria=NULL;
    1406         ParameterInputs* inputs=NULL;
     1614
     1615
     1616        /*inputs: */
     1617        bool onwater;
     1618        bool onbed;
     1619        bool shelf;
     1620
     1621        /*retrieve inputs :*/
     1622        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1623        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1624        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    14071625
    14081626        /*If on water, skip: */
    1409         if(this->properties.onwater)return;
     1627        if(onwater)return;
    14101628
    14111629        /*recover objects from hooks: */
     
    14131631        matpar=(Matpar*)hmatpar.delivers();
    14141632        numpar=(Numpar*)hnumpar.delivers();
    1415 
    1416         /*recover pointers: */
    1417         inputs=(ParameterInputs*)vinputs;
    14181633
    14191634        /* Get node coordinates and dof list: */
     
    14271642        heatcapacity=matpar->GetHeatCapacity();
    14281643        thermalconductivity=matpar->GetThermalConductivity();
    1429 
    1430         /*recover extra inputs from users, dt and velocity: */
    1431         found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
    1432         if(!found)ISSMERROR(" could not find velocity in inputs!");
    1433         found=inputs->Recover("dt",&dt);
    1434         if(!found)ISSMERROR(" could not find dt in inputs!");
    1435 
    1436         for(i=0;i<numgrids;i++){
    1437                 vx_list[i]=vxvyvz_list[i][0];
    1438                 vy_list[i]=vxvyvz_list[i][1];
    1439                 vz_list[i]=vxvyvz_list[i][2];
    1440         }
    1441 
     1644        dt=numpar->dt;
    14421645
    14431646        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     
    14951698
    14961699                        //Build the D matrix
    1497                         GetParameterValue(&u, &vx_list[0],gauss_coord);
    1498                         GetParameterValue(&v, &vy_list[0],gauss_coord);
    1499                         GetParameterValue(&w, &vz_list[0],gauss_coord);
     1700                        inputs->GetParameterValue(&u, gauss_coord,VxEnum);
     1701                        inputs->GetParameterValue(&v, gauss_coord,VyEnum);
     1702                        inputs->GetParameterValue(&w, gauss_coord,VzEnum);
    15001703
    15011704                        D_scalar=gauss_weight*Jdet;
     
    15771780
    15781781        //Ice/ocean heat exchange flux on ice shelf base
    1579         if(this->properties.onbed && this->properties.shelf){
     1782        if(onbed && shelf){
    15801783
    15811784                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1582                 tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
     1785                tria->CreateKMatrixThermal(Kgg, analysis_type,sub_analysis_type);
    15831786                delete tria;
    15841787        }
     
    15911794        if (analysis_type==ControlAnalysisEnum){
    15921795
    1593                 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     1796                CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type);
    15941797        }
    15951798        else if (analysis_type==DiagnosticAnalysisEnum){
     
    15971800                if (sub_analysis_type==HorizAnalysisEnum){
    15981801
    1599                         CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     1802                        CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type);
    16001803                }
    16011804                else if (sub_analysis_type==VertAnalysisEnum){
    16021805
    1603                         CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
     1806                        CreatePVectorDiagnosticVert( pg,analysis_type,sub_analysis_type);
    16041807                }
    16051808                else if (sub_analysis_type==StokesAnalysisEnum){
    16061809
    1607                         CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
     1810                        CreatePVectorDiagnosticStokes( pg,analysis_type,sub_analysis_type);
    16081811                }
    16091812                else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
     
    16111814        else if (analysis_type==SlopecomputeAnalysisEnum){
    16121815
    1613                 CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
     1816                CreatePVectorSlopeCompute( pg,analysis_type,sub_analysis_type);
    16141817        }
    16151818        else if (analysis_type==PrognosticAnalysisEnum){
    16161819
     1820                CreatePVectorPrognostic( pg,analysis_type,sub_analysis_type);
     1821        }
     1822        else if (analysis_type==BalancedthicknessAnalysisEnum){
     1823
    16171824                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    16181825        }
    1619         else if (analysis_type==BalancedthicknessAnalysisEnum){
     1826        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    16201827
    16211828                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    16221829        }
    1623         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    1624 
    1625                 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    1626         }
    16271830        else if (analysis_type==ThermalAnalysisEnum){
    16281831
    1629                 CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
     1832                CreatePVectorThermal( pg,analysis_type,sub_analysis_type);
    16301833        }
    16311834        else if (analysis_type==MeltingAnalysisEnum){
    16321835
    1633                 CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
     1836                CreatePVectorMelting( pg,analysis_type,sub_analysis_type);
    16341837        }
    16351838        else{
     
    17201923
    17211924        /*element vector at the gaussian points: */
    1722         double  pe_g[numdof];
     1925        double  pe_g[numdof]={0.0};
    17231926        double  pe_g_gaussian[numdof];
    1724 
    1725         ParameterInputs* inputs=NULL;
    17261927
    17271928        /*dynamic objects pointed to by hooks: */
     
    17321933        Tria* tria=NULL;
    17331934
     1935        /*inputs: */
     1936        bool onwater;
     1937        bool collapse;
     1938        bool onbed;
     1939
     1940        /*retrieve inputs :*/
     1941        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1942        inputs->GetParameterValue(&collapse,CollapseEnum);
     1943        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1944
    17341945        /*If on water, skip load: */
    1735         if(this->properties.onwater)return;
    1736 
    1737         /*recover pointers: */
    1738         inputs=(ParameterInputs*)vinputs;
     1946        if(onwater)return;
    17391947
    17401948        /*recover objects from hooks: */
     
    17471955          the load vector. */
    17481956
    1749         if ((this->properties.collapse==1) && (this->properties.onbed==0)){
     1957        if ((collapse==1) && (onbed==0)){
    17501958                /*This element should be collapsed, but this element is not on the bedrock, therefore all its
    17511959                 * dofs have already been frozen! Do nothing: */
    17521960                return;
    17531961        }
    1754         else if ((this->properties.collapse==1) && (this->properties.onbed==1)){
     1962        else if ((collapse==1) && (onbed==1)){
    17551963
    17561964                /*This element should be collapsed into a tria element at its base. Create this tria element,
    17571965                 *and use its CreatePVector functionality to return an elementary load vector: */
    17581966                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1759                 tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
     1967                tria->CreatePVector(pg, analysis_type,sub_analysis_type);
    17601968                delete tria;
    17611969                return;
     
    17681976                GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    17691977                GetDofList(&doflist[0],&numberofdofspernode);
    1770 
    1771                 /* Set pe_g to 0: */
    1772                 for(i=0;i<numdof;i++) pe_g[i]=0.0;
    17731978
    17741979                /*Get gaussian points and weights :*/
     
    17931998
    17941999                                /*Compute thickness at gaussian point: */
    1795                                 GetParameterValue(&thickness, &this->properties.h[0],gauss_coord);
     2000                                inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum);
    17962001
    17972002                                /*Compute slope at gaussian point: */
    1798                                 GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_coord);
     2003                                inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum);
    17992004
    18002005                                /* Get Jacobian determinant: */
     
    18202025                } //for (ig1=0; ig1<num_area_gauss; ig1++)
    18212026
    1822         } //else if ((collapse==1) && (this->properties.onbed==1))
     2027        } //else if ((collapse==1) && (onbed==1))
    18232028
    18242029        /*Add pe_g to global vector pg: */
     
    18572062        double             bed_normal[3];
    18582063        double         bed;
    1859         double         vxvyvz_list[numgrids][3];
    18602064
    18612065        /* gaussian points: */
     
    18982102        double     D_scalar;
    18992103        double     tBD[27][8];
    1900         double     P_terms[numdof];
    1901 
    1902         ParameterInputs* inputs=NULL;
     2104        double     P_terms[numdof]={0.0};
     2105
    19032106        Tria*            tria=NULL;
    19042107
     
    19092112        Numpar* numpar=NULL;
    19102113
     2114        /*inputs: */
     2115        bool onwater;
     2116        bool onbed;
     2117        bool shelf;
     2118
     2119        /*retrieve inputs :*/
     2120        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2121        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2122        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     2123
    19112124        /*If on water, skip load: */
    1912         if(this->properties.onwater)return;
    1913 
    1914         /*recover pointers: */
    1915         inputs=(ParameterInputs*)vinputs;
     2125        if(onwater)return;
    19162126
    19172127        /*recover objects from hooks: */
     
    19212131        numpar=(Numpar*)hnumpar.delivers();
    19222132
    1923         /* Set P_terms to 0: */
    1924         for(i=0;i<numdof;i++){
    1925                 P_terms[i]=0;
    1926         }       
    19272133
    19282134        /*recovre material parameters: */
     
    19302136        rho_ice=matpar->GetRhoIce();
    19312137        gravity=matpar->GetG();
    1932 
    1933         /*recover extra inputs from users, at current convergence iteration: */
    1934         inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
    19352138
    19362139        /* Get node coordinates and dof list: */
     
    19602163
    19612164                        /*Compute strain rate and viscosity: */
    1962                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     2165                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
    19632166                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    19642167
     
    20142217
    20152218        /*Deal with 2d friction at the bedrock interface: */
    2016         if ( (this->properties.onbed==1) && (this->properties.shelf==1)){
     2219        if ( (onbed==1) && (shelf==1)){
    20172220
    20182221                for(i=0;i<numgrids2d;i++){
     
    20412244
    20422245                        /* Get bed at gaussian point */
    2043                         GetParameterValue(&bed,&this->properties.b[0],gauss_coord);
     2246                        inputs->GetParameterValue(&bed, gauss_coord,BedEnum);
    20442247
    20452248                        /*Get L matrix : */
     
    20622265                        }
    20632266                }
    2064         } //if ( (this->properties.onbed==1) && (this->properties.shelf==1))
     2267        } //if ( (onbed==1) && (shelf==1))
    20652268
    20662269        /*Reduce the matrix */
     
    21182321
    21192322        /*element vector at the gaussian points: */
    2120         double  pe_g[numdof];
     2323        double  pe_g[numdof]={0.0};
    21212324        double  pe_g_gaussian[numdof];
    21222325        double l1l6[6];
     
    21262329
    21272330        /*input parameters for structural analysis (diagnostic): */
    2128         double vx_list[numgrids]={0,0,0,0,0,0};
    2129         double vy_list[numgrids]={0,0,0,0,0,0};
    21302331        double du[3];
    21312332        double dv[3];
     
    21372338        Node**  nodes=NULL;
    21382339
    2139         ParameterInputs* inputs=NULL;
     2340        /*inputs: */
     2341        bool onwater;
     2342        bool onbed;
    21402343
    21412344        /*recover objects from hooks: */
    21422345        nodes=(Node**)hnodes.deliverp();
    21432346
    2144         /*recover pointers: */
    2145         inputs=(ParameterInputs*)vinputs;
    2146 
     2347        /*retrieve inputs :*/
     2348        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2349        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2350
     2351       
    21472352        /*If on water, skip: */
    2148         if(this->properties.onwater)return;
     2353        if(onwater)return;
    21492354
    21502355        /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the
    21512356         *diagnostic base vertical stifness: */
    2152         if(this->properties.onbed){
     2357        if(onbed){
    21532358                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
    2154                 tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type,sub_analysis_type);
     2359                tria->CreatePVectorDiagnosticBaseVert(pg, analysis_type,sub_analysis_type);
    21552360                delete tria;
    21562361        }
     
    21612366        GetDofList(&doflist[0],&numberofdofspernode);
    21622367
    2163         /* Set pe_g to 0: */
    2164         for(i=0;i<numdof;i++) pe_g[i]=0.0;
    2165 
    2166         /* recover input parameters: */
    2167         if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))ISSMERROR(" cannot compute vertical velocity without horizontal velocity!");
    2168         inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
    2169 
    21702368        /*Get gaussian points and weights :*/
    21712369        order_area_gauss=2;
     
    21892387
    21902388                        /*Get velocity derivative, with respect to x and y: */
    2191                         GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_coord);
    2192                         GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_coord);
     2389
     2390                        inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum);
     2391                        inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VyEnum);
    21932392                        dudx=du[0];
    21942393                        dvdy=dv[1];
     
    22352434        Tria*  tria=NULL;
    22362435
     2436        /*inputs: */
     2437        bool onwater;
     2438        bool onbed;
     2439
     2440        /*retrieve inputs :*/
     2441        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2442        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2443       
    22372444        /*If on water, skip: */
    2238         if(this->properties.onwater)return;
     2445        if(onwater)return;
    22392446
    22402447        /*Is this element on the bed? :*/
    2241         if(!this->properties.onbed)return;
     2448        if(!onbed)return;
    22422449
    22432450        /*Spawn Tria element from the base of the Penta: */
    22442451        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2245         tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
     2452        tria->CreatePVector(pg, analysis_type,sub_analysis_type);
    22462453        delete tria;
    22472454        return;
     
    22552462        Tria*  tria=NULL;
    22562463
     2464        /*inputs: */
     2465        bool onwater;
     2466        bool onbed;
     2467
     2468        /*retrieve inputs :*/
     2469        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2470        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2471
    22572472        /*If on water, skip: */
    2258         if(this->properties.onwater)return;
     2473        if(onwater)return;
    22592474
    22602475        /*Is this element on the bed? :*/
    2261         if(!this->properties.onbed)return;
     2476        if(!onbed)return;
    22622477
    22632478        /*Spawn Tria element from the base of the Penta: */
    22642479        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2265         tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
     2480        tria->CreatePVector(pg, analysis_type,sub_analysis_type);
    22662481        delete tria;
    22672482        return;
     
    22992514
    23002515        double dt;
    2301         double vx_list[numgrids];
    2302         double vy_list[numgrids];
    2303         double vz_list[numgrids];
    2304         double vxvyvz_list[numgrids][3];
    23052516        double temperature_list[numgrids];
    23062517        double temperature;
     
    23402551        /*Collapsed formulation: */
    23412552        Tria*  tria=NULL;
    2342         ParameterInputs* inputs=NULL;
    23432553
    23442554        /*dynamic objects pointed to by hooks: */
     
    23462556        Matpar* matpar=NULL;
    23472557        Matice* matice=NULL;
     2558        Numpar* numpar=NULL;
     2559
     2560        /*inputs: */
     2561        bool onwater;
     2562        bool onbed;
     2563        bool shelf;
     2564
     2565        /*retrieve inputs :*/
     2566        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2567        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2568        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    23482569
    23492570        /*If on water, skip: */
    2350         if(this->properties.onwater)return;
    2351 
    2352         /*recover pointers: */
    2353         inputs=(ParameterInputs*)vinputs;
     2571        if(onwater)return;
    23542572
    23552573        /*recover objects from hooks: */
     
    23572575        matpar=(Matpar*)hmatpar.delivers();
    23582576        matice=(Matice*)hmatice.delivers();
     2577        numpar=(Numpar*)hnumpar.delivers();
    23592578
    23602579        /* Get node coordinates and dof list: */
     
    23692588        beta=matpar->GetBeta();
    23702589        meltingpoint=matpar->GetMeltingPoint();
    2371 
    2372         /*recover extra inputs from users, dt and velocity: */
    2373         found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
    2374         if(!found)ISSMERROR(" could not find velocity in inputs!");
    2375         found=inputs->Recover("dt",&dt);
    2376         if(!found)ISSMERROR(" could not find dt in inputs!");
    2377 
    2378         if(dt){
    2379                 found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes);
    2380                 if(!found)ISSMERROR(" could not find temperature in inputs!");
    2381         }
    2382 
    2383         for(i=0;i<numgrids;i++){
    2384                 vx_list[i]=vxvyvz_list[i][0];
    2385                 vy_list[i]=vxvyvz_list[i][1];
    2386                 vz_list[i]=vxvyvz_list[i][2];
    2387         }       
     2590        dt=numpar->dt;
    23882591
    23892592        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     
    24092612
    24102613                        /*Compute strain rate and viscosity: */
    2411                         GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
     2614                        inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
    24122615                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    24132616
     
    24332636                        /* Build transient now */
    24342637                        if(dt){
    2435                                 GetParameterValue(&temperature, &temperature_list[0],gauss_coord);
     2638                                inputs->GetParameterValue(&temperature, gauss_coord,TemperatureEnum);
    24362639                                scalar_transient=temperature*Jdet*gauss_weight;
    24372640                                for(i=0;i<numgrids;i++){
     
    24462649
    24472650        /* Ice/ocean heat exchange flux on ice shelf base */
    2448         if(this->properties.onbed && this->properties.shelf){
     2651        if(onbed && shelf){
    24492652
    24502653                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2451                 tria->CreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type);
     2654                tria->CreatePVectorThermalShelf(pg, analysis_type,sub_analysis_type);
    24522655                delete tria;
    24532656        }
    24542657
    24552658        /* Geothermal flux on ice sheet base and basal friction */
    2456         if(this->properties.onbed && !this->properties.shelf){
     2659        if(onbed && !shelf){
    24572660
    24582661                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2459                 tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type);
     2662                tria->CreatePVectorThermalSheet(pg, analysis_type,sub_analysis_type);
    24602663                delete tria;
    24612664        }
     
    24772680        Tria* tria=NULL;
    24782681
     2682        /*inputs: */
     2683        bool onwater;
     2684        bool collapse;
     2685        bool onsurface;
     2686        bool onbed;
     2687
     2688        /*retrieve inputs :*/
     2689        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2690        inputs->GetParameterValue(&collapse,CollapseEnum);
     2691        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2692        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     2693
    24792694        /*If on water, skip: */
    2480         if(this->properties.onwater)return;
     2695        if(onwater)return;
    24812696
    24822697        /*Bail out if this element if:
    24832698         * -> Non collapsed and not on the surface
    24842699         * -> collapsed (2d model) and not on bed) */
    2485         if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
     2700        if ((!collapse && !onsurface) || (collapse && !onbed)){
    24862701                return;
    24872702        }
    2488         else if (this->properties.collapse){
     2703        else if (collapse){
    24892704
    24902705                /*This element should be collapsed into a tria element at its base. Create this tria element,
    24912706                 * and compute Du*/
    24922707                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    2493                 tria->Du(du_g,inputs,analysis_type,sub_analysis_type);
     2708                tria->Du(du_g,analysis_type,sub_analysis_type);
    24942709                delete tria;
    24952710                return;
     
    24982713
    24992714                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    2500                 tria->Du(du_g,inputs,analysis_type,sub_analysis_type);
     2715                tria->Du(du_g,analysis_type,sub_analysis_type);
    25012716                delete tria;
    25022717                return;
    25032718        }
    2504 }
    2505 /*}}}*/
    2506 /*FUNCTION Enum {{{1*/
    2507 int Penta::Enum(void){
    2508 
    2509         return PentaEnum;
    2510 
    25112719}
    25122720/*}}}*/
     
    25272735        nodes=(Node**)hnodes.deliverp();
    25282736
     2737        /*inputs: */
     2738        bool collapse;
     2739        bool onbed;
     2740
     2741        /*retrieve inputs :*/
     2742        inputs->GetParameterValue(&collapse,CollapseEnum);
     2743        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2744
    25292745        /*Figure out if we should extrude for this element: */
    25302746        if (iscollapsed){
    25312747                /*From higher level, we are told to extrude only elements that have the collapse flag on: */
    2532                 if (this->properties.collapse)extrude=1;
     2748                if (collapse)extrude=1;
    25332749                else extrude=0;
    25342750        }
     
    25402756        /*Now, extrusion starts from the bed on, so double check this element is on
    25412757         * the bedrock: */
    2542         if(this->properties.onbed==0)extrude=0;
     2758        if(onbed==0)extrude=0;
    25432759
    25442760        /*Go on and extrude field: */
     
    28223038/*}}}*/
    28233039/*FUNCTION GetBedList {{{1*/
    2824 void Penta::GetBedList(double* bed_list){
    2825 
    2826         int i;
    2827         for(i=0;i<6;i++)bed_list[i]=this->properties.b[i];
     3040void Penta::GetBedList(double* bedlist){
     3041
     3042        const int numgrids=6;
     3043        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}};
     3044       
     3045        inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum);
    28283046
    28293047}
     
    36823900/*FUNCTION GetOnBed {{{1*/
    36833901bool Penta::GetOnBed(){
    3684         return this->properties.onbed;
     3902       
     3903        bool onbed;
     3904
     3905        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     3906
     3907        return onbed;
    36853908}
    36863909/*}}}*/
     
    37573980/*}}}*/
    37583981/*FUNCTION GetShelf {{{1*/
    3759 int   Penta::GetShelf(){
    3760         return this->properties.shelf;
     3982bool   Penta::GetShelf(){
     3983        bool onshelf;
     3984
     3985        /*retrieve inputs :*/
     3986        inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
     3987
     3988        return onshelf;
    37613989}
    37623990/*}}}*/
     
    38214049void Penta::GetThicknessList(double* thickness_list){
    38224050
    3823         int i;
    3824         for(i=0;i<6;i++)thickness_list[i]=this->properties.h[i];
     4051        const int numgrids=6;
     4052        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}};
     4053        inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum);
     4054
    38254055}
    38264056/*}}}*/
     
    38284058void  Penta::Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type){
    38294059
     4060        /*inputs: */
     4061        bool onwater;
     4062
     4063        /*retrieve inputs :*/
     4064        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4065
    38304066        /*If on water, skip grad (=0): */
    3831         if(this->properties.onwater)return;
     4067        if(onwater)return;
    38324068
    38334069        if (strcmp(control_type,"drag")==0){
    3834                 GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type);
     4070                GradjDrag( grad_g,analysis_type,sub_analysis_type);
    38354071        }
    38364072        else if (strcmp(control_type,"B")==0){
    3837                 GradjB( grad_g, inputs,analysis_type,sub_analysis_type);
     4073                GradjB( grad_g, analysis_type,sub_analysis_type);
    38384074        }
    38394075        else ISSMERROR("%s%s","control type not supported yet: ",control_type);
     
    38454081        Tria* tria=NULL;
    38464082
     4083        /*inputs: */
     4084        bool onwater;
     4085        bool onbed;
     4086        bool shelf;
     4087
     4088        /*retrieve inputs :*/
     4089        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4090        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     4091        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     4092
    38474093        /*If on water, skip: */
    3848         if(this->properties.onwater)return;
     4094        if(onwater)return;
    38494095
    38504096        /*If on shelf, skip: */
    3851         if(this->properties.shelf)return;
     4097        if(shelf)return;
    38524098
    38534099        /*Bail out if this element does not touch the bedrock: */
    3854         if (!this->properties.onbed) return;
     4100        if (!onbed) return;
    38554101
    38564102        if (sub_analysis_type==HorizAnalysisEnum){
     
    38584104                /*MacAyeal or Pattyn*/
    38594105                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3860                 tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type);
     4106                tria->GradjDrag( grad_g,analysis_type,sub_analysis_type);
    38614107                delete tria;
    38624108                return;
     
    38664112                /*Stokes*/
    38674113                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3868                 tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type);
     4114                tria->GradjDragStokes( grad_g,analysis_type,sub_analysis_type);
    38694115                delete tria;
    38704116                return;
     
    38784124        Tria* tria=NULL;
    38794125
     4126        /*inputs: */
     4127        bool onwater;
     4128        bool collapse;
     4129        bool onbed;
     4130
     4131        /*retrieve inputs :*/
     4132        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4133        inputs->GetParameterValue(&collapse,CollapseEnum);
     4134        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     4135
    38804136        /*If on water, skip: */
    3881         if(this->properties.onwater)return;
    3882 
    3883         if (this->properties.collapse){
     4137        if(onwater)return;
     4138
     4139        if (collapse){
    38844140                /*Bail out element if collapsed (2d) and not on bed*/
    3885                 if (!this->properties.onbed) return;
     4141                if (!onbed) return;
    38864142
    38874143                /*This element should be collapsed into a tria element at its base. Create this tria element,
    38884144                 * and compute gardj*/
    38894145                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    3890                 tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type);
     4146                tria->GradjB(grad_g,analysis_type,sub_analysis_type);
    38914147                delete tria;
    38924148                return;
     
    38954151                /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
    38964152                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    3897                 tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type);
     4153                tria->GradjB(grad_g,analysis_type,sub_analysis_type);
    38984154                delete tria;
    38994155                return;
     
    39124168        Tria* tria=NULL;
    39134169
     4170        /*inputs: */
     4171        bool onwater;
     4172        bool collapse;
     4173        bool onsurface;
     4174        bool onbed;
     4175
     4176        /*retrieve inputs :*/
     4177        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4178        inputs->GetParameterValue(&collapse,CollapseEnum);
     4179        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     4180        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     4181
    39144182        /*If on water, return 0: */
    3915         if(this->properties.onwater)return 0;
     4183        if(onwater)return 0;
    39164184
    39174185        /*Bail out if this element if:
    39184186         * -> Non collapsed and not on the surface
    39194187         * -> collapsed (2d model) and not on bed) */
    3920         if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
     4188        if ((!collapse && !onsurface) || (collapse && !onbed)){
    39214189                return 0;
    39224190        }
    3923         else if (this->properties.collapse){
     4191        else if (collapse){
    39244192
    39254193                /*This element should be collapsed into a tria element at its base. Create this tria element,
    39264194                 * and compute Misfit*/
    39274195                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    3928                 J=tria->Misfit(inputs,analysis_type,sub_analysis_type);
     4196                J=tria->Misfit(analysis_type,sub_analysis_type);
    39294197                delete tria;
    39304198                return J;
     
    39334201
    39344202                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    3935                 J=tria->Misfit(inputs,analysis_type,sub_analysis_type);
     4203                J=tria->Misfit(analysis_type,sub_analysis_type);
    39364204                delete tria;
    39374205                return J;
     
    40454313        Tria* tria=NULL;
    40464314
     4315        /*inputs: */
     4316        bool onwater;
     4317        bool collapse;
     4318        bool onsurface;
     4319        bool onbed;
     4320
     4321        /*retrieve inputs :*/
     4322        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4323        inputs->GetParameterValue(&collapse,CollapseEnum);
     4324        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     4325        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     4326
    40474327        /*If on water, return 0: */
    4048         if(this->properties.onwater)return 0;
     4328        if(onwater)return 0;
    40494329
    40504330        /*Bail out if this element if:
    40514331         * -> Non collapsed and not on the surface
    40524332         * -> collapsed (2d model) and not on bed) */
    4053         if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){
     4333        if ((!collapse && !onsurface) || (collapse && !onbed)){
    40544334                return 0;
    40554335        }
    4056         else if (this->properties.collapse){
     4336        else if (collapse){
    40574337
    40584338                /*This element should be collapsed into a tria element at its base. Create this tria element,
    40594339                 * and compute SurfaceArea*/
    40604340                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    4061                 S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type);
     4341                S=tria->SurfaceArea(analysis_type,sub_analysis_type);
    40624342                delete tria;
    40634343                return S;
     
    40664346
    40674347                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    4068                 S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type);
     4348                S=tria->SurfaceArea(analysis_type,sub_analysis_type);
    40694349                delete tria;
    40704350                return S;
  • issm/trunk/src/c/objects/Penta.h

    r3613 r3620  
    1616#include "./Element.h"
    1717#include "./Matpar.h"
    18 #include "./Numpar.h"
    1918#include "./Matice.h"
    2019#include "./Tria.h"
     
    4544                /*}}}*/
    4645                /*FUNCTION object management {{{1*/
    47                 void  Configure(void* loads,void* nodes,void* materials,void* parameters);
     46                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
    4847                Object* copy();
    4948                void  DeepEcho();
     
    5756                int   MyRank();
    5857                void*  SpawnTria(int g0, int g1, int g2);
    59                 void  UpdateFromDakota(void* inputs);
    60                 void  UpdateInputs(double* solution, int analysis_type, int sub_analysis_type);
    6158                void  SetClone(int* minranks);
    6259
     
    107104                void  CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type);
    108105                void  CreatePVectorPrognostic( Vec pg,  int analysis_type,int sub_analysis_type);
    109                 void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    110                 void  CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    111                 void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    112                 void  CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    113106
    114107                void  CreateKMatrixDiagnosticStokes( Mat Kgg,  int analysis_type,int sub_analysis_type);
     
    136129                void  GetPhi(double* phi, double*  epsilon, double viscosity);
    137130                double MassFlux(double* segment,double* ug);
     131
     132                /*updates: */
     133                void  UpdateFromDakota(void* inputs);
     134                void  UpdateInputs(double* solution, int analysis_type, int sub_analysis_type);
     135                void  UpdateInputsDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type);
     136                void  UpdateInputsSlopeCompute( double* solution,int analysis_type,int sub_analysis_type);
     137                void  UpdateInputsPrognostic( double* solution,int analysis_type,int sub_analysis_type);
     138                void  UpdateInputsPrognostic2(double* solution,int analysis_type,int sub_analysis_type);
     139                void  UpdateInputsBalancedthickness( double* solution,int analysis_type,int sub_analysis_type);
     140                void  UpdateInputsBalancedthickness2( double* solution,int analysis_type,int sub_analysis_type);
     141                void  UpdateInputsBalancedvelocities( double* solution,int analysis_type,int sub_analysis_type);
     142
    138143                /*}}}*/
    139144
  • issm/trunk/src/c/objects/Sing.cpp

    r3612 r3620  
    1111#include "stdio.h"
    1212#include "./Sing.h"
     13#include "./SingVertexInput.h"
    1314#include <string.h>
    1415#include "../EnumDefinitions/EnumDefinitions.h"
    1516#include "../shared/shared.h"
    1617#include "../DataSet/DataSet.h"
     18#include "../DataSet/Inputs.h"
    1719#include "../include/typedefs.h"
    1820#include "../include/macros.h"
     
    2123/*FUNCTION Sing::constructor {{{1*/
    2224Sing::Sing(){
     25        this->inputs=NULL;
    2326        return;
    2427}
    2528/*}}}*/
    2629/*FUNCTION Sing constructor {{{1*/
    27 Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id, ElementProperties* singproperties):
     30Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id):
    2831        hnodes(sing_node_ids,1),
    2932        hmatice(&sing_matice_id,1),
    3033        hmatpar(&sing_matpar_id,1),
    31         hnumpar(&sing_numpar_id,1),
    32         properties(singproperties)
     34        hnumpar(&sing_numpar_id,1)
    3335{
    3436
    3537        /*all the initialization has been done by the initializer, just fill in the id: */
    3638        this->id=sing_id;
     39        this->inputs=new Inputs();
    3740
    3841        return;
     
    4043/*}}}*/
    4144/*FUNCTION Sing other constructor {{{1*/
    42 Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, ElementProperties* sing_properties):
     45Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs):
    4346        hnodes(sing_hnodes),
    4447        hmatice(sing_hmatice),
    4548        hmatpar(sing_hmatpar),
    46         hnumpar(sing_hnumpar),
    47         properties(sing_properties)
     49        hnumpar(sing_hnumpar)
    4850{
    4951
    5052        /*all the initialization has been done by the initializer, just fill in the id: */
    5153        this->id=sing_id;
    52 
     54        if(sing_inputs){
     55                this->inputs=(Inputs*)sing_inputs->Copy();
     56        }
     57        else{
     58                this->inputs=new Inputs();
     59        }
    5360        return;
    5461}
     
    8087        this->hnumpar.Init(&sing_numpar_id,1);
    8188
    82         /*properties: */
    83         sing_h=iomodel->thickness[i];
    84         sing_k=iomodel->drag[i];
    85 
    86         this->properties.Init(1,&sing_h, NULL, NULL, &sing_k, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, UNDEF,UNDEF, UNDEF,UNDEF,UNDEF);
     89        //intialize inputs, and add as many inputs per element as requested:
     90        this->inputs=new Inputs();
     91
     92        if (iomodel->thickness) this->inputs->AddInput(new SingVertexInput(ThicknessEnum,iomodel->thickness[i]));
     93        if (iomodel->drag_coefficient) this->inputs->AddInput(new SingVertexInput(DragCoefficientEnum,iomodel->drag_coefficient[i]));
    8794
    8895}
     
    9097/*FUNCTION Sing::destructor {{{1*/
    9198Sing::~Sing(){
     99        delete inputs;
    92100        return;
    93101}
     
    136144        hmatpar.DeepEcho();
    137145        hnumpar.DeepEcho();
    138         properties.DeepEcho();
     146        printf("   inputs\n");
     147        inputs->DeepEcho();
    139148
    140149        return;
     
    160169        hnumpar.Demarshall(&marshalled_dataset);
    161170
    162         /*demarshall properties: */
    163         properties.Demarshall(&marshalled_dataset);
     171        /*demarshall inputs: */
     172        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);
    164173
    165174        /*return: */
     
    178187        hmatpar.Echo();
    179188        hnumpar.Echo();
    180         properties.Echo();
    181 
    182         return;
     189        printf("   inputs\n");
     190        inputs->Echo();
    183191}
    184192/*}}}*/
     
    188196        char* marshalled_dataset=NULL;
    189197        int   enum_type=0;
     198        char* marshalled_inputs=NULL;
     199        int   marshalled_inputs_size;
    190200
    191201        /*recover marshalled_dataset: */
     
    207217        hnumpar.Marshall(&marshalled_dataset);
    208218
    209         /*Marshall properties: */
    210         properties.Marshall(&marshalled_dataset);
     219        /*Marshall inputs: */
     220        marshalled_inputs_size=inputs->MarshallSize();
     221        marshalled_inputs=inputs->Marshall();
     222        memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
     223        marshalled_dataset+=marshalled_inputs_size;
     224       
     225        xfree((void**)&marshalled_inputs);
    211226
    212227        *pmarshalled_dataset=marshalled_dataset;
     
    222237                +hmatpar.MarshallSize()
    223238                +hnumpar.MarshallSize()
    224                 +properties.MarshallSize()
     239                +inputs->MarshallSize()
    225240                +sizeof(int); //sizeof(int) for enum type
    226241}
     
    243258void  Sing::ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type){
    244259
    245         int i;
    246         const int numgrids=1;
    247         int doflist[numgrids];
    248         double pressure[numgrids];
     260        int    dof;
     261        double pressure;
     262        double thickness;
    249263        double rho_ice,g;
    250264
    251265        /*dynamic objects pointed to by hooks: */
    252         Node**  nodes=NULL;
    253266        Matpar* matpar=NULL;
    254267        Matice* matice=NULL;
    255         Numpar* numpar=NULL;
    256268
    257269        /*Get dof list on which we will plug the pressure values: */
    258         GetDofList1(&doflist[0]);
     270        GetDofList1(&dof);
    259271
    260272        /*recover objects from hooks: */
    261         nodes=(Node**)hnodes.deliverp();
    262273        matpar=(Matpar*)hmatpar.delivers();
    263274        matice=(Matice*)hmatice.delivers();
    264         numpar=(Numpar*)hnumpar.delivers();
    265275
    266276        /*pressure is lithostatic: */
    267277        rho_ice=matpar->GetRhoIce();
    268278        g=matpar->GetG();
    269         pressure[0]=rho_ice*g*this->properties.h[0];
     279        inputs->GetParameterValue(&thickness,ThicknessEnum);
     280        pressure=rho_ice*g*thickness;
    270281       
    271282        /*plug local pressure values into global pressure vector: */
    272         VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
     283        VecSetValue(p_g,dof,pressure,INSERT_VALUES);
    273284
    274285}
     
    282293/*}}}*/
    283294/*FUNCTION Sing::CostFunction {{{1*/
    284 double Sing::CostFunction(void*, int,int){
     295double Sing::CostFunction( int,int){
    285296        ISSMERROR(" not supported yet!");
    286297}
     
    293304        if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){
    294305
    295                 CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
     306                CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type);
    296307
    297308        }
     
    313324        int    numberofdofspernode;
    314325
    315         ParameterInputs* inputs=NULL;
    316 
    317         /*recover pointers: */
    318         inputs=(ParameterInputs*)vinputs;
    319        
    320326        GetDofList(&doflist[0],&numberofdofspernode);
    321327
     
    330336        if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){
    331337       
    332                         CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
     338                        CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type);
    333339
    334340        }
     
    365371        Numpar* numpar=NULL;
    366372
    367         ParameterInputs* inputs=NULL;
    368 
    369         /*recover pointers: */
    370         inputs=(ParameterInputs*)vinputs;
    371 
    372373        /*recover objects from hooks: */
    373374        node=(Node**)hnodes.deliverp();
     
    376377        numpar=(Numpar*)hnumpar.delivers();
    377378
    378         found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,numgrids,(void**)&node[0]);
    379         if(!found)ISSMERROR(" surfaceslopex missing from inputs!");
    380         found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,numgrids,(void**)&node[0]);
    381         if(!found)ISSMERROR(" surfaceslopey missing from inputs!");
     379        inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum);
     380        inputs->GetParameterValue(&slope[1],SurfaceSlopeyEnum);
    382381
    383382        GetDofList(&doflist[0],&numberofdofspernode);
     
    391390        n=matice->GetN();
    392391        B=matice->GetB();
    393         thickness=this->properties.h[0];
     392        inputs->GetParameterValue(&thickness,ThicknessEnum);
    394393
    395394        ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
     
    408407/*}}}*/
    409408/*FUNCTION Sing::Du {{{1*/
    410 void  Sing::Du(_p_Vec*,void*,int,int){
     409void  Sing::Du(Vec,int,int){
    411410        ISSMERROR(" not supported yet!");
    412411}
     
    504503/*}}}*/
    505504/*FUNCTION Sing::GetShelf {{{1*/
    506 int   Sing::GetShelf(){
     505bool   Sing::GetShelf(){
    507506        ISSMERROR(" not supported yet!");
    508507}
     
    514513/*}}}*/
    515514/*FUNCTION Sing::Gradj {{{1*/
    516 void  Sing::Gradj(_p_Vec*, void*, int, int ,char*){
     515void  Sing::Gradj(Vec, int, int ,char*){
    517516        ISSMERROR(" not supported yet!");
    518517}
    519518/*}}}*/
    520519/*FUNCTION Sing::GradB {{{1*/
    521 void  Sing::GradjB(_p_Vec*, void*, int,int){
     520void  Sing::GradjB(Vec, int,int){
    522521        ISSMERROR(" not supported yet!");
    523522}
    524523/*}}}*/
    525524/*FUNCTION Sing::GradjDrag {{{1*/
    526 void  Sing::GradjDrag(_p_Vec*, void*, int,int){
     525void  Sing::GradjDrag(Vec, int,int){
    527526        ISSMERROR(" not supported yet!");
    528527}
     
    534533/*}}}*/
    535534/*FUNCTION Sing::Misfit {{{1*/
    536 double Sing::Misfit(void*, int,int){
     535double Sing::Misfit( int,int){
    537536        ISSMERROR(" not supported yet!");
    538537}
     
    545544/*}}}*/
    546545/*FUNCTION Sing::SurfaceArea {{{1*/
    547 double Sing::SurfaceArea(void*, int,int){
     546double Sing::SurfaceArea( int,int){
    548547        ISSMERROR(" not supported yet!");
    549548}
  • issm/trunk/src/c/objects/Sing.h

    r3612 r3620  
    1212#include "../ModelProcessorx/IoModel.h"
    1313#include "./Hook.h"
    14 
    15 class Hook;
     14#include "../DataSet/Inputs.h"
    1615
    1716class Sing: public Element{
     
    2625                Hook hmatpar; //hook to 1 matpar
    2726                Hook hnumpar; //hook to 1 numpar
    28 
    2927                Inputs* inputs;
    3028
     
    3937                /*}}}*/
    4038                /*object management: {{{1*/
    41                 void  Configure(void* loads,void* nodes,void* materials,void* parameters);
     39                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
    4240                Object* copy();
    4341                void  DeepEcho();
     
    7270                void  GetBedList(double*);
    7371                void  GetThicknessList(double* thickness_list);
    74                 void  Du(_p_Vec*,void*,int,int);
    75                 void  Gradj(_p_Vec*, void*, int, int,char*);
    76                 void  GradjDrag(_p_Vec*, void*, int,int);
    77                 void  GradjB(_p_Vec*, void*, int,int);
    78                 double Misfit(void*,int,int);
    79                 double SurfaceArea(void*,int,int);
    80                 double CostFunction(void*,int,int);
     72                void  Du(Vec,int,int);
     73                void  Gradj(Vec, int, int,char*);
     74                void  GradjDrag(Vec , int,int);
     75                void  GradjB(Vec, int,int);
     76                double Misfit(int,int);
     77                double SurfaceArea(int,int);
     78                double CostFunction(int,int);
    8179                double MassFlux(double* segment,double* ug);
    8280                /*}}}*/
  • issm/trunk/src/c/objects/Tria.cpp

    r3612 r3620  
    3030}
    3131/*}}}*/
    32 /*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id){{{1*/
    33 Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id):
     32/*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id){{{1*/
     33Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id):
    3434        hnodes(tria_node_ids,3),
    3535        hmatice(&tria_matice_id,1),
    36         hmatpar(&tria_matpar_id,1),
    37         hnumpar(&tria_numpar_id,1)
     36        hmatpar(&tria_matpar_id,1)
    3837{
    3938
    4039        /*all the initialization has been done by the initializer, just fill in the id: */
    4140        this->id=tria_id;
     41        this->parameters=NULL;
    4242        this->inputs=new Inputs();
    4343
    44         return;
    45 }
    46 /*}}}*/
    47 /*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, Inputs* tria_inputs) {{{1*/
    48 Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs):
     44}
     45/*}}}*/
     46/*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* tria_inputs) {{{1*/
     47Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* tria_parameters, Inputs* tria_inputs):
    4948        hnodes(tria_hnodes),
    5049        hmatice(tria_hmatice),
    51         hmatpar(tria_hmatpar),
    52         hnumpar(tria_hnumpar)
     50        hmatpar(tria_hmatpar)
    5351{
    5452
     
    6159                this->inputs=new Inputs();
    6260        }
    63 
    64         return;
     61        /*point parameters: */
     62        this->parameters=tria_parameters;
    6563}
    6664/*}}}*/
     
    7270        int tria_matice_id;
    7371        int tria_matpar_id;
    74         int tria_numpar_id;
    7572        double nodeinputs[3];
    7673
     
    9491        tria_matice_id=index+1; //refers to the corresponding ice material object
    9592        tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
    96         tria_numpar_id=1; //refers to numerical parameters object
    9793
    9894        this->hnodes.Init(tria_node_ids,3);
    9995        this->hmatice.Init(&tria_matice_id,1);
    10096        this->hmatpar.Init(&tria_matpar_id,1);
    101         this->hnumpar.Init(&tria_numpar_id,1);
    10297
    10398        //intialize inputs, and add as many inputs per element as requested:
     
    143138        if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
    144139
     140        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
     141        this->parameters=NULL;
     142
    145143
    146144}
     
    149147Tria::~Tria(){
    150148        delete inputs;
     149        this->parameters=NULL;
    151150        return;
    152151}
     
    155154/*Object management: */
    156155/*FUNCTION Tria::Configure {{{1*/
    157 void  Tria::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){
     156void  Tria::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){
    158157
    159158        int i;
    160 
    161         DataSet* loadsin=NULL;
    162         DataSet* nodesin=NULL;
    163         DataSet* materialsin=NULL;
    164         DataSet* parametersin=NULL;
    165 
    166         /*Recover pointers :*/
    167         loadsin=(DataSet*)ploadsin;
    168         nodesin=(DataSet*)pnodesin;
    169         materialsin=(DataSet*)pmaterialsin;
    170         parametersin=(DataSet*)pparametersin;
    171159
    172160        /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
     
    175163        hmatice.configure(materialsin);
    176164        hmatpar.configure(materialsin);
    177         hnumpar.configure(parametersin);
     165
     166        /*point parameters to real dataset: */
     167        this->parameters=parametersin;
    178168
    179169}
     
    182172Object* Tria::copy() {
    183173
    184         return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs);
     174        return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs);
    185175
    186176}
     
    205195        hmatice.Demarshall(&marshalled_dataset);
    206196        hmatpar.Demarshall(&marshalled_dataset);
    207         hnumpar.Demarshall(&marshalled_dataset);
    208197       
    209198        /*demarshall inputs: */
    210199        inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);
     200
     201        /*parameters: may not exist even yet, so let Configure handle it: */
     202        this->parameters=NULL;
    211203
    212204        /*return: */
     
    224216        hmatice.DeepEcho();
    225217        hmatpar.DeepEcho();
    226         hnumpar.DeepEcho();
     218        printf("   parameters\n");
     219        parameters->DeepEcho();
    227220        printf("   inputs\n");
    228221        inputs->DeepEcho();
    229 
     222       
    230223        return;
    231224}
     
    240233        hmatice.Echo();
    241234        hmatpar.Echo();
    242         hnumpar.Echo();
     235        printf("   parameters\n");
     236        parameters->Echo();
    243237        printf("   inputs\n");
    244238        inputs->Echo();
    245 
    246         return;
    247239}
    248240/*}}}*/
     
    271263        hmatice.Marshall(&marshalled_dataset);
    272264        hmatpar.Marshall(&marshalled_dataset);
    273         hnumpar.Marshall(&marshalled_dataset);
    274265
    275266        /*Marshall inputs: */
     
    279270        marshalled_dataset+=marshalled_inputs_size;
    280271
     272        /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/
     273
     274        xfree((void**)&marshalled_inputs);
     275
    281276        *pmarshalled_dataset=marshalled_dataset;
    282277        return;
     
    290285                +hmatice.MarshallSize()
    291286                +hmatpar.MarshallSize()
    292                 +hnumpar.MarshallSize()
    293287                +inputs->MarshallSize()
    294288                +sizeof(int); //sizeof(int) for enum type
     
    312306        Matpar* matpar=NULL;
    313307        Matice* matice=NULL;
    314         Numpar* numpar=NULL;
    315308
    316309        /*recover objects from hooks: */
     
    318311        matpar=(Matpar*)hmatpar.delivers();
    319312        matice=(Matice*)hmatice.delivers();
    320         numpar=(Numpar*)hnumpar.delivers();
    321313
    322314        /*Update internal data if inputs holds new values: */
     
    468460        Matpar* matpar=NULL;
    469461        Matice* matice=NULL;
    470         Numpar* numpar=NULL;
    471462
    472463        /*recover objects from hooks: */
     
    474465        matpar=(Matpar*)hmatpar.delivers();
    475466        matice=(Matice*)hmatice.delivers();
    476         numpar=(Numpar*)hnumpar.delivers();
    477467
    478468        /*plug local pressure values into global pressure vector: */
     
    497487        Matpar* matpar=NULL;
    498488        Matice* matice=NULL;
    499         Numpar* numpar=NULL;
    500489
    501490        /*recover objects from hooks: */
     
    503492        matpar=(Matpar*)hmatpar.delivers();
    504493        matice=(Matice*)hmatice.delivers();
    505         numpar=(Numpar*)hnumpar.delivers();
    506494
    507495        /*Get dof list on which we will plug the pressure values: */
     
    536524        Matpar* matpar=NULL;
    537525        Matice* matice=NULL;
    538         Numpar* numpar=NULL;
    539526
    540527        /*recover objects from hooks: */
     
    542529        matpar=(Matpar*)hmatpar.delivers();
    543530        matice=(Matice*)hmatice.delivers();
    544         numpar=(Numpar*)hnumpar.delivers();
    545531
    546532        /*plug local pressure values into global pressure vector: */
     
    583569        double  dk[NDOF2];
    584570        double  dB[NDOF2];
     571        char*   control_type=NULL;
     572        double  cm_noisedmp;
    585573
    586574        /* Jacobian: */
     
    591579        Matpar* matpar=NULL;
    592580        Matice* matice=NULL;
    593         Numpar* numpar=NULL;
    594581
    595582        /*inputs: */
     
    608595        matpar=(Matpar*)hmatpar.delivers();
    609596        matice=(Matice*)hmatice.delivers();
    610         numpar=(Numpar*)hnumpar.delivers();
     597
     598        /*retrieve some parameters: */
     599        this->parameters->FindParam(&control_type,"control_type");
     600        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
    611601
    612602        /* Get node coordinates and dof list: */
     
    631621
    632622                /*Add Tikhonov regularization term to misfit*/
    633                 if (strcmp(numpar->control_type,"drag")==0){
     623                if (strcmp(control_type,"drag")==0){
    634624                        if (shelf){
    635625                                inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
    636                                 Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
     626                                Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
    637627
    638628                        }
    639629                }
    640                 else if (strcmp(numpar->control_type,"B")==0){
     630                else if (strcmp(control_type,"B")==0){
    641631                        inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBEnum);
    642                         Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
     632                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
    643633                }
    644634                else{
    645                         ISSMERROR("%s%s","unsupported control type: ",numpar->control_type);
     635                        ISSMERROR("%s%s","unsupported control type: ",control_type);
    646636                }
    647637
     
    652642        xfree((void**)&third_gauss_area_coord);
    653643        xfree((void**)&gauss_weights);
     644        xfree((void**)&control_type);
    654645
    655646        /*Return: */
     
    763754        Matpar* matpar=NULL;
    764755        Matice* matice=NULL;
    765         Numpar* numpar=NULL;
     756
     757        /*parameters: */
     758        double artdiff;
    766759
    767760
     
    770763        matpar=(Matpar*)hmatpar.delivers();
    771764        matice=(Matice*)hmatice.delivers();
    772         numpar=(Numpar*)hnumpar.delivers();
    773765
    774766        /* Get node coordinates and dof list: */
     
    780772        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
    781773
     774        /*retrieve some parameters: */
     775        this->parameters->FindParam(&artdiff,"artdiff");
     776
    782777        //Create Artificial diffusivity once for all if requested
    783         if(numpar->artdiff){
     778        if(artdiff){
    784779                //Get the Jacobian determinant
    785780                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
     
    848843                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
    849844
    850                 if(numpar->artdiff){
     845                if(artdiff){
    851846
    852847                        /* Compute artificial diffusivity */
     
    918913        Matpar* matpar=NULL;
    919914        Matice* matice=NULL;
    920         Numpar* numpar=NULL;
    921915
    922916
     
    925919        matpar=(Matpar*)hmatpar.delivers();
    926920        matice=(Matice*)hmatice.delivers();
    927         numpar=(Numpar*)hnumpar.delivers();
    928921
    929922        /* Get node coordinates and dof list: */
     
    10341027        int     found=0;
    10351028
     1029        /*parameters: */
     1030        double artdiff;
     1031
    10361032        /*dynamic objects pointed to by hooks: */
    10371033        Node**  nodes=NULL;
    10381034        Matpar* matpar=NULL;
    10391035        Matice* matice=NULL;
    1040         Numpar* numpar=NULL;
    10411036
    10421037        /*recover objects from hooks: */
     
    10441039        matpar=(Matpar*)hmatpar.delivers();
    10451040        matice=(Matice*)hmatice.delivers();
    1046         numpar=(Numpar*)hnumpar.delivers();
    10471041
    10481042        /*Recover velocity: */
    10491043        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
    10501044        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
     1045
     1046        /*retrieve some parameters: */
     1047        this->parameters->FindParam(&artdiff,"artdiff");
    10511048
    10521049        /* Get node coordinates and dof list: */
     
    10751072
    10761073        //Create Artificial diffusivity once for all if requested
    1077         if(numpar->artdiff){
     1074        if(artdiff){
    10781075                //Get the Jacobian determinant
    10791076                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
     
    11311128                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
    11321129
    1133                 if(numpar->artdiff){
     1130                if(artdiff){
    11341131
    11351132                        /* Compute artificial diffusivity */
     
    11951192        double B[3][numdof];
    11961193        double Bprime[3][numdof];
    1197         double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
     1194        double D[3][3]={0.0};  // material matrix, simple scalar matrix.
    11981195        double D_scalar;
    11991196
     1197        /*parameters: */
     1198        double viscosity_overshoot;
     1199
    12001200        /* local element matrices: */
    1201         double Ke_gg[numdof][numdof]; //local element stiffness matrix
     1201        double Ke_gg[numdof][numdof]={0.0};
    12021202        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    12031203
     
    12121212        Matpar* matpar=NULL;
    12131213        Matice* matice=NULL;
    1214         Numpar* numpar=NULL;
    12151214
    12161215        /*inputs: */
     
    12201219        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    12211220        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     1221
     1222        /*retrieve some parameters: */
     1223        this->parameters->FindParam(&viscosity_overshoot,"viscosity_overshoot");
    12221224
    12231225        /*First, if we are on water, return empty matrix: */
     
    12281230        matpar=(Matpar*)hmatpar.delivers();
    12291231        matice=(Matice*)hmatice.delivers();
    1230         numpar=(Numpar*)hnumpar.delivers();
    12311232
    12321233        /* Get node coordinates and dof list: */
    12331234        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    12341235        GetDofList(&doflist[0],&numberofdofspernode);
    1235 
    1236         /* Set Ke_gg to 0: */
    1237         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    12381236
    12391237        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     
    12651263                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    12661264                        onto this scalar matrix, so that we win some computational time: */
    1267                 newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity);
     1265                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    12681266                D_scalar=newviscosity*thickness*gauss_weight*Jdet;
    12691267
     
    13331331
    13341332        /* local element matrices: */
    1335         double Ke_gg[numdof][numdof]; //local element stiffness matrix
     1333        double Ke_gg[numdof][numdof]={0.0};
    13361334        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    13371335       
     
    13651363        Matpar* matpar=NULL;
    13661364        Matice* matice=NULL;
    1367         Numpar* numpar=NULL;
    13681365
    13691366
     
    13721369        matpar=(Matpar*)hmatpar.delivers();
    13731370        matice=(Matice*)hmatice.delivers();
    1374         numpar=(Numpar*)hnumpar.delivers();
    13751371
    13761372        /*retrieve inputs :*/
     
    13811377        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    13821378        GetDofList(&doflist[0],&numberofdofspernode);
    1383 
    1384         /* Set Ke_gg to 0: */
    1385         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    13861379
    13871380        if (shelf){
     
    15221515
    15231516        /* local element matrices: */
    1524         double Ke_gg[numdof][numdof]; //local element stiffness matrix
     1517        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
    15251518        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    15261519
     
    15291522        Matpar* matpar=NULL;
    15301523        Matice* matice=NULL;
    1531         Numpar* numpar=NULL;
    15321524
    15331525        /*recover objects from hooks: */
     
    15351527        matpar=(Matpar*)hmatpar.delivers();
    15361528        matice=(Matice*)hmatice.delivers();
    1537         numpar=(Numpar*)hnumpar.delivers();
    15381529
    15391530        /* Get node coordinates and dof list: */
    15401531        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    15411532        GetDofList(&doflist[0],&numberofdofspernode);
    1542 
    1543         /* Set Ke_gg to 0: */
    1544         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    15451533
    15461534        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     
    16551643        Matpar* matpar=NULL;
    16561644        Matice* matice=NULL;
    1657         Numpar* numpar=NULL;
    16581645
    16591646        /*recover objects from hooks: */
     
    16611648        matpar=(Matpar*)hmatpar.delivers();
    16621649        matice=(Matice*)hmatice.delivers();
    1663         numpar=(Numpar*)hnumpar.delivers();
    16641650
    16651651        /*Recover constants of ice */
     
    17591745        double  K[2][2]={0.0};
    17601746        double  KDL[2][2]={0.0};
    1761         double  dt;
    17621747        int     dofs[2]={0,1};
    17631748        int     found;
     1749
     1750        /*parameters: */
     1751        double dt;
     1752        double artdiff;
    17641753
    17651754        /*dynamic objects pointed to by hooks: */
     
    17671756        Matpar* matpar=NULL;
    17681757        Matice* matice=NULL;
    1769         Numpar* numpar=NULL;
    17701758
    17711759        /*recover objects from hooks: */
     
    17731761        matpar=(Matpar*)hmatpar.delivers();
    17741762        matice=(Matice*)hmatice.delivers();
    1775         numpar=(Numpar*)hnumpar.delivers();
    1776 
    1777         /*recover dt: */
    1778         dt=numpar->dt;
     1763
     1764        /*retrieve some parameters: */
     1765        this->parameters->FindParam(&dt,"dt");
     1766        this->parameters->FindParam(&artdiff,"artdiff");
    17791767
    17801768        /* Get node coordinates and dof list: */
     
    17871775
    17881776        //Create Artificial diffusivity once for all if requested
    1789         if(numpar->artdiff==1){
     1777        if(artdiff==1){
    17901778                //Get the Jacobian determinant
    17911779                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
     
    18671855                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
    18681856
    1869                 if(numpar->artdiff==1){
     1857                if(artdiff==1){
    18701858
    18711859                        /* Compute artificial diffusivity */
     
    19341922        /*input parameters for structural analysis (diagnostic): */
    19351923        double  vx,vy;
    1936         double  dt;
    19371924        int     dofs[1]={0};
    19381925        int     found;
     1926
     1927        /*parameters: */
     1928        double dt;
    19391929
    19401930        /*dynamic objects pointed to by hooks: */
     
    19421932        Matpar* matpar=NULL;
    19431933        Matice* matice=NULL;
    1944         Numpar* numpar=NULL;
    19451934
    19461935        /*recover objects from hooks: */
     
    19481937        matpar=(Matpar*)hmatpar.delivers();
    19491938        matice=(Matice*)hmatice.delivers();
    1950         numpar=(Numpar*)hnumpar.delivers();
    1951 
    1952         /*recover dt: */
    1953         dt=numpar->dt;
     1939
     1940        /*retrieve some parameters: */
     1941        this->parameters->FindParam(&dt,"dt");
    19541942
    19551943        /* Get node coordinates and dof list: */
     
    20502038
    20512039        /* local element matrices: */
    2052         double Ke_gg[numdof][numdof]; //local element stiffness matrix
     2040        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
    20532041        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    20542042       
     
    20592047        Matpar* matpar=NULL;
    20602048        Matice* matice=NULL;
    2061         Numpar* numpar=NULL;
    20622049
    20632050        /*recover objects from hooks: */
     
    20652052        matpar=(Matpar*)hmatpar.delivers();
    20662053        matice=(Matice*)hmatice.delivers();
    2067         numpar=(Numpar*)hnumpar.delivers();
    20682054
    20692055        /* Get node coordinates and dof list: */
    20702056        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    20712057        GetDofList(&doflist[0],&numberofdofspernode);
    2072 
    2073         /* Set Ke_gg to 0: */
    2074         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    20752058
    20762059        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     
    21332116        double rho_ice;
    21342117        double heatcapacity;
    2135         double dt;
    21362118
    21372119        int     num_gauss,ig;
     
    21512133        double  D_scalar;
    21522134
     2135        /*parameters: */
     2136        double dt;
     2137
    21532138        /*dynamic objects pointed to by hooks: */
    21542139        Node**  nodes=NULL;
    21552140        Matpar* matpar=NULL;
    21562141        Matice* matice=NULL;
    2157         Numpar* numpar=NULL;
    21582142
    21592143        /*recover objects from hooks: */
     
    21612145        matpar=(Matpar*)hmatpar.delivers();
    21622146        matice=(Matice*)hmatice.delivers();
    2163         numpar=(Numpar*)hnumpar.delivers();
    2164 
    2165         /*recover extra inputs from users, dt: */
    2166         dt=numpar->dt;
     2147
     2148        /*retrieve some parameters: */
     2149        this->parameters->FindParam(&dt,"dt");
    21672150
    21682151        /* Get node coordinates and dof list: */
     
    23042287        Matpar* matpar=NULL;
    23052288        Matice* matice=NULL;
    2306         Numpar* numpar=NULL;
    23072289
    23082290        /*recover objects from hooks: */
     
    23102292        matpar=(Matpar*)hmatpar.delivers();
    23112293        matice=(Matice*)hmatice.delivers();
    2312         numpar=(Numpar*)hnumpar.delivers();
    23132294
    23142295        /* Get node coordinates and dof list: */
     
    23912372        Matpar* matpar=NULL;
    23922373        Matice* matice=NULL;
    2393         Numpar* numpar=NULL;
    23942374
    23952375        /*recover objects from hooks: */
     
    23972377        matpar=(Matpar*)hmatpar.delivers();
    23982378        matice=(Matice*)hmatice.delivers();
    2399         numpar=(Numpar*)hnumpar.delivers();
    24002379
    24012380        /* Get node coordinates and dof list: */
     
    24782457        Matpar* matpar=NULL;
    24792458        Matice* matice=NULL;
    2480         Numpar* numpar=NULL;
    24812459
    24822460        /*recover objects from hooks: */
     
    24842462        matpar=(Matpar*)hmatpar.delivers();
    24852463        matice=(Matice*)hmatice.delivers();
    2486         numpar=(Numpar*)hnumpar.delivers();
    2487 
    2488         /*recover objects from hooks: */
    2489         nodes=(Node**)hnodes.deliverp();
    2490         matpar=(Matpar*)hmatpar.delivers();
    2491         matice=(Matice*)hmatice.delivers();
    2492         numpar=(Numpar*)hnumpar.delivers();
    24932464
    24942465        /* Get node coordinates and dof list: */
     
    25622533
    25632534        /*element vector at the gaussian points: */
    2564         double  pe_g[numdof];
     2535        double  pe_g[numdof]={0.0};
    25652536        double  pe_g_gaussian[numdof];
    25662537
     
    25782549        Matpar* matpar=NULL;
    25792550        Matice* matice=NULL;
    2580         Numpar* numpar=NULL;
    25812551
    25822552        /*recover objects from hooks: */
     
    25842554        matpar=(Matpar*)hmatpar.delivers();
    25852555        matice=(Matice*)hmatice.delivers();
    2586         numpar=(Numpar*)hnumpar.delivers();
    25872556
    25882557        /* Get node coordinates and dof list: */
    25892558        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    25902559        GetDofList(&doflist[0],&numberofdofspernode);
    2591 
    2592         /* Set pe_g to 0: */
    2593         for(i=0;i<numdof;i++) pe_g[i]=0.0;
    25942560
    25952561        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     
    26802646
    26812647        /*element vector at the gaussian points: */
    2682         double  pe_g[numdof];
     2648        double  pe_g[numdof]={0.0};
    26832649        double  pe_g_gaussian[numdof];
    26842650
     
    26902656        Matpar* matpar=NULL;
    26912657        Matice* matice=NULL;
    2692         Numpar* numpar=NULL;
    26932658
    26942659        /*inputs: */
     
    27072672        matpar=(Matpar*)hmatpar.delivers();
    27082673        matice=(Matice*)hmatice.delivers();
    2709         numpar=(Numpar*)hnumpar.delivers();
    27102674
    27112675        /* Get node coordinates and dof list: */
    27122676        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    27132677        GetDofList(&doflist[0],&numberofdofspernode);
    2714 
    2715         /* Set pe_g to 0: */
    2716         for(i=0;i<numdof;i++) pe_g[i]=0.0;
    27172678
    27182679
     
    28122773        double  melting_g;
    28132774        double  thickness_g;
     2775
     2776        /*parameters: */
    28142777        double  dt;
    28152778
     
    28182781        Matpar* matpar=NULL;
    28192782        Matice* matice=NULL;
    2820         Numpar* numpar=NULL;
    28212783
    28222784        /*recover objects from hooks: */
     
    28242786        matpar=(Matpar*)hmatpar.delivers();
    28252787        matice=(Matice*)hmatice.delivers();
    2826         numpar=(Numpar*)hnumpar.delivers();
    2827 
    2828         /*recover dt: */
    2829         dt=numpar->dt;
     2788
     2789        /*retrieve some parameters: */
     2790        this->parameters->FindParam(&dt,"dt");
    28302791
    28312792        /* Get node coordinates and dof list: */
     
    29092870        Matpar* matpar=NULL;
    29102871        Matice* matice=NULL;
    2911         Numpar* numpar=NULL;
    29122872
    29132873        /*recover objects from hooks: */
     
    29152875        matpar=(Matpar*)hmatpar.delivers();
    29162876        matice=(Matice*)hmatice.delivers();
    2917         numpar=(Numpar*)hnumpar.delivers();
    2918 
    2919         /*recover dt: */
    2920         dt=numpar->dt;
    2921 
     2877
     2878        /*retrieve some parameters: */
     2879        this->parameters->FindParam(&dt,"dt");
    29222880
    29232881        /* Get node coordinates and dof list: */
     
    29932951
    29942952        /*element vector at the gaussian points: */
    2995         double  pe_g[numdof];
     2953        double  pe_g[numdof]={0.0};
    29962954        double  pe_g_gaussian[numdof];
    29972955        double  slope[2];
     
    30012959        Matpar* matpar=NULL;
    30022960        Matice* matice=NULL;
    3003         Numpar* numpar=NULL;
    30042961
    30052962        /*recover objects from hooks: */
     
    30072964        matpar=(Matpar*)hmatpar.delivers();
    30082965        matice=(Matice*)hmatice.delivers();
    3009         numpar=(Numpar*)hnumpar.delivers();
    30102966
    30112967        /* Get node coordinates and dof list: */
    30122968        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    30132969        GetDofList(&doflist[0],&numberofdofspernode);
    3014 
    3015         /* Set pe_g to 0: */
    3016         for(i=0;i<numdof;i++) pe_g[i]=0.0;
    30172970
    30182971
     
    31123065        Matpar* matpar=NULL;
    31133066        Matice* matice=NULL;
    3114         Numpar* numpar=NULL;
    31153067
    31163068        /*recover objects from hooks: */
     
    31183070        matpar=(Matpar*)hmatpar.delivers();
    31193071        matice=(Matice*)hmatice.delivers();
    3120         numpar=(Numpar*)hnumpar.delivers();
    31213072       
    31223073        /* Get node coordinates and dof list: */
     
    31323083        beta=matpar->GetBeta();
    31333084        meltingpoint=matpar->GetMeltingPoint();
    3134         dt=numpar->dt;
    3135 
     3085       
     3086        /*retrieve some solution parameters: */
     3087        this->parameters->FindParam(&dt,"dt");
    31363088
    31373089        /* Ice/ocean heat exchange flux on ice shelf base */
     
    32313183        Matpar* matpar=NULL;
    32323184        Matice* matice=NULL;
    3233         Numpar* numpar=NULL;
    32343185
    32353186        /*recover objects from hooks: */
     
    32373188        matpar=(Matpar*)hmatpar.delivers();
    32383189        matice=(Matice*)hmatice.delivers();
    3239         numpar=(Numpar*)hnumpar.delivers();
    32403190
    32413191        /*retrieve inputs :*/
     
    32493199        rho_ice=matpar->GetRhoIce();
    32503200        heatcapacity=matpar->GetHeatCapacity();
    3251         dt=numpar->dt;
     3201
     3202        /*retrieve some parameters: */
     3203        this->parameters->FindParam(&dt,"dt");
    32523204
    32533205
     
    33693321        double  obs_velocity_mag,velocity_mag;
    33703322        double  dux,duy;
     3323        double  meanvel, epsvel;
    33713324
    33723325        /*element vector : */
    3373         double  due_g[numdof]={0,0,0,0,0,0};
     3326        double  due_g[numdof]={0.0};
    33743327        double  due_g_gaussian[numdof];
    33753328
     
    33913344        Matpar* matpar=NULL;
    33923345        Matice* matice=NULL;
    3393         Numpar* numpar=NULL;
    33943346
    33953347        /*recover objects from hooks: */
     
    33973349        matpar=(Matpar*)hmatpar.delivers();
    33983350        matice=(Matice*)hmatice.delivers();
    3399         numpar=(Numpar*)hnumpar.delivers();
     3351
     3352        /*retrieve some parameters: */
     3353        this->parameters->FindParam(&meanvel,"meanvel");
     3354        this->parameters->FindParam(&epsvel,"epsvel");
    34003355
    34013356        /* Get node coordinates and dof list: */
     
    34603415                 */
    34613416                for (i=0;i<numgrids;i++){
    3462                         scalex=pow(numpar->meanvel/(obs_vx_list[i]+numpar->epsvel),2);
    3463                         scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),2);
     3417                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
     3418                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
    34643419                        if(obs_vx_list[i]==0)scalex=0;
    34653420                        if(obs_vy_list[i]==0)scaley=0;
     
    34823437                 */
    34833438                for (i=0;i<numgrids;i++){
    3484                         velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+numpar->epsvel; //epsvel to avoid velocity being nil.
    3485                         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.
    3486                         scale=-8*pow(numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
     3439                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
     3440                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
     3441                        scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    34873442                        dux_list[i]=scale*vx_list[i];
    34883443                        duy_list[i]=scale*vy_list[i];
     
    35013456                 */
    35023457                for (i=0;i<numgrids;i++){
    3503                         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);
     3458                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
    35043459                        dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
    35053460                        duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
     
    35183473                 */
    35193474                for (i=0;i<numgrids;i++){
    3520                         dux_list[i] = - pow(numpar->meanvel,(double)2)*(
    3521                                                 log((fabs(vx_list[i])+numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)) * 1/(vx_list[i]+numpar->epsvel));
    3522                         duy_list[i] = - pow(numpar->meanvel,(double)2)*(
    3523                                                 log((fabs(vy_list[i])+numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)) * 1/(vy_list[i]+numpar->epsvel));
     3475                        dux_list[i] = - pow(meanvel,(double)2)*(
     3476                                                log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
     3477                        duy_list[i] = - pow(meanvel,(double)2)*(
     3478                                                log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
    35243479                }
    35253480        }
     
    42284183
    42294184        /*element vector at the gaussian points: */
    4230         double  grade_g[numgrids];
     4185        double  grade_g[numgrids]={0.0};
    42314186        double  grade_g_gaussian[numgrids];
    42324187
     
    42524207        double  dB[NDOF2];
    42534208        double  B_gauss;
     4209       
     4210        /*parameters: */
     4211        double  cm_noisedmp;
     4212        double  cm_mindmp_slope;
     4213        double  cm_mindmp_value;
     4214        double  cm_maxdmp_value;
     4215        double  cm_maxdmp_slope;
    42544216
    42554217        /*dynamic objects pointed to by hooks: */
     
    42574219        Matpar* matpar=NULL;
    42584220        Matice* matice=NULL;
    4259         Numpar* numpar=NULL;
    42604221
    42614222        /*recover objects from hooks: */
     
    42634224        matpar=(Matpar*)hmatpar.delivers();
    42644225        matice=(Matice*)hmatice.delivers();
    4265         numpar=(Numpar*)hnumpar.delivers();
     4226
     4227        /*retrieve some parameters: */
     4228        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
     4229        this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value");
     4230        this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope");
     4231        this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value");
     4232        this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope");
    42664233
    42674234        /* Get node coordinates and dof list: */
    42684235        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    42694236        GetDofList1(&doflist1[0]);
    4270 
    4271         /* Set grade_g to 0: */
    4272         for(i=0;i<numgrids;i++) grade_g[i]=0.0;
    42734237
    42744238        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     
    43174281
    43184282                        //Add regularization term
    4319                         grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);
     4283                        grade_g_gaussian[i]-=cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);
    43204284
    43214285                        //min dampening
    4322                         if(B_gauss<numpar->cm_mindmp_value){
    4323                                 grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4286                        if(B_gauss<cm_mindmp_value){
     4287                                grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
    43244288                        }
    43254289
    43264290                        //max dampening
    4327                         if(B_gauss>numpar->cm_maxdmp_value){
    4328                                 grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4291                        if(B_gauss>cm_maxdmp_value){
     4292                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
    43294293                        }
    43304294
     
    43944358
    43954359        /*element vector at the gaussian points: */
    4396         double  grade_g[numgrids]={0,0,0};
     4360        double  grade_g[numgrids]={0.0};
    43974361        double  grade_g_gaussian[numgrids];
    43984362
     
    44084372        /*inputs: */
    44094373        bool shelf;
     4374
     4375        /*parameters: */
     4376        double  cm_noisedmp;
     4377        double  cm_mindmp_slope;
     4378        double  cm_mindmp_value;
     4379        double  cm_maxdmp_value;
     4380        double  cm_maxdmp_slope;
    44104381
    44114382        /*dynamic objects pointed to by hooks: */
     
    44134384        Matpar* matpar=NULL;
    44144385        Matice* matice=NULL;
    4415         Numpar* numpar=NULL;
    44164386
    44174387        /*recover objects from hooks: */
     
    44194389        matpar=(Matpar*)hmatpar.delivers();
    44204390        matice=(Matice*)hmatice.delivers();
    4421         numpar=(Numpar*)hnumpar.delivers();
    44224391
    44234392        /*retrieve inputs :*/
    44244393        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     4394
     4395        /*retrieve some parameters: */
     4396        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
     4397        this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value");
     4398        this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope");
     4399        this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value");
     4400        this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope");
    44254401
    44264402        /*Get out if shelf*/
     
    45204496
    45214497                        //noise dampening d/dki(1/2*(dk/dx)^2)
    4522                         grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     4498                        grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    45234499                       
    45244500                        //min dampening
    4525                         if(drag<numpar->cm_mindmp_value){
    4526                                 grade_g_gaussian[i]+=numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4501                        if(drag<cm_mindmp_value){
     4502                                grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
    45274503                        }
    45284504
    45294505                        //max dampening
    4530                         if(drag>numpar->cm_maxdmp_value){
    4531                                 grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4506                        if(drag>cm_maxdmp_value){
     4507                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
    45324508                        }
    45334509                }
     
    45644540        double vy_list[numgrids];
    45654541        double vz_list[numgrids];
    4566         double vxvyvz_list[numgrids][3];
    4567         double adjx_list[numgrids];
    4568         double adjy_list[numgrids];
    4569         double adjz_list[numgrids];
    4570         double adjxyz_list[numgrids][3];
    4571 
    45724542        double drag;
    4573         int    dofs1[1]={0};
    4574         int    dofs3[3]={0,1,2};
     4543        double  thickness_list[numgrids];
     4544        double  bed_list[numgrids];
     4545        double  dragcoefficient_list[numgrids];
     4546        double  drag_p,drag_q;
     4547        double alpha_complement_list[numgrids];
     4548        double alpha_complement;
    45754549
    45764550        /* gaussian points: */
     
    45824556        double  gauss_weight;
    45834557        double  gauss_l1l2l3[3];
     4558        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    45844559
    45854560        /* parameters: */
     
    45914566        double  dk[NDOF2];
    45924567
    4593         /*drag: */
    4594         double  pcoeff,qcoeff;
    4595         double alpha_complement_list[numgrids];
    4596         double alpha_complement;
    4597 
    45984568        /*element vector at the gaussian points: */
    4599         double  grade_g[numgrids];
     4569        double  grade_g[numgrids]={0.0};
    46004570        double  grade_g_gaussian[numgrids];
    46014571
     
    46124582        bool shelf;
    46134583        int  drag_type;
     4584
     4585        /*parameters: */
     4586        double  cm_noisedmp;
     4587        double  cm_mindmp_slope;
     4588        double  cm_mindmp_value;
     4589        double  cm_maxdmp_value;
     4590        double  cm_maxdmp_slope;
    46144591
    46154592        /*dynamic objects pointed to by hooks: */
     
    46174594        Matpar* matpar=NULL;
    46184595        Matice* matice=NULL;
    4619         Numpar* numpar=NULL;
    46204596
    46214597        /*retrieve inputs :*/
    46224598        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    46234599        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4600
     4601        /*retrieve some parameters: */
     4602        this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp");
     4603        this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value");
     4604        this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope");
     4605        this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value");
     4606        this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope");
    46244607
    46254608        /*Get out if shelf*/
     
    46304613        matpar=(Matpar*)hmatpar.delivers();
    46314614        matice=(Matice*)hmatice.delivers();
    4632         numpar=(Numpar*)hnumpar.delivers();
    46334615
    46344616        /* Get node coordinates and dof list: */
     
    46364618        GetDofList1(&doflist1[0]);
    46374619
    4638         /* Set grade_g to 0: */
    4639         for(i=0;i<numgrids;i++) grade_g[i]=0.0;
    4640 
    4641         /* recover input parameters: */
    4642         inputs->Recover("drag",&this->properties.k[0],1,dofs1,numgrids,(void**)nodes);
    4643         inputs->Recover("bed",&this->properties.b[0],1,dofs1,numgrids,(void**)nodes);
    4644         inputs->Recover("thickness",&this->properties.h[0],1,dofs1,numgrids,(void**)nodes);
    4645         if(!inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
    4646                 ISSMERROR("missing velocity input parameter");
    4647         }
    4648         if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){
    4649                 ISSMERROR("missing adjoint input parameter");
    4650         }
    4651 
    4652         /*Initialize parameter lists: */
    4653         for(i=0;i<numgrids;i++){
    4654                 vx_list[i]=vxvyvz_list[i][0];
    4655                 vy_list[i]=vxvyvz_list[i][1];
    4656                 vz_list[i]=vxvyvz_list[i][2];
    4657                 adjx_list[i]=adjxyz_list[i][0];
    4658                 adjy_list[i]=adjxyz_list[i][1];
    4659                 adjz_list[i]=adjxyz_list[i][2];
     4620        /*Build alpha_complement_list: */
     4621        if (drag_type==2){
     4622
     4623                /*Allocate friction object: */
     4624                Friction* friction=NewFriction();
     4625
     4626                inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
     4627                inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
     4628                inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],3,VzAverageEnum);
     4629                inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum);
     4630                inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum);
     4631                inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum);
     4632                inputs->GetParameterValue(&drag_p,DragPEnum);
     4633                inputs->GetParameterValue(&drag_q,DragQEnum);
     4634
     4635
     4636                /*Initialize all fields: */
     4637                friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
     4638                strcpy(friction->element_type,"2d");
     4639
     4640                friction->gravity=matpar->GetG();
     4641                friction->rho_ice=matpar->GetRhoIce();
     4642                friction->rho_water=matpar->GetRhoWater();
     4643                friction->K=&dragcoefficient_list[0];
     4644                friction->bed=&bed_list[0];
     4645                friction->thickness=&thickness_list[0];
     4646                friction->vx=&vx_list[0];
     4647                friction->vy=&vy_list[0];
     4648                friction->p=drag_p;
     4649                friction->q=drag_q;
     4650
     4651
     4652                if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");
     4653
     4654                /*Compute alpha complement: */
     4655                FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
     4656
     4657                /*Erase friction object: */
     4658                DeleteFriction(&friction);
     4659        }
     4660        else{
     4661                alpha_complement_list[0]=0;
     4662                alpha_complement_list[1]=0;
     4663                alpha_complement_list[2]=0;
    46604664        }
    46614665
     
    46714675                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    46724676
    4673                 /*Build alpha_complement_list: */
    4674                 if (drag_type==2){
    4675 
    4676                         /*Allocate friction object: */
    4677                         Friction* friction=NewFriction();
    4678 
    4679                         /*Initialize all fields: */
    4680                         friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    4681                         strcpy(friction->element_type,"2d");
    4682 
    4683                         friction->gravity=matpar->GetG();
    4684                         friction->rho_ice=matpar->GetRhoIce();
    4685                         friction->rho_water=matpar->GetRhoWater();
    4686                         friction->K=&this->properties.k[0];
    4687                         friction->bed=&this->properties.b[0];
    4688                         friction->thickness=&this->properties.h[0];
    4689                         friction->velocities=&vxvyvz_list[0][0];
    4690                         friction->p=this->properties.p;
    4691                         friction->q=this->properties.q;
    4692 
    4693                         if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");
    4694 
    4695                         /*Compute alpha complement: */
    4696                         FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
    4697 
    4698                         /*Erase friction object: */
    4699                         DeleteFriction(&friction);
    4700                 }
    4701                 else{
    4702                         alpha_complement_list[0]=0;
    4703                         alpha_complement_list[1]=0;
    4704                         alpha_complement_list[2]=0;
    4705                 }
    4706 
    47074677                /*Recover alpha_complement and k: */
    47084678                GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
    4709                 GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);
     4679                inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
    47104680
    47114681                /*recover lambda mu and xi: */
    4712                 GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
    4713                 GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
    4714                 GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);
     4682                inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum);
     4683                inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum);
     4684                inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum);
    47154685
    47164686                /*recover vx vy and vz: */
    4717                 GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
    4718                 GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
    4719                 GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);
     4687                inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum);
     4688                inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum);
     4689                inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum);
    47204690
    47214691                /*Get normal vecyor to the bed */
     
    47364706
    47374707                /*Get k derivative: dk/dx */
    4738                 GetParameterDerivativeValue(&dk[0], &this->properties.k[0],&xyz_list[0][0], gauss_l1l2l3);
     4708                inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
    47394709
    47404710                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     
    47484718
    47494719                        //Add regularization term
    4750                         grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     4720                        grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    47514721
    47524722                        //min dampening
    4753                         if(drag<numpar->cm_mindmp_value){
    4754                                 grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4723                        if(drag<cm_mindmp_value){
     4724                                grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
    47554725                        }
    47564726
    47574727                        //max dampening
    4758                         if(drag>numpar->cm_maxdmp_value){
    4759                                 grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4728                        if(drag>cm_maxdmp_value){
     4729                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
    47604730                        }
    47614731                }
     
    48394809
    48404810        /*get thickness and velocity at two segment extremities: */
    4841         GetParameterValue(&h1, &this->properties.h[0],gauss_1);
    4842         GetParameterValue(&h2, &this->properties.h[0],gauss_2);
    4843         GetParameterValue(&vx1, &vx_list[0],gauss_1);
    4844         GetParameterValue(&vy1, &vy_list[0],gauss_1);
    4845         GetParameterValue(&vx2, &vx_list[0],gauss_2);
    4846         GetParameterValue(&vy2, &vy_list[0],gauss_2);
    4847 
     4811        inputs->GetParameterValue(&h1, &gauss_1[0],ThicknessEnum);
     4812        inputs->GetParameterValue(&h2, &gauss_2[0],ThicknessEnum);
     4813        inputs->GetParameterValue(&vx1, &gauss_1[0],VxEnum);
     4814        inputs->GetParameterValue(&vx2, &gauss_2[0],VxEnum);
     4815        inputs->GetParameterValue(&vy1, &gauss_1[0],VyEnum);
     4816        inputs->GetParameterValue(&vy2, &gauss_2[0],VyEnum);
    48484817
    48494818        mass_flux= rho_ice*length*( 
     
    48714840
    48724841        /* grid data: */
    4873         double vxvy_list[numgrids][2];
    48744842        double vx_list[numgrids];
    48754843        double vy_list[numgrids];
    4876         double obs_vxvy_list[numgrids][2];
    48774844        double obs_vx_list[numgrids];
    48784845        double obs_vy_list[numgrids];
     
    48884855        double  gauss_weight;
    48894856        double  gauss_l1l2l3[3];
     4857        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    48904858
    48914859        /* parameters: */
     
    49064874        Matpar* matpar=NULL;
    49074875        Matice* matice=NULL;
    4908         Numpar* numpar=NULL;
    49094876
    49104877        /*inputs: */
    49114878        bool onwater;
     4879       
     4880        double  meanvel, epsvel;
    49124881
    49134882        /*retrieve inputs :*/
     
    49214890        matpar=(Matpar*)hmatpar.delivers();
    49224891        matice=(Matice*)hmatice.delivers();
    4923         numpar=(Numpar*)hnumpar.delivers();
    49244892
    49254893        /* Get node coordinates and dof list: */
     
    49274895
    49284896        /* Recover input data: */
    4929         if(!inputs->Recover("fit",&fit)) ISSMERROR(" missing fit input parameter");
    4930         if(fit==3 && !inputs->Recover("surfacearea",&S)){
    4931                 ISSMERROR("missing surface area input parameter");
    4932         }
    4933         if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
    4934                 ISSMERROR("missing velocity_obs input parameter");
    4935         }
    4936         if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){
    4937                 ISSMERROR("missing velocity input parameter");
    4938         }
    4939         if(!inputs->Recover("weights",&weights_list[0],1,dofs1,numgrids,(void**)nodes)){
    4940                 ISSMERROR("missing weights input parameter");
    4941         }
    4942 
    4943         /*Initialize velocities: */
    4944         for(i=0;i<numgrids;i++){
    4945                 obs_vx_list[i]=obs_vxvy_list[i][0];
    4946                 obs_vy_list[i]=obs_vxvy_list[i][1];
    4947                 vx_list[i]=vxvy_list[i][0];
    4948                 vy_list[i]=vxvy_list[i][1];
    4949         }
    4950 
     4897        inputs->GetParameterValue(&fit,FitEnum);
     4898        if(fit==3){
     4899                inputs->GetParameterValue(&S,SurfaceAreaEnum);
     4900        }
     4901        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     4902        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     4903        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     4904        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     4905        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     4906
     4907        /*retrieve some parameters: */
     4908        this->parameters->FindParam(&meanvel,"meanvel");
     4909        this->parameters->FindParam(&epsvel,"epsvel");
     4910       
    49514911        /* Compute Misfit at the 3 nodes
    49524912         * Here we integrate linearized functions:
     
    49794939                 */
    49804940                for (i=0;i<numgrids;i++){
    4981                         scalex=pow(numpar->meanvel/(obs_vx_list[i]+numpar->epsvel),(double)2);
    4982                         scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),(double)2);
     4941                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
     4942                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
    49834943                        if(obs_vx_list[i]==0)scalex=0;
    49844944                        if(obs_vy_list[i]==0)scaley=0;
     
    49954955                */
    49964956                for (i=0;i<numgrids;i++){
    4997                         velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid velocity being nil.
    4998                         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.
    4999                         misfit_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
     4957                        velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
     4958                        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.
     4959                        misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);
    50004960                }
    50014961        }
     
    50204980                 */
    50214981                for (i=0;i<numgrids;i++){
    5022                         misfit_list[i]=0.5*pow(numpar->meanvel,(double)2)*(
    5023                           pow(log((fabs(vx_list[i])+numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)),(double)2) +
    5024                           pow(log((fabs(vy_list[i])+numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)),(double)2) );
     4982                        misfit_list[i]=0.5*pow(meanvel,(double)2)*(
     4983                          pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
     4984                          pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) );
    50254985                }
    50264986        }
  • issm/trunk/src/c/objects/Tria.h

    r3612 r3620  
    3131                Hook hmatice; //hook to 1 matice
    3232                Hook hmatpar; //hook to 1 matpar
    33                 Hook hnumpar; //hook to 1 numpar
    3433
    35                 Inputs* inputs;
     34                Parameters* parameters; //pointer to solution parameters
     35                Inputs*  inputs;
    3636       
    3737        public:
     
    3939                /*FUNCTION constructors, destructors {{{1*/
    4040                Tria();
    41                 Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id);
    42                 Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs);
     41                Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id);
     42                Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* parameters, Inputs* tria_inputs);
    4343                Tria(int i, IoModel* iomodel);
    4444                ~Tria();
    4545                /*}}}*/
    4646                /*FUNCTION object management {{{1*/
    47                 void  Configure(void* loads,void* nodes,void* materials,void* parameters);
     47                void  Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters);
    4848                Object* copy();
    4949                void  DeepEcho();
Note: See TracChangeset for help on using the changeset viewer.