Changeset 18403


Ignore:
Timestamp:
08/15/14 15:35:51 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: simplified balancethickness2 to nu (got rid of nux,nuy) and separating cost functions

Location:
issm/trunk-jpl/src
Files:
6 added
2 deleted
17 edited

Legend:

Unmodified
Added
Removed

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

    r18394 r18403  
    5252        IssmDouble  NUMxUbar,NUMyUbar,DENUbar;
    5353        IssmDouble  vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
    54         IssmDouble  nux,nuy,phi,dphi[2];
     54        IssmDouble  nu,phi,dphi[2];
    5555        int        *responses = NULL;
    5656        IssmDouble *xyz_list  = NULL;
    … …  
    7373        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
    7474        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
    75         Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
    76         Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
     75        Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
    7776
    7877        /* Start  looping on the number of gaussian points: */
    … …  
    8786                vxobs_input->GetInputValue(&vxobs,gauss);
    8887                vyobs_input->GetInputValue(&vyobs,gauss);
    89                 nux_input->GetInputValue(&nux,gauss);
    90                 nuy_input->GetInputValue(&nuy,gauss);
     88                nu_input->GetInputValue(&nu,gauss);
    9189                potential_input->GetInputValue(&phi,gauss);
    9290                potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    9391                thicknessobs_input->GetInputValue(&hobs,gauss);
    9492
    95                 vxobsbar = nux*vxobs;
    96                 vyobsbar = nuy*vyobs;
    97 
    98                 vbarobs2 = (nux*nux*vxobs*vxobs + nuy*nuy*vyobs*vyobs);
     93                vxobsbar = nu*vxobs;
     94                vyobsbar = nu*vyobs;
     95
     96                vbarobs2 = (nu*nu*vxobs*vxobs + nu*nu*vyobs*vyobs);
    9997                vbarobs  = sqrt(vbarobs2);
    10098                hu2 = hobs*hobs*vbarobs2;
    … …  
    105103
    106104                        switch(responses[resp]){
     105                                case IrrotationalH2MisfitEnum:
     106                                        /*J = (H^2 - Hobs^2)^2*/
     107                                        for(i=0;i<numnodes;i++){
     108                                                NUMxH2 = 2.*dbasis[0*numnodes+i]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
     109                                                NUMyH2 = 2.*dbasis[1*numnodes+i]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
     110                                                DENH2 = vbarobs2*vbarobs2+1.e-14;
     111                                                pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight;
     112                                        }
     113                                        break;
     114                                case IrrotationalDirectionMisfitEnum:
     115                                        /*J = 1/2 (vbar ^ gard(phi))^2*/
     116                                        for(i=0;i<numnodes;i++){
     117                                                pe->values[i]+= weight*Jdet*gauss->weight*
     118                                                  nu*nu*(vyobs*dphi[0] - vxobs*dphi[1])*
     119                                                  (vyobs*dbasis[0*numnodes+i] - vxobs*dbasis[1*numnodes+i]);
     120                                        }
     121                                        break;
    107122                                case Balancethickness2MisfitEnum:
    108                                         /*J = (H^2 - Hobs^2)^2*/
    109                                         //for(i=0;i<numnodes;i++){
    110                                         //      NUMxH2 = 2.*dbasis[0*numnodes+i]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
    111                                         //      NUMyH2 = 2.*dbasis[1*numnodes+i]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
    112                                         //      DENH2 = vbarobs2*vbarobs2+1.e-14;
    113                                         //      pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight;
    114                                         //}
    115123                                        /*J = phi^2*/
    116124                                        //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; //OK
    117125                                        /*J = grad phi ^2*/
    118126                                        //for(i=0;i<numnodes;i++) pe->values[i]+= (dphi[0]*dbasis[0*numnodes+i] + dphi[1]*dbasis[1*numnodes+i])*weight*Jdet*gauss->weight; //OK
    119                                         /*J = (ubar - nux*uobs)^2*/
     127                                        /*J = (ubar - nu*uobs)^2*/
    120128                                        //for(i=0;i<numnodes;i++){
    121129                                        //      NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0*numnodes+i];
    … …  
    124132                                        //      pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
    125133                                        //}
    126                                         /*J = 1/2 (vbar ^ gard(phi))^2*/
    127                                         for(i=0;i<numnodes;i++){
    128                                                 pe->values[i]+= weight*Jdet*gauss->weight*
    129                                                   (nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*
    130                                                   (nuy*vyobs*dbasis[0*numnodes+i] - nux*vxobs*dbasis[1*numnodes+i]);
    131                                         }
    132 
     134                                        _error_("Not implemented yet");
    133135                                        break;
    134136                                default:
    … …  
    169171        /*Deal with first part (partial derivative a J with respect to k)*/
    170172        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    171                 case Balancethickness2MisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
     173                case Balancethickness2MisfitEnum:
     174                        /*Nothing, \partial J/\partial k = 0*/
     175                        break;
     176                case IrrotationalH2MisfitEnum:
     177                        switch(control_type){
     178                                case BalancethicknessApparentMassbalanceEnum:
     179                                        /*Nothing, \partial J/\partial k = 0*/
     180                                        break;
     181                                case BalancethicknessNuEnum:
     182                                        GradientJ1Nu(element,gradient,control_index);
     183                                        break;
     184                                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
     185                        }
     186                        break;
     187                case IrrotationalDirectionMisfitEnum:
     188                        switch(control_type){
     189                                case BalancethicknessApparentMassbalanceEnum:
     190                                        /*Nothing, \partial J/\partial k = 0*/
     191                                        break;
     192                                case BalancethicknessNuEnum:
     193                                        GradientJ2Nu(element,gradient,control_index); break; //Might need to be renamed ?
     194                                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
     195                        }
     196                        break;
    172197                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    173198        }
    … …  
    176201        switch(control_type){
    177202                case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break;
    178                 case BalancethicknessNuxEnum: GradientJNux(element,gradient,control_index); break;//should be in first part....
    179                 case BalancethicknessNuyEnum: GradientJNuy(element,gradient,control_index); break;//should be in first part....
     203                case BalancethicknessNuEnum: /*Nothing: state equations do not depend on k*/ break;
    180204                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
    181205        }
    … …  
    188212
    189213        /*Intermediaries*/
    190         IssmDouble Jdet,weight;
    191         IssmDouble lambda;
     214        IssmDouble lambda,Jdet;
    192215        IssmDouble *xyz_list= NULL;
    193216
    … …  
    203226        element->GetVerticesCoordinates(&xyz_list);
    204227        element->GradientIndexing(&vertexpidlist[0],control_index);
    205         Input* adjoint_input = element->GetInput(AdjointEnum);                            _assert_(adjoint_input);
    206         Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     228        Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input);
    207229
    208230        Gauss* gauss=element->NewGauss(2);
    … …  
    212234                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    213235                element->NodalFunctionsP1(basis,gauss);
    214                 weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
    215236                adjoint_input->GetInputValue(&lambda,gauss);
    216237
    217238                /*Build gradient vector (actually -dJ/da): */
    218239                for(int i=0;i<numvertices;i++){
    219                         ge[i]+= - weight*Jdet*gauss->weight*basis[i]*lambda;
     240                        ge[i]+= - Jdet*gauss->weight*basis[i]*lambda;
    220241                        _assert_(!xIsNan<IssmDouble>(ge[i]));
    221242                }
    … …  
    230251        delete gauss;
    231252}/*}}}*/
    232 void AdjointBalancethickness2Analysis::GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     253void AdjointBalancethickness2Analysis::GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    233254
    234255        /*Intermediaries*/
    235         IssmDouble Jdet,weight;
    236         IssmDouble vxobs,vyobs,thicknessobs;
    237         IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
     256        IssmDouble vxobs,vyobs,thicknessobs,weight;
     257        IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
    238258        IssmDouble *xyz_list= NULL;
    239259
    … …  
    252272        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
    253273        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
    254         Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
    255         Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
     274        Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
    256275        Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
    257276        Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
    … …  
    263282                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    264283                element->NodalFunctionsP1(basis,gauss);
    265                 weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
     284                weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum);
    266285                vxobs_input->GetInputValue(&vxobs,gauss);
    267286                vyobs_input->GetInputValue(&vyobs,gauss);
    268                 nux_input->GetInputValue(&nux,gauss);
    269                 nuy_input->GetInputValue(&nuy,gauss);
     287                nu_input->GetInputValue(&nu,gauss);
    270288                potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    271289                thicknessobs_input->GetInputValue(&thicknessobs,gauss);
    272290
    273291                dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
    274                 nuvobs2   = nux*vxobs*nux*vxobs + nuy*vyobs*nuy*vyobs;
     292                nuvobs2   = nu*nu*(vxobs*vxobs + vyobs*vyobs);
    275293
    276294                /*Build gradient vector (actually -dJ/da): */
    277295                for(int i=0;i<numvertices;i++){
    278                         //ge[i]+= - weight*Jdet*gauss->weight*(-vxobs)*dphi[1]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i]; //THIS IS FOR Q//V
    279                         ge[i]+= - weight*Jdet*gauss->weight*basis[i]*(
    280                                                 -2.*nux*vxobs*vxobs*(
    281                                                         dphinorm2*dphinorm2 - thicknessobs*thicknessobs*dphinorm2*nuvobs2
    282                                                         )/ ( (pow(nux,6)+pow(nuy,6))*pow(vxobs,6) + 3*(nux*vxobs*nuy*vyobs)*(nux*vxobs*nuy*vyobs)*nuvobs2)
     296                        ge[i]+= -weight*Jdet*gauss->weight*basis[i]*(
     297                                                -2./pow(nu,3) * dphinorm2/(vxobs*vxobs + vyobs*vyobs) * (dphinorm2/nuvobs2 - thicknessobs*thicknessobs)
    283298                                                );
    284299                        _assert_(!xIsNan<IssmDouble>(ge[i]));
    … …  
    294309        delete gauss;
    295310}/*}}}*/
    296 void AdjointBalancethickness2Analysis::GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     311void AdjointBalancethickness2Analysis::GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    297312
    298313        /*Intermediaries*/
    299         IssmDouble Jdet,weight;
    300         IssmDouble vxobs,vyobs,thicknessobs;
    301         IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
     314        IssmDouble vxobs,vyobs,thicknessobs,weight;
     315        IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
    302316        IssmDouble *xyz_list= NULL;
    303317
    … …  
    313327        element->GetVerticesCoordinates(&xyz_list);
    314328        element->GradientIndexing(&vertexpidlist[0],control_index);
    315         Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    316         Input* vxobs_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
    317         Input* vyobs_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
    318         Input* nux_input          = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
    319         Input* nuy_input          = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
    320         Input* potential_input    = element->GetInput(PotentialEnum);             _assert_(potential_input);
    321         Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
    322 
     329        Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     330        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
     331        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
     332        Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
     333        Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
     334        Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
    323335
    324336        Gauss* gauss=element->NewGauss(2);
    … …  
    328340                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    329341                element->NodalFunctionsP1(basis,gauss);
    330                 weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
     342                weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum);
    331343                vxobs_input->GetInputValue(&vxobs,gauss);
    332344                vyobs_input->GetInputValue(&vyobs,gauss);
    333                 nux_input->GetInputValue(&nux,gauss);
    334                 nuy_input->GetInputValue(&nuy,gauss);
     345                nu_input->GetInputValue(&nu,gauss);
    335346                potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    336347                thicknessobs_input->GetInputValue(&thicknessobs,gauss);
    337348
    338349                dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
    339                 nuvobs2   = nux*vxobs*nux*vxobs + nuy*vyobs*nuy*vyobs;
    340 
     350                nuvobs2   = nu*vxobs*nu*vxobs + nu*vyobs*nu*vyobs;
    341351
    342352                /*Build gradient vector (actually -dJ/da): */
    343353                for(int i=0;i<numvertices;i++){
    344                         //ge[i]+= - weight*Jdet*gauss->weight*(+vyobs)*dphi[0]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
    345                         ge[i]+= - weight*Jdet*gauss->weight*basis[i]*(
    346                                                 -2.*nuy*vyobs*vyobs*(
    347                                                         dphinorm2*dphinorm2 - thicknessobs*thicknessobs*dphinorm2*nuvobs2
    348                                                         )/ ( (pow(nux,6)+pow(nuy,6))*pow(vyobs,6) + 3*(nux*vxobs*nuy*vyobs)*(nux*vxobs*nuy*vyobs)*nuvobs2)
    349                                                 );
     354                        ge[i]+= - weight*Jdet*gauss->weight*nu*(vyobs*dphi[0] -vxobs*dphi[1])*(vyobs*dphi[0] -vxobs*dphi[1])*basis[i];
    350355                        _assert_(!xIsNan<IssmDouble>(ge[i]));
    351356                }

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

    r18390 r18403  
    2929                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    3030                void GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index);
    31                 void GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index);
    32                 void GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index);
     31                void GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index);
     32                void GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3333                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3434                void UpdateConstraints(FemModel* femmodel);

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

    r18060 r18403  
    199199
    200200        /*Intermediaries*/
    201         IssmDouble Jdet,weight;
    202         IssmDouble thickness,Dlambda[3],dp[3];
     201        IssmDouble thickness,Jdet,Dlambda[3],dp[3];
    203202        IssmDouble *xyz_list= NULL;
    204203
    … …  
    247246
    248247        /*Intermediaries*/
    249         IssmDouble Jdet,weight;
    250         IssmDouble thickness,Dlambda[3],dp[3];
     248        IssmDouble thickness,Jdet,Dlambda[3],dp[3];
    251249        IssmDouble *xyz_list= NULL;
    252250

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

    r18057 r18403  
    2727
    2828        iomodel->FetchDataToInput(elements,BalancethicknessApparentMassbalanceEnum);
    29         iomodel->FetchDataToInput(elements,BalancethicknessNuxEnum);
    30         iomodel->FetchDataToInput(elements,BalancethicknessNuyEnum);
     29        iomodel->FetchDataToInput(elements,BalancethicknessNuEnum);
    3130        iomodel->FetchDataToInput(elements,BalancethicknessVxObsEnum);
    3231        iomodel->FetchDataToInput(elements,BalancethicknessVyObsEnum);
    … …  
    203202        /*Intermediaries */
    204203        int         Hinterpolation;
    205         IssmDouble  vx,vy,vbar,nux,nuy,normdphi,dphi[2];
     204        IssmDouble  vx,vy,vbar,nu,normdphi,dphi[2];
    206205        IssmDouble* xyz_list = NULL;
    207206        int       * doflist  = NULL;
    … …  
    233232        Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input);
    234233        Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input);
    235         Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
    236         Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
     234        Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
    237235
    238236        switch(element->GetElementType()){
    … …  
    247245                vx_input->GetInputValue(&vx,gauss);
    248246                vy_input->GetInputValue(&vy,gauss);
    249                 nux_input->GetInputValue(&nux,gauss);
    250                 nuy_input->GetInputValue(&nuy,gauss);
     247                nu_input->GetInputValue(&nu,gauss);
    251248                potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    252249
    253                 vx = vx*nux; vy = vy*nuy;
     250                vx = vx*nu; vy = vy*nu;
    254251                vbar = sqrt(vx*vx + vy*vy) + 1.e-10;
    255252                normdphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);

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

    r18400 r18403  
    10541054                                name==BalancethicknessThickeningRateEnum ||
    10551055                                name==BalancethicknessApparentMassbalanceEnum ||
    1056                                 name==BalancethicknessNuxEnum ||
    1057                                 name==BalancethicknessNuyEnum ||
     1056                                name==BalancethicknessNuEnum ||
    10581057                                name==SigmaNNEnum ||
    10591058                                name==SurfaceSlopeXEnum ||

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

    r18293 r18403  
    10321032                                case ThicknessEnum:
    10331033                                case FrictionCoefficientEnum:
    1034                                 case BalancethicknessNuxEnum:
    1035                                 case BalancethicknessNuyEnum:
     1034                                case BalancethicknessNuEnum:
    10361035                                        if(iomodel->Data(control)){
    10371036                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(control)[tria_vertex_ids[j]-1];

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

    r18393 r18403  
    592592                                case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break;
    593593                                case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
    594                                 case Balancethickness2MisfitEnum:  Balancethickness2Misfitx(&double_result); break;
     594                                case Balancethickness2MisfitEnum:   Balancethickness2Misfitx(&double_result); break;
     595                                case IrrotationalH2MisfitEnum:      IrrotationalH2Misfitx(&double_result); break;
     596                                case IrrotationalDirectionMisfitEnum:  IrrotationalDirectionMisfitx(&double_result); break;
     597                                case IrrotationalVelMisfitEnum:        IrrotationalVelMisfitx(&double_result); break;
     598                                case IrrotationalAlongGradientNuEnum:  IrrotationalAlongGradientNux(&double_result); break;
     599                                case IrrotationalAcrossGradientNuEnum: IrrotationalAcrossGradientNux(&double_result); break;
    595600
    596601                           /*Vector */
    … …  
    13051310
    13061311        IssmDouble  weight,thicknessobs,thickness,potential,dpotential[2];
    1307         IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
     1312        IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nu;
    13081313        IssmDouble  Jdet;
    13091314        IssmDouble* xyz_list = NULL;
    … …  
    13261331                Input* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
    13271332                Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_input);
    1328                 Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
    1329                 Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
     1333                Input* nu_input       = element->GetInput(BalancethicknessNuEnum);   _assert_(nu_input);
    13301334
    13311335                /* Start  looping on the number of gaussian points: */
    … …  
    13481352                        vx_input->GetInputValue(&vxbar,gauss);
    13491353                        vy_input->GetInputValue(&vybar,gauss);
    1350                         nux_input->GetInputValue(&nux,gauss);
    1351                         nuy_input->GetInputValue(&nuy,gauss);
    1352 
    1353                         vxbarobs = nux*vxobs;
    1354                         vybarobs = nuy*vyobs;
     1354                        nu_input->GetInputValue(&nu,gauss);
     1355
     1356                        vxbarobs = nu*vxobs;
     1357                        vybarobs = nu*vyobs;
    13551358
    13561359                        /*J = (H^2 - Hobs^2)^2*/
    … …  
    13791382        *presponse=J;
    13801383
     1384}/*}}}*/
     1385void FemModel::IrrotationalH2Misfitx(IssmDouble* presponse){/*{{{*/
     1386
     1387        /*output: */
     1388        IssmDouble J=0.;
     1389        IssmDouble J_sum;
     1390
     1391        IssmDouble  weight,thicknessobs,thickness;
     1392        IssmDouble  Jdet;
     1393        IssmDouble* xyz_list = NULL;
     1394
     1395        /*Compute Misfit: */
     1396        for(int i=0;i<elements->Size();i++){
     1397                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1398
     1399                /*If on water, return 0: */
     1400                if(!element->IsIceInElement()) continue;
     1401
     1402                /* Get node coordinates*/
     1403                element->GetVerticesCoordinates(&xyz_list);
     1404                Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1405                Input* thickness_input   =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     1406                Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
     1407
     1408                /* Start  looping on the number of gaussian points: */
     1409                Gauss* gauss=element->NewGauss(2);
     1410                for(int ig=gauss->begin();ig<gauss->end();ig++){
     1411
     1412                        gauss->GaussPoint(ig);
     1413                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1414
     1415                        /*Get all parameters at gaussian point*/
     1416                        weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum);
     1417                        thickness_input->GetInputValue(&thickness,gauss);
     1418                        thicknessobs_input->GetInputValue(&thicknessobs,gauss);
     1419
     1420                        /*J = (H^2 - Hobs^2)^2*/
     1421                        J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
     1422                }
     1423
     1424                /*clean up and Return: */
     1425                xDelete<IssmDouble>(xyz_list);
     1426                delete gauss;
     1427        }
     1428
     1429        /*Sum all J from all cpus of the cluster:*/
     1430        ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1431        ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1432        J=J_sum;
     1433
     1434        /*Assign output pointers: */
     1435        *presponse=J;
     1436
     1437}/*}}}*/
     1438void FemModel::IrrotationalDirectionMisfitx(IssmDouble* presponse){/*{{{*/
     1439
     1440        /*output: */
     1441        IssmDouble J=0.;
     1442        IssmDouble J_sum;
     1443
     1444        IssmDouble  weight,thicknessobs,thickness,dpotential[2];
     1445        IssmDouble  vx,vy,vxobs,vyobs,nu;
     1446        IssmDouble  Jdet;
     1447        IssmDouble* xyz_list = NULL;
     1448
     1449        /*Compute Misfit: */
     1450        for(int i=0;i<elements->Size();i++){
     1451                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1452
     1453                /*If on water, return 0: */
     1454                if(!element->IsIceInElement()) continue;
     1455
     1456                /* Get node coordinates*/
     1457                element->GetVerticesCoordinates(&xyz_list);
     1458                Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1459                Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
     1460                Input* vxobs_input       =element->GetInput(BalancethicknessVxObsEnum);              _assert_(vxobs_input);
     1461                Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum);              _assert_(vyobs_input);
     1462                Input* nu_input          = element->GetInput(BalancethicknessNuEnum);                _assert_(nu_input);
     1463
     1464                /* Start  looping on the number of gaussian points: */
     1465                Gauss* gauss=element->NewGauss(2);
     1466                for(int ig=gauss->begin();ig<gauss->end();ig++){
     1467
     1468                        gauss->GaussPoint(ig);
     1469                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1470
     1471                        /*Get all parameters at gaussian point*/
     1472                        weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum);
     1473                        potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss);
     1474                        vxobs_input->GetInputValue(&vxobs,gauss);
     1475                        vyobs_input->GetInputValue(&vyobs,gauss);
     1476                        nu_input->GetInputValue(&nu,gauss);
     1477
     1478                        /*J = 1/2 (vbar ^ gard(phi))^2*/
     1479                        J += 0.5*(nu*vyobs*dpotential[0] - nu*vxobs*dpotential[1])*(nu*vyobs*dpotential[0] - nu*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
     1480                }
     1481
     1482                /*clean up and Return: */
     1483                xDelete<IssmDouble>(xyz_list);
     1484                delete gauss;
     1485        }
     1486
     1487        /*Sum all J from all cpus of the cluster:*/
     1488        ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1489        ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1490        J=J_sum;
     1491
     1492        /*Assign output pointers: */
     1493        *presponse=J;
     1494
     1495}/*}}}*/
     1496void FemModel::IrrotationalVelMisfitx(IssmDouble* presponse){/*{{{*/
     1497        _error_("not implemented yet");
     1498}/*}}}*/
     1499void FemModel::IrrotationalAlongGradientNux(IssmDouble* presponse){/*{{{*/
     1500        _error_("not implemented yet");
     1501}/*}}}*/
     1502void FemModel::IrrotationalAcrossGradientNux(IssmDouble* presponse){/*{{{*/
     1503        _error_("not implemented yet");
    13811504}/*}}}*/
    13821505void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/

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

    r18363 r18403  
    8181                void BalancethicknessMisfitx(IssmDouble* pV);
    8282                void Balancethickness2Misfitx(IssmDouble* pV);
     83                void IrrotationalH2Misfitx(IssmDouble* pV);
     84                void IrrotationalDirectionMisfitx(IssmDouble* pV);
     85                void IrrotationalVelMisfitx(IssmDouble* pV);
     86                void IrrotationalAlongGradientNux(IssmDouble* pV);
     87                void IrrotationalAcrossGradientNux(IssmDouble* pV);
    8388                #ifdef  _HAVE_DAKOTA_
    8489                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

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

    r18193 r18403  
    1212
    1313        int         solution_type,n;
     14        int         num_responses;
    1415        IssmDouble  j0,j;
    1516        IssmDouble  Ialpha,exponent,alpha;
     17        IssmDouble* jlist = NULL;
    1618        IssmDouble *G = NULL;
    1719        IssmDouble *X = NULL;
    … …  
    2527        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    2628        femmodel->parameters->SetParam(false,SaveResultsEnum);
     29        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    2730
    2831        /*Get initial guess*/
    … …  
    4548        if(VerboseControl()) _printf0_("   Compute Adjoint\n");
    4649        adjointcore(femmodel);
     50
    4751        if(VerboseControl()) _printf0_("   Compute Initial cost function\n");
    48         femmodel->CostFunctionx(&j0,NULL,NULL);
     52        femmodel->CostFunctionx(&j0,&jlist,NULL);
     53        _printf0_("Initial J(x+dk)      |  List of contributions\n");
     54        _printf0_("_____________________________________________\n");
     55        _printf0_("J(x) = "<<setw(12)<<setprecision(7)<<j0<<"  |  ");
     56        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<jlist[i]);
     57        _printf0_("\n");
     58        xDelete<IssmDouble>(jlist);
     59
    4960        if(VerboseControl()) _printf0_("   Compute Gradient\n");
    5061        Gradjx(&G,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    … …  
    95106        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,output,num,2,1,0));
    96107        #endif
     108        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
     109        femmodel->OutputControlsx(&femmodel->results);
    97110
    98111        /*Clean up and return*/

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

    r18003 r18403  
    4141                        case ThicknessEnum:
    4242                        case FrictionCoefficientEnum:
    43                         case BalancethicknessNuxEnum:
    44                         case BalancethicknessNuyEnum:
     43                        case BalancethicknessNuEnum:
    4544                        case BalancethicknessApparentMassbalanceEnum:
    4645                                iomodel->FetchData(1,control);

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

    r18385 r18403  
    301301        BalancethicknessApparentMassbalanceEnum,
    302302        Balancethickness2MisfitEnum,
    303         BalancethicknessNuxEnum,
    304         BalancethicknessNuyEnum,
     303        IrrotationalH2MisfitEnum,
     304        IrrotationalDirectionMisfitEnum,
     305        IrrotationalVelMisfitEnum,
     306        IrrotationalAlongGradientNuEnum,
     307        IrrotationalAcrossGradientNuEnum,
     308        BalancethicknessNuEnum,
    305309        BalancethicknessVxObsEnum,
    306310        BalancethicknessVyObsEnum,

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

    r18385 r18403  
    309309                case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance";
    310310                case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
    311                 case BalancethicknessNuxEnum : return "BalancethicknessNux";
    312                 case BalancethicknessNuyEnum : return "BalancethicknessNuy";
     311                case IrrotationalH2MisfitEnum : return "IrrotationalH2Misfit";
     312                case IrrotationalDirectionMisfitEnum : return "IrrotationalDirectionMisfit";
     313                case IrrotationalVelMisfitEnum : return "IrrotationalVelMisfit";
     314                case IrrotationalAlongGradientNuEnum : return "IrrotationalAlongGradientNu";
     315                case IrrotationalAcrossGradientNuEnum : return "IrrotationalAcrossGradientNu";
     316                case BalancethicknessNuEnum : return "BalancethicknessNu";
    313317                case BalancethicknessVxObsEnum : return "BalancethicknessVxObs";
    314318                case BalancethicknessVyObsEnum : return "BalancethicknessVyObs";

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

    r18385 r18403  
    315315              else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum;
    316316              else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum;
    317               else if (strcmp(name,"BalancethicknessNux")==0) return BalancethicknessNuxEnum;
    318               else if (strcmp(name,"BalancethicknessNuy")==0) return BalancethicknessNuyEnum;
     317              else if (strcmp(name,"IrrotationalH2Misfit")==0) return IrrotationalH2MisfitEnum;
     318              else if (strcmp(name,"IrrotationalDirectionMisfit")==0) return IrrotationalDirectionMisfitEnum;
     319              else if (strcmp(name,"IrrotationalVelMisfit")==0) return IrrotationalVelMisfitEnum;
     320              else if (strcmp(name,"IrrotationalAlongGradientNu")==0) return IrrotationalAlongGradientNuEnum;
     321              else if (strcmp(name,"IrrotationalAcrossGradientNu")==0) return IrrotationalAcrossGradientNuEnum;
     322              else if (strcmp(name,"BalancethicknessNu")==0) return BalancethicknessNuEnum;
    319323              else if (strcmp(name,"BalancethicknessVxObs")==0) return BalancethicknessVxObsEnum;
    320324              else if (strcmp(name,"BalancethicknessVyObs")==0) return BalancethicknessVyObsEnum;
    … …  
    379383              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
    380384              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    381               else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    382389              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    383390              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    384391              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
     392              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
    389393              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    390394              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    … …  
    502506              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    503507              else if (strcmp(name,"Fill")==0) return FillEnum;
    504               else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    505512              else if (strcmp(name,"Friction")==0) return FrictionEnum;
    506513              else if (strcmp(name,"Internal")==0) return InternalEnum;
    507514              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
     515              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    512516              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    513517              else if (strcmp(name,"Pressure")==0) return PressureEnum;
    … …  
    625629              else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
    626630              else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
    627               else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
    628635              else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
    629636              else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
    630637              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
     638              else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    635639              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    636640              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;

  • issm/trunk-jpl/src/m/classes/balancethickness2.m

    r18169 r18403  
    88                spcpotential         = NaN;
    99                apparent_massbalance = NaN;
    10                 nux                  = NaN;
    11                 nuy                  = NaN;
     10                nu                  = NaN;
    1211                vx_obs               = NaN;
    1312                vy_obs               = NaN;
    … …  
    3231                        md = checkfield(md,'fieldname','balancethickness.spcpotential','size',[md.mesh.numberofvertices 1]);
    3332                        md = checkfield(md,'fieldname','balancethickness.apparent_massbalance','size',[md.mesh.numberofvertices 1],'NaN',1);
    34                         md = checkfield(md,'fieldname','balancethickness.nux','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1);
    35                         md = checkfield(md,'fieldname','balancethickness.nuy','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1);
     33                        md = checkfield(md,'fieldname','balancethickness.nu','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1);
    3634                        md = checkfield(md,'fieldname','balancethickness.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
    3735                        md = checkfield(md,'fieldname','balancethickness.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
    … …  
    4341                        fielddisplay(obj,'spcpotential','potential constraints (NaN means no constraint)');
    4442                        fielddisplay(obj,'apparent_massbalance','Apparent mass balance [m/yr]');
    45                         fielddisplay(obj,'nux','vx_bar = nux vx_s (in ]0 1])');
    46                         fielddisplay(obj,'nuy','vy_bar = nuy vy_s (in ]0 1])');
     43                        fielddisplay(obj,'nu','v_bar = nu v_s (in ]0 1])');
    4744                        fielddisplay(obj,'vx_obs','observed vx');
    4845                        fielddisplay(obj,'vy_obs','observed vy');
    … …  
    5350                        WriteData(fid,'object',obj,'class','Balancethickness','fieldname','spcpotential','format','DoubleMat','mattype',1);
    5451                        WriteData(fid,'object',obj,'class','Balancethickness','fieldname','apparent_massbalance','format','DoubleMat','mattype',1);
    55                         WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nux','format','DoubleMat','mattype',1);
    56                         WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nuy','format','DoubleMat','mattype',1);
    57                         WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vx_obs','format','DoubleMat','mattype',1);
    58                         WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vy_obs','format','DoubleMat','mattype',1);
     52                        WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nu','format','DoubleMat','mattype',1);
     53                        if(numel(obj.vx_obs)==md.mesh.numberofelements),
     54                                mattype=2;
     55                        else
     56                                mattype=1;
     57                        end
     58                        WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vx_obs','format','DoubleMat','mattype',mattype);
     59                        WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vy_obs','format','DoubleMat','mattype',mattype);
    5960                        WriteData(fid,'object',obj,'class','Balancethickness','fieldname','thickness_obs','format','DoubleMat','mattype',1);
    6061                end % }}}

  • issm/trunk-jpl/src/m/classes/inversion.m

    r18216 r18403  
    132132                        md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
    133133                        md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
    134                                 {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessNux' 'BalancethicknessNuy' 'BalancethicknessApparentMassbalance'});
     134                                {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness',...
     135                                'BalancethicknessNux' 'BalancethicknessNuy' 'BalancethicknessApparentMassbalance'});
    135136                        md = checkfield(md,'fieldname','inversion.nsteps','numel',1,'>=',0);
    136137                        md = checkfield(md,'fieldname','inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);

  • issm/trunk-jpl/src/m/classes/inversionvalidation.m

    r18076 r18403  
    5353                        md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
    5454                        md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
    55                                 {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessApparentMassbalance'});
    56                         md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:507]);
     55                                {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness',...
     56                                'BalancethicknessNu' 'BalancethicknessApparentMassbalance'});
     57                        md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506 601:604]);
    5758                        md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
    5859                        md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
    … …  
    139140                        pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
    140141                        pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
    141                         pos=find(obj.cost_functions==507); data(pos)=Balancethickness2MisfitEnum();
     142                        pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum();
     143                        pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum();
     144                        pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum();
     145                        pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum();
     146                        pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum();
    142147                        WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
    143148                        WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r18385 r18403  
    301301def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0]
    302302def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0]
    303 def BalancethicknessNuxEnum(): return StringToEnum("BalancethicknessNux")[0]
    304 def BalancethicknessNuyEnum(): return StringToEnum("BalancethicknessNuy")[0]
     303def IrrotationalH2MisfitEnum(): return StringToEnum("IrrotationalH2Misfit")[0]
     304def IrrotationalDirectionMisfitEnum(): return StringToEnum("IrrotationalDirectionMisfit")[0]
     305def IrrotationalVelMisfitEnum(): return StringToEnum("IrrotationalVelMisfit")[0]
     306def IrrotationalAlongGradientNuEnum(): return StringToEnum("IrrotationalAlongGradientNu")[0]
     307def IrrotationalAcrossGradientNuEnum(): return StringToEnum("IrrotationalAcrossGradientNu")[0]
     308def BalancethicknessNuEnum(): return StringToEnum("BalancethicknessNu")[0]
    305309def BalancethicknessVxObsEnum(): return StringToEnum("BalancethicknessVxObs")[0]
    306310def BalancethicknessVyObsEnum(): return StringToEnum("BalancethicknessVyObs")[0]
Note: See TracChangeset for help on using the changeset viewer.