Changeset 17641


Ignore:
Timestamp:
04/03/14 17:30:29 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added GL Migration in 3d, testing required

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

Legend:

Unmodified
Added
Removed

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

    r17637 r17641  
    448448int Penta::FiniteElement(void){
    449449        return this->element_type;
     450}
     451/*}}}*/
     452/*FUNCTION Penta::FSContactMigration{{{*/
     453void Penta::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){
     454
     455        if(!IsOnBase()) return;
     456
     457        /*Intermediaries*/
     458        IssmDouble* xyz_list = NULL;
     459        IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
     460        IssmDouble  bed_normal[3];
     461        IssmDouble  epsilon[6]; /* epsilon=[exx eyy ezz exy exz eyz];*/
     462        IssmDouble  surface=0,value=0;
     463        bool grounded;
     464
     465        /* Get node coordinates and dof list: */
     466        GetVerticesCoordinates(&xyz_list);
     467
     468        /*Retrieve all inputs we will be needing: */
     469        Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
     470        Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
     471        Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
     472        Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
     473        Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
     474        Input* vz_input       = inputs->GetInput(VzEnum);       _assert_(vz_input);
     475
     476        /*Create gauss point in the middle of the basal edge*/
     477        Gauss* gauss=NewGaussBase(1);
     478        gauss->GaussPoint(0);
     479
     480        if(!IsFloating()){
     481                /*Check for basal force only if grounded and touching GL*/
     482                if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
     483                        this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     484                        this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     485                        pressure_input->GetInputValue(&pressure, gauss);
     486                        base_input->GetInputValue(&base, gauss); _assert_(base<0.);
     487
     488                        /*Compute Stress*/
     489                        IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
     490                        IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
     491                        IssmDouble sigma_zz=2.*viscosity*epsilon[2]-pressure;
     492                        IssmDouble sigma_xy=2.*viscosity*epsilon[3];
     493                        IssmDouble sigma_xz=2.*viscosity*epsilon[4];
     494                        IssmDouble sigma_yz=2.*viscosity*epsilon[5];
     495
     496                        /*Get normal vector to the bed */
     497                        NormalBase(&bed_normal[0],xyz_list);
     498
     499                        /*basalforce*/
     500                        sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + sigma_zz*bed_normal[2]*bed_normal[2]
     501                          + 2.*sigma_xy*bed_normal[0]*bed_normal[1] + 2.*sigma_xz*bed_normal[0]*bed_normal[2] + 2.*sigma_yz*bed_normal[1]*bed_normal[2];
     502
     503                        /*Compute water pressure*/
     504                        IssmDouble rho_ice   = matpar->GetRhoIce();
     505                        IssmDouble rho_water = matpar->GetRhoWater();
     506                        IssmDouble gravity   = matpar->GetG();
     507                        water_pressure=gravity*rho_water*base;
     508
     509                        /*Compare basal stress to water pressure and determine whether it should ground*/
     510                        if (sigma_nn<water_pressure) grounded=true;
     511                        else                         grounded=false;
     512                }
     513                else{
     514                        grounded=true;
     515                }
     516        }
     517        else{
     518                /*Check for basal elevation if floating*/
     519                base_input->GetInputValue(&base, gauss);
     520                bed_input->GetInputValue(&bed, gauss);
     521                if (base<bed) grounded=true;
     522                else          grounded=false;
     523        }
     524        for(int i=0;i<NUMVERTICES;i++){
     525                if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
     526                else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
     527        }
     528
     529        /*clean up*/
     530        delete gauss;
     531        xDelete<IssmDouble>(xyz_list);
    450532}
    451533/*}}}*/

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

    r17585 r17641  
    6060                void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    6161                int    FiniteElement(void);
    62                 void   FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){_error_("not implemented yet");};
     62                void   FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
    6363                void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
    6464                void   Delta18oParameterization(void);

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

    r17638 r17641  
    13441344
    13451345        if(!IsFloating()){
    1346 
     1346                /*Check for basal force only if grounded and touching GL*/
    13471347                if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
    1348                 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    1349                 this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
    1350                 pressure_input->GetInputValue(&pressure, gauss);
    1351                 base_input->GetInputValue(&base, gauss); _assert_(base<0.);
    1352 
    1353                 /*Compute Stress*/
    1354                 IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
    1355                 IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
    1356                 IssmDouble sigma_xy=2.*viscosity*epsilon[2];
    1357 
    1358                 /*Get normal vector to the bed */
    1359                 NormalBase(&bed_normal[0],xyz_list);
    1360 
    1361                 /*basalforce*/
    1362                 sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
    1363 
    1364                 /*Compute water pressure*/
    1365                 IssmDouble rho_ice   = matpar->GetRhoIce();
    1366                 IssmDouble rho_water = matpar->GetRhoWater();
    1367                 IssmDouble gravity   = matpar->GetG();
    1368                 water_pressure=gravity*rho_water*base;
    1369 
    1370                 /*Compare basal stress to water pressure and determine whether it should ground*/
    1371                 if (sigma_nn<water_pressure) grounded=true;
    1372                 else                         grounded=false;
     1348                        this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     1349                        this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
     1350                        pressure_input->GetInputValue(&pressure, gauss);
     1351                        base_input->GetInputValue(&base, gauss); _assert_(base<0.);
     1352
     1353                        /*Compute Stress*/
     1354                        IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
     1355                        IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
     1356                        IssmDouble sigma_xy=2.*viscosity*epsilon[2];
     1357
     1358                        /*Get normal vector to the bed */
     1359                        NormalBase(&bed_normal[0],xyz_list);
     1360
     1361                        /*basalforce*/
     1362                        sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
     1363
     1364                        /*Compute water pressure*/
     1365                        IssmDouble rho_ice   = matpar->GetRhoIce();
     1366                        IssmDouble rho_water = matpar->GetRhoWater();
     1367                        IssmDouble gravity   = matpar->GetG();
     1368                        water_pressure=gravity*rho_water*base;
     1369
     1370                        /*Compare basal stress to water pressure and determine whether it should ground*/
     1371                        if (sigma_nn<water_pressure) grounded=true;
     1372                        else                         grounded=false;
    13731373                }
    13741374                else{
    … …  
    13771377        }
    13781378        else{
     1379                /*Check for basal elevation if floating*/
    13791380                base_input->GetInputValue(&base, gauss);
    13801381                bed_input->GetInputValue(&bed, gauss);

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

    r17579 r17641  
    9494                }
    9595                else{
    96                         _error_("not supported");
     96                        _error_("not supported (vertexfloating="<<serial_vertexfloating[i]<<" vertexgrounded="<<serial_vertexgrounded[i]<<")");
    9797                }
    9898        }
Note: See TracChangeset for help on using the changeset viewer.