Changeset 15386


Ignore:
Timestamp:
07/01/13 23:32:16 (12 years ago)
Author:
Eric.Larour
Message:

CHG: moved code around, to decrease the scope, and hopefully increase AD mode computations, which have too many actives.

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r15384 r15386  
    7575                virtual void   GetVectorFromResults(Vector<IssmDouble>* vector,int id,int enum_in,int interp)=0;
    7676                virtual void   InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max)=0;
    77                 virtual void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
    7877                virtual int*   GetHorizontalNeighboorSids(void)=0;
    7978                virtual IssmDouble TimeAdapt()=0;
    80                 virtual void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding)=0;
    81                 virtual void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
    8279                virtual void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0;
    8380                virtual void   Delta18oParameterization(void)=0;
    8481                virtual void   SmbGradients()=0;
    85                 virtual int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0;
    8682                virtual void UpdateConstraints(void)=0;
    8783                virtual void   ResetCoordinateSystem()=0;
    … …  
    137133                virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
    138134                #endif
     135               
     136                #ifdef _HAVE_GROUNDINGLINE_
     137                virtual void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding)=0;
     138                virtual void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
     139                virtual int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0;
     140                #endif
     141
     142                #ifdef _HAVE_WRAPPERS_
     143                virtual void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
     144                #endif
     145
    139146};
    140147#endif

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

    r15377 r15386  
    142142
    143143/*Other*/
    144 /*FUNCTION Penta::AverageOntoPartition {{{*/
    145 void  Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
    146         _error_("Not supported yet!");
    147 }
    148 /*}}}*/
    149144/*FUNCTION Penta::BedNormal {{{*/
    150145void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){
    … …  
    23182313
    23192314}/*}}}*/
    2320 /*FUNCTION Penta::MigrateGroundingLine{{{*/
    2321 void  Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
    2322 
    2323         int     i,migration_style;
    2324         bool    floatingelement = false;
    2325         bool    groundedelement = false;
    2326         IssmDouble  bed_hydro,yts,gl_melting_rate;
    2327         IssmDouble  rho_water,rho_ice,density;
    2328         IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];
    2329         IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
    2330 
    2331         if(!IsOnBed()) return;
    2332 
    2333         /*Recover info at the vertices: */
    2334         parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    2335         parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
    2336         parameters->FindParam(&yts,ConstantsYtsEnum);
    2337         GetInputListOnVertices(&h[0],ThicknessEnum);
    2338         GetInputListOnVertices(&s[0],SurfaceEnum);
    2339         GetInputListOnVertices(&b[0],BedEnum);
    2340         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2341         if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
    2342         rho_water=matpar->GetRhoWater();
    2343         rho_ice=matpar->GetRhoIce();
    2344         density=rho_ice/rho_water;
    2345 
    2346         /*go through vertices, and update inputs, considering them to be PentaVertex type: */
    2347         for(i=0;i<NUMVERTICES;i++){
    2348                 /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
    2349                 if(reCast<bool,IssmDouble>(old_floating_ice[nodes[i]->Sid()])){
    2350                         if(b[i]<=ba[i]){
    2351                                 b[i]=ba[i];
    2352                                 s[i]=b[i]+h[i];
    2353                                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
    2354                                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
    2355                         }
    2356                 }
    2357                 /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
    2358                 /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
    2359                 else{
    2360                         bed_hydro=-density*h[i];
    2361                         if (bed_hydro>ba[i]){
    2362                                 /*Unground only if the element is connected to the ice shelf*/
    2363                                 if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
    2364                                         s[i]=(1-density)*h[i];
    2365                                         b[i]=-density*h[i];
    2366                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2367                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2368                                 }
    2369                                 else if(migration_style==SoftMigrationEnum && reCast<int,IssmDouble>(sheet_ungrounding[nodes[i]->Sid()])){
    2370                                         s[i]=(1-density)*h[i];
    2371                                         b[i]=-density*h[i];
    2372                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2373                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2374                                 }
    2375                                 else{
    2376                                         if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
    2377                                 }
    2378                         }
    2379                 }
    2380         }
    2381 
    2382         /*SubelementMigrationEnum: if one grounded, all grounded*/
    2383         if(migration_style==SubelementMigrationEnum){
    2384                 for(i=0;i<NUMVERTICES;i++){
    2385                         if(nodes[i]->IsGrounded()){
    2386                                 groundedelement=true;
    2387                                 break;
    2388                         }
    2389                 }
    2390                 floatingelement=!groundedelement;
    2391         }
    2392         else{
    2393                 for(i=0;i<NUMVERTICES;i++){
    2394                         if(nodes[i]->IsFloating()){
    2395                                 floatingelement=true;
    2396                                 break;
    2397                         }
    2398                 }
    2399         }
    2400 
    2401    /*Add basal melting rate if element just ungrounded*/
    2402         if(!this->IsFloating() && floatingelement==true){
    2403                 for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
    2404                 this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
    2405         }
    2406 
    2407         /*Update inputs*/
    2408         this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
    2409         this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
    2410    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
    2411 
    2412         /*Recalculate phi*/
    2413         if(migration_style==SubelementMigrationEnum){
    2414                 for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
    2415                 this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
    2416                 this->InputExtrude(GLlevelsetEnum,ElementEnum);
    2417         }
    2418 
    2419         /*Extrude inputs*/
    2420         this->InputExtrude(SurfaceEnum,ElementEnum);
    2421         this->InputExtrude(BedEnum,ElementEnum);
    2422         this->InputExtrude(MaskElementonfloatingiceEnum,ElementEnum);
    2423         this->InputExtrude(MaskVertexonfloatingiceEnum,NodeEnum);
    2424         this->InputExtrude(MaskVertexongroundediceEnum,NodeEnum);
    2425 }
    2426 /*}}}*/
    24272315/*FUNCTION Penta::MinEdgeLength{{{*/
    24282316IssmDouble Penta::MinEdgeLength(IssmDouble xyz_list[6][3]){
    … …  
    25812469        /*clean-up*/
    25822470        delete gauss;
    2583 }
    2584 /*}}}*/
    2585 /*FUNCTION Penta::PotentialUngrounding{{{*/
    2586 void  Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
    2587 
    2588         int     i;
    2589         IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
    2590         IssmDouble  bed_hydro;
    2591         IssmDouble  rho_water,rho_ice,density;
    2592 
    2593         /*material parameters: */
    2594         rho_water=matpar->GetRhoWater();
    2595         rho_ice=matpar->GetRhoIce();
    2596         density=rho_ice/rho_water;
    2597         GetInputListOnVertices(&h[0],ThicknessEnum);
    2598         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2599 
    2600         /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
    2601         for(i=0;i<NUMVERTICES;i++){
    2602                 /*Find if grounded vertices want to start floating*/
    2603                 if (!nodes[i]->IsFloating()){
    2604                         bed_hydro=-density*h[i];
    2605                         if (bed_hydro>ba[i]){
    2606                                 /*Vertex that could potentially unground, flag it*/
    2607                                 potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
    2608                         }
    2609                 }
    2610         }
    26112471}
    26122472/*}}}*/
    … …  
    32063066}
    32073067/*}}}*/
    3208 /*FUNCTION Penta::UpdatePotentialUngrounding{{{*/
    3209 int Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
    3210 
    3211         int i;
    3212         int nflipped=0;
    3213 
    3214         /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
    3215         for(i=0;i<NUMVERTICES;i++){
    3216                 if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
    3217                         vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
    3218 
    3219                         /*If node was not on ice shelf, we flipped*/
    3220                         if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
    3221                                 nflipped++;
    3222                         }
    3223                 }
    3224         }
    3225         return nflipped;
    3226 }
    3227 /*}}}*/
    32283068/*FUNCTION Penta::UpdateConstraints{{{*/
    32293069void Penta::UpdateConstraints(void){
    … …  
    94259265        xDelete<int>(doflist);
    94269266}
    9427 /*{{{*/
     9267/*}}}*/
    94289268/*FUNCTION Penta::HydrologyEPLGetActive {{{*/
    94299269void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
    … …  
    95169356/*}}}*/
    95179357#endif
     9358
     9359#ifdef _HAVE_GROUNDINGLINE_
     9360/*FUNCTION Penta::MigrateGroundingLine{{{*/
     9361void  Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
     9362
     9363        int     i,migration_style;
     9364        bool    floatingelement = false;
     9365        bool    groundedelement = false;
     9366        IssmDouble  bed_hydro,yts,gl_melting_rate;
     9367        IssmDouble  rho_water,rho_ice,density;
     9368        IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];
     9369        IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
     9370
     9371        if(!IsOnBed()) return;
     9372
     9373        /*Recover info at the vertices: */
     9374        parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
     9375        parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
     9376        parameters->FindParam(&yts,ConstantsYtsEnum);
     9377        GetInputListOnVertices(&h[0],ThicknessEnum);
     9378        GetInputListOnVertices(&s[0],SurfaceEnum);
     9379        GetInputListOnVertices(&b[0],BedEnum);
     9380        GetInputListOnVertices(&ba[0],BathymetryEnum);
     9381        if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
     9382        rho_water=matpar->GetRhoWater();
     9383        rho_ice=matpar->GetRhoIce();
     9384        density=rho_ice/rho_water;
     9385
     9386        /*go through vertices, and update inputs, considering them to be PentaVertex type: */
     9387        for(i=0;i<NUMVERTICES;i++){
     9388                /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
     9389                if(reCast<bool,IssmDouble>(old_floating_ice[nodes[i]->Sid()])){
     9390                        if(b[i]<=ba[i]){
     9391                                b[i]=ba[i];
     9392                                s[i]=b[i]+h[i];
     9393                                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
     9394                                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
     9395                        }
     9396                }
     9397                /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
     9398                /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
     9399                else{
     9400                        bed_hydro=-density*h[i];
     9401                        if (bed_hydro>ba[i]){
     9402                                /*Unground only if the element is connected to the ice shelf*/
     9403                                if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
     9404                                        s[i]=(1-density)*h[i];
     9405                                        b[i]=-density*h[i];
     9406                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     9407                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     9408                                }
     9409                                else if(migration_style==SoftMigrationEnum && reCast<int,IssmDouble>(sheet_ungrounding[nodes[i]->Sid()])){
     9410                                        s[i]=(1-density)*h[i];
     9411                                        b[i]=-density*h[i];
     9412                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     9413                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     9414                                }
     9415                                else{
     9416                                        if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
     9417                                }
     9418                        }
     9419                }
     9420        }
     9421
     9422        /*SubelementMigrationEnum: if one grounded, all grounded*/
     9423        if(migration_style==SubelementMigrationEnum){
     9424                for(i=0;i<NUMVERTICES;i++){
     9425                        if(nodes[i]->IsGrounded()){
     9426                                groundedelement=true;
     9427                                break;
     9428                        }
     9429                }
     9430                floatingelement=!groundedelement;
     9431        }
     9432        else{
     9433                for(i=0;i<NUMVERTICES;i++){
     9434                        if(nodes[i]->IsFloating()){
     9435                                floatingelement=true;
     9436                                break;
     9437                        }
     9438                }
     9439        }
     9440
     9441   /*Add basal melting rate if element just ungrounded*/
     9442        if(!this->IsFloating() && floatingelement==true){
     9443                for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
     9444                this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
     9445        }
     9446
     9447        /*Update inputs*/
     9448        this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
     9449        this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
     9450   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
     9451
     9452        /*Recalculate phi*/
     9453        if(migration_style==SubelementMigrationEnum){
     9454                for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
     9455                this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
     9456                this->InputExtrude(GLlevelsetEnum,ElementEnum);
     9457        }
     9458
     9459        /*Extrude inputs*/
     9460        this->InputExtrude(SurfaceEnum,ElementEnum);
     9461        this->InputExtrude(BedEnum,ElementEnum);
     9462        this->InputExtrude(MaskElementonfloatingiceEnum,ElementEnum);
     9463        this->InputExtrude(MaskVertexonfloatingiceEnum,NodeEnum);
     9464        this->InputExtrude(MaskVertexongroundediceEnum,NodeEnum);
     9465}
     9466/*}}}*/
     9467#endif
     9468
     9469#ifdef _HAVE_WRAPPERS_
     9470/*FUNCTION Penta::PotentialUngrounding{{{*/
     9471void  Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
     9472
     9473        int     i;
     9474        IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
     9475        IssmDouble  bed_hydro;
     9476        IssmDouble  rho_water,rho_ice,density;
     9477
     9478        /*material parameters: */
     9479        rho_water=matpar->GetRhoWater();
     9480        rho_ice=matpar->GetRhoIce();
     9481        density=rho_ice/rho_water;
     9482        GetInputListOnVertices(&h[0],ThicknessEnum);
     9483        GetInputListOnVertices(&ba[0],BathymetryEnum);
     9484
     9485        /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
     9486        for(i=0;i<NUMVERTICES;i++){
     9487                /*Find if grounded vertices want to start floating*/
     9488                if (!nodes[i]->IsFloating()){
     9489                        bed_hydro=-density*h[i];
     9490                        if (bed_hydro>ba[i]){
     9491                                /*Vertex that could potentially unground, flag it*/
     9492                                potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
     9493                        }
     9494                }
     9495        }
     9496}
     9497/*}}}*/
     9498/*FUNCTION Penta::AverageOntoPartition {{{*/
     9499void  Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
     9500        _error_("Not supported yet!");
     9501}
     9502/*}}}*/
     9503/*FUNCTION Penta::UpdatePotentialUngrounding{{{*/
     9504int Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
     9505
     9506        int i;
     9507        int nflipped=0;
     9508
     9509        /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
     9510        for(i=0;i<NUMVERTICES;i++){
     9511                if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
     9512                        vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
     9513
     9514                        /*If node was not on ice shelf, we flipped*/
     9515                        if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
     9516                                nflipped++;
     9517                        }
     9518                }
     9519        }
     9520        return nflipped;
     9521}
     9522/*}}}*/
     9523#endif

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

    r15384 r15386  
    7575                /*}}}*/
    7676                /*Element virtual functions definitions: {{{*/
    77                 void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    7877                void   BasalFrictionCreateInput(void);
    7978                void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
    … …  
    106105
    107106                void   InputToResult(int enum_type,int step,IssmDouble time);
    108                 void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding);
    109                 void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    110107                void   RequestedOutput(int output_enum,int step,IssmDouble time);
    111108                void   ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results);
    … …  
    117114                IssmDouble SurfaceArea(void);
    118115                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    119                 int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
    120116                int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
    121117                IssmDouble TimeAdapt();
    … …  
    173169                void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
    174170                #endif
     171               
     172                #ifdef _HAVE_GROUNDINGLINE_
     173                void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding);
     174                void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
     175                int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
     176                #endif
     177
     178                #ifdef _HAVE_WRAPPERS_
     179                void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
     180                #endif
     181
    175182                /*}}}*/
    176183                /*Penta specific routines:{{{*/

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

    r15377 r15386  
    126126
    127127/*Other*/
    128 /*FUNCTION Tria::AverageOntoPartition {{{*/
    129 void  Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
    130 
    131         bool       already = false;
    132         int        i,j;
    133         int        partition[NUMVERTICES];
    134         int        offsetsid[NUMVERTICES];
    135         int        offsetdof[NUMVERTICES];
    136         IssmDouble area;
    137         IssmDouble mean;
    138 
    139         /*First, get the area: */
    140         area=this->GetArea();
    141 
    142         /*Figure out the average for this element: */
    143         this->GetVertexSidList(&offsetsid[0]);
    144         this->GetVertexPidList(&offsetdof[0]);
    145         mean=0;
    146         for(i=0;i<NUMVERTICES;i++){
    147                 partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
    148                 mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
    149         }
    150 
    151         /*Add contribution: */
    152         for(i=0;i<NUMVERTICES;i++){
    153                 already=false;
    154                 for(j=0;j<i;j++){
    155                         if (partition[i]==partition[j]){
    156                                 already=true;
    157                                 break;
    158                         }
    159                 }
    160                 if(!already){
    161                         partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
    162                         partition_areas->SetValue(partition[i],area,ADD_VAL);
    163                 };
    164         }
    165 }
    166 /*}}}*/
    167128/*FUNCTION Tria::SetwiseNodeConnectivity{{{*/
    168129void Tria::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int set1_enum,int set2_enum){
    … …  
    282243}
    283244/*}}}*/
    284 /*FUNCTION Tria::CreateKMatrixMelting {{{*/
    285 ElementMatrix* Tria::CreateKMatrixMelting(void){
    286 
    287         /*Constants*/
    288         const int  numdof=NUMVERTICES*NDOF1;
    289 
    290         /*Intermediaries */
    291         IssmDouble heatcapacity,latentheat;
    292         IssmDouble Jdet,D_scalar;
    293         IssmDouble xyz_list[NUMVERTICES][3];
    294         IssmDouble basis[3];
    295         GaussTria *gauss=NULL;
    296 
    297         /*Initialize Element matrix*/
    298         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    299 
    300         /*Retrieve all inputs and parameters*/
    301         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    302         latentheat=matpar->GetLatentHeat();
    303         heatcapacity=matpar->GetHeatCapacity();
    304 
    305         /* Start looping on the number of gauss  (nodes on the bedrock) */
    306         gauss=new GaussTria(2);
    307         for(int ig=gauss->begin();ig<gauss->end();ig++){
    308 
    309                 gauss->GaussPoint(ig);
    310 
    311                 GetNodalFunctions(&basis[0],gauss);
    312                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);
    313 
    314                 D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
    315 
    316                 TripleMultiply(&basis[0],numdof,1,0,
    317                                         &D_scalar,1,1,0,
    318                                         &basis[0],1,numdof,0,
    319                                         &Ke->values[0],1);
    320         }
    321 
    322         /*Clean up and return*/
    323         delete gauss;
    324         return Ke;
    325 }
    326 /*}}}*/
    327 /*FUNCTION Tria::CreateKMatrixPrognostic {{{*/
    328 ElementMatrix* Tria::CreateKMatrixPrognostic(void){
    329 
    330         switch(GetElementType()){
    331                 case P1Enum:
    332                         return CreateKMatrixPrognostic_CG();
    333                 case P1DGEnum:
    334                         return CreateKMatrixPrognostic_DG();
    335                 default:
    336                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    337         }
    338 }
    339 /*}}}*/
    340 /*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{*/
    341 ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){
    342 
    343         /*Constants*/
    344         const int    numdof=NDOF1*NUMVERTICES;
    345 
    346         /*Intermediaries */
    347         int        stabilization;
    348         int        dim;
    349         IssmDouble Jdettria,DL_scalar,dt,h;
    350         IssmDouble vel,vx,vy,dvxdx,dvydy;
    351         IssmDouble dvx[2],dvy[2];
    352         IssmDouble v_gauss[2]={0.0};
    353         IssmDouble xyz_list[NUMVERTICES][3];
    354         IssmDouble basis[NUMVERTICES];
    355         IssmDouble B[2][NUMVERTICES];
    356         IssmDouble Bprime[2][NUMVERTICES];
    357         IssmDouble K[2][2]                        ={0.0};
    358         IssmDouble KDL[2][2]                      ={0.0};
    359         IssmDouble DL[2][2]                        ={0.0};
    360         IssmDouble DLprime[2][2]                   ={0.0};
    361 
    362         /*Initialize Element matrix*/
    363         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    364 
    365         /*Retrieve all inputs and parameters*/
    366         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    367         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    368         this->parameters->FindParam(&dim,MeshDimensionEnum);
    369         this->parameters->FindParam(&stabilization,PrognosticStabilizationEnum);
    370         Input* vxaverage_input=NULL;
    371         Input* vyaverage_input=NULL;
    372         if(dim==2){
    373                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    374                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    375         }
    376         else{
    377                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    378                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    379         }
    380         h=sqrt(2*this->GetArea());
    381 
    382         /* Start  looping on the number of gaussian points: */
    383         GaussTria *gauss=new GaussTria(2);
    384         for(int ig=gauss->begin();ig<gauss->end();ig++){
    385 
    386                 gauss->GaussPoint(ig);
    387 
    388                 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    389                 GetNodalFunctions(&basis[0],gauss);
    390 
    391                 vxaverage_input->GetInputValue(&vx,gauss);
    392                 vyaverage_input->GetInputValue(&vy,gauss);
    393                 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    394                 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    395 
    396                 DL_scalar=gauss->weight*Jdettria;
    397 
    398                 TripleMultiply(&basis[0],1,numdof,1,
    399                                         &DL_scalar,1,1,0,
    400                                         &basis[0],1,numdof,0,
    401                                         &Ke->values[0],1);
    402 
    403                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    404                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    405 
    406                 dvxdx=dvx[0];
    407                 dvydy=dvy[1];
    408                 DL_scalar=dt*gauss->weight*Jdettria;
    409 
    410                 DL[0][0]=DL_scalar*dvxdx;
    411                 DL[1][1]=DL_scalar*dvydy;
    412                 DLprime[0][0]=DL_scalar*vx;
    413                 DLprime[1][1]=DL_scalar*vy;
    414 
    415                 TripleMultiply( &B[0][0],2,numdof,1,
    416                                         &DL[0][0],2,2,0,
    417                                         &B[0][0],2,numdof,0,
    418                                         &Ke->values[0],1);
    419 
    420                 TripleMultiply( &B[0][0],2,numdof,1,
    421                                         &DLprime[0][0],2,2,0,
    422                                         &Bprime[0][0],2,numdof,0,
    423                                         &Ke->values[0],1);
    424 
    425                 if(stabilization==2){
    426                         /*Streamline upwinding*/
    427                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    428                         K[0][0]=h/(2*vel)*vx*vx;
    429                         K[1][0]=h/(2*vel)*vy*vx;
    430                         K[0][1]=h/(2*vel)*vx*vy;
    431                         K[1][1]=h/(2*vel)*vy*vy;
    432                 }
    433                 else if(stabilization==1){
    434                         /*MacAyeal*/
    435                         vxaverage_input->GetInputAverage(&vx);
    436                         vyaverage_input->GetInputAverage(&vy);
    437                         K[0][0]=h/2.0*fabs(vx);
    438                         K[0][1]=0.;
    439                         K[1][0]=0.;
    440                         K[1][1]=h/2.0*fabs(vy);
    441                 }
    442                 if(stabilization==1 || stabilization==2){
    443                         KDL[0][0]=DL_scalar*K[0][0];
    444                         KDL[1][0]=DL_scalar*K[1][0];
    445                         KDL[0][1]=DL_scalar*K[0][1];
    446                         KDL[1][1]=DL_scalar*K[1][1];
    447                         TripleMultiply( &Bprime[0][0],2,numdof,1,
    448                                                 &KDL[0][0],2,2,0,
    449                                                 &Bprime[0][0],2,numdof,0,
    450                                                 &Ke->values[0],1);
    451                 }
    452         }
    453 
    454         /*Clean up and return*/
    455         delete gauss;
    456         return Ke;
    457 }
    458 /*}}}*/
    459 /*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{*/
    460 ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
    461 
    462         /*Constants*/
    463         const int    numdof=NDOF1*NUMVERTICES;
    464 
    465         /*Intermediaries */
    466         int        dim;
    467         IssmDouble xyz_list[NUMVERTICES][3];
    468         IssmDouble Jdettria,dt,vx,vy;
    469         IssmDouble basis[NUMVERTICES];
    470         IssmDouble B[2][NUMVERTICES];
    471         IssmDouble Bprime[2][NUMVERTICES];
    472         IssmDouble DL[2][2]={0.0};
    473         IssmDouble DLprime[2][2]={0.0};
    474         IssmDouble DL_scalar;
    475         GaussTria  *gauss=NULL;
    476 
    477         /*Initialize Element matrix*/
    478         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    479 
    480         /*Retrieve all inputs and parameters*/
    481         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    482         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    483         this->parameters->FindParam(&dim,MeshDimensionEnum);
    484         Input* vxaverage_input=NULL;
    485         Input* vyaverage_input=NULL;
    486         if(dim==2){
    487                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    488                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    489         }
    490         else{
    491                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    492                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    493         }
    494 
    495         /* Start  looping on the number of gaussian points: */
    496         gauss=new GaussTria(2);
    497         for(int ig=gauss->begin();ig<gauss->end();ig++){
    498 
    499                 gauss->GaussPoint(ig);
    500 
    501                 vxaverage_input->GetInputValue(&vx,gauss);
    502                 vyaverage_input->GetInputValue(&vy,gauss);
    503 
    504                 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    505                 GetNodalFunctions(&basis[0],gauss);
    506 
    507                 DL_scalar=gauss->weight*Jdettria;
    508 
    509                 TripleMultiply(&basis[0],1,numdof,1,
    510                                         &DL_scalar,1,1,0,
    511                                         &basis[0],1,numdof,0,
    512                                         &Ke->values[0],1);
    513 
    514                 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    515                 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    516                 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
    517 
    518                 DL_scalar=-dt*gauss->weight*Jdettria;
    519 
    520                 DLprime[0][0]=DL_scalar*vx;
    521                 DLprime[1][1]=DL_scalar*vy;
    522 
    523                 TripleMultiply( &B[0][0],2,numdof,1,
    524                                         &DLprime[0][0],2,2,0,
    525                                         &Bprime[0][0],2,numdof,0,
    526                                         &Ke->values[0],1);
    527         }
    528 
    529         /*Clean up and return*/
    530         delete gauss;
    531         return Ke;
    532 }
    533 /*}}}*/
    534245/*FUNCTION Tria::CreateMassMatrix {{{*/
    535246ElementMatrix* Tria::CreateMassMatrix(void){
    … …  
    641352                delete pe;
    642353        }
    643 }
    644 /*}}}*/
    645 /*FUNCTION Tria::CreatePVectorPrognostic{{{*/
    646 ElementVector* Tria::CreatePVectorPrognostic(void){
    647 
    648         switch(GetElementType()){
    649                 case P1Enum:
    650                         return CreatePVectorPrognostic_CG();
    651                 case P1DGEnum:
    652                         return CreatePVectorPrognostic_DG();
    653                 default:
    654                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    655         }
    656 }
    657 /*}}}*/
    658 /*FUNCTION Tria::CreatePVectorPrognostic_CG {{{*/
    659 ElementVector* Tria::CreatePVectorPrognostic_CG(void){
    660 
    661         /*Constants*/
    662         const int    numdof=NDOF1*NUMVERTICES;
    663 
    664         /*Intermediaries */
    665         IssmDouble Jdettria,dt;
    666         IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
    667         IssmDouble xyz_list[NUMVERTICES][3];
    668         IssmDouble basis[NUMVERTICES];
    669         GaussTria* gauss=NULL;
    670 
    671         /*Initialize Element vector*/
    672         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    673 
    674         /*Retrieve all inputs and parameters*/
    675         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    676         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    677         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    678         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    679         Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
    680         Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
    681 
    682         /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    683         /* Start  looping on the number of gaussian points: */
    684         gauss=new GaussTria(2);
    685         for(int ig=gauss->begin();ig<gauss->end();ig++){
    686 
    687                 gauss->GaussPoint(ig);
    688 
    689                 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    690                 GetNodalFunctions(&basis[0],gauss);
    691 
    692                 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    693                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
    694                 thickness_input->GetInputValue(&thickness_g,gauss);
    695                 if(basal_melting_correction_input)
    696                  basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);
    697                 else
    698                  basal_melting_correction_g=0.;
    699 
    700                 for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
    701         }
    702 
    703         /*Clean up and return*/
    704         delete gauss;
    705         return pe;
    706 }
    707 /*}}}*/
    708 /*FUNCTION Tria::CreatePVectorPrognostic_DG {{{*/
    709 ElementVector* Tria::CreatePVectorPrognostic_DG(void){
    710 
    711         /*Constants*/
    712         const int    numdof=NDOF1*NUMVERTICES;
    713 
    714         /*Intermediaries */
    715         IssmDouble Jdettria,dt;
    716         IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
    717         IssmDouble xyz_list[NUMVERTICES][3];
    718         IssmDouble basis[NUMVERTICES];
    719         GaussTria* gauss=NULL;
    720 
    721         /*Initialize Element vector*/
    722         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    723 
    724         /*Retrieve all inputs and parameters*/
    725         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    726         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    727         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    728         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    729         Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
    730 
    731         /* Start  looping on the number of gaussian points: */
    732         gauss=new GaussTria(2);
    733         for(int ig=gauss->begin();ig<gauss->end();ig++){
    734 
    735                 gauss->GaussPoint(ig);
    736 
    737                 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    738                 GetNodalFunctions(&basis[0],gauss);
    739 
    740                 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    741                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
    742                 thickness_input->GetInputValue(&thickness_g,gauss);
    743 
    744                 for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
    745         }
    746 
    747         /*Clean up and return*/
    748         delete gauss;
    749         return pe;
    750354}
    751355/*}}}*/
    … …  
    22721876
    22731877}/*}}}*/
    2274 /*FUNCTION Tria::MigrateGroundingLine{{{*/
    2275 void  Tria::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
    2276 
    2277         int     i,migration_style;
    2278         bool    floatingelement = false;
    2279         bool    groundedelement = false;
    2280         IssmDouble  bed_hydro,yts,gl_melting_rate;
    2281         IssmDouble  rho_water,rho_ice,density;
    2282         IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];;
    2283         IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
    2284 
    2285         /*Recover info at the vertices: */
    2286         parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    2287         parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
    2288         parameters->FindParam(&yts,ConstantsYtsEnum);
    2289         GetInputListOnVertices(&h[0],ThicknessEnum);
    2290         GetInputListOnVertices(&s[0],SurfaceEnum);
    2291         GetInputListOnVertices(&b[0],BedEnum);
    2292         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2293         if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
    2294         rho_water=matpar->GetRhoWater();
    2295         rho_ice=matpar->GetRhoIce();
    2296         density=rho_ice/rho_water;
    2297 
    2298         /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    2299         for(i=0;i<NUMVERTICES;i++){
    2300                 /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
    2301                 if(reCast<bool>(old_floating_ice[nodes[i]->Sid()])){
    2302                         if(b[i]<=ba[i]){
    2303                                 b[i]=ba[i];
    2304                                 s[i]=b[i]+h[i];
    2305                                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
    2306                                 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
    2307                         }
    2308                 }
    2309                 /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
    2310                 /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
    2311                 else{
    2312                         bed_hydro=-density*h[i];
    2313                         if (bed_hydro>ba[i]){
    2314                                 /*Unground only if the element is connected to the ice shelf*/
    2315                                 if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
    2316                                         s[i]=(1-density)*h[i];
    2317                                         b[i]=-density*h[i];
    2318                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2319                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2320                                 }
    2321                                 else if(migration_style==SoftMigrationEnum && reCast<bool>(sheet_ungrounding[nodes[i]->Sid()])){
    2322                                         s[i]=(1-density)*h[i];
    2323                                         b[i]=-density*h[i];
    2324                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
    2325                                         nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
    2326                                 }
    2327                                 else{
    2328                                         if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
    2329                                 }
    2330                         }
    2331                 }
    2332         }
    2333 
    2334         /*SubelementMigrationEnum: if one grounded, all grounded*/
    2335         if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
    2336                 for(i=0;i<NUMVERTICES;i++){
    2337                         if(nodes[i]->IsGrounded()){
    2338                                 groundedelement=true;
    2339                                 break;
    2340                         }
    2341                 }
    2342                 floatingelement=!groundedelement;
    2343         }
    2344         else{
    2345                 /*Otherwise: if one floating, all floating*/
    2346                 for(i=0;i<NUMVERTICES;i++){
    2347                         if(nodes[i]->IsFloating()){
    2348                                 floatingelement=true;
    2349                                 break;
    2350                         }
    2351                 }
    2352         }
    2353 
    2354    /*Add basal melting rate if element just ungrounded*/
    2355         if(!this->IsFloating() && floatingelement==true){
    2356                 for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
    2357                 this->inputs->AddInput(new TriaInput(BasalforcingsMeltingRateEnum,&melting[0],P1Enum));
    2358         }
    2359 
    2360         /*Update inputs*/
    2361    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
    2362         this->inputs->AddInput(new TriaInput(SurfaceEnum,&s[0],P1Enum));
    2363         this->inputs->AddInput(new TriaInput(BedEnum,&b[0],P1Enum));
    2364 
    2365         /*Recalculate phi*/
    2366         if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
    2367                 for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
    2368                 this->inputs->AddInput(new TriaInput(GLlevelsetEnum,&phi[0],P1Enum));
    2369         }
    2370 }
    2371 /*}}}*/
    23721878/*FUNCTION Tria::NodalValue {{{*/
    23731879int    Tria::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){
    … …  
    24521958        *pnumvertices=NUMVERTICES;
    24531959        *pnumnodes=numnodes;
    2454 }
    2455 /*}}}*/
    2456 /*FUNCTION Tria::PotentialUngrounding{{{*/
    2457 void  Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
    2458 
    2459         int     i;
    2460         IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
    2461         IssmDouble  bed_hydro;
    2462         IssmDouble  rho_water,rho_ice,density;
    2463 
    2464         /*material parameters: */
    2465         rho_water=matpar->GetRhoWater();
    2466         rho_ice=matpar->GetRhoIce();
    2467         density=rho_ice/rho_water;
    2468         GetInputListOnVertices(&h[0],ThicknessEnum);
    2469         GetInputListOnVertices(&ba[0],BathymetryEnum);
    2470 
    2471         /*go through vertices, and figure out which ones are grounded and want to unground: */
    2472         for(i=0;i<NUMVERTICES;i++){
    2473                 /*Find if grounded vertices want to start floating*/
    2474                 if (!nodes[i]->IsFloating()){
    2475                         bed_hydro=-density*h[i];
    2476                         if (bed_hydro>ba[i]){
    2477                                 /*Vertex that could potentially unground, flag it*/
    2478                                 potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
    2479                         }
    2480                 }
    2481         }
    24821960}
    24831961/*}}}*/
    … …  
    28342312        if(IsOnWater()) return;
    28352313
    2836 }
    2837 /*}}}*/
    2838 /*FUNCTION Tria::UpdatePotentialUngrounding{{{*/
    2839 int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
    2840 
    2841         int i;
    2842         int nflipped=0;
    2843 
    2844         /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
    2845         for(i=0;i<3;i++){
    2846                 if (reCast<bool>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
    2847                         vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
    2848 
    2849                         /*If node was not on ice shelf, we flipped*/
    2850                         if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
    2851                                 nflipped++;
    2852                         }
    2853                 }
    2854         }
    2855         return nflipped;
    28562314}
    28572315/*}}}*/
    … …  
    58525310#endif
    58535311
     5312#ifdef _HAVE_THERMAL_
     5313/*FUNCTION Tria::CreateKMatrixMelting {{{*/
     5314ElementMatrix* Tria::CreateKMatrixMelting(void){
     5315
     5316        /*Constants*/
     5317        const int  numdof=NUMVERTICES*NDOF1;
     5318
     5319        /*Intermediaries */
     5320        IssmDouble heatcapacity,latentheat;
     5321        IssmDouble Jdet,D_scalar;
     5322        IssmDouble xyz_list[NUMVERTICES][3];
     5323        IssmDouble basis[3];
     5324        GaussTria *gauss=NULL;
     5325
     5326        /*Initialize Element matrix*/
     5327        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     5328
     5329        /*Retrieve all inputs and parameters*/
     5330        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5331        latentheat=matpar->GetLatentHeat();
     5332        heatcapacity=matpar->GetHeatCapacity();
     5333
     5334        /* Start looping on the number of gauss  (nodes on the bedrock) */
     5335        gauss=new GaussTria(2);
     5336        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5337
     5338                gauss->GaussPoint(ig);
     5339
     5340                GetNodalFunctions(&basis[0],gauss);
     5341                GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);
     5342
     5343                D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
     5344
     5345                TripleMultiply(&basis[0],numdof,1,0,
     5346                                        &D_scalar,1,1,0,
     5347                                        &basis[0],1,numdof,0,
     5348                                        &Ke->values[0],1);
     5349        }
     5350
     5351        /*Clean up and return*/
     5352        delete gauss;
     5353        return Ke;
     5354}
     5355/*}}}*/
     5356#endif
     5357
    58545358#ifdef _HAVE_HYDROLOGY_
    58555359/*FUNCTION Tria::AllActive{{{*/
    … …  
    66326136#endif
    66336137
     6138#ifdef _HAVE_PROGNOSTIC_
     6139/*FUNCTION Tria::CreateKMatrixPrognostic {{{*/
     6140ElementMatrix* Tria::CreateKMatrixPrognostic(void){
     6141
     6142        switch(GetElementType()){
     6143                case P1Enum:
     6144                        return CreateKMatrixPrognostic_CG();
     6145                case P1DGEnum:
     6146                        return CreateKMatrixPrognostic_DG();
     6147                default:
     6148                        _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
     6149        }
     6150}
     6151/*}}}*/
     6152/*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{*/
     6153ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){
     6154
     6155        /*Constants*/
     6156        const int    numdof=NDOF1*NUMVERTICES;
     6157
     6158        /*Intermediaries */
     6159        int        stabilization;
     6160        int        dim;
     6161        IssmDouble Jdettria,DL_scalar,dt,h;
     6162        IssmDouble vel,vx,vy,dvxdx,dvydy;
     6163        IssmDouble dvx[2],dvy[2];
     6164        IssmDouble v_gauss[2]={0.0};
     6165        IssmDouble xyz_list[NUMVERTICES][3];
     6166        IssmDouble basis[NUMVERTICES];
     6167        IssmDouble B[2][NUMVERTICES];
     6168        IssmDouble Bprime[2][NUMVERTICES];
     6169        IssmDouble K[2][2]                        ={0.0};
     6170        IssmDouble KDL[2][2]                      ={0.0};
     6171        IssmDouble DL[2][2]                        ={0.0};
     6172        IssmDouble DLprime[2][2]                   ={0.0};
     6173
     6174        /*Initialize Element matrix*/
     6175        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     6176
     6177        /*Retrieve all inputs and parameters*/
     6178        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6179        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6180        this->parameters->FindParam(&dim,MeshDimensionEnum);
     6181        this->parameters->FindParam(&stabilization,PrognosticStabilizationEnum);
     6182        Input* vxaverage_input=NULL;
     6183        Input* vyaverage_input=NULL;
     6184        if(dim==2){
     6185                vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
     6186                vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
     6187        }
     6188        else{
     6189                vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
     6190                vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     6191        }
     6192        h=sqrt(2*this->GetArea());
     6193
     6194        /* Start  looping on the number of gaussian points: */
     6195        GaussTria *gauss=new GaussTria(2);
     6196        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6197
     6198                gauss->GaussPoint(ig);
     6199
     6200                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6201                GetNodalFunctions(&basis[0],gauss);
     6202
     6203                vxaverage_input->GetInputValue(&vx,gauss);
     6204                vyaverage_input->GetInputValue(&vy,gauss);
     6205                vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     6206                vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     6207
     6208                DL_scalar=gauss->weight*Jdettria;
     6209
     6210                TripleMultiply(&basis[0],1,numdof,1,
     6211                                        &DL_scalar,1,1,0,
     6212                                        &basis[0],1,numdof,0,
     6213                                        &Ke->values[0],1);
     6214
     6215                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     6216                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     6217
     6218                dvxdx=dvx[0];
     6219                dvydy=dvy[1];
     6220                DL_scalar=dt*gauss->weight*Jdettria;
     6221
     6222                DL[0][0]=DL_scalar*dvxdx;
     6223                DL[1][1]=DL_scalar*dvydy;
     6224                DLprime[0][0]=DL_scalar*vx;
     6225                DLprime[1][1]=DL_scalar*vy;
     6226
     6227                TripleMultiply( &B[0][0],2,numdof,1,
     6228                                        &DL[0][0],2,2,0,
     6229                                        &B[0][0],2,numdof,0,
     6230                                        &Ke->values[0],1);
     6231
     6232                TripleMultiply( &B[0][0],2,numdof,1,
     6233                                        &DLprime[0][0],2,2,0,
     6234                                        &Bprime[0][0],2,numdof,0,
     6235                                        &Ke->values[0],1);
     6236
     6237                if(stabilization==2){
     6238                        /*Streamline upwinding*/
     6239                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
     6240                        K[0][0]=h/(2*vel)*vx*vx;
     6241                        K[1][0]=h/(2*vel)*vy*vx;
     6242                        K[0][1]=h/(2*vel)*vx*vy;
     6243                        K[1][1]=h/(2*vel)*vy*vy;
     6244                }
     6245                else if(stabilization==1){
     6246                        /*MacAyeal*/
     6247                        vxaverage_input->GetInputAverage(&vx);
     6248                        vyaverage_input->GetInputAverage(&vy);
     6249                        K[0][0]=h/2.0*fabs(vx);
     6250                        K[0][1]=0.;
     6251                        K[1][0]=0.;
     6252                        K[1][1]=h/2.0*fabs(vy);
     6253                }
     6254                if(stabilization==1 || stabilization==2){
     6255                        KDL[0][0]=DL_scalar*K[0][0];
     6256                        KDL[1][0]=DL_scalar*K[1][0];
     6257                        KDL[0][1]=DL_scalar*K[0][1];
     6258                        KDL[1][1]=DL_scalar*K[1][1];
     6259                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     6260                                                &KDL[0][0],2,2,0,
     6261                                                &Bprime[0][0],2,numdof,0,
     6262                                                &Ke->values[0],1);
     6263                }
     6264        }
     6265
     6266        /*Clean up and return*/
     6267        delete gauss;
     6268        return Ke;
     6269}
     6270/*}}}*/
     6271/*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{*/
     6272ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
     6273
     6274        /*Constants*/
     6275        const int    numdof=NDOF1*NUMVERTICES;
     6276
     6277        /*Intermediaries */
     6278        int        dim;
     6279        IssmDouble xyz_list[NUMVERTICES][3];
     6280        IssmDouble Jdettria,dt,vx,vy;
     6281        IssmDouble basis[NUMVERTICES];
     6282        IssmDouble B[2][NUMVERTICES];
     6283        IssmDouble Bprime[2][NUMVERTICES];
     6284        IssmDouble DL[2][2]={0.0};
     6285        IssmDouble DLprime[2][2]={0.0};
     6286        IssmDouble DL_scalar;
     6287        GaussTria  *gauss=NULL;
     6288
     6289        /*Initialize Element matrix*/
     6290        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     6291
     6292        /*Retrieve all inputs and parameters*/
     6293        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6294        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6295        this->parameters->FindParam(&dim,MeshDimensionEnum);
     6296        Input* vxaverage_input=NULL;
     6297        Input* vyaverage_input=NULL;
     6298        if(dim==2){
     6299                vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
     6300                vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
     6301        }
     6302        else{
     6303                vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
     6304                vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     6305        }
     6306
     6307        /* Start  looping on the number of gaussian points: */
     6308        gauss=new GaussTria(2);
     6309        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6310
     6311                gauss->GaussPoint(ig);
     6312
     6313                vxaverage_input->GetInputValue(&vx,gauss);
     6314                vyaverage_input->GetInputValue(&vy,gauss);
     6315
     6316                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6317                GetNodalFunctions(&basis[0],gauss);
     6318
     6319                DL_scalar=gauss->weight*Jdettria;
     6320
     6321                TripleMultiply(&basis[0],1,numdof,1,
     6322                                        &DL_scalar,1,1,0,
     6323                                        &basis[0],1,numdof,0,
     6324                                        &Ke->values[0],1);
     6325
     6326                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
     6327                GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     6328                GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
     6329
     6330                DL_scalar=-dt*gauss->weight*Jdettria;
     6331
     6332                DLprime[0][0]=DL_scalar*vx;
     6333                DLprime[1][1]=DL_scalar*vy;
     6334
     6335                TripleMultiply( &B[0][0],2,numdof,1,
     6336                                        &DLprime[0][0],2,2,0,
     6337                                        &Bprime[0][0],2,numdof,0,
     6338                                        &Ke->values[0],1);
     6339        }
     6340
     6341        /*Clean up and return*/
     6342        delete gauss;
     6343        return Ke;
     6344}
     6345/*}}}*/
     6346/*FUNCTION Tria::CreatePVectorPrognostic{{{*/
     6347ElementVector* Tria::CreatePVectorPrognostic(void){
     6348
     6349        switch(GetElementType()){
     6350                case P1Enum:
     6351                        return CreatePVectorPrognostic_CG();
     6352                case P1DGEnum:
     6353                        return CreatePVectorPrognostic_DG();
     6354                default:
     6355                        _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
     6356        }
     6357}
     6358/*}}}*/
     6359/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{*/
     6360ElementVector* Tria::CreatePVectorPrognostic_CG(void){
     6361
     6362        /*Constants*/
     6363        const int    numdof=NDOF1*NUMVERTICES;
     6364
     6365        /*Intermediaries */
     6366        IssmDouble Jdettria,dt;
     6367        IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
     6368        IssmDouble xyz_list[NUMVERTICES][3];
     6369        IssmDouble basis[NUMVERTICES];
     6370        GaussTria* gauss=NULL;
     6371
     6372        /*Initialize Element vector*/
     6373        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     6374
     6375        /*Retrieve all inputs and parameters*/
     6376        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6377        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6378        Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
     6379        Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
     6380        Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
     6381        Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
     6382
     6383        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
     6384        /* Start  looping on the number of gaussian points: */
     6385        gauss=new GaussTria(2);
     6386        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6387
     6388                gauss->GaussPoint(ig);
     6389
     6390                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6391                GetNodalFunctions(&basis[0],gauss);
     6392
     6393                surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
     6394                basal_melting_input->GetInputValue(&basal_melting_g,gauss);
     6395                thickness_input->GetInputValue(&thickness_g,gauss);
     6396                if(basal_melting_correction_input)
     6397                 basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);
     6398                else
     6399                 basal_melting_correction_g=0.;
     6400
     6401                for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
     6402        }
     6403
     6404        /*Clean up and return*/
     6405        delete gauss;
     6406        return pe;
     6407}
     6408/*}}}*/
     6409/*FUNCTION Tria::CreatePVectorPrognostic_DG {{{*/
     6410ElementVector* Tria::CreatePVectorPrognostic_DG(void){
     6411
     6412        /*Constants*/
     6413        const int    numdof=NDOF1*NUMVERTICES;
     6414
     6415        /*Intermediaries */
     6416        IssmDouble Jdettria,dt;
     6417        IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
     6418        IssmDouble xyz_list[NUMVERTICES][3];
     6419        IssmDouble basis[NUMVERTICES];
     6420        GaussTria* gauss=NULL;
     6421
     6422        /*Initialize Element vector*/
     6423        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     6424
     6425        /*Retrieve all inputs and parameters*/
     6426        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     6427        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6428        Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
     6429        Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
     6430        Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
     6431
     6432        /* Start  looping on the number of gaussian points: */
     6433        gauss=new GaussTria(2);
     6434        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6435
     6436                gauss->GaussPoint(ig);
     6437
     6438                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
     6439                GetNodalFunctions(&basis[0],gauss);
     6440
     6441                surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
     6442                basal_melting_input->GetInputValue(&basal_melting_g,gauss);
     6443                thickness_input->GetInputValue(&thickness_g,gauss);
     6444
     6445                for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
     6446        }
     6447
     6448        /*Clean up and return*/
     6449        delete gauss;
     6450        return pe;
     6451}
     6452/*}}}*/
     6453#endif
     6454
    66346455#ifdef _HAVE_DAKOTA_
    66356456/*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/
    … …  
    70586879/*}}}*/
    70596880#endif
     6881
     6882#ifdef _HAVE_GROUNDINGLINE_
     6883/*FUNCTION Tria::MigrateGroundingLine{{{*/
     6884void  Tria::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
     6885
     6886        int     i,migration_style;
     6887        bool    floatingelement = false;
     6888        bool    groundedelement = false;
     6889        IssmDouble  bed_hydro,yts,gl_melting_rate;
     6890        IssmDouble  rho_water,rho_ice,density;
     6891        IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];;
     6892        IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
     6893
     6894        /*Recover info at the vertices: */
     6895        parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
     6896        parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
     6897        parameters->FindParam(&yts,ConstantsYtsEnum);
     6898        GetInputListOnVertices(&h[0],ThicknessEnum);
     6899        GetInputListOnVertices(&s[0],SurfaceEnum);
     6900        GetInputListOnVertices(&b[0],BedEnum);
     6901        GetInputListOnVertices(&ba[0],BathymetryEnum);
     6902        if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
     6903        rho_water=matpar->GetRhoWater();
     6904        rho_ice=matpar->GetRhoIce();
     6905        density=rho_ice/rho_water;
     6906
     6907        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
     6908        for(i=0;i<NUMVERTICES;i++){
     6909                /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
     6910                if(reCast<bool>(old_floating_ice[nodes[i]->Sid()])){
     6911                        if(b[i]<=ba[i]){
     6912                                b[i]=ba[i];
     6913                                s[i]=b[i]+h[i];
     6914                                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
     6915                                nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
     6916                        }
     6917                }
     6918                /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
     6919                /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
     6920                else{
     6921                        bed_hydro=-density*h[i];
     6922                        if (bed_hydro>ba[i]){
     6923                                /*Unground only if the element is connected to the ice shelf*/
     6924                                if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
     6925                                        s[i]=(1-density)*h[i];
     6926                                        b[i]=-density*h[i];
     6927                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     6928                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     6929                                }
     6930                                else if(migration_style==SoftMigrationEnum && reCast<bool>(sheet_ungrounding[nodes[i]->Sid()])){
     6931                                        s[i]=(1-density)*h[i];
     6932                                        b[i]=-density*h[i];
     6933                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
     6934                                        nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
     6935                                }
     6936                                else{
     6937                                        if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
     6938                                }
     6939                        }
     6940                }
     6941        }
     6942
     6943        /*SubelementMigrationEnum: if one grounded, all grounded*/
     6944        if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
     6945                for(i=0;i<NUMVERTICES;i++){
     6946                        if(nodes[i]->IsGrounded()){
     6947                                groundedelement=true;
     6948                                break;
     6949                        }
     6950                }
     6951                floatingelement=!groundedelement;
     6952        }
     6953        else{
     6954                /*Otherwise: if one floating, all floating*/
     6955                for(i=0;i<NUMVERTICES;i++){
     6956                        if(nodes[i]->IsFloating()){
     6957                                floatingelement=true;
     6958                                break;
     6959                        }
     6960                }
     6961        }
     6962
     6963   /*Add basal melting rate if element just ungrounded*/
     6964        if(!this->IsFloating() && floatingelement==true){
     6965                for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
     6966                this->inputs->AddInput(new TriaInput(BasalforcingsMeltingRateEnum,&melting[0],P1Enum));
     6967        }
     6968
     6969        /*Update inputs*/
     6970   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
     6971        this->inputs->AddInput(new TriaInput(SurfaceEnum,&s[0],P1Enum));
     6972        this->inputs->AddInput(new TriaInput(BedEnum,&b[0],P1Enum));
     6973
     6974        /*Recalculate phi*/
     6975        if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
     6976                for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
     6977                this->inputs->AddInput(new TriaInput(GLlevelsetEnum,&phi[0],P1Enum));
     6978        }
     6979}
     6980/*}}}*/
     6981/*FUNCTION Tria::PotentialUngrounding{{{*/
     6982void  Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
     6983
     6984        int     i;
     6985        IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
     6986        IssmDouble  bed_hydro;
     6987        IssmDouble  rho_water,rho_ice,density;
     6988
     6989        /*material parameters: */
     6990        rho_water=matpar->GetRhoWater();
     6991        rho_ice=matpar->GetRhoIce();
     6992        density=rho_ice/rho_water;
     6993        GetInputListOnVertices(&h[0],ThicknessEnum);
     6994        GetInputListOnVertices(&ba[0],BathymetryEnum);
     6995
     6996        /*go through vertices, and figure out which ones are grounded and want to unground: */
     6997        for(i=0;i<NUMVERTICES;i++){
     6998                /*Find if grounded vertices want to start floating*/
     6999                if (!nodes[i]->IsFloating()){
     7000                        bed_hydro=-density*h[i];
     7001                        if (bed_hydro>ba[i]){
     7002                                /*Vertex that could potentially unground, flag it*/
     7003                                potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
     7004                        }
     7005                }
     7006        }
     7007}
     7008/*}}}*/
     7009/*FUNCTION Tria::UpdatePotentialUngrounding{{{*/
     7010int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
     7011
     7012        int i;
     7013        int nflipped=0;
     7014
     7015        /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
     7016        for(i=0;i<3;i++){
     7017                if (reCast<bool>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
     7018                        vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
     7019
     7020                        /*If node was not on ice shelf, we flipped*/
     7021                        if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
     7022                                nflipped++;
     7023                        }
     7024                }
     7025        }
     7026        return nflipped;
     7027}
     7028/*}}}*/
     7029#endif
     7030
     7031#ifdef _HAVE_WRAPPERS_
     7032/*FUNCTION Tria::AverageOntoPartition {{{*/
     7033void  Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
     7034
     7035        bool       already = false;
     7036        int        i,j;
     7037        int        partition[NUMVERTICES];
     7038        int        offsetsid[NUMVERTICES];
     7039        int        offsetdof[NUMVERTICES];
     7040        IssmDouble area;
     7041        IssmDouble mean;
     7042
     7043        /*First, get the area: */
     7044        area=this->GetArea();
     7045
     7046        /*Figure out the average for this element: */
     7047        this->GetVertexSidList(&offsetsid[0]);
     7048        this->GetVertexPidList(&offsetdof[0]);
     7049        mean=0;
     7050        for(i=0;i<NUMVERTICES;i++){
     7051                partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
     7052                mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
     7053        }
     7054
     7055        /*Add contribution: */
     7056        for(i=0;i<NUMVERTICES;i++){
     7057                already=false;
     7058                for(j=0;j<i;j++){
     7059                        if (partition[i]==partition[j]){
     7060                                already=true;
     7061                                break;
     7062                        }
     7063                }
     7064                if(!already){
     7065                        partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
     7066                        partition_areas->SetValue(partition[i],area,ADD_VAL);
     7067                };
     7068        }
     7069}
     7070/*}}}*/
     7071#endif

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

    r15384 r15386  
    7272                /*}}}*/
    7373                /*Element virtual functions definitions: {{{*/
    74                 void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    7574                void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
    7675                void   ComputeStrainRate(Vector<IssmDouble>* eps);
    … …  
    105104                void   DeleteResults(void);
    106105                void   MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
    107                 void   MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding);
    108106                int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
    109                 void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    110107                void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
    111108                void   RequestedOutput(int output_enum,int step,IssmDouble time);
    … …  
    117114                IssmDouble SurfaceArea(void);
    118115                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
    119                 int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
    120116                IssmDouble TimeAdapt();
    121117                int*   GetHorizontalNeighboorSids(void);
    … …  
    175171                IssmDouble SurfaceAverageVelMisfit(int weight_index);
    176172                void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
     173                #endif
     174                #ifdef _HAVE_GROUNDINGLINE_
     175                void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
     176                void   MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding);
     177                int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
     178                #endif
     179
     180                #ifdef _HAVE_WRAPPERS_
     181                void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    177182                #endif
    178183
Note: See TracChangeset for help on using the changeset viewer.