source: issm/oecreview/Archive/18296-19100/ISSM-18757-18758.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 28.4 KB

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

     
    1919}
    2020/*}}}*/
    2121void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    22         int  finiteelement;
    2322
    2423        /*Finite element type*/
    25         finiteelement = P1Enum;
     24        int finiteelement = P1Enum;
    2625
    2726        /*Update elements: */
    2827        int counter=0;
    … …  
    3736        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    3837        iomodel->FetchDataToInput(elements,VxEnum);
    3938        iomodel->FetchDataToInput(elements,VyEnum);
    40         iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
     39
     40        /*Get calving parameters*/
     41        int calvinglaw;
     42        iomodel->Constant(&calvinglaw,CalvingLawEnum);
     43        switch(calvinglaw){
     44                case DefaultCalvingEnum:
     45                        iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
     46                        break;
     47                case CalvingLevermannEnum:
     48                        iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
     49                        break;
     50                default:
     51                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     52        }
    4153}
    4254/*}}}*/
    4355void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    … …  
    99111        Element* basalelement = element->SpawnBasalElement();
    100112
    101113        /*Intermediaries */
    102         int  dim, domaintype;
    103         bool iscalvingrate;
     114        int  stabilization=2;
     115        int  dim, domaintype, calvinglaw;
     116        bool iscalving;
    104117        int i, row, col;
    105118        IssmDouble kappa;
    106119        IssmDouble Jdet, dt, D_scalar;
    107120        IssmDouble h,hx,hy,hz;
    108121        IssmDouble vel;
    109 //      IssmDouble norm_dlsf;
     122        IssmDouble norm_dlsf, calvingrate;
    110123        IssmDouble* xyz_list = NULL;
    111124
    112125        /*Get problem dimension and whether there is calving or not*/
    113         basalelement->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
     126        basalelement->FindParam(&iscalving,TransientIscalvingEnum);
    114127        basalelement->FindParam(&domaintype,DomainTypeEnum);
     128        basalelement->FindParam(&calvinglaw,CalvingLawEnum);
    115129        switch(domaintype){
    116130                case Domain2DverticalEnum:   dim = 1; break;
    117131                case Domain2DhorizontalEnum: dim = 2; break;
    … …  
    130144        IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
    131145        IssmDouble*    v        = xNew<IssmDouble>(dim);
    132146        IssmDouble*    w        = xNew<IssmDouble>(dim);
    133         IssmDouble*    c        = xNew<IssmDouble>(dim);
    134         //IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
     147        IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
     148        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
    135149
    136150        /*Retrieve all inputs and parameters*/
    137151        basalelement->GetVerticesCoordinates(&xyz_list);
    138152        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    139         Input* vx_input=NULL;
    140         Input* vy_input=NULL;
    141         Input* calvingratex_input=NULL;
    142         Input* calvingratey_input=NULL;
     153        Input* vx_input           = NULL;
     154        Input* vy_input           = NULL;
     155        Input* calvingratex_input = NULL;
     156        Input* calvingratey_input = NULL;
     157        Input* lsf_slopex_input   = NULL;
     158        Input* lsf_slopey_input   = NULL;
     159        Input* calvingrate_input  = NULL;
     160
     161        /*Load velocities*/
    143162        if(domaintype==Domain2DhorizontalEnum){
    144163                vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    145164                vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
    146                 if(iscalvingrate){
    147                         calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    148                         calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    149                 }
    150165        }
    151166        else{
    152167                if(dim==1){
    153168                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    154                         if(iscalvingrate){
    155                                 calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    156                         }
    157169                }
    158170                if(dim==2){
    159171                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
    160172                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
    161                         if(iscalvingrate){
    162                                 calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    163                                 calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    164                         }
    165173                }
    166174        }
    167175
    168 //      Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    169 //      Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    170         //Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
    171        
     176        /*Load calving inputs*/
     177        if(iscalving){
     178                switch(calvinglaw){
     179                        case DefaultCalvingEnum:
     180                                lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     181                                if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     182                                calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     183                                break;
     184                        case CalvingLevermannEnum:
     185                                if(domaintype==Domain2DhorizontalEnum){
     186                                        calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     187                                        calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     188                                }
     189                                else{
     190                                        if(dim==1){
     191                                                calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     192                                        }
     193                                        if(dim==2){
     194                                                calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     195                                                calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     196                                        }
     197                                }
     198                                break;
     199                        default:
     200                                _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     201                }
     202        }
     203
    172204        /* Start  looping on the number of gaussian points: */
    173205        Gauss* gauss=basalelement->NewGauss(2);
    174206        for(int ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    192224                GetBprime(Bprime,basalelement,xyz_list,gauss);
    193225                vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity
    194226                vy_input->GetInputValue(&v[1],gauss);
    195                 if(iscalvingrate){
    196                         calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
    197                         calvingratey_input->GetInputValue(&c[1],gauss);
    198                         for(i=0;i<dim;i++) w[i]=v[i]-c[i];
     227
     228                /*Get calving speed*/
     229                if(iscalving){
     230                        switch(calvinglaw){
     231                                case DefaultCalvingEnum:
     232                                        lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     233                                        if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     234                                        calvingrate_input->GetInputValue(&calvingrate,gauss);
     235
     236                                        norm_dlsf=0.;
     237                                        for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     238                                        norm_dlsf=sqrt(norm_dlsf);
     239
     240                                        if(norm_dlsf>1.e-10)
     241                                         for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
     242                                        else
     243                                         for(i=0;i<dim;i++) c[i]=0.;
     244                                        break;
     245                                case CalvingLevermannEnum:
     246                                        calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
     247                                        if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
     248                                        break;
     249                                default:
     250                                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     251                        }
    199252                }
    200                 else{
    201                         for(i=0;i<dim;i++) w[i]=v[i];
    202                 }
    203                 //lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    204                 //lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    205                 //calvingrate_input->GetInputValue(&calvingrate,gauss);
    206253
    207                 //norm_dlsf=0.;
    208                 //for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    209                 //norm_dlsf=sqrt(norm_dlsf);
     254                /*Levelset speed is ice velocity - calving rate*/
     255                for(i=0;i<dim;i++) w[i]=v[i]-c[i];
    210256
    211                 //if(norm_dlsf>1.e-10)
    212                 //      for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
    213                 //else
    214                 //for(i=0;i<dim;i++) c[i]=0.;
    215                
    216 
    217                 for(row=0;row<dim;row++)
    218                         for(col=0;col<dim;col++)
     257                /*Compute D*/
     258                for(row=0;row<dim;row++){
     259                        for(col=0;col<dim;col++){
    219260                                if(row==col)
    220                                         D[row*dim+col]=D_scalar*w[row];
     261                                 D[row*dim+col]=D_scalar*w[row];
    221262                                else
    222                                    D[row*dim+col]=0.;
     263                                 D[row*dim+col]=0.;
     264                        }
     265                }
    223266
    224267                TripleMultiply(B,dim,numnodes,1,
    225268                                        D,dim,dim,0,
    … …  
    227270                                        &Ke->values[0],1);
    228271
    229272                /* Stabilization */
    230                 int stabilization=2;
    231273                vel=0.;
    232                 for(i=0;i<dim;i++) vel+=w[i]*w[i];
     274                for(i=0;i<dim;i++) vel+=v[i]*v[i];
    233275                vel=sqrt(vel)+1.e-14;
    234276                switch(stabilization){
    235277                        case 0:
    … …  
    255297                        case 2:
    256298                                /* Streamline Upwinding */
    257299                                basalelement->ElementSizes(&hx,&hy,&hz);
    258                                 h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
     300                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
    259301                                for(row=0;row<dim;row++)
    260302                                        for(col=0;col<dim;col++)
    261                                                 D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
     303                                                D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
    262304
    263305                                TripleMultiply(Bprime,dim,numnodes,1,
    264306                                                        D,dim,dim,0,
    … …  
    279321        xDelete<IssmDouble>(v);
    280322        xDelete<IssmDouble>(w);
    281323        xDelete<IssmDouble>(c);
    282         //xDelete<IssmDouble>(dlsf);
     324        xDelete<IssmDouble>(dlsf);
    283325        delete gauss;
    284326        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    285327        return Ke;

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

     
    208208        DamageEnum,
    209209        NewDamageEnum,
    210210        StressIntensityFactorEnum,
     211        CalvingLawEnum,
     212        CalvingCalvingrateEnum,
     213        CalvingLevermannEnum,
     214        DefaultCalvingEnum,
     215        CalvingRequestedOutputsEnum,
     216        CalvinglevermannCoeffEnum,
    211217        CalvingratexEnum,
    212218        CalvingrateyEnum,
    213219        CalvingratexAverageEnum,
    … …  
    249255        Domain3DEnum,
    250256        MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
    251257        MasstransportHydrostaticAdjustmentEnum,
    252         MasstransportIscalvingrateEnum,
    253         MasstransportLevermannCalvingCoeffEnum,
    254258        MasstransportIsfreesurfaceEnum,
    255259        MasstransportMinThicknessEnum,
    256260        MasstransportPenaltyFactorEnum,
    257261        MasstransportSpcthicknessEnum,
    258         MasstransportCalvingrateEnum,
    259262        MasstransportStabilizationEnum,
    260263        MasstransportVertexPairingEnum,
    261264        MasstransportNumRequestedOutputsEnum,
    … …  
    313316        TransientIsgiaEnum,
    314317        TransientIsdamageevolutionEnum,
    315318        TransientIshydrologyEnum,
     319        TransientIscalvingEnum,
    316320        TransientNumRequestedOutputsEnum,
    317321        TransientRequestedOutputsEnum,
    318322        PotentialEnum,

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

     
    216216                case DamageEnum : return "Damage";
    217217                case NewDamageEnum : return "NewDamage";
    218218                case StressIntensityFactorEnum : return "StressIntensityFactor";
     219                case CalvingLawEnum : return "CalvingLaw";
     220                case CalvingCalvingrateEnum : return "CalvingCalvingrate";
     221                case CalvingLevermannEnum : return "CalvingLevermann";
     222                case DefaultCalvingEnum : return "DefaultCalving";
     223                case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs";
     224                case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
    219225                case CalvingratexEnum : return "Calvingratex";
    220226                case CalvingrateyEnum : return "Calvingratey";
    221227                case CalvingratexAverageEnum : return "CalvingratexAverage";
    … …  
    257263                case Domain3DEnum : return "Domain3D";
    258264                case MiscellaneousNameEnum : return "MiscellaneousName";
    259265                case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";
    260                 case MasstransportIscalvingrateEnum : return "MasstransportIscalvingrate";
    261                 case MasstransportLevermannCalvingCoeffEnum : return "MasstransportLevermannCalvingCoeff";
    262266                case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface";
    263267                case MasstransportMinThicknessEnum : return "MasstransportMinThickness";
    264268                case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor";
    265269                case MasstransportSpcthicknessEnum : return "MasstransportSpcthickness";
    266                 case MasstransportCalvingrateEnum : return "MasstransportCalvingrate";
    267270                case MasstransportStabilizationEnum : return "MasstransportStabilization";
    268271                case MasstransportVertexPairingEnum : return "MasstransportVertexPairing";
    269272                case MasstransportNumRequestedOutputsEnum : return "MasstransportNumRequestedOutputs";
    … …  
    321324                case TransientIsgiaEnum : return "TransientIsgia";
    322325                case TransientIsdamageevolutionEnum : return "TransientIsdamageevolution";
    323326                case TransientIshydrologyEnum : return "TransientIshydrology";
     327                case TransientIscalvingEnum : return "TransientIscalving";
    324328                case TransientNumRequestedOutputsEnum : return "TransientNumRequestedOutputs";
    325329                case TransientRequestedOutputsEnum : return "TransientRequestedOutputs";
    326330                case PotentialEnum : return "Potential";

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

     
    219219              else if (strcmp(name,"Damage")==0) return DamageEnum;
    220220              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    221221              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
     222              else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
     223              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
     224              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
     225              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
     226              else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum;
     227              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    222228              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    223229              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
    224230              else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
    … …  
    253259              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    254260              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    255261              else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
    256               else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
     262         else stage=3;
     263   }
     264   if(stage==3){
     265              if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
    257266              else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
    258267              else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
    259268              else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum;
    260269              else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
    261270              else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
    266               else if (strcmp(name,"MasstransportIscalvingrate")==0) return MasstransportIscalvingrateEnum;
    267               else if (strcmp(name,"MasstransportLevermannCalvingCoeff")==0) return MasstransportLevermannCalvingCoeffEnum;
     271              else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
    268272              else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
    269273              else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
    270274              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
    271275              else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
    272               else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
    273276              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    274277              else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
    275278              else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
    … …  
    327330              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    328331              else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
    329332              else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
     333              else if (strcmp(name,"TransientIscalving")==0) return TransientIscalvingEnum;
    330334              else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
    331335              else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
    332336              else if (strcmp(name,"Potential")==0) return PotentialEnum;
    … …  
    378382              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    379383              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    380384              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    381               else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
    382389              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    383390              else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
    384391              else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
     392              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
    389393              else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
    390394              else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
    391395              else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
    … …  
    501505              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    502506              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    503507              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    504               else if (strcmp(name,"Seg")==0) return SegEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"Seg")==0) return SegEnum;
    505512              else if (strcmp(name,"SegInput")==0) return SegInputEnum;
    506513              else if (strcmp(name,"Tria")==0) return TriaEnum;
    507514              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"Tetra")==0) return TetraEnum;
     515              else if (strcmp(name,"Tetra")==0) return TetraEnum;
    512516              else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
    513517              else if (strcmp(name,"Penta")==0) return PentaEnum;
    514518              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    … …  
    624628              else if (strcmp(name,"P2")==0) return P2Enum;
    625629              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    626630              else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
    627               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    628635              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    629636              else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
    630637              else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"P1P1")==0) return P1P1Enum;
     638              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    635639              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    636640              else if (strcmp(name,"MINI")==0) return MINIEnum;
    637641              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    … …  
    747751              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    748752              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
    749753              else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
    750               else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
    751758              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    752759              else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
    753760              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
     761              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
    758762              else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
    759763              else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
    760764              else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;

  • ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

     
    6666        parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
    6767        parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
    6868        parameters->AddObject(iomodel->CopyConstantObject(InversionTypeEnum));
     69        parameters->AddObject(iomodel->CopyConstantObject(CalvingLawEnum));
    6970        if(solution_type==SeaiceSolutionEnum){
    7071        }
    7172        else{
    … …  
    8283                parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
    8384                parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
    8485                parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
    85                 parameters->AddObject(iomodel->CopyConstantObject(MasstransportIscalvingrateEnum));
    86                 parameters->AddObject(iomodel->CopyConstantObject(MasstransportLevermannCalvingCoeffEnum));
     86                parameters->AddObject(iomodel->CopyConstantObject(TransientIscalvingEnum));
    8787                parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
    8888                parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
    8989

  • ../trunk-jpl/src/c/cores/transient_core.cpp

     
    2020
    2121        /*parameters: */
    2222        IssmDouble starttime,finaltime,dt,yts;
    23         bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalvingrate;
     23        bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving;
    2424        bool       save_results,dakota_analysis;
    2525        bool       time_adapt=false;
    2626        int        output_frequency;
    27         int        domaintype,groundingline_migration;
     27        int        domaintype,groundingline_migration,calvinglaw;
    2828        int        numoutputs         = 0;
    2929        Analysis  *analysis = NULL;
    3030        char**     requested_outputs = NULL;
    … …  
    5252        femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
    5353        femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
    5454        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
    55         femmodel->parameters->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
     55        femmodel->parameters->FindParam(&iscalving,TransientIscalvingEnum);
     56        femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
    5657        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    5758        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    5859        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
    … …  
    101102                }
    102103
    103104                if(islevelset){
    104                         if(iscalvingrate){
     105                        if(iscalving && calvinglaw==CalvingLevermannEnum){
    105106                                if(VerboseSolution()) _printf0_("   computing calving rate\n");
    106107                                femmodel->StrainRateparallelx();
    107108                                femmodel->StrainRateperpendicularx();

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

     
    438438        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
    439439        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
    440440        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);                                 _assert_(strainperpendicular_input);
    441         this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
     441        Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
    442442
    443443        /* Start looping on the number of vertices: */
    444444        gauss=new GaussTria();
    … …  
    451451                vel=vx*vx+vy*vy;
    452452                strainparallel_input->GetInputValue(&strainparallel,gauss);
    453453                strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
     454                levermanncoeff_input->GetInputValue(&propcoeff,gauss);
    454455
    455456                /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
    456457                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
    … …  
    464465        /*Add input*/
    465466        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
    466467        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
    467         this->inputs->AddInput(new TriaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
     468        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
    468469
    469470        /*Clean up and return*/
    470471        delete gauss;

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

     
    366366        Input* vx_input=inputs->GetInput(VxEnum);                                                                                                                                               _assert_(vx_input);
    367367        Input* vy_input=inputs->GetInput(VyEnum);                                                                                                                                               _assert_(vy_input);
    368368        Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                                                           _assert_(strainparallel_input);
    369         Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);             _assert_(strainperpendicular_input);
    370         this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
     369        Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);              _assert_(strainperpendicular_input);
     370        Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
    371371
    372372        /* Start looping on the number of vertices: */
    373373        gauss=new GaussPenta();
    … …  
    380380                vel=vx*vx+vy*vy;
    381381                strainparallel_input->GetInputValue(&strainparallel,gauss);
    382382                strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
     383                levermanncoeff_input->GetInputValue(&propcoeff,gauss);
    383384
    384385                /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
    385386                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;   
    … …  
    393394        /*Add input*/
    394395        this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
    395396        this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
    396         this->inputs->AddInput(new PentaInput(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
     397        this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
    397398
    398399        /*Clean up and return*/
    399400        delete gauss;

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

     
    10911091                                name==LevelsetfunctionSlopeXEnum ||
    10921092                                name==LevelsetfunctionSlopeYEnum ||
    10931093                                name==LevelsetfunctionPicardEnum ||
    1094                                 name==MasstransportCalvingrateEnum ||
     1094                                //name==CalvingCalvingrateEnum ||
    10951095                                name==GradientEnum ||
    10961096                                name==OldGradientEnum  ||
    10971097                                name==ConvergedEnum ||
    … …  
    12111211                                break;
    12121212                        case CalvingratexEnum:
    12131213                        case CalvingrateyEnum:
    1214                         case MasstransportCalvingrateEnum:
     1214                        case CalvingCalvingrateEnum:
    12151215                                this->StrainRateparallel();
    12161216                                this->StrainRateperpendicular();
    12171217                                this->CalvingRateLevermann();
Note: See TracBrowser for help on using the repository browser.