Ignore:
Timestamp:
12/08/20 08:45:53 (4 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 25834

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        2323issm_ocean
        2424issm_dakota
         25issm_post

  • issm/trunk/src/c/analyses/AdjointHorizAnalysis.cpp

    r24686 r25836  
    44#include "../shared/shared.h"
    55#include "../modules/modules.h"
    6 #include "../classes/Inputs2/DatasetInput2.h"
     6#include "../classes/Inputs/DatasetInput.h"
    77
    88/*Model processing*/
    … …  
    1919        _error_("not implemented");
    2020}/*}}}*/
    21 void AdjointHorizAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     21void AdjointHorizAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    2222           _error_("not implemented yet");
    2323}/*}}}*/
    … …  
    8787        /*Retrieve all inputs and parameters*/
    8888        element->GetVerticesCoordinates(&xyz_list);
    89         Input2* vx_input = element->GetInput2(VxEnum);_assert_(vx_input);
    90         Input2* vy_input = element->GetInput2(VyEnum);_assert_(vy_input);
    91         Input2* vz_input = NULL;
     89        Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
     90        Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
     91        Input* vz_input = NULL;
    9292        if(dim==3){
    93                 vz_input = element->GetInput2(VzEnum);
     93                vz_input = element->GetInput(VzEnum);
    9494        }
    9595        else{
    … …  
    102102        /* Start  looping on the number of gaussian points: */
    103103        Gauss* gauss=element->NewGauss(5);
    104         for(int ig=gauss->begin();ig<gauss->end();ig++){
    105                 gauss->GaussPoint(ig);
     104        while(gauss->next()){
    106105
    107106                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    171170        /*Retrieve all inputs and parameters*/
    172171        element->GetVerticesCoordinates(&xyz_list);
    173         Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
    174         Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
     172        Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
     173        Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
    175174
    176175        /*Allocate dbasis*/
    … …  
    179178        /* Start  looping on the number of gaussian points: */
    180179        Gauss* gauss=element->NewGauss(5);
    181         for(int ig=gauss->begin();ig<gauss->end();ig++){
    182                 gauss->GaussPoint(ig);
     180        while(gauss->next()){
    183181
    184182                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    249247                case Domain3DEnum:
    250248                        if(!element->IsOnBase()) return NULL;
    251                         basalelement = element->SpawnBasalElement();
     249                        basalelement = element->SpawnBasalElement(true);
    252250                        break;
    253251                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    … …  
    276274        /*Retrieve all inputs and parameters*/
    277275        basalelement->GetVerticesCoordinates(&xyz_list);
    278         Input2* vx_input        = basalelement->GetInput2(VxEnum);       _assert_(vx_input);
    279         Input2* vy_input        = basalelement->GetInput2(VyEnum);       _assert_(vy_input);
    280         Input2* thickness_input = basalelement->GetInput2(ThicknessEnum); _assert_(thickness_input);
     276        Input* vx_input        = basalelement->GetInput(VxEnum);       _assert_(vx_input);
     277        Input* vy_input        = basalelement->GetInput(VyEnum);       _assert_(vy_input);
     278        Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input);
    281279
    282280        /*Allocate dbasis*/
    … …  
    285283        /* Start  looping on the number of gaussian points: */
    286284        Gauss* gauss=basalelement->NewGauss(2);
    287         for(int ig=gauss->begin();ig<gauss->end();ig++){
    288                 gauss->GaussPoint(ig);
     285        while(gauss->next()){
    289286
    290287                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    375372        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    376373        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    377         DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    378         Input2* vx_input      = element->GetInput2(VxEnum);                                 _assert_(vx_input);
    379         Input2* vxobs_input   = element->GetInput2(InversionVxObsEnum);                    _assert_(vxobs_input);
    380         Input2* vy_input    = NULL;
    381         Input2* vyobs_input = NULL;
     374        DatasetInput* weights_input = element->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     375        Input* vx_input      = element->GetInput(VxEnum);             _assert_(vx_input);
     376        Input* vxobs_input   = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
     377        Input* vy_input    = NULL;
     378        Input* vyobs_input = NULL;
    382379        if(domaintype!=Domain2DverticalEnum){
    383                 vy_input      = element->GetInput2(VyEnum);                                 _assert_(vy_input);
    384                 vyobs_input   = element->GetInput2(InversionVyObsEnum);                    _assert_(vyobs_input);
     380                vy_input      = element->GetInput(VyEnum);             _assert_(vy_input);
     381                vyobs_input   = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
    385382        }
    386383        IssmDouble epsvel  = 2.220446049250313e-16;
    … …  
    388385
    389386        /*Get Surface if required by one response*/
     387        Input* S_input = NULL;
    390388        for(int resp=0;resp<num_responses;resp++){
    391389                if(responses[resp]==SurfaceAverageVelMisfitEnum){
    392                         element->GetInputValue(&S,SurfaceAreaEnum); break;
     390                        S_input = element->GetInput(SurfaceAreaEnum);  _assert_(S_input); break;
    393391                }
    394392        }
    … …  
    396394        /* Start  looping on the number of gaussian points: */
    397395        Gauss* gauss=element->NewGaussTop(4);
    398         for(int ig=gauss->begin();ig<gauss->end();ig++){
    399                 gauss->GaussPoint(ig);
     396        while(gauss->next()){
    400397
    401398                element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
    … …  
    507504                                         *        du      S  2 sqrt(...)           obs
    508505                                         */
     506                                        S_input->GetInputValue(&S,gauss);
    509507                                        for(i=0;i<vnumnodes;i++){
    510508                                                if (domaintype!=Domain2DverticalEnum){
    … …  
    612610        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    613611        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    614         DatasetInput2* weights_input = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    615         Input2* vx_input      = element->GetInput2(VxEnum);                                 _assert_(vx_input);
    616         Input2* vxobs_input   = element->GetInput2(InversionVxObsEnum);                     _assert_(vxobs_input);
    617         Input2* vy_input=NULL;
    618         Input2* vyobs_input=NULL;
     612        DatasetInput* weights_input = element->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     613        Input* vx_input      = element->GetInput(VxEnum);                                 _assert_(vx_input);
     614        Input* vxobs_input   = element->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
     615        Input* vy_input=NULL;
     616        Input* vyobs_input=NULL;
    619617        if(domaintype!=Domain2DverticalEnum){
    620                 vy_input      = element->GetInput2(VyEnum);                                 _assert_(vy_input);
    621                 vyobs_input   = element->GetInput2(InversionVyObsEnum);                     _assert_(vyobs_input);
     618                vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
     619                vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
    622620        }
    623621        IssmDouble epsvel  = 2.220446049250313e-16;
    … …  
    625623
    626624        /*Get Surface if required by one response*/
     625        Input* S_input = NULL;
    627626        for(int resp=0;resp<num_responses;resp++){
    628627                if(responses[resp]==SurfaceAverageVelMisfitEnum){
    629                         element->GetInputValue(&S,SurfaceAreaEnum); break;
     628                        S_input = element->GetInput(SurfaceAreaEnum);  _assert_(S_input); break;
    630629                }
    631630        }
    … …  
    633632        /* Start  looping on the number of gaussian points: */
    634633        Gauss* gauss=element->NewGaussTop(4);
    635         for(int ig=gauss->begin();ig<gauss->end();ig++){
    636                 gauss->GaussPoint(ig);
     634        while(gauss->next()){
    637635
    638636                element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
    … …  
    742740                                         *        du      S  2 sqrt(...)           obs
    743741                                         */
     742                                        S_input->GetInputValue(&S,gauss);
    744743                                        for(i=0;i<numnodes;i++){
    745744                                                if(domaintype!=Domain2DverticalEnum){
    … …  
    839838                case Domain3DEnum:
    840839                        if(!element->IsOnBase()) return NULL;
    841                         basalelement = element->SpawnBasalElement();
     840                        basalelement = element->SpawnBasalElement(true);
    842841                        break;
    843842                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    … …  
    863862        basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    864863        basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    865         DatasetInput2* weights_input = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    866         Input2* vx_input      = basalelement->GetInput2(VxEnum);                                               _assert_(vx_input);
    867         Input2* vxobs_input   = basalelement->GetInput2(InversionVxObsEnum);                                   _assert_(vxobs_input);
    868         Input2* vy_input=NULL;
    869         Input2* vyobs_input=NULL;
     864        DatasetInput* weights_input = basalelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     865        Input* vx_input      = basalelement->GetInput(VxEnum);                                               _assert_(vx_input);
     866        Input* vxobs_input   = basalelement->GetInput(InversionVxObsEnum);                                   _assert_(vxobs_input);
     867        Input* vy_input=NULL;
     868        Input* vyobs_input=NULL;
    870869        if(domaintype!=Domain2DverticalEnum){
    871                 vy_input      = basalelement->GetInput2(VyEnum);              _assert_(vy_input);
    872                 vyobs_input   = basalelement->GetInput2(InversionVyObsEnum);  _assert_(vyobs_input);
     870                vy_input      = basalelement->GetInput(VyEnum);              _assert_(vy_input);
     871                vyobs_input   = basalelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
    873872        }
    874873        IssmDouble epsvel  = 2.220446049250313e-16;
    … …  
    876875
    877876        /*Get Surface if required by one response*/
     877        Input* S_input = NULL;
    878878        for(int resp=0;resp<num_responses;resp++){
    879879                if(responses[resp]==SurfaceAverageVelMisfitEnum){
    880                         basalelement->GetInputValue(&S,SurfaceAreaEnum); break;
     880                        S_input = element->GetInput(SurfaceAreaEnum);  _assert_(S_input); break;
    881881                }
    882882        }
    … …  
    884884        /* Start  looping on the number of gaussian points: */
    885885        Gauss* gauss=basalelement->NewGauss(4);
    886         for(int ig=gauss->begin();ig<gauss->end();ig++){
    887                 gauss->GaussPoint(ig);
     886        while(gauss->next()){
    888887
    889888                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    993992                                         *        du      S  2 sqrt(...)           obs
    994993                                         */
     994                                        S_input->GetInputValue(&S,gauss);
    995995                                        for(i=0;i<numnodes;i++){
    996996                                                if(domaintype!=Domain2DverticalEnum){
    … …  
    10721072           _error_("not implemented yet");
    10731073}/*}}}*/
    1074 void           AdjointHorizAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     1074void           AdjointHorizAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    10751075        /*The gradient of the cost function is calculated in 2 parts.
    10761076         *
    … …  
    10971097        if(control_type!=MaterialsRheologyBbarEnum &&
    10981098                control_type!=FrictionCoefficientEnum   &&
    1099                 control_type!=FrictionAsEnum   &&
     1099                control_type!=FrictionCEnum             &&
     1100                control_type!=FrictionAsEnum            &&
    11001101                control_type!=DamageDbarEnum            &&
    11011102                control_type!=MaterialsRheologyBEnum){
    … …  
    11101111                case SurfaceLogVxVyMisfitEnum:    /*Nothing, \partial J/\partial k = 0*/ break;
    11111112                case SurfaceAverageVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
    1112                 case DragCoefficientAbsGradientEnum: GradientJDragGradient(element,gradient,control_index); break;
    1113                 case RheologyBbarAbsGradientEnum:    GradientJBbarGradient(element,gradient,control_index); break;
    1114                 case RheologyBAbsGradientEnum:       GradientJBGradient(element,gradient,control_index);    break;
    1115                 case RheologyBInitialguessMisfitEnum:  GradientJBinitial(element,gradient,control_index);    break;
     1113                case DragCoefficientAbsGradientEnum:   GradientJDragGradient(element,gradient,control_interp,control_index); break;
     1114                case RheologyBbarAbsGradientEnum:      GradientJBbarGradient(element,gradient,control_interp,control_index); break;
     1115                case RheologyBAbsGradientEnum:         GradientJBGradient(element,gradient,control_interp,control_index);    break;
     1116                case RheologyBInitialguessMisfitEnum:  GradientJBinitial(element,gradient,control_interp,control_index);     break;
    11161117                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    11171118        }
    … …  
    11201121        switch(control_type){
    11211122                case FrictionCoefficientEnum:
     1123                case FrictionCEnum:
    11221124                        switch(approximation){
    1123                                 case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break;
    1124                                 case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_index); break;
    1125                                 case HOApproximationEnum:  GradientJDragHO( element,gradient,control_index); break;
    1126                                 case FSApproximationEnum:  GradientJDragFS( element,gradient,control_index); break;
     1125                                case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_interp,control_index); break;
     1126                                case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_interp,control_index); break;
     1127                                case HOApproximationEnum:  GradientJDragHO( element,gradient,control_interp,control_index); break;
     1128                                case FSApproximationEnum:  GradientJDragFS( element,gradient,control_interp,control_index); break;
    11271129                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11281130                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11311133                case FrictionAsEnum:
    11321134                        switch(approximation){
    1133                                 case SSAApproximationEnum: GradientJDragHydroSSA(element,gradient,control_index); break;
    1134                                 case L1L2ApproximationEnum:GradientJDragHydroL1L2(element,gradient,control_index); break;
    1135                                 case HOApproximationEnum:  GradientJDragHydroHO( element,gradient,control_index); break;
    1136                                 case FSApproximationEnum:  GradientJDragHydroFS( element,gradient,control_index); break;
     1135                                case SSAApproximationEnum: GradientJDragHydroSSA(element,gradient,control_interp,control_index); break;
     1136                                case L1L2ApproximationEnum:GradientJDragHydroL1L2(element,gradient,control_interp,control_index); break;
     1137                                case HOApproximationEnum:  GradientJDragHydroHO( element,gradient,control_interp,control_index); break;
     1138                                case FSApproximationEnum:  GradientJDragHydroFS( element,gradient,control_interp,control_index); break;
    11371139                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11381140                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11411143                case MaterialsRheologyBbarEnum:
    11421144                        switch(approximation){
    1143                                 case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_index); break;
    1144                                 case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_index); break;
    1145                                 case HOApproximationEnum:  GradientJBbarHO( element,gradient,control_index); break;
    1146                                 case FSApproximationEnum:  GradientJBbarFS( element,gradient,control_index); break;
     1145                                case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_interp,control_index); break;
     1146                                case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_interp,control_index); break;
     1147                                case HOApproximationEnum:  GradientJBbarHO( element,gradient,control_interp,control_index); break;
     1148                                case FSApproximationEnum:  GradientJBbarFS( element,gradient,control_interp,control_index); break;
    11471149                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11481150                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11511153                case MaterialsRheologyBEnum:
    11521154                        switch(approximation){
    1153                                 case SSAApproximationEnum: GradientJBSSA(element,gradient,control_index); break;
    1154                                 case HOApproximationEnum:  GradientJBHO( element,gradient,control_index); break;
    1155                                 case FSApproximationEnum:  GradientJBFS( element,gradient,control_index); break;
     1155                                case SSAApproximationEnum: GradientJBSSA(element,gradient,control_interp,control_index); break;
     1156                                case HOApproximationEnum:  GradientJBHO( element,gradient,control_interp,control_index); break;
     1157                                case FSApproximationEnum:  GradientJBFS( element,gradient,control_interp,control_index); break;
    11561158                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11571159                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11601162                case DamageDbarEnum:
    11611163                        switch(approximation){
    1162                                 case SSAApproximationEnum: GradientJDSSA(element,gradient,control_index); break;
     1164                                case SSAApproximationEnum: GradientJDSSA(element,gradient,control_interp,control_index); break;
    11631165                                case NoneApproximationEnum: /*Gradient is 0*/                 break;
    11641166                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11721174
    11731175}/*}}}*/
    1174 void           AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1176void           AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    11751177        /*WARNING: We use SSA as an estimate for now*/
    1176         this->GradientJBbarSSA(element,gradient,control_index);
    1177 }/*}}}*/
    1178 void           AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1178        this->GradientJBbarSSA(element,gradient,control_interp,control_index);
     1179}/*}}}*/
     1180void           AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    11791181
    11801182        /*Intermediaries*/
    … …  
    11941196                case Domain3DEnum:
    11951197                        if(!element->IsOnBase()) return;
    1196                         basalelement = element->SpawnBasalElement();
     1198                        basalelement = element->SpawnBasalElement(true);
    11971199                        break;
    11981200                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    … …  
    12151217        basalelement->GetVerticesCoordinates(&xyz_list);
    12161218        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1217         Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
    1218         DatasetInput2* weights_input   = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1219        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
     1220        DatasetInput* weights_input   = basalelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    12191221
    12201222        /* Start  looping on the number of gaussian points: */
    12211223        Gauss* gauss=basalelement->NewGauss(2);
    1222         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1223                 gauss->GaussPoint(ig);
     1224        while(gauss->next()){
    12241225
    12251226                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    12511252        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    12521253}/*}}}*/
    1253 void           AdjointHorizAnalysis::GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1254void           AdjointHorizAnalysis::GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    12541255
    12551256           /*Same as SSA*/
    1256            return this->GradientJBbarSSA(element,gradient,control_index);
    1257 }/*}}}*/
    1258 void           AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1257           return this->GradientJBbarSSA(element,gradient,control_interp,control_index);
     1258}/*}}}*/
     1259void           AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    12591260
    12601261        /*WARNING: We use SSA as an estimate for now*/
    1261         this->GradientJBbarSSA(element,gradient,control_index);
    1262 }/*}}}*/
    1263 void           AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1262        this->GradientJBbarSSA(element,gradient,control_interp,control_index);
     1263}/*}}}*/
     1264void           AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    12641265
    12651266        /*Intermediaries*/
    … …  
    12811282                case Domain3DEnum:
    12821283                        if(!element->IsOnBase()) return;
    1283                         basalelement = element->SpawnBasalElement();
     1284                        basalelement = element->SpawnBasalElement(true);
    12841285                        dim          = 2;
    12851286                        break;
    … …  
    13041305        basalelement->GetVerticesCoordinates(&xyz_list);
    13051306        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1306         Input2* thickness_input = basalelement->GetInput2(ThicknessEnum);             _assert_(thickness_input);
    1307         Input2* vx_input        = basalelement->GetInput2(VxEnum);                    _assert_(vx_input);
    1308         Input2* vy_input        = basalelement->GetInput2(VyEnum);                    _assert_(vy_input);
    1309         Input2* adjointx_input  = basalelement->GetInput2(AdjointxEnum);              _assert_(adjointx_input);
    1310         Input2* adjointy_input  = basalelement->GetInput2(AdjointyEnum);              _assert_(adjointy_input);
    1311         Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
     1307        Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
     1308        Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
     1309        Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
     1310        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
     1311        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
     1312        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    13121313
    13131314        /* Start  looping on the number of gaussian points: */
    13141315        Gauss* gauss=basalelement->NewGauss(4);
    1315         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1316                 gauss->GaussPoint(ig);
     1316        while(gauss->next()){
    13171317
    13181318                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    13281328
    13291329                /*Build gradient vector (actually -dJ/dB): */
    1330                 for(int i=0;i<numvertices;i++){
    1331                         ge[i]+=-dmudB*thickness*(
     1330                if(control_interp==P1Enum){
     1331                        for(int i=0;i<numvertices;i++){
     1332                                ge[i]+=-dmudB*thickness*(
     1333                                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     1334                                                        )*Jdet*gauss->weight*basis[i];
     1335                                _assert_(!xIsNan<IssmDouble>(ge[i]));
     1336                        }
     1337                }
     1338                else if(control_interp==P0Enum){
     1339                        ge[0]+=-dmudB*thickness*(
    13321340                                                (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
    1333                                                 )*Jdet*gauss->weight*basis[i];
    1334                         _assert_(!xIsNan<IssmDouble>(ge[i]));
    1335                 }
    1336         }
    1337         gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1341                                                )*Jdet*gauss->weight;
     1342                        _assert_(!xIsNan<IssmDouble>(ge[0]));
     1343                }
     1344                else{
     1345                        _error_("not supported");
     1346                }
     1347        }
     1348        if(control_interp==P1Enum){
     1349                gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1350        }
     1351        else if(control_interp==P0Enum){
     1352                gradient->SetValue(vertexpidlist[0],ge[0],ADD_VAL);
     1353        }
     1354        else{
     1355                _error_("not supported");
     1356        }
    13381357
    13391358        /*Clean up and return*/
    … …  
    13451364        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    13461365}/*}}}*/
    1347 void           AdjointHorizAnalysis::GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1366void           AdjointHorizAnalysis::GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    13481367        /*WARNING: We use HO as an estimate for now*/
    1349         this->GradientJBHO(element,gradient,control_index);
    1350 }/*}}}*/
    1351 void           AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1368        this->GradientJBHO(element,gradient,control_interp,control_index);
     1369}/*}}}*/
     1370void           AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    13521371
    13531372        /*Intermediaries*/
    … …  
    13821401        element->GetVerticesCoordinates(&xyz_list);
    13831402        element->GradientIndexing(&vertexpidlist[0],control_index);
    1384         Input2* rheology_input = element->GetInput2(MaterialsRheologyBEnum);              _assert_(rheology_input);
    1385         DatasetInput2* weights_input   = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1403        Input* rheology_input = element->GetInput(MaterialsRheologyBEnum);              _assert_(rheology_input);
     1404        DatasetInput* weights_input   = element->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    13861405        /* Start  looping on the number of gaussian points: */
    13871406        Gauss* gauss=element->NewGauss(2);
    1388         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1389                 gauss->GaussPoint(ig);
     1407        while(gauss->next()){
    13901408                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    13911409                element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
    … …  
    14151433        delete gauss;
    14161434}/*}}}*/
    1417 void           AdjointHorizAnalysis::GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1435void           AdjointHorizAnalysis::GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    14181436
    14191437        /*Intermediaries*/
    … …  
    14481466        element->GetVerticesCoordinates(&xyz_list);
    14491467        element->GradientIndexing(&vertexpidlist[0],control_index);
    1450         Input2* rheology_input  = element->GetInput2(MaterialsRheologyBbarEnum);              _assert_(rheology_input);
    1451         Input2* rheology0_input = element->GetInput2(RheologyBInitialguessEnum);              _assert_(rheology0_input);
    1452         DatasetInput2* weights_input   = element->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1468        Input* rheology_input  = element->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheology_input);
     1469        Input* rheology0_input = element->GetInput(RheologyBInitialguessEnum);              _assert_(rheology0_input);
     1470        DatasetInput* weights_input   = element->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    14531471
    14541472        /* Start  looping on the number of gaussian points: */
    14551473        Gauss* gauss=element->NewGauss(2);
    1456         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1457                 gauss->GaussPoint(ig);
     1474        while(gauss->next()){
    14581475                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    14591476                element->NodalFunctionsP1(basis,gauss);
    … …  
    14791496        delete gauss;
    14801497}/*}}}*/
    1481 void           AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1498void           AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    14821499        /*Intermediaries*/
    14831500        int      domaintype,dim;
    … …  
    15031520        element->GetVerticesCoordinates(&xyz_list);
    15041521        element->GradientIndexing(&vertexpidlist[0],control_index);
    1505         Input2* thickness_input = element->GetInput2(ThicknessEnum);             _assert_(thickness_input);
    1506         Input2* vx_input        = element->GetInput2(VxEnum);                    _assert_(vx_input);
    1507         Input2* vy_input        = NULL;
    1508         Input2* adjointx_input  = element->GetInput2(AdjointxEnum);              _assert_(adjointx_input);
    1509         Input2* adjointy_input  = NULL;
    1510         Input2* rheologyb_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(rheologyb_input);
     1522        Input* thickness_input = element->GetInput(ThicknessEnum);             _assert_(thickness_input);
     1523        Input* vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
     1524        Input* vy_input        = NULL;
     1525        Input* adjointx_input  = element->GetInput(AdjointxEnum);              _assert_(adjointx_input);
     1526        Input* adjointy_input  = NULL;
     1527        Input* rheologyb_input = element->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input);
    15111528        if(domaintype!=Domain2DverticalEnum){
    1512                 vy_input        = element->GetInput2(VyEnum);                   _assert_(vy_input);
    1513                 adjointy_input  = element->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
     1529                vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     1530                adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
    15141531        }
    15151532        /* Start  looping on the number of gaussian points: */
    15161533        Gauss* gauss=element->NewGauss(4);
    1517         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1518                 gauss->GaussPoint(ig);
     1534        while(gauss->next()){
    15191535
    15201536                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    15551571        delete gauss;
    15561572}/*}}}*/
    1557 void           AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1573void           AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    15581574
    15591575        /*Intermediaries*/
    … …  
    15981614        basalelement->GetVerticesCoordinates(&xyz_list);
    15991615        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1600         Input2* thickness_input = basalelement->GetInput2(ThicknessEnum);             _assert_(thickness_input);
    1601         Input2* vx_input        = basalelement->GetInput2(VxEnum);                    _assert_(vx_input);
    1602         Input2* vy_input        = basalelement->GetInput2(VyEnum);                    _assert_(vy_input);
    1603         Input2* adjointx_input  = basalelement->GetInput2(AdjointxEnum);              _assert_(adjointx_input);
    1604         Input2* adjointy_input  = basalelement->GetInput2(AdjointyEnum);              _assert_(adjointy_input);
    1605         Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBEnum); _assert_(rheologyb_input);
     1616        Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
     1617        Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
     1618        Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
     1619        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
     1620        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
     1621        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input);
    16061622
    16071623        /* Start  looping on the number of gaussian points: */
    16081624        Gauss* gauss=basalelement->NewGauss(4);
    1609         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1610                 gauss->GaussPoint(ig);
     1625        while(gauss->next()){
    16111626
    16121627                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    16391654        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    16401655}/*}}}*/
    1641 void           AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1656void           AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    16421657
    16431658        /*return if floating (gradient is 0)*/
    … …  
    16461661        /*Intermediaries*/
    16471662        int      domaintype,dim;
     1663        int frictionlaw;
    16481664        Element* basalelement;
    16491665
    … …  
    16841700        basalelement->GetVerticesCoordinates(&xyz_list);
    16851701        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1686         Input2* dragcoefficient_input = basalelement->GetInput2(FrictionCoefficientEnum);                _assert_(dragcoefficient_input);
    1687         DatasetInput2* weights_input         = basalelement->GetDatasetInput2(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1702
     1703        /* get the friction law: if 2-Weertman, 11-Schoof, use a special name for the coefficient*/
     1704        element->FindParam(&frictionlaw, FrictionLawEnum);
     1705        Input* dragcoefficient_input;
     1706        switch(frictionlaw) {
     1707                case 2:
     1708                case 11:
     1709                        dragcoefficient_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoefficient_input);
     1710                        break;
     1711                default:
     1712                        dragcoefficient_input = basalelement->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
     1713        }
     1714
     1715        DatasetInput* weights_input         = basalelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    16881716
    16891717        /* Start  looping on the number of gaussian points: */
    16901718        Gauss* gauss=basalelement->NewGauss(2);
    1691         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1692                 gauss->GaussPoint(ig);
     1719        while(gauss->next()){
    16931720
    16941721                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    17211748
    17221749}/*}}}*/
    1723 void           AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1750void           AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    17241751
    17251752        /*return if floating or not on bed (gradient is 0)*/
    … …  
    17531780        element->GetVerticesCoordinatesBase(&xyz_list_base);
    17541781        element->GradientIndexing(&vertexpidlist[0],control_index);
    1755         Input2* vx_input        = element->GetInput2(VxEnum);                   _assert_(vx_input);
    1756         Input2* vy_input        = element->GetInput2(VyEnum);                   _assert_(vy_input);
    1757         Input2* adjointx_input  = element->GetInput2(AdjointxEnum);             _assert_(adjointx_input);
    1758         Input2* adjointy_input  = element->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
    1759         Input2* vz_input        = NULL;
    1760         Input2* adjointz_input  = NULL;
     1782        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
     1783        Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     1784        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     1785        Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
     1786        Input* vz_input        = NULL;
     1787        Input* adjointz_input  = NULL;
    17611788        if(domaintype!=Domain2DverticalEnum){
    1762                 vz_input        = element->GetInput2(VzEnum);                   _assert_(vy_input);
    1763                 adjointz_input  = element->GetInput2(AdjointzEnum);             _assert_(adjointz_input);
    1764         }
    1765         Input2* dragcoeff_input = element->GetInput2(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
     1789                vz_input        = element->GetInput(VzEnum);                   _assert_(vy_input);
     1790                adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
     1791        }
     1792        Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
    17661793
    17671794        /* Start  looping on the number of gaussian points: */
    17681795        Gauss* gauss=element->NewGaussBase(4);
    1769         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1770                 gauss->GaussPoint(ig);
     1796        while(gauss->next()){
    17711797
    17721798                adjointx_input->GetInputValue(&lambda, gauss);
    … …  
    18171843        delete friction;
    18181844}/*}}}*/
    1819 void           AdjointHorizAnalysis::GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1845void           AdjointHorizAnalysis::GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    18201846
    18211847        /*Same as SSA*/
    1822         return this->GradientJDragSSA(element,gradient,control_index);
    1823 }/*}}}*/
    1824 void           AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1848        return this->GradientJDragSSA(element,gradient,control_interp,control_index);
     1849}/*}}}*/
     1850void           AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    18251851
    18261852        /*return if floating or not on bed (gradient is 0)*/
    … …  
    18521878        element->GetVerticesCoordinatesBase(&xyz_list_base);
    18531879        element->GradientIndexing(&vertexpidlist[0],control_index);
    1854         Input2* vx_input        = element->GetInput2(VxEnum);                   _assert_(vx_input);
    1855         Input2* vy_input        = NULL;
    1856         Input2* adjointx_input  = element->GetInput2(AdjointxEnum);             _assert_(adjointx_input);
    1857         Input2* adjointy_input  = NULL;
    1858         Input2* dragcoeff_input = element->GetInput2(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
     1880        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
     1881        Input* vy_input        = NULL;
     1882        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     1883        Input* adjointy_input  = NULL;
     1884        Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
    18591885        if(domaintype!=Domain2DverticalEnum){
    1860                 vy_input        = element->GetInput2(VyEnum);                   _assert_(vy_input);
    1861                 adjointy_input  = element->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
     1886                vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     1887                adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
    18621888        }
    18631889        /* Start  looping on the number of gaussian points: */
    18641890        Gauss* gauss=element->NewGaussBase(4);
    1865         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1866                 gauss->GaussPoint(ig);
     1891        while(gauss->next()){
    18671892
    18681893                adjointx_input->GetInputValue(&lambda, gauss);
    … …  
    18961921        delete friction;
    18971922}/*}}}*/
    1898 void           AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1923void           AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    18991924
    19001925        /*return if floating (gradient is 0)*/
    … …  
    19331958        /*Fetch number of vertices for this finite element*/
    19341959        int numvertices = basalelement->GetNumberOfVertices();
     1960        int frictionlaw;
    19351961
    19361962        /*Initialize some vectors*/
    … …  
    19451971        basalelement->GetVerticesCoordinates(&xyz_list);
    19461972        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    1947         Input2* vx_input        = basalelement->GetInput2(VxEnum);                   _assert_(vx_input);
    1948         Input2* vy_input        = basalelement->GetInput2(VyEnum);                   _assert_(vy_input);
    1949         Input2* adjointx_input  = basalelement->GetInput2(AdjointxEnum);             _assert_(adjointx_input);
    1950         Input2* adjointy_input  = basalelement->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
    1951         Input2* dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
     1973        Input* vx_input        = basalelement->GetInput(VxEnum);                   _assert_(vx_input);
     1974        Input* vy_input        = basalelement->GetInput(VyEnum);                   _assert_(vy_input);
     1975        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     1976        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);             _assert_(adjointy_input);
     1977
     1978        /* get the friction law: 1- Budd, 11-Schoof*/
     1979        element->FindParam(&frictionlaw, FrictionLawEnum);
     1980        Input* dragcoeff_input;
     1981        switch(frictionlaw) {
     1982                case 1:
     1983                        dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input);
     1984                        break;
     1985                case 2:
     1986                case 11:
     1987                        dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input);
     1988                        break;
     1989                default:
     1990                        _error_("Friction law "<< frictionlaw <<" not supported in the inversion.");
     1991        }
    19521992
    19531993        /* Start  looping on the number of gaussian points: */
    19541994        Gauss* gauss=basalelement->NewGauss(4);
    1955         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1956                 gauss->GaussPoint(ig);
     1995        while(gauss->next()){
    19571996
    19581997                adjointx_input->GetInputValue(&lambda, gauss);
    … …  
    19842023        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    19852024}/*}}}*/
    1986 
    1987 void AdjointHorizAnalysis::GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2025void           AdjointHorizAnalysis::GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    19882026
    19892027        /*return if floating or not on bed (gradient is 0)*/
    … …  
    20172055        element->GetVerticesCoordinatesBase(&xyz_list_base);
    20182056        element->GradientIndexing(&vertexpidlist[0],control_index);
    2019         Input2* vx_input        = element->GetInput2(VxEnum);                   _assert_(vx_input);
    2020         Input2* vy_input        = element->GetInput2(VyEnum);                   _assert_(vy_input);
    2021         Input2* adjointx_input  = element->GetInput2(AdjointxEnum);             _assert_(adjointx_input);
    2022         Input2* adjointy_input  = element->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
    2023         Input2* vz_input        = NULL;
    2024         Input2* adjointz_input  = NULL;
     2057        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
     2058        Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     2059        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     2060        Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
     2061        Input* vz_input        = NULL;
     2062        Input* adjointz_input  = NULL;
    20252063        if(domaintype!=Domain2DverticalEnum){
    2026                 vz_input        = element->GetInput2(VzEnum);                   _assert_(vy_input);
    2027                 adjointz_input  = element->GetInput2(AdjointzEnum);             _assert_(adjointz_input);
     2064                vz_input        = element->GetInput(VzEnum);                   _assert_(vy_input);
     2065                adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
    20282066        }
    20292067
    20302068        /* Start  looping on the number of gaussian points: */
    20312069        Gauss* gauss=element->NewGaussBase(4);
    2032         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2033                 gauss->GaussPoint(ig);
     2070        while(gauss->next()){
    20342071
    20352072                adjointx_input->GetInputValue(&lambda, gauss);
    … …  
    20792116        delete friction;
    20802117}/*}}}*/
    2081 void           AdjointHorizAnalysis::GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2118void           AdjointHorizAnalysis::GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    20822119
    20832120        /*Same as SSA*/
    2084         return this->GradientJDragSSA(element,gradient,control_index);
    2085 }/*}}}*/
    2086 void           AdjointHorizAnalysis::GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2121        return this->GradientJDragSSA(element,gradient,control_interp,control_index);
     2122}/*}}}*/
     2123void           AdjointHorizAnalysis::GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    20872124
    20882125        /*return if floating or not on bed (gradient is 0)*/
    … …  
    21142151        element->GetVerticesCoordinatesBase(&xyz_list_base);
    21152152        element->GradientIndexing(&vertexpidlist[0],control_index);
    2116         Input2* vx_input        = element->GetInput2(VxEnum);                   _assert_(vx_input);
    2117         Input2* vy_input        = NULL;
    2118         Input2* adjointx_input  = element->GetInput2(AdjointxEnum);             _assert_(adjointx_input);
    2119         Input2* adjointy_input  = NULL;
     2153        Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
     2154        Input* vy_input        = NULL;
     2155        Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
     2156        Input* adjointy_input  = NULL;
    21202157        if(domaintype!=Domain2DverticalEnum){
    2121                 vy_input        = element->GetInput2(VyEnum);                   _assert_(vy_input);
    2122                 adjointy_input  = element->GetInput2(AdjointyEnum);             _assert_(adjointy_input);
     2158                vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
     2159                adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
    21232160        }
    21242161        /* Start  looping on the number of gaussian points: */
    21252162        Gauss* gauss=element->NewGaussBase(4);
    2126         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2127                 gauss->GaussPoint(ig);
     2163        while(gauss->next()){
    21282164
    21292165                adjointx_input->GetInputValue(&lambda, gauss);
    … …  
    21562192        delete friction;
    21572193}/*}}}*/
    2158 
    2159 void AdjointHorizAnalysis::GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2194void           AdjointHorizAnalysis::GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    21602195
    21612196        /*return if floating (gradient is 0)*/
    … …  
    22072242        basalelement->GetVerticesCoordinates(&xyz_list);
    22082243        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    2209         Input2* vx_input        = basalelement->GetInput2(VxEnum);          _assert_(vx_input);
    2210         Input2* vy_input        = basalelement->GetInput2(VyEnum);          _assert_(vy_input);
    2211         Input2* adjointx_input  = basalelement->GetInput2(AdjointxEnum);    _assert_(adjointx_input);
    2212         Input2* adjointy_input  = basalelement->GetInput2(AdjointyEnum);    _assert_(adjointy_input);
     2244        Input* vx_input        = basalelement->GetInput(VxEnum);          _assert_(vx_input);
     2245        Input* vy_input        = basalelement->GetInput(VyEnum);          _assert_(vy_input);
     2246        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);    _assert_(adjointx_input);
     2247        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);    _assert_(adjointy_input);
    22132248
    22142249        IssmDouble  q_exp;
    … …  
    22232258
    22242259        /*Recover parameters: */
    2225         Input2* qinput = basalelement->GetInput2(FrictionQEnum);
    2226         Input2* cinput = basalelement->GetInput2(FrictionCEnum);
    2227         Input2* Asinput = basalelement->GetInput2(FrictionAsEnum);
    2228         Input2* nInput =basalelement->GetInput2(MaterialsRheologyNEnum);
    2229         Input2* Ninput = basalelement->GetInput2(FrictionEffectivePressureEnum);       
     2260        Input* qinput = basalelement->GetInput(FrictionQEnum);
     2261        Input* cinput = basalelement->GetInput(FrictionCEnum);
     2262        Input* Asinput = basalelement->GetInput(FrictionAsEnum);
     2263        Input* nInput =basalelement->GetInput(MaterialsRheologyNEnum);
     2264        Input* Ninput = basalelement->GetInput(FrictionEffectivePressureEnum); 
    22302265        /* Start  looping on the number of gaussian points: */
    22312266        Gauss* gauss=basalelement->NewGauss(4);
    2232         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2233                 gauss->GaussPoint(ig);
     2267        while(gauss->next()){
    22342268
    22352269                adjointx_input->GetInputValue(&lambda, gauss);
    … …  
    22822316        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    22832317}/*}}}*/
    2284 
    2285 void           AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2318void           AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    22862319
    22872320        /*Intermediaries*/
    … …  
    23032336                case Domain3DEnum:
    23042337                        if(!element->IsOnBase()) return;
    2305                         basalelement = element->SpawnBasalElement();
     2338                        basalelement = element->SpawnBasalElement(true);
    23062339                        dim          = 2;
    23072340                        break;
    … …  
    23262359        basalelement->GetVerticesCoordinates(&xyz_list);
    23272360        basalelement->GradientIndexing(&vertexpidlist[0],control_index);
    2328         Input2* thickness_input = basalelement->GetInput2(ThicknessEnum);             _assert_(thickness_input);
    2329         Input2* vx_input        = basalelement->GetInput2(VxEnum);                    _assert_(vx_input);
    2330         Input2* vy_input        = basalelement->GetInput2(VyEnum);                    _assert_(vy_input);
    2331         Input2* adjointx_input  = basalelement->GetInput2(AdjointxEnum);              _assert_(adjointx_input);
    2332         Input2* adjointy_input  = basalelement->GetInput2(AdjointyEnum);              _assert_(adjointy_input);
    2333         Input2* rheologyb_input = basalelement->GetInput2(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
     2361        Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
     2362        Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
     2363        Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
     2364        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
     2365        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
     2366        Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    23342367
    23352368        /* Start  looping on the number of gaussian points: */
    23362369        Gauss* gauss=basalelement->NewGauss(4);
    2337         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2338                 gauss->GaussPoint(ig);
     2370        while(gauss->next()){
    23392371
    23402372                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    24512483
    24522484        /*Add vx and vy as inputs to the tria element: */
    2453         element->AddInput2(AdjointxEnum,lambdax,element->VelocityInterpolation());
    2454         element->AddInput2(AdjointyEnum,lambday,element->VelocityInterpolation());
    2455         if(domaintype!=Domain2DverticalEnum) element->AddInput2(AdjointzEnum,lambdaz,element->VelocityInterpolation());
     2485        element->AddInput(AdjointxEnum,lambdax,element->VelocityInterpolation());
     2486        element->AddInput(AdjointyEnum,lambday,element->VelocityInterpolation());
     2487        if(domaintype!=Domain2DverticalEnum) element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation());
    24562488
    24572489        element->FindParam(&fe_FS,FlowequationFeFSEnum);
    24582490        if(fe_FS!=LATaylorHoodEnum && fe_FS!=LACrouzeixRaviartEnum)     
    2459          element->AddInput2(AdjointpEnum,lambdap,element->PressureInterpolation());     
     2491         element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation());     
    24602492
    24612493        /*Free ressources:*/
    … …  
    24962528        for(i=0;i<numnodes;i++){
    24972529                if(domaintype!=Domain2DverticalEnum){
    2498                         lambdax[i]=values[i*NDOF2+0];
    2499                         lambday[i]=values[i*NDOF2+1];
     2530                        lambdax[i]=values[i*2+0];
     2531                        lambday[i]=values[i*2+1];
    25002532                }
    25012533                else {lambdax[i]=values[i];lambday[i]=0;}
    … …  
    25082540
    25092541        /*Add vx and vy as inputs to the tria element: */
    2510         element->AddInput2(AdjointxEnum,lambdax,element->GetElementType());
    2511         element->AddInput2(AdjointyEnum,lambday,element->GetElementType());
     2542        element->AddInput(AdjointxEnum,lambdax,element->GetElementType());
     2543        element->AddInput(AdjointyEnum,lambday,element->GetElementType());
    25122544
    25132545        /*Free ressources:*/
Note: See TracChangeset for help on using the changeset viewer.