Changeset 18476


Ignore:
Timestamp:
09/02/14 09:23:04 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removed irrotational balancethickness2, working on SIA model

Location:
issm/trunk-jpl/src
Files:
2 added
10 deleted
18 edited

Legend:

Unmodified
Added
Removed

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

    r18403 r18476  
    4545}/*}}}*/
    4646ElementVector* AdjointBalancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
     47        _error_("Not implemented yet");
    4748
    48         /*Intermediaries */
    49         int         num_responses,i;
    50         IssmDouble  hobs,hu2,weight,Jdet;
    51         IssmDouble  NUMxH2,NUMyH2,DENH2;
    52         IssmDouble  NUMxUbar,NUMyUbar,DENUbar;
    53         IssmDouble  vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
    54         IssmDouble  nu,phi,dphi[2];
    55         int        *responses = NULL;
    56         IssmDouble *xyz_list  = NULL;
    57 
    58         /*Fetch number of nodes and dof for this finite element*/
    59         int numnodes = element->GetNumberOfNodes();
    60 
    61         /*Initialize Element vector and vectors*/
    62         ElementVector* pe     = element->NewElementVector(SSAApproximationEnum);
    63         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    64         IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    65 
    66         /*Retrieve all inputs and parameters*/
    67         element->GetVerticesCoordinates(&xyz_list);
    68         element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    69         element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    70         Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
    71         Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    72         Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
    73         Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
    74         Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
    75         Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
    76 
    77         /* Start  looping on the number of gaussian points: */
    78         Gauss* gauss=element->NewGauss(2);
    79         for(int ig=gauss->begin();ig<gauss->end();ig++){
    80                 gauss->GaussPoint(ig);
    81 
    82                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    83                 element->NodalFunctions(basis,gauss);
    84                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    85 
    86                 vxobs_input->GetInputValue(&vxobs,gauss);
    87                 vyobs_input->GetInputValue(&vyobs,gauss);
    88                 nu_input->GetInputValue(&nu,gauss);
    89                 potential_input->GetInputValue(&phi,gauss);
    90                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    91                 thicknessobs_input->GetInputValue(&hobs,gauss);
    92 
    93                 vxobsbar = nu*vxobs;
    94                 vyobsbar = nu*vyobs;
    95 
    96                 vbarobs2 = (nu*nu*vxobs*vxobs + nu*nu*vyobs*vyobs);
    97                 vbarobs  = sqrt(vbarobs2);
    98                 hu2 = hobs*hobs*vbarobs2;
    99 
    100                 /*Loop over all requested responses*/
    101                 for(int resp=0;resp<num_responses;resp++){
    102                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
    103 
    104                         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;
    122                                 case Balancethickness2MisfitEnum:
    123                                         /*J = phi^2*/
    124                                         //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; //OK
    125                                         /*J = grad phi ^2*/
    126                                         //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
    127                                         /*J = (ubar - nu*uobs)^2*/
    128                                         //for(i=0;i<numnodes;i++){
    129                                         //      NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0*numnodes+i];
    130                                         //      NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1*numnodes+i];
    131                                         //      DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
    132                                         //      pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
    133                                         //}
    134                                         _error_("Not implemented yet");
    135                                         break;
    136                                 default:
    137                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    138                         }
    139                 }
    140         }
    141 
    142         /*Clean up and return*/
    143         xDelete<int>(responses);
    144         xDelete<IssmDouble>(xyz_list);
    145         xDelete<IssmDouble>(basis);
    146         xDelete<IssmDouble>(dbasis);
    147         delete gauss;
    148         return pe;
    14949}/*}}}*/
    15050void AdjointBalancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    17474                        /*Nothing, \partial J/\partial k = 0*/
    17575                        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;
    19776                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    19877        }
    … …  
    20079        /*Deal with second term*/
    20180        switch(control_type){
    202                 case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break;
    203                 case BalancethicknessNuEnum: /*Nothing: state equations do not depend on k*/ break;
    20481                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
    20582        }
    … …  
    251128        delete gauss;
    252129}/*}}}*/
    253 void AdjointBalancethickness2Analysis::GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    254 
    255         /*Intermediaries*/
    256         IssmDouble vxobs,vyobs,thicknessobs,weight;
    257         IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
    258         IssmDouble *xyz_list= NULL;
    259 
    260         /*Fetch number of vertices for this finite element*/
    261         int numvertices = element->GetNumberOfVertices();
    262 
    263         /*Initialize some vectors*/
    264         IssmDouble* basis         = xNew<IssmDouble>(numvertices);
    265         IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
    266         int*        vertexpidlist = xNew<int>(numvertices);
    267 
    268         /*Retrieve all inputs we will be needing: */
    269         element->GetVerticesCoordinates(&xyz_list);
    270         element->GradientIndexing(&vertexpidlist[0],control_index);
    271         Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    272         Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
    273         Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
    274         Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
    275         Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
    276         Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
    277 
    278         Gauss* gauss=element->NewGauss(2);
    279         for(int ig=gauss->begin();ig<gauss->end();ig++){
    280                 gauss->GaussPoint(ig);
    281 
    282                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    283                 element->NodalFunctionsP1(basis,gauss);
    284                 weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum);
    285                 vxobs_input->GetInputValue(&vxobs,gauss);
    286                 vyobs_input->GetInputValue(&vyobs,gauss);
    287                 nu_input->GetInputValue(&nu,gauss);
    288                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    289                 thicknessobs_input->GetInputValue(&thicknessobs,gauss);
    290 
    291                 dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
    292                 nuvobs2   = nu*nu*(vxobs*vxobs + vyobs*vyobs);
    293 
    294                 /*Build gradient vector (actually -dJ/da): */
    295                 for(int i=0;i<numvertices;i++){
    296                         ge[i]+= -weight*Jdet*gauss->weight*basis[i]*(
    297                                                 -2./pow(nu,3) * dphinorm2/(vxobs*vxobs + vyobs*vyobs) * (dphinorm2/nuvobs2 - thicknessobs*thicknessobs)
    298                                                 );
    299                         _assert_(!xIsNan<IssmDouble>(ge[i]));
    300                 }
    301         }
    302         gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
    303 
    304         /*Clean up and return*/
    305         xDelete<IssmDouble>(ge);
    306         xDelete<IssmDouble>(xyz_list);
    307         xDelete<IssmDouble>(basis);
    308         xDelete<int>(vertexpidlist);
    309         delete gauss;
    310 }/*}}}*/
    311 void AdjointBalancethickness2Analysis::GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    312 
    313         /*Intermediaries*/
    314         IssmDouble vxobs,vyobs,thicknessobs,weight;
    315         IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
    316         IssmDouble *xyz_list= NULL;
    317 
    318         /*Fetch number of vertices for this finite element*/
    319         int numvertices = element->GetNumberOfVertices();
    320 
    321         /*Initialize some vectors*/
    322         IssmDouble* basis         = xNew<IssmDouble>(numvertices);
    323         IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
    324         int*        vertexpidlist = xNew<int>(numvertices);
    325 
    326         /*Retrieve all inputs we will be needing: */
    327         element->GetVerticesCoordinates(&xyz_list);
    328         element->GradientIndexing(&vertexpidlist[0],control_index);
    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);
    335 
    336         Gauss* gauss=element->NewGauss(2);
    337         for(int ig=gauss->begin();ig<gauss->end();ig++){
    338                 gauss->GaussPoint(ig);
    339 
    340                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    341                 element->NodalFunctionsP1(basis,gauss);
    342                 weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum);
    343                 vxobs_input->GetInputValue(&vxobs,gauss);
    344                 vyobs_input->GetInputValue(&vyobs,gauss);
    345                 nu_input->GetInputValue(&nu,gauss);
    346                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    347                 thicknessobs_input->GetInputValue(&thicknessobs,gauss);
    348 
    349                 dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
    350                 nuvobs2   = nu*vxobs*nu*vxobs + nu*vyobs*nu*vyobs;
    351 
    352                 /*Build gradient vector (actually -dJ/da): */
    353                 for(int i=0;i<numvertices;i++){
    354                         ge[i]+= - weight*Jdet*gauss->weight*nu*(vyobs*dphi[0] -vxobs*dphi[1])*(vyobs*dphi[0] -vxobs*dphi[1])*basis[i];
    355                         _assert_(!xIsNan<IssmDouble>(ge[i]));
    356                 }
    357         }
    358         gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
    359 
    360         /*Clean up and return*/
    361         xDelete<IssmDouble>(ge);
    362         xDelete<IssmDouble>(xyz_list);
    363         xDelete<IssmDouble>(basis);
    364         xDelete<int>(vertexpidlist);
    365         delete gauss;
    366 }/*}}}*/
    367130void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    368131        element->InputUpdateFromSolutionOneDof(solution,AdjointEnum);

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

    r18403 r18476  
    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 GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index);
    32                 void GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3331                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3432                void UpdateConstraints(FemModel* femmodel);

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

    r18450 r18476  
    2626        }
    2727
    28         iomodel->FetchDataToInput(elements,BalancethicknessApparentMassbalanceEnum);
    29         iomodel->FetchDataToInput(elements,BalancethicknessNuEnum);
    30         iomodel->FetchDataToInput(elements,BalancethicknessVxObsEnum);
    31         iomodel->FetchDataToInput(elements,BalancethicknessVyObsEnum);
    32         iomodel->FetchDataToInput(elements,BalancethicknessThicknessObsEnum);
    33         iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
     28        iomodel->FetchDataToInput(elements,ThicknessEnum);
     29        iomodel->FetchDataToInput(elements,SurfaceEnum);
     30        iomodel->FetchDataToInput(elements,BaseEnum);
    3431        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     32        iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
     33        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
     34        iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
    3535}/*}}}*/
    3636void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    … …  
    4343
    4444        int finiteelement = P1Enum;
    45         IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcpotentialEnum,Balancethickness2AnalysisEnum,finiteelement);
     45        IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,Balancethickness2AnalysisEnum,finiteelement);
    4646
    4747}/*}}}*/
    … …  
    6464
    6565        /*Intermediaries */
    66         IssmDouble Jdet,D_scalar;
     66        IssmDouble  Jdet,D;
    6767        IssmDouble* xyz_list = NULL;
    6868
    … …  
    7474        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    7575
     76        /*Create input D*/
     77        this->CreateDiffusionCoefficient(element);
     78
    7679        /*Retrieve all inputs and parameters*/
    7780        element->GetVerticesCoordinates(&xyz_list);
     81        Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input);
    7882
    7983        /* Start  looping on the number of gaussian points: */
    … …  
    8387                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    8488                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     89                D_input->GetInputValue(&D,gauss);
    8590
    8691                for(int i=0;i<numnodes;i++){
    8792                        for(int j=0;j<numnodes;j++){
    88                                 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
     93                                Ke->values[i*numnodes+j] += D*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
    8994                        }
    9095                }
    … …  
    99104ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
    100105
    101         /*compute all load vectors for this element*/
    102         ElementVector* pe1=CreatePVectorVolume(element);
    103         ElementVector* pe2=CreatePVectorBoundary(element);
    104         ElementVector* pe =new ElementVector(pe1,pe2);
    105 
    106         /*clean-up and return*/
    107         delete pe1;
    108         delete pe2;
    109         return pe;
    110 }/*}}}*/
    111 ElementVector* Balancethickness2Analysis::CreatePVectorVolume(Element* element){/*{{{*/
    112 
    113106        /*Intermediaries */
    114         IssmDouble  adot,Jdet;
     107        IssmDouble  dhdt,mb,ms,dD[2],db[2],Jdet;
    115108        IssmDouble* xyz_list = NULL;
    116109
    … …  
    124117        /*Retrieve all inputs and parameters*/
    125118        element->GetVerticesCoordinates(&xyz_list);
    126         Input* adot_input   = element->GetInput(BalancethicknessApparentMassbalanceEnum); _assert_(adot_input);
     119        Input* ms_input   = element->GetInput(SurfaceforcingsMassBalanceEnum);                _assert_(ms_input);
     120        Input* mb_input   = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);       _assert_(mb_input);
     121        Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum);            _assert_(dhdt_input);
     122        Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);       _assert_(D_input);
     123        Input* bed_input  = element->GetInput(BaseEnum);                                      _assert_(bed_input);
    127124
    128125        /* Start  looping on the number of gaussian points: */
    … …  
    133130                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    134131                element->NodalFunctions(basis,gauss);
    135                 adot_input->GetInputValue(&adot,gauss);
    136 
    137                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*adot*basis[i];
     132
     133                ms_input->GetInputValue(&ms,gauss);
     134                mb_input->GetInputValue(&mb,gauss);
     135                dhdt_input->GetInputValue(&dhdt,gauss);
     136                bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss);
     137                D_input->GetInputDerivativeValue(&dD[0],xyz_list,gauss);
     138                db[0]=0.;
     139                db[1]=0.;
     140
     141                /*Since grad(b) is constant div(D grad(b) ) = grad(D).grad(b)*/
     142                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt +dD[0]*db[0]+dD[1]*db[1])*basis[i];
    138143        }
    139144
    … …  
    144149        return pe;
    145150}/*}}}*/
    146 ElementVector* Balancethickness2Analysis::CreatePVectorBoundary(Element* element){/*{{{*/
    147 
    148         /*If no front, return NULL*/
    149         if(!element->IsFaceOnBoundary()) return NULL;
    150 
    151         /*Intermediaries*/
    152         IssmDouble  Jdet,thickness,vx,vy;
    153         IssmDouble *xyz_list = NULL;
    154         IssmDouble *xyz_list_front = NULL;
    155         IssmDouble  normal[2];
    156 
    157         /*Fetch number of nodes for this finite element*/
    158         int numnodes = element->GetNumberOfNodes();
    159 
    160         /*Initialize Element vector and other vectors*/
    161         ElementVector* pe    = element->NewElementVector();
    162         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    163 
    164         /*Retrieve all inputs and parameters*/
    165         Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum); _assert_(thickness_input);
    166         Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum);        _assert_(vx_input);
    167         Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum);        _assert_(vy_input);
    168 
    169         element->GetVerticesCoordinates(&xyz_list);
    170         element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    171         element->NormalSection(&normal[0],xyz_list_front);
    172 
    173         /*Start looping on Gaussian points*/
    174         Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
    175         for(int ig=gauss->begin();ig<gauss->end();ig++){
    176 
    177                 gauss->GaussPoint(ig);
    178                 thickness_input->GetInputValue(&thickness,gauss);
    179                 vx_input->GetInputValue(&vx,gauss);
    180                 vy_input->GetInputValue(&vy,gauss);
    181                 element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    182                 element->NodalFunctions(basis,gauss);
    183 
    184                 for(int i=0;i<numnodes;i++) pe->values[i] += Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i];
    185         }
    186 
    187         /*Clean up and return*/
    188         xDelete<IssmDouble>(xyz_list);
    189         xDelete<IssmDouble>(xyz_list_front);
    190         xDelete<IssmDouble>(basis);
    191         delete gauss;
    192         return pe;
    193 }/*}}}*/
    194151void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    195152           _error_("not implemented yet");
    … …  
    199156}/*}}}*/
    200157void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    201 
    202         /*Intermediaries */
    203         int         Hinterpolation;
    204         IssmDouble  vx,vy,vbar,nu,normdphi,dphi[2],H;
    205         IssmDouble* xyz_list = NULL;
    206         int       * doflist  = NULL;
    207 
    208         /*Fetch number of nodes and dof for this finite element*/
    209         int numnodes    = element->GetNumberOfNodes();
    210         int numvertices = element->GetNumberOfVertices();
    211 
    212         /*Fetch dof list and allocate solution vector*/
    213         element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    214         IssmDouble* values         = xNew<IssmDouble>(numnodes);
    215         IssmDouble* thickness_list = xNew<IssmDouble>(numvertices);
    216         IssmDouble* vx_list        = xNew<IssmDouble>(numvertices);
    217         IssmDouble* vy_list        = xNew<IssmDouble>(numvertices);
    218 
    219         /*Use the dof list to index into the solution vector: */
    220         for(int i=0;i<numnodes;i++){
    221                 values[i]=solution[doflist[i]];
    222 
    223                 /*Check solution*/
    224                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    225         }
    226 
    227         element->AddInput(PotentialEnum,values,element->GetElementType());
    228 
    229         /*Retrieve all inputs and parameters*/
    230         element->GetVerticesCoordinates(&xyz_list);
    231         Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
    232         Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input);
    233         Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input);
    234         Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
    235         Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum);             _assert_(thickness_input);
    236 
    237         switch(element->GetElementType()){
    238                 case P1Enum: Hinterpolation = P0Enum; break;
    239                 default:     _error_("not implemented");
    240         }
    241 
    242         Gauss* gauss=element->NewGauss();
    243         for (int iv=0;iv<1;iv++){
    244                 gauss->GaussNode(Hinterpolation,iv);//P0 Only for now
    245 
    246                 vx_input->GetInputValue(&vx,gauss);
    247                 vy_input->GetInputValue(&vy,gauss);
    248                 nu_input->GetInputValue(&nu,gauss);
    249                 thickness_input->GetInputValue(&H,gauss);
    250                 potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
    251 
    252                 vx = vx*nu; vy = vy*nu;
    253                 vbar = sqrt(vx*vx + vy*vy) + 1.e-10;
    254                 normdphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
    255 
    256                 thickness_list[iv] = normdphi/vbar;
    257                 vx_list[iv]        = -1./H * dphi[0];
    258                 vy_list[iv]        = -1./H * dphi[1];
    259         }
    260         element->AddInput(ThicknessEnum,thickness_list,Hinterpolation);
    261         element->AddInput(VxEnum,vx_list,Hinterpolation);
    262         element->AddInput(VyEnum,vy_list,Hinterpolation);
    263 
    264         /*Clean up and return*/
    265         delete gauss;
    266         xDelete<int>(doflist);
    267         xDelete<IssmDouble>(xyz_list);
    268         xDelete<IssmDouble>(values);
    269         xDelete<IssmDouble>(thickness_list);
    270         xDelete<IssmDouble>(vx_list);
    271         xDelete<IssmDouble>(vy_list);
     158        element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
    272159}/*}}}*/
    273160void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    … …  
    275162        return;
    276163}/*}}}*/
     164
     165/*Specifics*/
     166void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/
     167
     168        /*Intermediaries */
     169        IssmDouble       omega,h,mu0,ds[2],Cmu,B;
     170        const int        n = 3.;
     171        const IssmDouble Hstar = 500.;
     172        const IssmDouble Lstar = 500.e+3;
     173        IssmDouble *xyz_list  = NULL;
     174
     175        /*Fetch number of vertices and allocate output*/
     176        int  numnodes = element->GetNumberOfNodes();
     177        IssmDouble* D      = xNew<IssmDouble>(numnodes);
     178        IssmDouble* Dgradb = xNew<IssmDouble>(numnodes);
     179
     180        /*retrieve what we need: */
     181        element->GetVerticesCoordinates(&xyz_list);
     182        Input* thickness_input  = element->GetInput(ThicknessEnum);          _assert_(thickness_input);
     183        Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
     184        Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
     185        IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
     186
     187        /*Calculate damage evolution source term: */
     188        Gauss* gauss=element->NewGauss();
     189        for (int i=0;i<numnodes;i++){
     190                gauss->GaussNode(element->GetElementType(),i);
     191               
     192                B_input->GetInputValue(&B,gauss);
     193                thickness_input->GetInputValue(&h,gauss);
     194                surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
     195
     196                mu0   = pow(2.,(1.-3*n)/(2.*n)) * B;
     197                omega = pow(rhog,3) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
     198                Cmu   = 0.;
     199
     200                D[i] = omega*(ds[0]*ds[0]+ds[1]*ds[1])*pow(h,4)*(1./5.*h - Cmu);
     201        }
     202
     203        /*Add input*/
     204        element->AddInput(BalancethicknessDiffusionCoefficientEnum,D,element->GetElementType());
     205       
     206        /*Clean up and return*/
     207        xDelete<IssmDouble>(D);
     208        delete gauss;
     209}/*}}}*/

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

    r18057 r18476  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
    28                 ElementVector* CreatePVectorVolume(Element* element);
    29                 ElementVector* CreatePVectorBoundary(Element* element);
    3028                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3129                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    3230                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3331                void UpdateConstraints(FemModel* femmodel);
     32
     33                /*Specifics*/
     34                void CreateDiffusionCoefficient(Element* element);
    3435};
    3536#endif

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

    r18470 r18476  
    10531053                                name==BedEnum ||
    10541054                                name==BalancethicknessThickeningRateEnum ||
    1055                                 name==BalancethicknessApparentMassbalanceEnum ||
    1056                                 name==BalancethicknessNuEnum ||
    10571055                                name==SigmaNNEnum ||
    10581056                                name==SurfaceSlopeXEnum ||

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

    r18474 r18476  
    10321032                                case ThicknessEnum:
    10331033                                case FrictionCoefficientEnum:
    1034                                 case BalancethicknessNuEnum:
    10351034                                        if(iomodel->Data(control)){
    10361035                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(control)[tria_vertex_ids[j]-1];

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

    r18452 r18476  
    593593                                case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
    594594                                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;
    600595
    601596                           /*Vector */
    … …  
    13201315                /*If on water, return 0: */
    13211316                if(!element->IsIceInElement()) continue;
    1322 
    1323                 /* Get node coordinates*/
    1324                 element->GetVerticesCoordinates(&xyz_list);
    1325                 Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    1326                 Input* thickness_input   =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    1327                 Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
    1328                 Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
    1329                 Input* vxobs_input       =element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
    1330                 Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
    1331                 Input* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
    1332                 Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_input);
    1333                 Input* nu_input       = element->GetInput(BalancethicknessNuEnum);   _assert_(nu_input);
    1334 
    1335                 /* Start  looping on the number of gaussian points: */
    1336                 Gauss* gauss=element->NewGauss(2);
    1337                 for(int ig=gauss->begin();ig<gauss->end();ig++){
    1338 
    1339                         gauss->GaussPoint(ig);
    1340 
    1341                         /* Get Jacobian determinant: */
    1342                         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1343 
    1344                         /*Get all parameters at gaussian point*/
    1345                         weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
    1346                         thickness_input->GetInputValue(&thickness,gauss);
    1347                         thicknessobs_input->GetInputValue(&thicknessobs,gauss);
    1348                         potential_input->GetInputValue(&potential,gauss);
    1349                         potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss);
    1350                         vxobs_input->GetInputValue(&vxobs,gauss);
    1351                         vyobs_input->GetInputValue(&vyobs,gauss);
    1352                         vx_input->GetInputValue(&vxbar,gauss);
    1353                         vy_input->GetInputValue(&vybar,gauss);
    1354                         nu_input->GetInputValue(&nu,gauss);
    1355 
    1356                         vxbarobs = nu*vxobs;
    1357                         vybarobs = nu*vyobs;
    1358 
    1359                         /*J = (H^2 - Hobs^2)^2*/
    1360                         J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
    1361                         /*J = phi^2*/
    1362                         //J +=.5*potential*potential*weight*Jdet*gauss->weight;// OK
    1363                         /*J = grad phi^2*/
    1364                         //J +=.5*(dpotential[0]*dpotential[0] + dpotential[1]*dpotential[1])*weight*Jdet*gauss->weight;
    1365                         /*J = (ubar - nux*uobs)^2*/
    1366                         //J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
    1367                         /*J = 1/2 (vbar ^ gard(phi))^2*/
    1368                         //J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
    1369                 }
    1370 
    1371                 /*clean up and Return: */
    1372                 xDelete<IssmDouble>(xyz_list);
    1373                 delete gauss;
     1317                 _error_("Not implemented");
     1318
    13741319        }
    13751320
    … …  
    13821327        *presponse=J;
    13831328
    1384 }/*}}}*/
    1385 void 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 }/*}}}*/
    1438 void 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 }/*}}}*/
    1496 void FemModel::IrrotationalVelMisfitx(IssmDouble* presponse){/*{{{*/
    1497         _error_("not implemented yet");
    1498 }/*}}}*/
    1499 void FemModel::IrrotationalAlongGradientNux(IssmDouble* presponse){/*{{{*/
    1500         _error_("not implemented yet");
    1501 }/*}}}*/
    1502 void FemModel::IrrotationalAcrossGradientNux(IssmDouble* presponse){/*{{{*/
    1503         _error_("not implemented yet");
    15041329}/*}}}*/
    15051330void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/

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

    r18403 r18476  
    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);
    8883                #ifdef  _HAVE_DAKOTA_
    8984                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

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

    r17850 r18476  
    2626        if(save_results){
    2727                if(VerboseSolution()) _printf0_("   saving results\n");
    28                 int outputs[4] = {ThicknessEnum,PotentialEnum,VxEnum,VyEnum};
    29                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
     28                int outputs[1] = {ThicknessEnum};
     29                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
    3030        }
    3131

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

    r18403 r18476  
    4141                        case ThicknessEnum:
    4242                        case FrictionCoefficientEnum:
    43                         case BalancethicknessNuEnum:
    4443                        case BalancethicknessApparentMassbalanceEnum:
    4544                                iomodel->FetchData(1,control);

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

    r18474 r18476  
    301301        BalancethicknessApparentMassbalanceEnum,
    302302        Balancethickness2MisfitEnum,
    303         IrrotationalH2MisfitEnum,
    304         IrrotationalDirectionMisfitEnum,
    305         IrrotationalVelMisfitEnum,
    306         IrrotationalAlongGradientNuEnum,
    307         IrrotationalAcrossGradientNuEnum,
    308         BalancethicknessNuEnum,
    309         BalancethicknessVxObsEnum,
    310         BalancethicknessVyObsEnum,
    311         BalancethicknessThicknessObsEnum,
     303        BalancethicknessDiffusionCoefficientEnum,
    312304        /*}}}*/
    313305        /*Surfaceforcings{{{*/

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

    r18474 r18476  
    309309                case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance";
    310310                case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
    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";
    317                 case BalancethicknessVxObsEnum : return "BalancethicknessVxObs";
    318                 case BalancethicknessVyObsEnum : return "BalancethicknessVyObs";
    319                 case BalancethicknessThicknessObsEnum : return "BalancethicknessThicknessObs";
     311                case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient";
    320312                case SurfaceforcingsEnum : return "Surfaceforcings";
    321313                case SMBEnum : return "SMB";

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

    r18474 r18476  
    315315              else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum;
    316316              else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum;
    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;
    323               else if (strcmp(name,"BalancethicknessVxObs")==0) return BalancethicknessVxObsEnum;
    324               else if (strcmp(name,"BalancethicknessVyObs")==0) return BalancethicknessVyObsEnum;
    325               else if (strcmp(name,"BalancethicknessThicknessObs")==0) return BalancethicknessThicknessObsEnum;
     317              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    326318              else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum;
    327319              else if (strcmp(name,"SMB")==0) return SMBEnum;
    … …  
    383375              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
    384376              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
     377              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    389378              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    390379              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    … …  
    394383              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    395384              else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
    396               else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    397389              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    398390              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    … …  
    506498              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    507499              else if (strcmp(name,"Fill")==0) return FillEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
     500              else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    512501              else if (strcmp(name,"Friction")==0) return FrictionEnum;
    513502              else if (strcmp(name,"Internal")==0) return InternalEnum;
    … …  
    517506              else if (strcmp(name,"Pressure")==0) return PressureEnum;
    518507              else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
    519               else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
    520512              else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
    521513              else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
    … …  
    629621              else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
    630622              else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
     623              else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
    635624              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
    636625              else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
    … …  
    640629              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    641630              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    642               else if (strcmp(name,"MinVx")==0) return MinVxEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"MinVx")==0) return MinVxEnum;
    643635              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    644636              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;

  • issm/trunk-jpl/src/dox/issm.dox

    r18053 r18476  
    4747</th>
    4848<tr>
    49 <th  bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td  bgcolor=#FFFFFF style="text-align:right;">453</td><td  bgcolor=#FFFFFF style="text-align:right;">15861</td><td  bgcolor=#FFFFFF style="text-align:right;">18188</td><td  bgcolor=#FFFFFF style="text-align:right;">68849</td><td  bgcolor=#FFFFFF style="text-align:right;">102898</td>
     49<th  bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td  bgcolor=#FFFFFF style="text-align:right;">445</td><td  bgcolor=#FFFFFF style="text-align:right;">15967</td><td  bgcolor=#FFFFFF style="text-align:right;">16145</td><td  bgcolor=#FFFFFF style="text-align:right;">70397</td><td  bgcolor=#FFFFFF style="text-align:right;">102509</td>
    5050</tr>
    5151<tr>
    52 <th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">1378</td><td  bgcolor=#C6E2FF style="text-align:right;">8341</td><td  bgcolor=#C6E2FF style="text-align:right;">16259</td><td  bgcolor=#C6E2FF style="text-align:right;">38645</td><td  bgcolor=#C6E2FF style="text-align:right;">63245</td>
     52<th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">1422</td><td  bgcolor=#C6E2FF style="text-align:right;">8497</td><td  bgcolor=#C6E2FF style="text-align:right;">16638</td><td  bgcolor=#C6E2FF style="text-align:right;">39538</td><td  bgcolor=#C6E2FF style="text-align:right;">64673</td>
    5353</tr>
    5454<tr>
    55 <th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">421</td><td  bgcolor=#FFFFFF style="text-align:right;">3504</td><td  bgcolor=#FFFFFF style="text-align:right;">3753</td><td  bgcolor=#FFFFFF style="text-align:right;">15489</td><td  bgcolor=#FFFFFF style="text-align:right;">22746</td>
     55<th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">411</td><td  bgcolor=#FFFFFF style="text-align:right;">3443</td><td  bgcolor=#FFFFFF style="text-align:right;">3528</td><td  bgcolor=#FFFFFF style="text-align:right;">15243</td><td  bgcolor=#FFFFFF style="text-align:right;">22214</td>
    5656</tr>
    5757<tr>
    58 <th  bgcolor=#C6E2FF style="text-align:left;"> m4 </th><td  bgcolor=#C6E2FF style="text-align:right;">9</td><td  bgcolor=#C6E2FF style="text-align:right;">1588</td><td  bgcolor=#C6E2FF style="text-align:right;">151</td><td  bgcolor=#C6E2FF style="text-align:right;">11565</td><td  bgcolor=#C6E2FF style="text-align:right;">13304</td>
     58<th  bgcolor=#C6E2FF style="text-align:left;"> m4 </th><td  bgcolor=#C6E2FF style="text-align:right;">8</td><td  bgcolor=#C6E2FF style="text-align:right;">1036</td><td  bgcolor=#C6E2FF style="text-align:right;">149</td><td  bgcolor=#C6E2FF style="text-align:right;">9756</td><td  bgcolor=#C6E2FF style="text-align:right;">10941</td>
    5959</tr>
    6060<tr>
    61 <th  bgcolor=#FFFFFF style="text-align:left;"> Python </th><td  bgcolor=#FFFFFF style="text-align:right;">144</td><td  bgcolor=#FFFFFF style="text-align:right;">2330</td><td  bgcolor=#FFFFFF style="text-align:right;">2559</td><td  bgcolor=#FFFFFF style="text-align:right;">9348</td><td  bgcolor=#FFFFFF style="text-align:right;">14237</td>
    62 </tr>
    63 <tr>
    64 <th  bgcolor=#C6E2FF style="text-align:left;"> XML </th><td  bgcolor=#C6E2FF style="text-align:right;">4</td><td  bgcolor=#C6E2FF style="text-align:right;">160</td><td  bgcolor=#C6E2FF style="text-align:right;">92</td><td  bgcolor=#C6E2FF style="text-align:right;">3075</td><td  bgcolor=#C6E2FF style="text-align:right;">3327</td>
    65 </tr>
    66 <tr>
    67 <th  bgcolor=#FFFFFF style="text-align:left;"> Java </th><td  bgcolor=#FFFFFF style="text-align:right;">18</td><td  bgcolor=#FFFFFF style="text-align:right;">719</td><td  bgcolor=#FFFFFF style="text-align:right;">891</td><td  bgcolor=#FFFFFF style="text-align:right;">2321</td><td  bgcolor=#FFFFFF style="text-align:right;">3931</td>
     61<th  bgcolor=#FFFFFF style="text-align:left;"> Python </th><td  bgcolor=#FFFFFF style="text-align:right;">151</td><td  bgcolor=#FFFFFF style="text-align:right;">2421</td><td  bgcolor=#FFFFFF style="text-align:right;">2623</td><td  bgcolor=#FFFFFF style="text-align:right;">9733</td><td  bgcolor=#FFFFFF style="text-align:right;">14777</td>
    6862</tr>
    6963<tr>
    … …  
    7165</tr>
    7266<tr>
    73 <th  bgcolor=#FFFFFF style="text-align:left;"> Bourne  Shell </th><td  bgcolor=#FFFFFF style="text-align:right;">3</td><td  bgcolor=#FFFFFF style="text-align:right;">61</td><td  bgcolor=#FFFFFF style="text-align:right;">88</td><td  bgcolor=#FFFFFF style="text-align:right;">266</td><td  bgcolor=#FFFFFF style="text-align:right;">415</td>
     67<th  bgcolor=#FFFFFF style="text-align:left;"> Bourne  Shell </th><td  bgcolor=#FFFFFF style="text-align:right;">2</td><td  bgcolor=#FFFFFF style="text-align:right;">59</td><td  bgcolor=#FFFFFF style="text-align:right;">84</td><td  bgcolor=#FFFFFF style="text-align:right;">262</td><td  bgcolor=#FFFFFF style="text-align:right;">405</td>
    7468</tr>
    7569<tr>
    76 <th  bgcolor=#C6E2FF style="text-align:left;"> XSD </th><td  bgcolor=#C6E2FF style="text-align:right;">2</td><td  bgcolor=#C6E2FF style="text-align:right;">30</td><td  bgcolor=#C6E2FF style="text-align:right;">36</td><td  bgcolor=#C6E2FF style="text-align:right;">229</td><td  bgcolor=#C6E2FF style="text-align:right;">295</td>
    77 </tr>
    78 <tr>
    79 <th  bgcolor=#FFFFFF style="text-align:left;"> Perl </th><td  bgcolor=#FFFFFF style="text-align:right;">1</td><td  bgcolor=#FFFFFF style="text-align:right;">6</td><td  bgcolor=#FFFFFF style="text-align:right;">9</td><td  bgcolor=#FFFFFF style="text-align:right;">196</td><td  bgcolor=#FFFFFF style="text-align:right;">211</td>
    80 </tr>
    81 <tr>
    82 <th  bgcolor=#C6E2FF style="text-align:left;"> Ant </th><td  bgcolor=#C6E2FF style="text-align:right;">1</td><td  bgcolor=#C6E2FF style="text-align:right;">16</td><td  bgcolor=#C6E2FF style="text-align:right;">7</td><td  bgcolor=#C6E2FF style="text-align:right;">103</td><td  bgcolor=#C6E2FF style="text-align:right;">126</td>
    83 </tr>
    84 <tr>
    85 <th  bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td  bgcolor=#FFFFFF style="text-align:right;">2441</td><td  bgcolor=#FFFFFF style="text-align:right;">32620</td><td  bgcolor=#FFFFFF style="text-align:right;">42335</td><td  bgcolor=#FFFFFF style="text-align:right;">150451</td><td  bgcolor=#FFFFFF style="text-align:right;">225406</td>
     70<th  bgcolor=#C6E2FF style="text-align:left;"> SUM: </th><td  bgcolor=#C6E2FF style="text-align:right;">2446</td><td  bgcolor=#C6E2FF style="text-align:right;">31427</td><td  bgcolor=#C6E2FF style="text-align:right;">39469</td><td  bgcolor=#C6E2FF style="text-align:right;">145294</td><td  bgcolor=#C6E2FF style="text-align:right;">216190</td>
    8671</tr>
    8772</table>

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

    r18405 r18476  
    239239                        pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
    240240                        pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
    241                         pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum();
    242                         pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum();
    243                         pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum();
    244                         pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum();
    245                         pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum();
    246241                        WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
    247242                        WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');

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

    r18403 r18476  
    140140                        pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
    141141                        pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
    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();
    147142                        WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
    148143                        WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');

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

    r18405 r18476  
    170170                        pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
    171171                        pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
    172                         pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum();
    173                         pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum();
    174                         pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum();
    175                         pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum();
    176                         pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum();
    177172                        WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
    178173                        WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');

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

    r18403 r18476  
    301301def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0]
    302302def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0]
    303 def IrrotationalH2MisfitEnum(): return StringToEnum("IrrotationalH2Misfit")[0]
    304 def IrrotationalDirectionMisfitEnum(): return StringToEnum("IrrotationalDirectionMisfit")[0]
    305 def IrrotationalVelMisfitEnum(): return StringToEnum("IrrotationalVelMisfit")[0]
    306 def IrrotationalAlongGradientNuEnum(): return StringToEnum("IrrotationalAlongGradientNu")[0]
    307 def IrrotationalAcrossGradientNuEnum(): return StringToEnum("IrrotationalAcrossGradientNu")[0]
    308 def BalancethicknessNuEnum(): return StringToEnum("BalancethicknessNu")[0]
    309 def BalancethicknessVxObsEnum(): return StringToEnum("BalancethicknessVxObs")[0]
    310 def BalancethicknessVyObsEnum(): return StringToEnum("BalancethicknessVyObs")[0]
    311 def BalancethicknessThicknessObsEnum(): return StringToEnum("BalancethicknessThicknessObs")[0]
     303def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0]
    312304def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0]
    313305def SMBEnum(): return StringToEnum("SMB")[0]
    … …  
    590582def OneLayerP4zEnum(): return StringToEnum("OneLayerP4z")[0]
    591583def CrouzeixRaviartEnum(): return StringToEnum("CrouzeixRaviart")[0]
     584def LACrouzeixRaviartEnum(): return StringToEnum("LACrouzeixRaviart")[0]
    592585def SaveResultsEnum(): return StringToEnum("SaveResults")[0]
    593586def BoolExternalResultEnum(): return StringToEnum("BoolExternalResult")[0]
Note: See TracChangeset for help on using the changeset viewer.