Changeset 25317


Ignore:
Timestamp:
07/31/20 09:21:13 (5 years ago)
Author:
Mathieu Morlighem
Message:

NEW: enable P0 controls

Location:
issm/trunk-jpl/src/c
Files:
100 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp

    r24335 r25317  
    116116        _error_("not implemented yet");
    117117}/*}}}*/
    118 void           AdjointBalancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     118void           AdjointBalancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    119119        /*The gradient of the cost function is calculated in 2 parts.
    120120         *

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           GradientJdHdt(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3131                void           GradientJOmega(Element* element,Vector<IssmDouble>* gradient,int control_index);

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp

    r24335 r25317  
    156156        _error_("not implemented yet");
    157157}/*}}}*/
    158 void           AdjointBalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     158void           AdjointBalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    159159        /*The gradient of the cost function is calculated in 2 parts.
    160160         *
    … …  
    194194        /*Deal with second term*/
    195195        switch(control_type){
    196                 case BalancethicknessSpcthicknessEnum:   GradientJDirichlet(element,gradient,control_index); break;
    197                 case BalancethicknessThickeningRateEnum: GradientJDhDt(element,gradient,control_index); break;
    198                 case VxEnum:                             GradientJVx(  element,gradient,control_index); break;
    199                 case VyEnum:                             GradientJVy(  element,gradient,control_index); break;
     196                case BalancethicknessSpcthicknessEnum:   GradientJDirichlet(element,gradient,control_interp,control_index); break;
     197                case BalancethicknessThickeningRateEnum: GradientJDhDt(element,gradient,control_interp,control_index); break;
     198                case VxEnum:                             GradientJVx(  element,gradient,control_interp,control_index); break;
     199                case VyEnum:                             GradientJVy(  element,gradient,control_interp,control_index); break;
    200200                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
    201201        }
    … …  
    205205
    206206}/*}}}*/
    207 void           AdjointBalancethicknessAnalysis::GradientJDirichlet(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     207void           AdjointBalancethicknessAnalysis::GradientJDirichlet(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    208208
    209209        /*Fetch number of vertices for this finite element*/
    … …  
    242242        delete Ke;
    243243}/*}}}*/
    244 void           AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     244void           AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    245245
    246246        /*Fetch number of vertices for this finite element*/
    … …  
    266266        xDelete<int>(vertexpidlist);
    267267}/*}}}*/
    268 void           AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     268void           AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    269269
    270270        /*Intermediaries*/
    … …  
    313313        delete gauss;
    314314}/*}}}*/
    315 void           AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     315void           AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    316316
    317317        /*Intermediaries*/

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    30                 void           GradientJDirichlet(Element* element,Vector<IssmDouble>* gradient,int control_index);
    31                 void           GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index);
    32                 void           GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index);
    33                 void           GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
     30                void           GradientJDirichlet(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     31                void           GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     32                void           GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     33                void           GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
    3434                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3535                void           UpdateConstraints(FemModel* femmodel);

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

    r25240 r25317  
    10781078           _error_("not implemented yet");
    10791079}/*}}}*/
    1080 void           AdjointHorizAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     1080void           AdjointHorizAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    10811081        /*The gradient of the cost function is calculated in 2 parts.
    10821082         *
    … …  
    11031103        if(control_type!=MaterialsRheologyBbarEnum &&
    11041104                control_type!=FrictionCoefficientEnum   &&
    1105                 control_type!=FrictionCEnum   &&
    1106                 control_type!=FrictionAsEnum   &&
     1105                control_type!=FrictionCEnum             &&
     1106                control_type!=FrictionAsEnum            &&
    11071107                control_type!=DamageDbarEnum            &&
    11081108                control_type!=MaterialsRheologyBEnum){
    … …  
    11171117                case SurfaceLogVxVyMisfitEnum:    /*Nothing, \partial J/\partial k = 0*/ break;
    11181118                case SurfaceAverageVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
    1119                 case DragCoefficientAbsGradientEnum: GradientJDragGradient(element,gradient,control_index); break;
    1120                 case RheologyBbarAbsGradientEnum:    GradientJBbarGradient(element,gradient,control_index); break;
    1121                 case RheologyBAbsGradientEnum:       GradientJBGradient(element,gradient,control_index);    break;
    1122                 case RheologyBInitialguessMisfitEnum:  GradientJBinitial(element,gradient,control_index);    break;
     1119                case DragCoefficientAbsGradientEnum:   GradientJDragGradient(element,gradient,control_interp,control_index); break;
     1120                case RheologyBbarAbsGradientEnum:      GradientJBbarGradient(element,gradient,control_interp,control_index); break;
     1121                case RheologyBAbsGradientEnum:         GradientJBGradient(element,gradient,control_interp,control_index);    break;
     1122                case RheologyBInitialguessMisfitEnum:  GradientJBinitial(element,gradient,control_interp,control_index);     break;
    11231123                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    11241124        }
    … …  
    11291129                case FrictionCEnum:
    11301130                        switch(approximation){
    1131                                 case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_index); break;
    1132                                 case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_index); break;
    1133                                 case HOApproximationEnum:  GradientJDragHO( element,gradient,control_index); break;
    1134                                 case FSApproximationEnum:  GradientJDragFS( element,gradient,control_index); break;
     1131                                case SSAApproximationEnum: GradientJDragSSA(element,gradient,control_interp,control_index); break;
     1132                                case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_interp,control_index); break;
     1133                                case HOApproximationEnum:  GradientJDragHO( element,gradient,control_interp,control_index); break;
     1134                                case FSApproximationEnum:  GradientJDragFS( element,gradient,control_interp,control_index); break;
    11351135                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11361136                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11391139                case FrictionAsEnum:
    11401140                        switch(approximation){
    1141                                 case SSAApproximationEnum: GradientJDragHydroSSA(element,gradient,control_index); break;
    1142                                 case L1L2ApproximationEnum:GradientJDragHydroL1L2(element,gradient,control_index); break;
    1143                                 case HOApproximationEnum:  GradientJDragHydroHO( element,gradient,control_index); break;
    1144                                 case FSApproximationEnum:  GradientJDragHydroFS( element,gradient,control_index); break;
     1141                                case SSAApproximationEnum: GradientJDragHydroSSA(element,gradient,control_interp,control_index); break;
     1142                                case L1L2ApproximationEnum:GradientJDragHydroL1L2(element,gradient,control_interp,control_index); break;
     1143                                case HOApproximationEnum:  GradientJDragHydroHO( element,gradient,control_interp,control_index); break;
     1144                                case FSApproximationEnum:  GradientJDragHydroFS( element,gradient,control_interp,control_index); break;
    11451145                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11461146                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11491149                case MaterialsRheologyBbarEnum:
    11501150                        switch(approximation){
    1151                                 case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_index); break;
    1152                                 case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_index); break;
    1153                                 case HOApproximationEnum:  GradientJBbarHO( element,gradient,control_index); break;
    1154                                 case FSApproximationEnum:  GradientJBbarFS( element,gradient,control_index); break;
     1151                                case SSAApproximationEnum: GradientJBbarSSA(element,gradient,control_interp,control_index); break;
     1152                                case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_interp,control_index); break;
     1153                                case HOApproximationEnum:  GradientJBbarHO( element,gradient,control_interp,control_index); break;
     1154                                case FSApproximationEnum:  GradientJBbarFS( element,gradient,control_interp,control_index); break;
    11551155                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11561156                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11591159                case MaterialsRheologyBEnum:
    11601160                        switch(approximation){
    1161                                 case SSAApproximationEnum: GradientJBSSA(element,gradient,control_index); break;
    1162                                 case HOApproximationEnum:  GradientJBHO( element,gradient,control_index); break;
    1163                                 case FSApproximationEnum:  GradientJBFS( element,gradient,control_index); break;
     1161                                case SSAApproximationEnum: GradientJBSSA(element,gradient,control_interp,control_index); break;
     1162                                case HOApproximationEnum:  GradientJBHO( element,gradient,control_interp,control_index); break;
     1163                                case FSApproximationEnum:  GradientJBFS( element,gradient,control_interp,control_index); break;
    11641164                                case NoneApproximationEnum: /*Gradient is 0*/                    break;
    11651165                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11681168                case DamageDbarEnum:
    11691169                        switch(approximation){
    1170                                 case SSAApproximationEnum: GradientJDSSA(element,gradient,control_index); break;
     1170                                case SSAApproximationEnum: GradientJDSSA(element,gradient,control_interp,control_index); break;
    11711171                                case NoneApproximationEnum: /*Gradient is 0*/                 break;
    11721172                                default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    11801180
    11811181}/*}}}*/
    1182 void           AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1182void           AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    11831183        /*WARNING: We use SSA as an estimate for now*/
    1184         this->GradientJBbarSSA(element,gradient,control_index);
    1185 }/*}}}*/
    1186 void           AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1184        this->GradientJBbarSSA(element,gradient,control_interp,control_index);
     1185}/*}}}*/
     1186void           AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    11871187
    11881188        /*Intermediaries*/
    … …  
    12591259        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    12601260}/*}}}*/
    1261 void           AdjointHorizAnalysis::GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1261void           AdjointHorizAnalysis::GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    12621262
    12631263           /*Same as SSA*/
    1264            return this->GradientJBbarSSA(element,gradient,control_index);
    1265 }/*}}}*/
    1266 void           AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1264           return this->GradientJBbarSSA(element,gradient,control_interp,control_index);
     1265}/*}}}*/
     1266void           AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    12671267
    12681268        /*WARNING: We use SSA as an estimate for now*/
    1269         this->GradientJBbarSSA(element,gradient,control_index);
    1270 }/*}}}*/
    1271 void           AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1269        this->GradientJBbarSSA(element,gradient,control_interp,control_index);
     1270}/*}}}*/
     1271void           AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    12721272
    12731273        /*Intermediaries*/
    … …  
    13361336
    13371337                /*Build gradient vector (actually -dJ/dB): */
    1338                 for(int i=0;i<numvertices;i++){
    1339                         ge[i]+=-dmudB*thickness*(
     1338                if(control_interp==P1Enum){
     1339                        for(int i=0;i<numvertices;i++){
     1340                                ge[i]+=-dmudB*thickness*(
     1341                                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     1342                                                        )*Jdet*gauss->weight*basis[i];
     1343                                _assert_(!xIsNan<IssmDouble>(ge[i]));
     1344                        }
     1345                }
     1346                else if(control_interp==P0Enum){
     1347                        ge[0]+=-dmudB*thickness*(
    13401348                                                (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
    1341                                                 )*Jdet*gauss->weight*basis[i];
    1342                         _assert_(!xIsNan<IssmDouble>(ge[i]));
    1343                 }
    1344         }
    1345         gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1349                                                )*Jdet*gauss->weight;
     1350                        _assert_(!xIsNan<IssmDouble>(ge[0]));
     1351                }
     1352                else{
     1353                        _error_("not supported");
     1354                }
     1355        }
     1356        if(control_interp==P1Enum){
     1357                gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     1358        }
     1359        else if(control_interp==P0Enum){
     1360                gradient->SetValue(vertexpidlist[0],ge[0],ADD_VAL);
     1361        }
     1362        else{
     1363                _error_("not supported");
     1364        }
    13461365
    13471366        /*Clean up and return*/
    … …  
    13531372        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    13541373}/*}}}*/
    1355 void           AdjointHorizAnalysis::GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1374void           AdjointHorizAnalysis::GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    13561375        /*WARNING: We use HO as an estimate for now*/
    1357         this->GradientJBHO(element,gradient,control_index);
    1358 }/*}}}*/
    1359 void           AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1376        this->GradientJBHO(element,gradient,control_interp,control_index);
     1377}/*}}}*/
     1378void           AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    13601379
    13611380        /*Intermediaries*/
    … …  
    14231442        delete gauss;
    14241443}/*}}}*/
    1425 void           AdjointHorizAnalysis::GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1444void           AdjointHorizAnalysis::GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    14261445
    14271446        /*Intermediaries*/
    … …  
    14871506        delete gauss;
    14881507}/*}}}*/
    1489 void           AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1508void           AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    14901509        /*Intermediaries*/
    14911510        int      domaintype,dim;
    … …  
    15631582        delete gauss;
    15641583}/*}}}*/
    1565 void           AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1584void           AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    15661585
    15671586        /*Intermediaries*/
    … …  
    16471666        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    16481667}/*}}}*/
    1649 void           AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1668void           AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    16501669
    16511670        /*return if floating (gradient is 0)*/
    … …  
    17411760
    17421761}/*}}}*/
    1743 void           AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1762void           AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    17441763
    17451764        /*return if floating or not on bed (gradient is 0)*/
    … …  
    18371856        delete friction;
    18381857}/*}}}*/
    1839 void           AdjointHorizAnalysis::GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1858void           AdjointHorizAnalysis::GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    18401859
    18411860        /*Same as SSA*/
    1842         return this->GradientJDragSSA(element,gradient,control_index);
    1843 }/*}}}*/
    1844 void           AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1861        return this->GradientJDragSSA(element,gradient,control_interp,control_index);
     1862}/*}}}*/
     1863void           AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    18451864
    18461865        /*return if floating or not on bed (gradient is 0)*/
    … …  
    19161935        delete friction;
    19171936}/*}}}*/
    1918 void           AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     1937void           AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    19191938
    19201939        /*return if floating (gradient is 0)*/
    … …  
    20182037        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    20192038}/*}}}*/
    2020 
    2021 void AdjointHorizAnalysis::GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2039void           AdjointHorizAnalysis::GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    20222040
    20232041        /*return if floating or not on bed (gradient is 0)*/
    … …  
    21132131        delete friction;
    21142132}/*}}}*/
    2115 void           AdjointHorizAnalysis::GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2133void           AdjointHorizAnalysis::GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    21162134
    21172135        /*Same as SSA*/
    2118         return this->GradientJDragSSA(element,gradient,control_index);
    2119 }/*}}}*/
    2120 void           AdjointHorizAnalysis::GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2136        return this->GradientJDragSSA(element,gradient,control_interp,control_index);
     2137}/*}}}*/
     2138void           AdjointHorizAnalysis::GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    21212139
    21222140        /*return if floating or not on bed (gradient is 0)*/
    … …  
    21902208        delete friction;
    21912209}/*}}}*/
    2192 
    2193 void AdjointHorizAnalysis::GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2210void           AdjointHorizAnalysis::GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    21942211
    21952212        /*return if floating (gradient is 0)*/
    … …  
    23162333        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    23172334}/*}}}*/
    2318 
    2319 void           AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     2335void           AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/
    23202336
    23212337        /*Intermediaries*/

  • issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h

    r24335 r25317  
    3535                ElementVector* CreatePVectorSSA(Element* element);
    3636                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    37                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    38                 void           GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
    39                 void           GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_index);
    40                 void           GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_index);
    41                 void           GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
    42                 void           GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
    43                 void           GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
    44                 void           GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
    45                 void           GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index);
    46                 void           GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
    47                 void           GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
    48                 void           GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
    49                 void           GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index);
    50                 void           GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
    51                 void           GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
    52                 void           GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
    53                 void           GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_index);
    54                 void           GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_index);
    55                 void           GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_index);
    56                 void           GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
    57                 void           GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index);
     37                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
     38                void           GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     39                void           GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     40                void           GradientJBinitial(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     41                void           GradientJBbarL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     42                void           GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     43                void           GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     44                void           GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     45                void           GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     46                void           GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     47                void           GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     48                void           GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     49                void           GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     50                void           GradientJDragL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     51                void           GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     52                void           GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     53                void           GradientJDragHydroFS(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     54                void           GradientJDragHydroL1L2(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     55                void           GradientJDragHydroHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     56                void           GradientJDragHydroSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
     57                void           GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index);
    5858                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    5959                void           InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);

  • issm/trunk-jpl/src/c/analyses/Analysis.h

    r24335 r25317  
    4949                virtual ElementVector* CreatePVector(Element* element)=0;
    5050                virtual void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element)=0;
    51                 virtual void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index)=0;
     51                virtual void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index)=0;
    5252                virtual void           InputUpdateFromSolution(IssmDouble* solution,Element* element)=0;
    5353                virtual void           UpdateConstraints(FemModel* femmodel)=0;

  • issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp

    r25118 r25317  
    174174                element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
    175175}/*}}}*/
    176 void           Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     176void           Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    177177        _error_("Not implemented yet");
    178178}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp

    r25262 r25317  
    424424           _error_("not implemented yet");
    425425}/*}}}*/
    426 void           BalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     426void           BalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    427427        /* WARNING: this gradient is valid for Soft balance thickness only */
    428428

  • issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h

    r25262 r25317  
    3131                ElementVector* CreatePVectorDG(Element* element);
    3232                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    33                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     33                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3434                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3535                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.cpp

    r24335 r25317  
    4545           _error_("not implemented yet");
    4646}/*}}}*/
    47 void BalancethicknessSoftAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     47void BalancethicknessSoftAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    4848        _error_("Not implemented yet");
    4949}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/BalancethicknessSoftAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp

    r25118 r25317  
    239239           _error_("not implemented yet");
    240240}/*}}}*/
    241 void           BalancevelocityAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     241void           BalancevelocityAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    242242        _error_("Not implemented yet");
    243243}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

    r25263 r25317  
    677677        element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum);
    678678}/*}}}*/
    679 void           DamageEvolutionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     679void           DamageEvolutionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    680680        _error_("Not implemented yet");
    681681}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h

    r25263 r25317  
    3131                ElementVector* CreatePVector(Element* element);
    3232                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    33                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     33                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3434                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3535                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp

    r25265 r25317  
    140140           _error_("not implemented yet");
    141141}/*}}}*/
    142 void           DepthAverageAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     142void           DepthAverageAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    143143        _error_("Not implemented yet");
    144144}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h

    r25265 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r25118 r25317  
    15211521        return (1.-waterfraction)*kappa_i + waterfraction*kappa_w;
    15221522}/*}}}*/
    1523 void           EnthalpyAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     1523void           EnthalpyAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    15241524        _error_("Not implemented yet");
    15251525}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h

    r24939 r25317  
    4848                static int        GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate);
    4949                static IssmDouble GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure);
    50                 void              GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     50                void              GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    5151                void              InputUpdateFromSolution(IssmDouble* solution,Element* element);
    5252                static void       PostProcessing(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/EnumToAnalysis.h

    r24335 r25317  
    2222                ElementVector* CreatePVector(Element* element);
    2323                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    24                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     24                void GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    2525                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    2626                void UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp

    r25146 r25317  
    205205           _error_("not implemented yet");
    206206}/*}}}*/
    207 void           EsaAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     207void           EsaAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    208208        _error_("Not implemented yet");
    209209}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/EsaAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp

    r24433 r25317  
    211211        _error_("not implemented yet");
    212212}/*}}}*/
    213 void           ExtrapolationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     213void           ExtrapolationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    214214        _error_("Not implemented yet");
    215215}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h

    r24431 r25317  
    2727        ElementVector* CreatePVector(Element* element);
    2828        void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29         void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29        void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030        void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131        int                             GetExtrapolationCase(Element* element);

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp

    r24335 r25317  
    9595           _error_("not implemented yet");
    9696}/*}}}*/
    97 void           ExtrudeFromBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     97void           ExtrudeFromBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    9898        _error_("Not implemented yet");
    9999}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp

    r24335 r25317  
    9595           _error_("not implemented yet");
    9696}/*}}}*/
    97 void           ExtrudeFromTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     97void           ExtrudeFromTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    9898        _error_("Not implemented yet");
    9999}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp

    r25264 r25317  
    348348           _error_("not implemented yet");
    349349}/*}}}*/
    350 void           FreeSurfaceBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     350void           FreeSurfaceBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    351351        _error_("Not implemented yet");
    352352}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h

    r25264 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp

    r25264 r25317  
    322322           _error_("not implemented yet");
    323323}/*}}}*/
    324 void           FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     324void           FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    325325        _error_("Not implemented yet");
    326326}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h

    r25264 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp

    r25262 r25317  
    186186        _error_("not implemented yet");
    187187}/*}}}*/
    188 void           GLheightadvectionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     188void           GLheightadvectionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    189189        _error_("Not implemented yet");
    190190}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.h

    r25262 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/GiaAnalysis.cpp

    r25080 r25317  
    106106           _error_("not implemented yet");
    107107}/*}}}*/
    108 void           GiaAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     108void           GiaAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    109109        _error_("Not implemented yet");
    110110}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/GiaAnalysis.h

    r25066 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

    r25308 r25317  
    384384        element->GetSolutionFromInputsOneDof(solution,EplHeadSubstepEnum);
    385385}/*}}}*/
    386 void HydrologyDCEfficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     386void HydrologyDCEfficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    387387        _error_("Not implemented yet");
    388388}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h

    r25308 r25317  
    3030                ElementVector* CreatePVector(Element* element);
    3131                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    32                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     32                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index);
    3333                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3434                void UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

    r25308 r25317  
    450450        element->GetSolutionFromInputsOneDof(solution,SedimentHeadSubstepEnum);
    451451}/*}}}*/
    452 void HydrologyDCInefficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     452void HydrologyDCInefficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    453453        _error_("Not implemented yet");
    454454}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h

    r25227 r25317  
    2828                ElementVector* CreatePVector(Element* element);
    2929                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    30                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index);
    3131                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

    r24861 r25317  
    430430
    431431}/*}}}*/
    432 void           HydrologyGlaDSAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     432void           HydrologyGlaDSAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    433433        _error_("Not implemented yet");
    434434}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp

    r24861 r25317  
    7979        _error_("not implemented");
    8080}/*}}}*/
    81 void           HydrologyPismAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     81void           HydrologyPismAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    8282        _error_("Not implemented yet");
    8383}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp

    r24861 r25317  
    343343        element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum);
    344344}/*}}}*/
    345 void           HydrologyShaktiAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     345void           HydrologyShaktiAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    346346        _error_("Not implemented yet");
    347347}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

    r25266 r25317  
    281281        element->GetSolutionFromInputsOneDof(solution,WatercolumnEnum);
    282282}/*}}}*/
    283 void           HydrologyShreveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     283void           HydrologyShreveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    284284        _error_("Not implemented yet");
    285285}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h

    r25266 r25317  
    2828                ElementVector* CreatePVector(Element* element);
    2929                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    30                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3131                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp

    r25263 r25317  
    198198           _error_("not implemented yet");
    199199}/*}}}*/
    200 void           L2ProjectionBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     200void           L2ProjectionBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    201201        _error_("Not implemented yet");
    202202}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp

    r25228 r25317  
    223223           _error_("not implemented yet");
    224224}/*}}}*/
    225 void           L2ProjectionEPLAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     225void           L2ProjectionEPLAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    226226        _error_("Not implemented yet");
    227227}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r24861 r25317  
    621621        _error_("not implemented yet");
    622622}/*}}}*/
    623 void           LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     623void           LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    624624        _error_("Not implemented yet");
    625625}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h

    r24429 r25317  
    2828                IssmDouble     GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
    2929                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    30                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3131                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp

    r24335 r25317  
    6767           _error_("not implemented yet");
    6868}/*}}}*/
    69 void           LoveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     69void           LoveAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    7070        _error_("Not implemented yet");
    7171}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/LoveAnalysis.h

    r24335 r25317  
    2929                void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    3030                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    31                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     31                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3232                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3333                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r25266 r25317  
    790790        element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
    791791}/*}}}*/
    792 void           MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     792void           MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    793793        _error_("Not implemented yet");
    794794}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h

    r24730 r25317  
    3131                ElementVector* CreatePVectorDG(Element* element);
    3232                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    33                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     33                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3434                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3535                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp

    r25261 r25317  
    125125           _error_("not implemented yet");
    126126}/*}}}*/
    127 void           MeltingAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     127void           MeltingAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    128128        _error_("Not implemented yet");
    129129}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/MeltingAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

    r25275 r25317  
    413413           _error_("not implemented yet");
    414414}/*}}}*/
    415 void           SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     415void           SealevelriseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    416416        _error_("Not implemented yet");
    417417}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r25209 r25317  
    507507           _error_("not implemented yet");
    508508}/*}}}*/
    509 void           SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     509void           SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    510510        _error_("Not implemented yet");
    511511}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp

    r24335 r25317  
    206206           _error_("not implemented yet");
    207207}/*}}}*/
    208 void           SmoothAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     208void           SmoothAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    209209        _error_("Not implemented yet");
    210210}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/SmoothAnalysis.h

    r24335 r25317  
    2727                ElementVector* CreatePVector(Element* element);
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     29                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3030                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3131                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r25239 r25317  
    12371237        xDelete<int>(doflist);
    12381238}/*}}}*/
    1239 void           StressbalanceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     1239void           StressbalanceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    12401240        _error_("Not implemented yet");
    12411241}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r24335 r25317  
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2929                void           GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
    30                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3131                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp

    r24933 r25317  
    556556        xDelete<IssmDouble>(values);
    557557}/*}}}*/
    558 void           StressbalanceSIAAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     558void           StressbalanceSIAAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    559559        _error_("Not implemented yet");
    560560}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h

    r24335 r25317  
    3131                ElementVector* CreatePVector3D(Element* element);
    3232                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    33                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     33                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3434                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3535                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

    r25118 r25317  
    505505        element->GetSolutionFromInputsOneDof(solution,VzEnum);
    506506}/*}}}*/
    507 void           StressbalanceVerticalAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     507void           StressbalanceVerticalAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    508508        _error_("Not implemented yet");
    509509}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h

    r24335 r25317  
    3333                ElementVector* CreatePVectorVolume(Element* element);
    3434                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    35                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     35                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3636                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3737                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

    r25265 r25317  
    742742        element->GetSolutionFromInputsOneDof(solution,TemperatureEnum);
    743743}/*}}}*/
    744 void           ThermalAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     744void           ThermalAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    745745        _error_("Not implemented yet");
    746746}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h

    r25261 r25317  
    3232                ElementVector* CreatePVectorVolume(Element* element);
    3333                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    34                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     34                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3535                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3636                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp

    r24429 r25317  
    182182        _error_("not implemented yet");
    183183}/*}}}*/
    184 void           UzawaPressureAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     184void           UzawaPressureAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    185185        _error_("Not implemented yet");
    186186}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.h

    r24335 r25317  
    2828                void           GetM(IssmDouble* M,Element* element,Gauss* gauss);
    2929                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    30                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void           GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index);
    3131                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void           UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r25242 r25317  
    15091509        int num_controls;
    15101510        bool isautodiff;
     1511        int* N=NULL;
     1512        int* M=NULL;
     1513        int* interp=NULL;
    15111514        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    1512         parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
     1515        parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     1516        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     1517        parameters->FindParam(&interp,NULL,ControlInputInterpolationEnum);
    15131518
    15141519        /*Get number of vertices*/
    15151520        const int NUM_VERTICES = this->GetNumberOfVertices();
    1516         if(isautodiff){
    1517                 int* N=NULL;
    1518                 int* M=NULL;
    1519                 int start = 0;
    1520                 parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    1521                 parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    1522                 if(control_index>0) {
    1523                         for(int n=0;n<control_index;n++){
    1524                                 start+=N[n]*M[n];
    1525                         }
    1526                 }
    1527 
     1521
     1522        /*Get starting index of gradient for this control*/
     1523        int start = 0;
     1524        for(int n=0;n<control_index;n++) start+=N[n]*M[n];
     1525
     1526        /*Create index*/
     1527        if(interp[control_index]==P1Enum){
    15281528                for(int n=0;n<N[control_index];n++){
    15291529                        for(int i=0;i<NUM_VERTICES;i++){
    1530                                 indexing[i+n*NUM_VERTICES]=this->vertices[i]->Sid() + start + n*M[control_index];
     1530                                indexing[i+n*NUM_VERTICES]= start + n*M[control_index] + this->vertices[i]->Sid();
    15311531                        }
    15321532                }
    15331533        }
     1534        else if(interp[control_index]==P0Enum){
     1535                for(int n=0;n<N[control_index];n++){
     1536                                indexing[n]= start + n*M[control_index] + this->Sid();
     1537                }
     1538        }
    15341539        else{
    1535                 int M;
    1536                 parameters->FindParam(&M,ControlInputSizeMEnum);
    1537                 /*get gradient indices*/
    1538                 for(int i=0;i<NUM_VERTICES;i++){
    1539                         indexing[i]=this->vertices[i]->Sid() + (control_index)*M;
    1540                 }
    1541         }
     1540                _error_("interpolation not supported");
     1541        }
     1542
     1543        xDelete<int>(M);
     1544        xDelete<int>(N);
     1545        xDelete<int>(interp);
    15421546}
    15431547/*}}}*/
    … …  
    17761780                }
    17771781                this->AddControlInput(input_enum,inputs2,iomodel,values,values_min,values_max,P1Enum,id);
     1782        }
     1783        else if(M==iomodel->numberofelements && N==1){
     1784                values[0]     = vector[this->Sid()];
     1785                values_min[0] = scale*min_vector[this->Sid()];
     1786                values_max[0] = scale*max_vector[this->Sid()];
     1787                this->AddControlInput(input_enum,inputs2,iomodel,values,values_min,values_max,P0Enum,id);
    17781788        }
    17791789

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r25256 r25317  
    231231                virtual void       ComputeEsaStrainAndVorticity(void)=0;
    232232                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs2* inputs2in)=0;
    233                 virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M)=0;
     233                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp)=0;
    234234                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
    235                 virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
     235                virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp)=0;
    236236                virtual void       CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){_error_("not implemented yet");};
    237237                virtual void       CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet "<<this->ObjectEnum());};
    … …  
    335335                virtual void       SetElementInput(Inputs2* inputs2,int enum_in,IssmDouble values){_error_("not implemented yet");};
    336336                virtual void       SetElementInput(Inputs2* inputs2,int numindices,int* indices,IssmDouble* values,int enum_in){_error_("not implemented yet");};
    337                 virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M)=0;
     337                virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int M,int N)=0;
    338338                virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0;
    339339                virtual void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r25256 r25317  
    899899}
    900900/*}}}*/
    901 void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){/*{{{*/
    902 
    903         if(enum_type==MaterialsRheologyBbarEnum) enum_type = MaterialsRheologyBEnum;
    904         if(enum_type==DamageDbarEnum)            enum_type = DamageDEnum;
    905 
    906         _error_("not implemented");
    907         int    vertexpidlist[NUMVERTICES];
    908         IssmDouble grad_list[NUMVERTICES];
    909         Input2* grad_input=NULL;
    910         Input2* input=NULL;
    911 
    912         if(enum_type==MaterialsRheologyBbarEnum){
    913                 input=this->GetInput2(MaterialsRheologyBEnum);
    914         }
    915         else if(enum_type==DamageDbarEnum){
    916                 input=this->GetInput2(DamageDEnum);
    917         }
    918         else{
    919                 input=this->GetInput2(enum_type);
    920         }
    921         if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
    922         if(input->ObjectEnum()!=ControlInput2Enum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
    923 
    924         GradientIndexing(&vertexpidlist[0],control_index);
    925 
    926         //for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
    927         //grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
    928         //((ControlInput*)input)->SetGradient(grad_input);
    929         _error_("not implemented");
     901void       Penta::ControlInputSetGradient(IssmDouble* gradient,int control_enum,int control_index,int offset,int M,int N,int interp){/*{{{*/
     902
     903        IssmDouble  values[NUMVERTICES];
     904        int         lidlist[NUMVERTICES];
     905        int         idlist[NUMVERTICES];
     906
     907        if(control_enum==MaterialsRheologyBbarEnum) control_enum = MaterialsRheologyBEnum;
     908        if(control_enum==DamageDbarEnum)            control_enum = DamageDEnum;
     909
     910        ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"gradient");   _assert_(input);
     911        this->GetVerticesLidList(&lidlist[0]);
     912        GradientIndexing(&idlist[0],control_index);
     913
     914        /*Get values on vertices*/
     915        if(input->ObjectEnum()==PentaInput2Enum && input->GetInputInterpolationType()==P1Enum){
     916                _assert_(N==1);
     917                for(int i=0;i<NUMVERTICES;i++){
     918                        values[i] = gradient[idlist[i]];
     919                }
     920                input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
     921        }
     922        else if(input->ObjectEnum()==PentaInput2Enum && input->GetInputInterpolationType()==P0Enum){
     923                _assert_(N==1);
     924                input->SetInput(P0Enum,this->lid,gradient[idlist[0]]);
     925        }
     926        else if(input->ObjectEnum()==TransientInput2Enum){
     927                for(int n=0;n<N;n++){
     928                        _error_("not implemented");
     929                        //Input* new_input = new PentaInput(control_enum,gradient,P1Enum);
     930                        //controlinput->SetInput(new_input,n);
     931                        //controlinput->Configure(parameters);
     932                }
     933        }
     934        else _error_("Type not supported");
    930935
    931936}
    … …  
    946951        this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0]);
    947952}/*}}}*/
    948 void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
     953void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp){/*{{{*/
    949954
    950955        int         sidlist[NUMVERTICES];
    … …  
    955960        IssmDouble  value,gradient;
    956961
     962        /*Get relevant inputs*/
    957963        if(control_enum==MaterialsRheologyBbarEnum) control_enum = MaterialsRheologyBEnum;
    958964        if(control_enum==DamageDbarEnum)            control_enum = DamageDEnum;
    959 
    960         this->GetVerticesConnectivityList(&connectivity[0]);
    961         this->GetVerticesSidList(&sidlist[0]);
    962         this->GetVerticesLidList(&lidlist[0]);
    963 
    964965        ElementInput2* control_value    = this->inputs2->GetControlInput2Data(control_enum,"value");    _assert_(control_value);
    965966        ElementInput2* control_gradient = this->inputs2->GetControlInput2Data(control_enum,"gradient"); _assert_(control_gradient);
    966         control_value->Serve(NUMVERTICES,&lidlist[0]);
    967         control_gradient->Serve(NUMVERTICES,&lidlist[0]);
    968 
    969         GaussPenta* gauss=new GaussPenta();
    970         for (int iv=0;iv<NUMVERTICES;iv++){
    971                 gauss->GaussVertex(iv);
    972 
    973                 control_value->GetInputValue(&value,gauss);
    974                 control_gradient->GetInputValue(&gradient,gauss);
    975 
    976                 values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
    977                 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
    978         }
    979         delete gauss;
    980 
    981         vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
    982         vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
     967
     968        if(control_interp==P1Enum){
     969                _assert_(control_value->GetInputInterpolationType()==P1Enum);
     970                _assert_(control_gradient->GetInputInterpolationType()==P1Enum);
     971
     972                this->GetVerticesConnectivityList(&connectivity[0]);
     973                this->GetVerticesSidList(&sidlist[0]);
     974                this->GetVerticesLidList(&lidlist[0]);
     975
     976                control_value->Serve(NUMVERTICES,&lidlist[0]);
     977                control_gradient->Serve(NUMVERTICES,&lidlist[0]);
     978
     979                GaussPenta* gauss=new GaussPenta();
     980                for (int iv=0;iv<NUMVERTICES;iv++){
     981                        gauss->GaussVertex(iv);
     982
     983                        control_value->GetInputValue(&value,gauss);
     984                        control_gradient->GetInputValue(&gradient,gauss);
     985
     986                        values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
     987                        gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
     988                }
     989                delete gauss;
     990
     991                vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
     992                vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
     993        }
     994        else if(control_interp==P0Enum){
     995                _assert_(control_value->GetInputInterpolationType()==P0Enum);
     996                _assert_(control_gradient->GetInputInterpolationType()==P0Enum);
     997
     998                control_value->Serve(1,&this->lid);
     999                control_gradient->Serve(1,&this->lid);
     1000
     1001                vector_control->SetValue(this->sid,control_value->element_values[0],ADD_VAL);
     1002                vector_gradient->SetValue(this->sid,control_gradient->element_values[0],ADD_VAL);
     1003        }
     1004        else{
     1005                _error_("not supported");
     1006        }
     1007
    9831008}/*}}}*/
    9841009void       Penta::CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){/*{{{*/
    … …  
    18221847void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    18231848
    1824         /*Get out if this is not an element input*/
    1825         if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
    1826 
    1827         /*Prepare index list*/
    1828         int idlist[NUMVERTICES];
    1829         GradientIndexing(&idlist[0],control_index);
    1830 
    1831         /*Get input (either in element or material)*/
     1849        _error_("NOT NEEDED ANYMORE");
     1850        /*Get input*/
    18321851        if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
    18331852        ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);
    18341853
    1835         /*Intermediaries*/
    1836         int numindices;
    1837         int indices[NUMVERTICES];
    1838 
    1839         /*Check interpolation*/
    1840         int interpolation = input->GetInterpolation();
    1841         switch(interpolation){
    1842                 case P1Enum:
    1843                         numindices = NUMVERTICES;
    1844                         for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
    1845                         input->Serve(numindices,&indices[0]);
    1846                         break;
    1847                 default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
    1848         }
    1849 
    1850         /* Start looping on the number of vertices: */
    1851         IssmDouble values[NUMVERTICES];
    1852         Gauss*gauss=this->NewGauss();
    1853         for(int iv=0;iv<NUMVERTICES;iv++){
    1854                 gauss->GaussVertex(iv);
    1855                 input->GetInputValue(&values[iv],gauss);
    1856         }
    1857         delete gauss;
    1858 
    1859         vector->SetValues(NUMVERTICES,idlist,&values[0],INS_VAL);
    1860 }
    1861 /*}}}*/
    1862 void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset){/*{{{*/
    1863 
    1864         int* idlist = NULL;
    1865         IssmDouble* values = NULL;
    1866 
    1867         /*Get input*/
    1868         ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);
     1854        /*Lid list once for all*/
     1855        int lidlist[NUMVERTICES];
     1856        for(int i=0;i<NUMVERTICES;i++) lidlist[i] = vertices[i]->lid;
    18691857
    18701858        /*Check what input we are dealing with*/
    18711859        switch(input->ObjectEnum()){
    18721860                case PentaInput2Enum:
    1873                                   {
    1874                                         PentaInput2* pentainput = xDynamicCast<PentaInput2*>(input);
    1875                                         if(pentainput->GetInputInterpolationType()!=P1Enum) _error_("not supported yet");
    1876 
    1877                                         /*Create list of indices and values for global vector*/
    1878                                         idlist = xNew<int>(NUMVERTICES);
    1879                                         values = xNew<IssmDouble>(NUMVERTICES);
    1880                                         GradientIndexing(&idlist[0],control_index);
    1881                                         for(int i=0;i<NUMVERTICES;i++){
    1882                                                 values[i] = pentainput->element_values[i];
     1861                          {
     1862                                IssmDouble values[NUMVERTICES];
     1863                                int        idlist[NUMVERTICES];
     1864
     1865                                PentaInput2* triainput = xDynamicCast<PentaInput2*>(input);
     1866
     1867                                /*Create list of indices and values for global vector*/
     1868                                GradientIndexing(&idlist[0],control_index);
     1869
     1870                                if(triainput->GetInputInterpolationType()==P1Enum){
     1871                                        input->Serve(NUMVERTICES,&lidlist[0]);
     1872                                        for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
     1873                                        vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     1874                                }
     1875                                else if(triainput->GetInputInterpolationType()==P0Enum){
     1876                                        input->Serve(1,&this->lid);
     1877                                        vector->SetValue(idlist[0],triainput->element_values[0],INS_VAL);
     1878                                }
     1879                                else{
     1880                                        _error_("not supported yet");
     1881                                }
     1882                                break;
     1883                          }
     1884
     1885                case TransientInputEnum:
     1886                                {
     1887                                        TransientInput2* transientinput = xDynamicCast<TransientInput2*>(input);
     1888                                        int  N = transientinput->numtimesteps;
     1889                                        int* M = NULL;
     1890                                        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     1891                                        int* idlist = xNew<int>(NUMVERTICES*N);
     1892                                        IssmDouble* values = xNew<IssmDouble>(NUMVERTICES*N);
     1893                                        for(int t=0;t<transientinput->numtimesteps;t++) {
     1894                                                IssmDouble time = transientinput->GetTimeByOffset(t);
     1895                                                _error_("not implemented");
     1896                                                //PentaInput* timeinput = xDynamicCast<PentaInput*>(transientinput->GetTimeInput(time));
     1897                                                //if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
     1898                                                //input->Serve(NUMVERTICES,&lidlist[0]);
     1899                                                ///*Create list of indices and values for global vector*/
     1900                                                //for(int i=0;i<NUMVERTICES;i++){
     1901                                                //              idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
     1902                                                //              values[N*i+t] = timeinput->values[i];
     1903                                                //}
    18831904                                        }
     1905                                        vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
     1906                                        xDelete<int>(M);
     1907                                        xDelete<int>(idlist);
     1908                                        xDelete<IssmDouble>(values);
     1909                                        break;
     1910                                }
     1911                default: _error_("input "<<EnumToStringx(input->ObjectEnum())<<" not supported yet");
     1912        }
     1913}
     1914/*}}}*/
     1915void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset){/*{{{*/
     1916
     1917        /*Get input*/
     1918        if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
     1919        ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);
     1920
     1921        /*Lid list once for all*/
     1922        int lidlist[NUMVERTICES];
     1923        for(int i=0;i<NUMVERTICES;i++) lidlist[i] = vertices[i]->lid;
     1924
     1925        /*Check what input we are dealing with*/
     1926        switch(input->ObjectEnum()){
     1927                case PentaInput2Enum:
     1928                          {
     1929                                IssmDouble values[NUMVERTICES];
     1930                                int        idlist[NUMVERTICES];
     1931
     1932                                PentaInput2* triainput = xDynamicCast<PentaInput2*>(input);
     1933
     1934                                /*Create list of indices and values for global vector*/
     1935                                GradientIndexing(&idlist[0],control_index);
     1936
     1937                                if(triainput->GetInputInterpolationType()==P1Enum){
     1938                                        input->Serve(NUMVERTICES,&lidlist[0]);
     1939                                        for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
    18841940                                        vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     1941                                }
     1942                                else if(triainput->GetInputInterpolationType()==P0Enum){
     1943                                        input->Serve(1,&this->lid);
     1944                                        vector->SetValue(idlist[0],triainput->element_values[0],INS_VAL);
     1945                                }
     1946                                else{
     1947                                        _error_("not supported yet");
     1948                                }
     1949                                break;
     1950                          }
     1951
     1952                case TransientInputEnum:
     1953                                {
     1954                                        TransientInput2* transientinput = xDynamicCast<TransientInput2*>(input);
     1955                                        int  N = transientinput->numtimesteps;
     1956                                        int* M = NULL;
     1957                                        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     1958                                        int* idlist = xNew<int>(NUMVERTICES*N);
     1959                                        IssmDouble* values = xNew<IssmDouble>(NUMVERTICES*N);
     1960                                        for(int t=0;t<transientinput->numtimesteps;t++) {
     1961                                                IssmDouble time = transientinput->GetTimeByOffset(t);
     1962                                                _error_("not implemented");
     1963                                                //PentaInput* timeinput = xDynamicCast<PentaInput*>(transientinput->GetTimeInput(time));
     1964                                                //if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
     1965                                                //input->Serve(NUMVERTICES,&lidlist[0]);
     1966                                                ///*Create list of indices and values for global vector*/
     1967                                                //for(int i=0;i<NUMVERTICES;i++){
     1968                                                //              idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
     1969                                                //              values[N*i+t] = timeinput->values[i];
     1970                                                //}
     1971                                        }
     1972                                        vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
     1973                                        xDelete<int>(M);
     1974                                        xDelete<int>(idlist);
     1975                                        xDelete<IssmDouble>(values);
    18851976                                        break;
    1886                                   }
    1887 
    1888                                 case TransientInputEnum:
    1889                                   {
    1890                                         _error_("not implemented (see Tria)");
    1891                                         break;
    1892                                   }
    1893                                 default: _error_("input "<<input->ObjectEnum()<<" not supported yet");
    1894 
    1895                 }
    1896                         xDelete<int>(idlist);
    1897                         xDelete<IssmDouble>(values);
    1898         }
    1899 /*}}}*/
     1977                                }
     1978                default: _error_("input "<<EnumToStringx(input->ObjectEnum())<<" not supported yet");
     1979        }
     1980}/*}}}*/
    19001981void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
    19011982
    … …  
    33033384}
    33043385/*}}}*/
    3305 void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N, int M){/*{{{*/
     3386void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int M, int N){/*{{{*/
    33063387
    33073388        IssmDouble  values[NUMVERTICES];
    3308         int         vertexpidlist[NUMVERTICES],control_init;
     3389        int         lidlist[NUMVERTICES];
     3390        int         idlist[NUMVERTICES],control_init;
    33093391
    33103392        /*Specific case for depth averaged quantities*/
    … …  
    33193401        }
    33203402
     3403        /*Get Domain type*/
     3404        int domaintype;
     3405        parameters->FindParam(&domaintype,DomainTypeEnum);
     3406
     3407        /*Specific case for depth averaged quantities*/
     3408        if(domaintype==Domain2DverticalEnum){
     3409                if(control_enum==MaterialsRheologyBbarEnum){
     3410                        control_enum=MaterialsRheologyBEnum;
     3411                        if(!IsOnBase()) return;
     3412                }
     3413                if(control_enum==DamageDbarEnum){
     3414                        control_enum=DamageDEnum;
     3415                        if(!IsOnBase()) return;
     3416                }
     3417        }
     3418
    33213419        /*Get out if this is not an element input*/
    33223420        if(!IsInputEnum(control_enum)) return;
    33233421
    33243422        /*Prepare index list*/
    3325         GradientIndexing(&vertexpidlist[0],control_index);
     3423        ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"value");   _assert_(input);
     3424        this->GetVerticesLidList(&lidlist[0]);
     3425        GradientIndexing(&idlist[0],control_index);
    33263426
    33273427        /*Get values on vertices*/
    3328         for(int i=0;i<NUMVERTICES;i++){
    3329                 values[i]=vector[vertexpidlist[i]];
    3330         }
    3331         _error_("not implemented");
    3332         //Input* new_input = new PentaInput(control_enum,values,P1Enum);
    3333         Input2* input=(Input2*)this->GetInput2(control_enum);   _assert_(input);
    3334         if(input->ObjectEnum()!=ControlInputEnum){
    3335                 _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    3336         }
    3337 
    3338         //((ControlInput*)input)->SetInput(new_input);
    3339 
     3428        if(input->ObjectEnum()==PentaInput2Enum && input->GetInputInterpolationType()==P1Enum){
     3429                _assert_(N==1);
     3430                for(int i=0;i<NUMVERTICES;i++){
     3431                        values[i] = vector[idlist[i]];
     3432                }
     3433                input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
     3434        }
     3435        else if(input->ObjectEnum()==PentaInput2Enum && input->GetInputInterpolationType()==P0Enum){
     3436                _assert_(N==1);
     3437                input->SetInput(P0Enum,this->lid,vector[idlist[0]]);
     3438        }
     3439        else if(input->ObjectEnum()==TransientInput2Enum){
     3440                for(int n=0;n<N;n++){
     3441                        _error_("not implemented");
     3442                        //Input* new_input = new PentaInput(control_enum,values,P1Enum);
     3443                        //controlinput->SetInput(new_input,n);
     3444                        //controlinput->Configure(parameters);
     3445                }
     3446        }
     3447        else _error_("Type not supported");
     3448
     3449        /*Extrude depending on the control*/
    33403450        if(control_init==MaterialsRheologyBbarEnum){
    3341                 this->InputExtrude(control_enum,-1);
     3451                this->ControlInputExtrude(control_enum,-1);
    33423452        }
    33433453        if(control_init==DamageDbarEnum){
    3344                 this->InputExtrude(control_enum,-1);
     3454                this->ControlInputExtrude(control_enum,-1);
    33453455        }
    33463456}
    33473457/*}}}*/
    33483458void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
     3459
     3460        _error_("not needed anymore?");
    33493461
    33503462        IssmDouble  values[NUMVERTICES];

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r25256 r25317  
    6363                void           ComputeStressTensor();
    6464                void           Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs2* inputs2in);
    65                 void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M);
     65                void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp);
    6666                void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    67                 void           ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
     67                void           ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp);
    6868                void                            CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
    6969                ElementMatrix* CreateBasalMassMatrix(void);
    … …  
    167167                void           SetElementInput(Inputs2* inputs2,int enum_in,IssmDouble values);
    168168                void           SetElementInput(Inputs2* inputs2,int numindices,int* indices,IssmDouble* values,int enum_in);
    169                 void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M);
     169                void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int M,int N);
    170170                void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
    171171                void           SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r25256 r25317  
    5050                void        ComputeStressTensor(){_error_("not implemented yet");};
    5151                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs2* inputs2in){_error_("not implemented yet");};
    52                 void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");};
     52                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp){_error_("not implemented yet");};
    5353                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    54                 void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
     54                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp){_error_("not implemented yet");};
    5555                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    5656                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    … …  
    134134                void        ResetFSBasalBoundaryCondition(void){_error_("not implemented yet");};
    135135                void        ResetHooks(){_error_("not implemented yet");};
    136                 void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M){_error_("not implemented yet");};
     136                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int M,int N){_error_("not implemented yet");};
    137137                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");};
    138138                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r25256 r25317  
    4848                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};
    4949                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs2* inputs2in);
    50                 void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");};
     50                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp){_error_("not implemented yet");};
    5151                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    52                 void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
     52                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp){_error_("not implemented yet");};
    5353                IssmDouble  DragCoefficientAbsGradient(void){_error_("not implemented yet");};
    5454                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    … …  
    140140                void        ResetHooks();
    141141                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
    142                 void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){_error_("not implemented yet");};
     142                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int M, int N){_error_("not implemented yet");};
    143143                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");};
    144144                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r25293 r25317  
    222222        /*Call inputs method*/
    223223        switch(interpolation_enum){
     224                case P0Enum:
     225                        inputs2->SetTriaControlInput(input_enum,TriaInput2Enum,interpolation_enum,id,1,&this->lid,values,values_min,values_max);
     226                        break;
    224227                case P1Enum:
    225228                        inputs2->SetTriaControlInput(input_enum,TriaInput2Enum,interpolation_enum,id,NUMVERTICES,vertexlids,values,values_min,values_max);
    … …  
    10971100        this->inputs2=inputs2in;
    10981101}/*}}}*/
    1099 void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M){/*{{{*/
    1100 
     1102void       Tria::ControlInputSetGradient(IssmDouble* gradient,int control_enum,int control_index,int offset,int M,int N,int interp){/*{{{*/
     1103
     1104        IssmDouble  values[NUMVERTICES];
     1105        int         lidlist[NUMVERTICES];
    11011106        int         idlist[NUMVERTICES];
    1102         int         vertexlids[NUMVERTICES];
    1103         int         gradidlist[NUMVERTICES];
    1104         IssmDouble  grad_list[NUMVERTICES];
    1105 
    1106         GradientIndexing(&gradidlist[0],control_index);
    1107         for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;
    1108 
    1109         if(N==1){
    1110                 this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0]);
    1111         }
    1112         else{
     1107
     1108        ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"gradient");   _assert_(input);
     1109        this->GetVerticesLidList(&lidlist[0]);
     1110        GradientIndexing(&idlist[0],control_index);
     1111
     1112        /*Get values on vertices*/
     1113        if(input->ObjectEnum()==TriaInput2Enum && input->GetInputInterpolationType()==P1Enum){
     1114                _assert_(N==1);
     1115                for(int i=0;i<NUMVERTICES;i++){
     1116                        values[i] = gradient[idlist[i]];
     1117                }
     1118                input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
     1119        }
     1120        else if(input->ObjectEnum()==TriaInput2Enum && input->GetInputInterpolationType()==P0Enum){
     1121                _assert_(N==1);
     1122                input->SetInput(P0Enum,this->lid,gradient[idlist[0]]);
     1123        }
     1124        else if(input->ObjectEnum()==TransientInput2Enum){
    11131125                for(int n=0;n<N;n++){
    1114                         for(int i=0;i<NUMVERTICES;i++){
    1115                                 idlist[i] = offset + this->vertices[i]->Sid()+n*M;
    1116                                 grad_list[i]=gradient[idlist[i]];
    1117                         }
    1118                         this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0],n);
    1119                 }
    1120         }
     1126                        _error_("not implemented");
     1127                        //Input* new_input = new TriaInput(control_enum,gradient,P1Enum);
     1128                        //controlinput->SetInput(new_input,n);
     1129                        //controlinput->Configure(parameters);
     1130                }
     1131        }
     1132        else _error_("Type not supported");
    11211133
    11221134}/*}}}*/
    11231135void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
     1136        _error_("NOT NEEDED ANYMORE (remove)");
    11241137
    11251138        int        idlist[NUMVERTICES];
    … …  
    11341147
    11351148}/*}}}*/
    1136 void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
     1149void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp){/*{{{*/
    11371150
    11381151        int         sidlist[NUMVERTICES];
    … …  
    11431156        IssmDouble  value,gradient;
    11441157
    1145         this->GetVerticesConnectivityList(&connectivity[0]);
    1146         this->GetVerticesSidList(&sidlist[0]);
    1147         this->GetVerticesLidList(&lidlist[0]);
    1148 
     1158        /*Get relevant inputs*/
    11491159        ElementInput2* control_value    = this->inputs2->GetControlInput2Data(control_enum,"value");    _assert_(control_value);
    11501160        ElementInput2* control_gradient = this->inputs2->GetControlInput2Data(control_enum,"gradient"); _assert_(control_gradient);
    1151         control_value->Serve(NUMVERTICES,&lidlist[0]);
    1152         control_gradient->Serve(NUMVERTICES,&lidlist[0]);
    1153 
    1154         GaussTria* gauss=new GaussTria();
    1155         for (int iv=0;iv<NUMVERTICES;iv++){
    1156                 gauss->GaussVertex(iv);
    1157 
    1158                 control_value->GetInputValue(&value,gauss);
    1159                 control_gradient->GetInputValue(&gradient,gauss);
    1160 
    1161                 values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
    1162                 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
    1163         }
    1164         delete gauss;
    1165 
    1166         vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
    1167         vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
     1161
     1162        if(control_interp==P1Enum){
     1163                _assert_(control_value->GetInputInterpolationType()==P1Enum);
     1164                _assert_(control_gradient->GetInputInterpolationType()==P1Enum);
     1165
     1166                this->GetVerticesConnectivityList(&connectivity[0]);
     1167                this->GetVerticesSidList(&sidlist[0]);
     1168                this->GetVerticesLidList(&lidlist[0]);
     1169
     1170                control_value->Serve(NUMVERTICES,&lidlist[0]);
     1171                control_gradient->Serve(NUMVERTICES,&lidlist[0]);
     1172
     1173                GaussTria* gauss=new GaussTria();
     1174                for (int iv=0;iv<NUMVERTICES;iv++){
     1175                        gauss->GaussVertex(iv);
     1176
     1177                        control_value->GetInputValue(&value,gauss);
     1178                        control_gradient->GetInputValue(&gradient,gauss);
     1179
     1180                        values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
     1181                        gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
     1182                }
     1183                delete gauss;
     1184
     1185                vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
     1186                vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
     1187        }
     1188        else if(control_interp==P0Enum){
     1189                _assert_(control_value->GetInputInterpolationType()==P0Enum);
     1190                _assert_(control_gradient->GetInputInterpolationType()==P0Enum);
     1191
     1192                control_value->Serve(1,&this->lid);
     1193                control_gradient->Serve(1,&this->lid);
     1194
     1195                vector_control->SetValue(this->sid,control_value->element_values[0],ADD_VAL);
     1196                vector_gradient->SetValue(this->sid,control_gradient->element_values[0],ADD_VAL);
     1197        }
     1198        else{
     1199                _error_("not supported");
     1200        }
    11681201
    11691202}/*}}}*/
    … …  
    23902423void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    23912424
     2425        _error_("DO WE NEED THIS FUNCTION?? REALLY?? see below");
     2426
    23922427        /*Get out if this is not an element input*/
    23932428        if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
    … …  
    24072442        int interpolation = input->GetInterpolation();
    24082443        switch(interpolation){
     2444                case P0Enum:
     2445                        input->Serve(1,&this->lid);
     2446                        break;
    24092447                case P1Enum:
    24102448                        numindices = NUMVERTICES;
    … …  
    24492487
    24502488                                TriaInput2* triainput = xDynamicCast<TriaInput2*>(input);
    2451                                 if(triainput->GetInputInterpolationType()!=P1Enum) _error_("not supported yet");
    2452                                 input->Serve(NUMVERTICES,&lidlist[0]);
    24532489
    24542490                                /*Create list of indices and values for global vector*/
    24552491                                GradientIndexing(&idlist[0],control_index);
    2456                                 for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
    2457                                 vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     2492
     2493                                if(triainput->GetInputInterpolationType()==P1Enum){
     2494                                        input->Serve(NUMVERTICES,&lidlist[0]);
     2495                                        for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
     2496                                        vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     2497                                }
     2498                                else if(triainput->GetInputInterpolationType()==P0Enum){
     2499                                        input->Serve(1,&this->lid);
     2500                                        vector->SetValue(idlist[0],triainput->element_values[0],INS_VAL);
     2501                                }
     2502                                else{
     2503                                        _error_("not supported yet");
     2504                                }
    24582505                                break;
    24592506                          }
    … …  
    24622509                                {
    24632510                                        TransientInput2* transientinput = xDynamicCast<TransientInput2*>(input);
    2464                                         int N = transientinput->numtimesteps;
    2465                                         int* M=NULL;
     2511                                        int  N = transientinput->numtimesteps;
     2512                                        int* M = NULL;
    24662513                                        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    24672514                                        int* idlist = xNew<int>(NUMVERTICES*N);
    … …  
    39073954}
    39083955/*}}}*/
    3909 void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){/*{{{*/
     3956void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int M, int N){/*{{{*/
    39103957
    39113958        IssmDouble  values[NUMVERTICES];
    39123959        int         lidlist[NUMVERTICES];
    3913         int         idlist[NUMVERTICES],control_init;
     3960        int         idlist[NUMVERTICES];
    39143961
    39153962        /*Get Domain type*/
    … …  
    39183965
    39193966        /*Specific case for depth averaged quantities*/
    3920         control_init=control_enum;
     3967        int control_init=control_enum;
    39213968        if(domaintype==Domain2DverticalEnum){
    39223969                if(control_enum==MaterialsRheologyBbarEnum){
    … …  
    39333980        if(!IsInputEnum(control_enum)) return;
    39343981
     3982        ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"value");   _assert_(input);
    39353983        this->GetVerticesLidList(&lidlist[0]);
    3936         ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"value");   _assert_(input);
     3984        GradientIndexing(&idlist[0],control_index);
    39373985
    39383986        /*Get values on vertices*/
    3939         for(int n=0;n<N;n++){
     3987        if(input->ObjectEnum()==TriaInput2Enum && input->GetInputInterpolationType()==P1Enum){
     3988                _assert_(N==1);
    39403989                for(int i=0;i<NUMVERTICES;i++){
    3941                         idlist[i]=offset + this->vertices[i]->Sid()+n*M;
    3942                         values[i]=vector[idlist[i]];
    3943                 }
    3944                 if(input->ObjectEnum()==TriaInput2Enum){
    3945                         input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
    3946                 }
    3947                 else if(input->ObjectEnum()==TransientInput2Enum){
     3990                        values[i] = vector[idlist[i]];
     3991                }
     3992                input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
     3993        }
     3994        else if(input->ObjectEnum()==TriaInput2Enum && input->GetInputInterpolationType()==P0Enum){
     3995                _assert_(N==1);
     3996                input->SetInput(P0Enum,this->lid,vector[idlist[0]]);
     3997        }
     3998        else if(input->ObjectEnum()==TransientInput2Enum){
     3999                for(int n=0;n<N;n++){
    39484000                        _error_("not implemented");
    39494001                        //Input* new_input = new TriaInput(control_enum,values,P1Enum);
    … …  
    39514003                        //controlinput->Configure(parameters);
    39524004                }
    3953                 else _error_("Type not supported");
    3954         }
     4005        }
     4006        else _error_("Type not supported");
    39554007}
    39564008/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r25256 r25317  
    6666                void        ComputeSurfaceNormalVelocity();
    6767                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs2* inputs2in);
    68                 void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M);
     68                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp);
    6969                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    70                 void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
     70                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum,int control_interp);
    7171                void        CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
    7272                int         EdgeOnBaseIndex();
    … …  
    130130                void        SetElementInput(Inputs2* inputs2,int enum_in,IssmDouble values);
    131131                void        SetElementInput(Inputs2* inputs2,int numindices,int* indices,IssmDouble* values,int enum_in);
    132                 void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M);
     132                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int M,int N);
    133133                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
    134134                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r25139 r25317  
    21742174        int         num_controls,step;
    21752175        IssmDouble  time;
    2176         int        *control_type = NULL;
     2176        int        *control_type   = NULL;
     2177        int        *control_interp = NULL;
     2178        int        *M = NULL;
     2179        int        *N = NULL;
    21772180
    21782181        /*recover results*/
    … …  
    21832186        this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    21842187        this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     2188        this->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     2189        this->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     2190        this->parameters->FindParam(&control_interp,NULL,ControlInputInterpolationEnum);
    21852191        this->parameters->FindParam(&step,StepEnum);
    21862192        this->parameters->FindParam(&time,TimeEnum);
    … …  
    21992205
    22002206                /*Allocate vector*/
    2201                 Vector<IssmPDouble> *vector_control  = new Vector<IssmPDouble>(this->vertices->NumberOfVertices());
    2202                 Vector<IssmPDouble> *vector_gradient = new Vector<IssmPDouble>(this->vertices->NumberOfVertices());
     2207                Vector<IssmPDouble> *vector_control  = new Vector<IssmPDouble>(M[i]*N[i]);
     2208                Vector<IssmPDouble> *vector_gradient = new Vector<IssmPDouble>(M[i]*N[i]);
    22032209
    22042210                /*Fill in vector*/
    22052211                for(int j=0;j<elements->Size();j++){
    22062212                        Element* element=(Element*)elements->GetObjectByOffset(j);
    2207                         element->ControlToVectors(vector_control,vector_gradient,control_enum);
     2213                        element->ControlToVectors(vector_control,vector_gradient,control_enum,control_interp[i]);
    22082214                }
    22092215                vector_control->Assemble();
    … …  
    22162222        /*Clean up and return*/
    22172223        xDelete<int>(control_type);
     2224        xDelete<int>(control_interp);
     2225        xDelete<int>(M);
     2226        xDelete<int>(N);
    22182227}
    22192228/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs2/ControlInput2.cpp

    r24335 r25317  
    2929        this->control_id  = id;
    3030        this->layout_enum = input_layout_enum;
    31 
    32         _assert_(interp==P1Enum);
    3331
    3432        switch(this->layout_enum){

  • issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp

    r25293 r25317  
    912912        /*Create it if necessary*/
    913913        if(this->inputs[id]){
     914                int* tenp = xNew<int>(3);
    914915                if(this->inputs[id]->ObjectEnum()!=PentaInput2Enum) _error_("cannot add Element values to a "<<EnumToStringx(this->inputs[id]->ObjectEnum()));
    915916        }

  • issm/trunk-jpl/src/c/classes/Params/Parameters.cpp

    r24335 r25317  
    386386
    387387        /*cast to "passive"*/
     388        #ifdef _HAVE_AD_
    388389        *pscalar=reCast<IssmPDouble>(intermediary);
     390        #else
     391        *pscalar=intermediary;
     392        #endif
    389393}
    390394/*}}}*/
    … …  
    401405
    402406        /*Make output passive*/
     407        #ifdef _HAVE_AD_
    403408        IssmPDouble* output = xNew<IssmPDouble>(n);
    404409        for(int i=0;i<n;i++) output[i] = reCast<IssmPDouble>(vector[i]);
     410        xDelete<IssmDouble>(vector);
     411        if(pvec) *pvec = output;
     412        #else
     413        if(pvec) *pvec = vector;
     414        #endif
    405415
    406416        /*assign output pointers*/
    407         if(pvec) *pvec = output;
    408417        if(pM)   *pM   = n;
    409418}/*}}}*/

  • issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp

    r23350 r25317  
    4646        double *X  = NULL;
    4747        double *G  = NULL;
     48
     49        /*Get control sizes*/
     50        int* M = NULL;
     51        int* N = NULL;
     52        femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     53        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    4854
    4955        /*Recover some parameters*/
    … …  
    7379        long      io           = 6;         /*Channel number for the output*/
    7480
    75         /*Optimization criterions*/
     81        /*Optimization criterions (need to be cast to long for m1qn3)*/
    7682        long niter = long(maxsteps); /*Maximum number of iterations*/
    77         long nsim  = long(maxiter);/*Maximum number of function calls*/
     83        long nsim  = long(maxiter);  /*Maximum number of function calls*/
    7884
    7985        /*Get initial guess*/
    8086        GetPassiveVectorFromControlInputsx(&X,&intn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    81         _assert_(intn==numberofvertices*num_controls);
    8287
    8388        /*Get problem dimension and initialize gradient and initial guess*/
    … …  
    8691
    8792        /*Scale control for M1QN3*/
     93        int offset = 0;
    8894        for(int c=0;c<num_controls;c++){
    89                 for(int i=0;i<numberofvertices;i++){
    90                         int index = numberofvertices*c+i;
    91                         X[index] = X[index]/scaling_factors[c];
     95                for(int i=0;i<M[c]*N[c];i++){
     96                        int index = offset+i;
     97                        X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
    9298                }
     99                offset += M[c]*N[c];
    93100        }
    94101
    … …  
    142149        GetPassiveVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    143150        GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     151
     152        offset = 0;
    144153        for(int c=0;c<num_controls;c++){
    145                 for(int i=0;i<numberofvertices;i++){
    146                         int index = numberofvertices*c+i;
     154                for(int i=0;i<M[c]*N[c];i++){
     155                        int index = offset+i;
    147156                        X[index] = X[index]*scaling_factors[c];
    148157                        if(X[index]>XU[index]) X[index]=XU[index];
    149158                        if(X[index]<XL[index]) X[index]=XL[index];
    150159                }
     160                offset += M[c]*N[c];
    151161        }
    152162
    … …  
    179189
    180190        /*Clean-up and return*/
     191        xDelete<int>(M);
     192        xDelete<int>(N);
    181193        xDelete<double>(G);
    182194        xDelete<double>(X);
    … …  
    201213        int num_responses,num_controls,numberofvertices,solution_type;
    202214        double* scaling_factors = NULL;
     215        int* M = NULL;
     216        int* N = NULL;
    203217        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     218        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     219        femmodel->parameters->FindParam(&M,NULL,ControlInputSizeNEnum);
    204220        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    205221        femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    … …  
    212228        GetPassiveVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    213229        GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     230
     231        int offset = 0;
    214232        for(int c=0;c<num_controls;c++){
    215                 for(int i=0;i<numberofvertices;i++){
    216                         int index = numberofvertices*c+i;
     233                for(int i=0;i<M[c]*N[c];i++){
     234                        int index = offset+i;
    217235                        X[index] = X[index]*scaling_factors[c];
    218236                        if(X[index]>XU[index]) X[index]=XU[index];
    219237                        if(X[index]<XL[index]) X[index]=XL[index];
    220238                }
    221         }
     239                offset += M[c]*N[c];
     240        }
     241
    222242        #ifdef _HAVE_AD_
    223243        IssmDouble* aX=xNew<IssmDouble>(*n);
    … …  
    277297        /*Constrain Gradient*/
    278298        IssmDouble  Gnorm = 0.;
     299        offset = 0;
    279300        for(int c=0;c<num_controls;c++){
    280                 for(int i=0;i<numberofvertices;i++){
    281                         int index = numberofvertices*c+i;
     301                for(int i=0;i<M[c]*N[c];i++){
     302                        int index = offset+i;
    282303                        if(X[index]>=XU[index]) G[index]=0.;
    283304                        if(X[index]<=XL[index]) G[index]=0.;
    … …  
    286307                        Gnorm += G[index]*G[index];
    287308                }
     309                offset += M[c]*N[c];
    288310        }
    289311        Gnorm = sqrt(Gnorm);
    … …  
    296318        /*Clean-up and return*/
    297319        *Jlisti = (*Jlisti) +1;
     320        xDelete<int>(M);
     321        xDelete<int>(N);
    298322        xDelete<double>(XU);
    299323        xDelete<double>(XL);

  • issm/trunk-jpl/src/c/modules/ControlInputSetGradientx/ControlInputSetGradientx.cpp

    r22740 r25317  
    99void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* gradient){
    1010
    11         bool isautodiff;
    12         parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    13         if(isautodiff){
    14                 /*Intermediaries*/
    15                 int  num_controls;
    16                 int *control_type = NULL;
    17                 int* M_all = NULL;
    18                 int* N_all = NULL;
     11        /*Intermediaries*/
     12        int  num_controls;
     13        int *control_type = NULL;
     14        int* M_all = NULL;
     15        int* N_all = NULL;
     16        int* interp_all = NULL;
    1917
    20                 /*Retrieve some parameters*/
    21                 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    22                 parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    23                 parameters->FindParam(&M_all,NULL,ControlInputSizeMEnum);
    24                 parameters->FindParam(&N_all,NULL,ControlInputSizeNEnum);
     18        /*Retrieve some parameters*/
     19        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     20        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     21        parameters->FindParam(&M_all,NULL,ControlInputSizeMEnum);
     22        parameters->FindParam(&N_all,NULL,ControlInputSizeNEnum);
     23        parameters->FindParam(&interp_all,NULL,ControlInputInterpolationEnum);
    2524
    26                 int offset = 0;
    27                 for(int i=0;i<num_controls;i++){
    28                         for(int j=0;j<elements->Size();j++){
    29                                 Element* element=(Element*)elements->GetObjectByOffset(j);
    30                                 element->ControlInputSetGradient(gradient,control_type[i],i,offset,N_all[i],M_all[i]);
    31                         }
    32                         offset+=M_all[i]*N_all[i];
     25        int offset = 0;
     26        for(int i=0;i<num_controls;i++){
     27                for(int j=0;j<elements->Size();j++){
     28                        Element* element=(Element*)elements->GetObjectByOffset(j);
     29                        element->ControlInputSetGradient(gradient,control_type[i],i,offset,M_all[i],N_all[i],interp_all[i]);
    3330                }
     31                offset+=M_all[i]*N_all[i];
     32        }
    3433
    35                 /*Clean up and return*/
    36                 xDelete<int>(control_type);
    37         }
    38         else{
    39                 int  num_controls;
    40                 int *control_type = NULL;
    41 
    42                 /*Retrieve some parameters*/
    43                 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    44                 parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    45 
    46                 int offset = 0;
    47                 for(int i=0;i<num_controls;i++){
    48                         for(int j=0;j<elements->Size();j++){
    49                                 Element* element=(Element*)elements->GetObjectByOffset(j);
    50                                 element->ControlInputSetGradient(gradient,control_type[i],i);
    51                         }
    52                 }
    53 
    54                 /*Clean up and return*/
    55                 xDelete<int>(control_type);
    56         }
     34        /*Clean up and return*/
     35        xDelete<int>(control_type);
     36        xDelete<int>(M_all);
     37        xDelete<int>(N_all);
     38        xDelete<int>(interp_all);
    5739
    5840}
    … …  
    6042
    6143        /*Serialize gradient*/
    62         IssmDouble* serial_gradient=NULL;
    63         serial_gradient=gradient->ToMPISerial();
     44        IssmDouble* serial_gradient=gradient->ToMPISerial();
    6445
    6546        ControlInputSetGradientx(elements,nodes,vertices, loads, materials, parameters,serial_gradient);

  • issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp

    r24335 r25317  
    99void GetVectorFromControlInputsx(Vector<IssmDouble>** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,const char* data){/*{{{*/
    1010
    11         bool isautodiff;
    12         parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    13         if(isautodiff){
    14                 int*  N = NULL;
    15                 int*  M = NULL;
    16                 int  num_controls;
    17                 int* control_type = NULL;
    18                 Vector<IssmDouble>*  vector=NULL;
     11        int  num_controls;
     12        int* N = NULL;
     13        int* M = NULL;
     14        int* control_type = NULL;
    1915
    20                 /*Retrieve some parameters*/
    21                 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    22                 parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    23                 parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    24                 parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     16        /*Retrieve some parameters*/
     17        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     18        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     19        parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     20        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    2521
    26                 /*1. Get vector size*/
    27                 int size = 0;
    28                 for(int i=0;i<num_controls;i++) size+=M[i]*N[i];
     22        /*1. Get vector size*/
     23        int size = 0;
     24        for(int i=0;i<num_controls;i++) size+=M[i]*N[i];
    2925
    30                 /*2. Allocate vector*/
    31                 vector=new Vector<IssmDouble>(size);
     26        /*2. Allocate vector*/
     27        Vector<IssmDouble>* vector=new Vector<IssmDouble>(size);
    3228
    33                 /*3. Populate vector*/
    34                 int offset = 0;
    35                 for(int i=0;i<num_controls;i++){
    36                         for(int j=0;j<elements->Size();j++){
    37                                 Element* element=(Element*)elements->GetObjectByOffset(j);
    38                                 element->GetVectorFromControlInputs(vector,control_type[i],i,data,offset);
    39                         }
    40                         offset += M[i]*N[i];
     29        /*3. Populate vector*/
     30        int offset = 0;
     31        for(int i=0;i<num_controls;i++){
     32                for(int j=0;j<elements->Size();j++){
     33                        Element* element=(Element*)elements->GetObjectByOffset(j);
     34                        element->GetVectorFromControlInputs(vector,control_type[i],i,data,offset);
    4135                }
     36                offset += M[i]*N[i];
     37        }
     38        vector->Assemble();
    4239
    43                 vector->Assemble();
    44 
    45                 /*Assign output pointers:*/
    46                 xDelete<int>(control_type);
    47                 xDelete<int>(M);
    48                 xDelete<int>(N);
    49 
    50                 *pvector=vector;
    51         }
    52         else{
    53                 int  num_controls;
    54                 int* control_type = NULL;
    55                 Vector<IssmDouble>*  vector=NULL;
    56 
    57                 /*Retrieve some parameters*/
    58                 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    59                 parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    60 
    61                 /*2. Allocate vector*/
    62                 vector=new Vector<IssmDouble>(num_controls*vertices->NumberOfVertices());
    63 
    64                 /*3. Populate vector*/
    65                 int offset = 0;
    66                 for(int i=0;i<num_controls;i++){
    67                         for(int j=0;j<elements->Size();j++){
    68                                 Element* element=(Element*)elements->GetObjectByOffset(j);
    69                                 element->GetVectorFromControlInputs(vector,control_type[i],i,data);
    70                         }
    71                 }
    72                 vector->Assemble();
    73 
    74                 /*Assign output pointers:*/
    75                 xDelete<int>(control_type);
    76                 *pvector=vector;
    77         }
    78 
     40        /*Assign output pointers:*/
     41        xDelete<int>(control_type);
     42        xDelete<int>(M);
     43        xDelete<int>(N);
     44        *pvector=vector;
    7945}/*}}}*/
    8046void GetVectorFromControlInputsx( IssmDouble** pvector,int *pN, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, const char* data){/*{{{*/

  • issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp

    r18258 r25317  
    99void Gradjx(Vector<IssmDouble>** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
    1010
    11         int          numberofvertices;
    12         int          num_controls,analysisenum;
    13         IssmDouble   norm_inf;
    14         IssmDouble  *norm_list     = NULL;
    15         int     *control_type  = NULL;
     11        int         numberofvertices;
     12        int         num_controls,analysisenum;
     13        IssmDouble  norm_inf;
     14        IssmDouble *norm_list      = NULL;
     15        int        *control_type   = NULL;
     16        int        *control_interp = NULL;
     17        int        *M = NULL;
     18        int        *N = NULL;
    1619        Vector<IssmDouble>  *gradient      = NULL;
    1720        Vector<IssmDouble> **gradient_list = NULL;
    … …  
    2124
    2225        /*retrieve some parameters: */
    23         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
     26        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);    _assert_(num_controls);
    2427        parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     28   parameters->FindParam(&control_interp,NULL,ControlInputInterpolationEnum);
     29        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     30        parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    2531        numberofvertices=vertices->NumberOfVertices();
    2632
    … …  
    3137        /*Allocate gradient_list */
    3238        gradient_list = xNew<Vector<IssmDouble>*>(num_controls);
    33         norm_list = xNew<IssmDouble>(num_controls);
    34         for(int i=0;i<num_controls;i++){
    35                 gradient_list[i]=new Vector<IssmDouble>(num_controls*numberofvertices);
    36         }
    37         gradient=new Vector<IssmDouble>(num_controls*numberofvertices);
     39        norm_list     = xNew<IssmDouble>(num_controls);
     40        int totalsize = 0;
     41        for(int i=0;i<num_controls;i++) totalsize += M[i]*N[i];
     42        for(int i=0;i<num_controls;i++) gradient_list[i]=new Vector<IssmDouble>(totalsize);
    3843
    3944        /*Compute all gradient_list*/
    … …  
    4146                for(int j=0;j<elements->Size();j++){
    4247                        Element* element=(Element*)elements->GetObjectByOffset(j);
    43                         analysis->GradientJ(gradient_list[i],element,control_type[i],i);
     48                        analysis->GradientJ(gradient_list[i],element,control_type[i],control_interp[i],i);
    4449                }
    4550                gradient_list[i]->Assemble();
    … …  
    4853
    4954        /*Add all gradient_list together*/
     55        gradient=new Vector<IssmDouble>(totalsize);
    5056        for(int i=0;i<num_controls;i++){
    5157                gradient->AXPY(gradient_list[i],1.0);
    … …  
    6268        xDelete<Vector<IssmDouble>*>(gradient_list);
    6369        xDelete<int>(control_type);
     70        xDelete<int>(control_interp);
     71        xDelete<int>(M);
     72        xDelete<int>(N);
    6473        if(pnorm_list){
    6574                *pnorm_list=norm_list;
    … …  
    7382void Gradjx(IssmDouble** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
    7483
    75         /*output: */
    76         IssmDouble* gradient=NULL;
     84        /*Get gradient: */
     85        Vector<IssmDouble>* vec_gradient=NULL;
     86        Gradjx(&vec_gradient,pnorm_list,elements,nodes, vertices,loads,materials,parameters);
    7787
    78         /*intermediary: */
    79         Vector<IssmDouble>* vec_gradient=NULL;
     88        /*Serialize*/
     89        IssmDouble* gradient=vec_gradient->ToMPISerial();
    8090
    81         Gradjx(&vec_gradient,pnorm_list,elements,nodes, vertices,loads,materials,parameters);
    82         gradient=vec_gradient->ToMPISerial();
    83 
    84         /*Free ressources:*/
     91        /*Free ressources and assign output pointer*/
    8592        delete vec_gradient;
    86 
    87         /*Assign output pointers:*/
    8893        *pgradient=gradient;
    8994}

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r25201 r25317  
    110110        }
    111111
     112        int* M_all = xNew<int>(num_controls);
     113        int* N_all = xNew<int>(num_controls);
     114        int* Interp_all = xNew<int>(num_controls);
     115
     116        int offset = 0;
    112117        for(int i=0;i<num_controls;i++){
    113118                control = control_enums[i];
    … …  
    133138                                _error_("Control " << EnumToStringx(control) << " not implemented yet");
    134139                }
    135                 if(M!=iomodel->numberofvertices && N!=1) _error_("not supported yet");
     140
     141                /*Transient independents not supported outside of AD*/
     142                if(N!=1) _error_("Transient controls not supported yet");
     143                N_all[i] = N;
     144                M_all[i] = M;
     145
     146                if(M==iomodel->numberofvertices){
     147                        Interp_all[i] = P1Enum;
     148                }
     149                else if(M==iomodel->numberofelements){
     150                        Interp_all[i] = P0Enum;
     151                }
     152                else{
     153                        _error_("Control size not supported");
     154                }
    136155
    137156                /*Special case if 3d*/
    … …  
    143162                for(int j=0;j<elements->Size();j++){
    144163                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    145                         element->ControlInputCreate(independent,&independents_min[i*iomodel->numberofvertices],&independents_max[i*iomodel->numberofvertices],inputs2,iomodel,M,N,scale,control,i+1);
     164                        element->ControlInputCreate(independent,&independents_min[offset],&independents_max[offset],inputs2,iomodel,M,N,scale,control,i+1);
    146165                }
    147166                xDelete<IssmDouble>(independent);
    148         }
     167
     168                offset += M*N;
     169        }
     170        parameters->AddObject(new IntVecParam(ControlInputSizeMEnum,M_all,num_controls));
     171        parameters->AddObject(new IntVecParam(ControlInputSizeNEnum,N_all,num_controls));
     172        parameters->AddObject(new IntVecParam(ControlInputInterpolationEnum,Interp_all,num_controls));
     173        xDelete<int>(M_all);
     174        xDelete<int>(N_all);
     175        xDelete<int>(Interp_all);
    149176        xDelete<IssmDouble>(independents_min);
    150177        xDelete<IssmDouble>(independents_max);
    151 
    152178
    153179        /*Free data: */

  • issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp

    r23293 r25317  
    99void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
    1010
    11         bool isautodiff;
    12         femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    13         if(isautodiff){
    14                 int  num_controls;
    15                 int* control_type = NULL;
    16                 int* M = NULL;
    17                 int* N = NULL;
     11        int  num_controls;
     12        int* control_type = NULL;
     13        int* M = NULL;
     14        int* N = NULL;
    1815
    19                 /*Retrieve some parameters*/
    20                 femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    21                 femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    22                 femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    23                 femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     16        /*Retrieve some parameters*/
     17        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     18        femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     19        femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     20        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    2421
    25                 int offset = 0;
    26                 for(int i=0;i<num_controls;i++){
    27                         for(int j=0;j<femmodel->elements->Size();j++){
    28                                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
    29                                 element->SetControlInputsFromVector(vector,control_type[i],i,offset,N[i],M[i]);
    30                         }
     22        int offset = 0;
     23        for(int i=0;i<num_controls;i++){
     24                for(int j=0;j<femmodel->elements->Size();j++){
     25                        Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
     26                        element->SetControlInputsFromVector(vector,control_type[i],i,offset,M[i],N[i]);
     27                }
    3128                offset += M[i]*N[i];
    32                 }
     29        }
    3330
    34                 xDelete<int>(control_type);
    35                 xDelete<int>(M);
    36                 xDelete<int>(N);
    37         }
    38         else{
    39 
    40                 int  num_controls;
    41                 int* control_type = NULL;
    42                 femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    43                 femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    44                 int offset = 0;
    45                 for(int i=0;i<num_controls;i++){
    46                         for(int j=0;j<femmodel->elements->Size();j++){
    47                                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
    48                                 element->SetControlInputsFromVector(vector,control_type[i],i);
    49                         }
    50                 }
    51                 xDelete<int>(control_type);
    52         }
     31        xDelete<int>(control_type);
     32        xDelete<int>(M);
     33        xDelete<int>(N);
    5334}
    5435

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r25308 r25317  
    111111syn keyword cConstant ControlInputSizeMEnum
    112112syn keyword cConstant ControlInputSizeNEnum
     113syn keyword cConstant ControlInputInterpolationEnum
    113114syn keyword cConstant DamageC1Enum
    114115syn keyword cConstant DamageC2Enum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r25308 r25317  
    105105        ControlInputSizeMEnum,
    106106        ControlInputSizeNEnum,
     107        ControlInputInterpolationEnum,
    107108        DamageC1Enum,
    108109        DamageC2Enum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r25308 r25317  
    113113                case ControlInputSizeMEnum : return "ControlInputSizeM";
    114114                case ControlInputSizeNEnum : return "ControlInputSizeN";
     115                case ControlInputInterpolationEnum : return "ControlInputInterpolation";
    115116                case DamageC1Enum : return "DamageC1";
    116117                case DamageC2Enum : return "DamageC2";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r25308 r25317  
    113113              else if (strcmp(name,"ControlInputSizeM")==0) return ControlInputSizeMEnum;
    114114              else if (strcmp(name,"ControlInputSizeN")==0) return ControlInputSizeNEnum;
     115              else if (strcmp(name,"ControlInputInterpolation")==0) return ControlInputInterpolationEnum;
    115116              else if (strcmp(name,"DamageC1")==0) return DamageC1Enum;
    116117              else if (strcmp(name,"DamageC2")==0) return DamageC2Enum;
    … …  
    136137              else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum;
    137138              else if (strcmp(name,"EarthId")==0) return EarthIdEnum;
    138               else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum;
     142              if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
     143              else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum;
    143144              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
    144145              else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
    … …  
    259260              else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
    260261              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    261               else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
     265              if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
     266              else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
    266267              else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum;
    267268              else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
    … …  
    382383              else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
    383384              else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
    384               else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum;
     388              if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
     389              else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum;
    389390              else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum;
    390391              else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum;
    … …  
    505506              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    506507              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    507               else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"BasalStressy")==0) return BasalStressyEnum;
     511              if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
     512              else if (strcmp(name,"BasalStressy")==0) return BasalStressyEnum;
    512513              else if (strcmp(name,"BasalStress")==0) return BasalStressEnum;
    513514              else if (strcmp(name,"Base")==0) return BaseEnum;
    … …  
    628629              else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
    629630              else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum;
    630               else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"HydrologyGapHeight")==0) return HydrologyGapHeightEnum;
     634              if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
     635              else if (strcmp(name,"HydrologyGapHeight")==0) return HydrologyGapHeightEnum;
    635636              else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
    636637              else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    … …  
    751752              else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
    752753              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    753               else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbEC")==0) return SmbECEnum;
     757              if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
     758              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
    758759              else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    759760              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    … …  
    874875              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    875876              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    876               else if (strcmp(name,"Vz")==0) return VzEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"VzFS")==0) return VzFSEnum;
     880              if (strcmp(name,"Vz")==0) return VzEnum;
     881              else if (strcmp(name,"VzFS")==0) return VzFSEnum;
    881882              else if (strcmp(name,"VzHO")==0) return VzHOEnum;
    882883              else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    … …  
    997998              else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
    998999              else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
    999               else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
     1003              if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
     1004              else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
    10041005              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    10051006              else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum;
    … …  
    11201121              else if (strcmp(name,"Gset")==0) return GsetEnum;
    11211122              else if (strcmp(name,"Gsl")==0) return GslEnum;
    1122               else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
     1126              if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
     1127              else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
    11271128              else if (strcmp(name,"Hook")==0) return HookEnum;
    11281129              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    … …  
    12431244              else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
    12441245              else if (strcmp(name,"P0DG")==0) return P0DGEnum;
    1245               else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"P1P1")==0) return P1P1Enum;
     1249              if (strcmp(name,"P1DG")==0) return P1DGEnum;
     1250              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    12501251              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    12511252              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    … …  
    13661367              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    13671368              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    1368               else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
     1372              if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
     1373              else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
    13731374              else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
    13741375              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
Note: See TracChangeset for help on using the changeset viewer.