Changeset 17633

Timestamp:

04/03/14 13:51:57 (11 years ago)

Author:

hongjuy

Message:

NEW: first version of contact GL criterion

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Element.cpp (modified) (1 diff)
• Tria.cpp (modified) (2 diffs)

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

@@ -670 +670 @@
                 else shelf=false;
         }
-        else if(migration_style==NoneEnum || migration_style==AgressiveMigrationEnum || migration_style==SoftMigrationEnum){ //Floating if all nodes are floating
+        else if(migration_style==NoneEnum || migration_style==AgressiveMigrationEnum || migration_style==SoftMigrationEnum || migration_style==ContactEnum){ //Floating if all nodes are floating
                 if(this->inputs->Min(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
                 else shelf=true;

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

@@ -1323 +1323 @@
         /*Intermediaries*/
         IssmDouble* xyz_list = NULL;
-
-        /*Figure out wether element is grounded or floating*/
-        bool grounded = true;
+        IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
+        IssmDouble  bed_normal[2];
+        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+        IssmDouble  surface=0,value=0;
+        bool grounded;
+
+        IssmDouble* basis = xNew<IssmDouble>(NUMVERTICES);
+
 
         /* Get node coordinates and dof list: */
-        ::GetVerticesCoordinates(xyz_list,vertices,NUMVERTICES);
+        GetVerticesCoordinates(&xyz_list);
 
         /*Retrieve all inputs we will be needing: */
-        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
-        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
-
-        _error_("STOP");
-
+        Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
+        Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+        Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
+        Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
+        Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
+
+        Gauss* gauss=NewGaussBase(1);
+        gauss->GaussPoint(0);
+
+        if(!IsFloating()){ 
+
+                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+                this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
+                pressure_input->GetInputValue(&pressure, gauss);
+                base_input->GetInputValue(&base, gauss); _assert_(base<0.);
+
+                /*Compute Stress*/
+                IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
+                IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
+                IssmDouble sigma_xy=2.*viscosity*epsilon[2];
+
+                /*Get normal vector to the bed */
+                NormalBase(&bed_normal[0],xyz_list);
+
+                /*basalforce*/
+                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];
+
+                /*Compute water pressure*/
+                IssmDouble rho_ice   = matpar->GetRhoIce();
+                IssmDouble rho_water = matpar->GetRhoWater();
+                IssmDouble gravity   = matpar->GetG();
+                water_pressure=gravity*rho_water*base;
+
+                /*Compare basal stress to water pressure and determine whether it should ground*/
+                if (sigma_nn<water_pressure) grounded=true;
+                else                         grounded=false;
+        }
+        else{
+                base_input->GetInputValue(&base, gauss);
+                bed_input->GetInputValue(&bed, gauss);
+                if (base<bed) grounded=true;
+                else          grounded=false;
+        }
         for(int i=0;i<NUMVERTICES;i++){
                 if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
@@ -1343 +1385 @@
 
         /*clean up*/
+        delete gauss;
         xDelete<IssmDouble>(xyz_list);
-
 }
 /*}}}*/