Changeset 14359

Timestamp:

02/21/13 08:54:44 (12 years ago)

Author:

seroussi

Message:

CHG: started work on grounding line

Location:

issm/trunk-jpl/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/classes/objects/Elements/Tria.cpp (modified) (8 diffs)
• c/classes/objects/Elements/Tria.h (modified) (1 diff)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (2 diffs)
• m/enum/EnumDefinitions.py (modified) (2 diffs)
• m/enum/MaximumNumberOfEnums.m (modified) (1 diff)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h ¶

@@ -490 +490 @@
         SoftMigrationEnum,
         SubelementMigrationEnum,
+        GLlevelsetEnum,
         /*}}}*/
         /*Solver{{{1*/

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

@@ -1113 +1113 @@
         return &this->horizontalneighborsids[0];
 
+}
+/*}}}*/
+/*FUNCTION Tria::GetGroundedPortion{{{*/
+IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){
+        /*Computeportion of the element that is grounded*/ 
+
+        bool mainlyfloating=false;
+        IssmDouble phi,s1,s2,area_init,area_grounded;
+        IssmDouble gl[3];
+        IssmDouble xyz_bis[3][3];
+        GaussTria *gauss    = NULL;
+
+        /*Recover parameters and values*/
+        GetInputListOnVertices(&gl[0],GLlevelsetEnum);
+        
+        /*Check that not all nodes are grounded or floating*/
+        if(gl[0]>0 && gl[1]>0 && gl[2]>0) _error_("Error. Three nodes of the element are grounded, should not compute the fraction of grounded element");
+        if(gl[0]<0 && gl[1]<0 && gl[2]<0) _error_("Error. Three nodes of the element are floating, should not compute the fraction of grounded element");
+
+        /*Figure out if two nodes are floating or grounded*/
+        if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=true;
+
+        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                /*Coordinates of point 2: same as initial point 2*/
+                xyz_bis[2][0]=*(xyz_list+3*2+0);
+                xyz_bis[2][1]=*(xyz_list+3*2+1);
+                xyz_bis[2][2]=*(xyz_list+3*2+2);
+
+                /*Portion of the segments*/
+                s1=gl[2]/(gl[2]-gl[1]);
+                s2=gl[2]/(gl[2]-gl[0]);
+                
+                /*New point 1*/
+                xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
+                xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
+                xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
+
+                /*New point 0*/
+                xyz_bis[0][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
+                xyz_bis[0][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
+                xyz_bis[0][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
+        }
+        else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                /*Coordinates of point 0: same as initial point 2*/
+                xyz_bis[0][0]=*(xyz_list+3*0+0);
+                xyz_bis[0][1]=*(xyz_list+3*0+1);
+                xyz_bis[0][2]=*(xyz_list+3*0+2);
+
+                /*Portion of the segments*/
+                s1=gl[0]/(gl[0]-gl[1]);
+                s2=gl[0]/(gl[0]-gl[2]);
+                
+                /*New point 1*/
+                xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
+                xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
+                xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
+
+                /*New point 2*/
+                xyz_bis[2][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
+                xyz_bis[2][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
+                xyz_bis[2][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
+        }
+        else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                /*Coordinates of point 1: same as initial point 2*/
+                xyz_bis[1][0]=*(xyz_list+3*1+0);
+                xyz_bis[1][1]=*(xyz_list+3*1+1);
+                xyz_bis[1][2]=*(xyz_list+3*1+2);
+
+                /*Portion of the segments*/
+                s1=gl[1]/(gl[1]-gl[0]);
+                s2=gl[1]/(gl[1]-gl[2]);
+                
+                /*New point 0*/
+                xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
+                xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
+                xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
+
+                /*New point 2*/
+                xyz_bis[2][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
+                xyz_bis[2][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
+                xyz_bis[2][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
+        }
+
+        /*Compute fraction of grounded element*/
+        GetJacobianDeterminant2d(&area_init, xyz_list,gauss);
+        GetJacobianDeterminant2d(&area_grounded, &xyz_bis[0][0],gauss);
+        if(mainlyfloating==true) area_grounded=area_init-area_grounded;
+        phi=area_grounded/area_init;
+        if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
+        return phi;
 }
 /*}}}*/
@@ -2035 +2125 @@
         int     i,migration_style;
         bool    elementonshelf = false;
+        bool    elementonsheet = false;
         IssmDouble  bed_hydro,yts,gl_melting_rate;
         IssmDouble  rho_water,rho_ice,density;
-        IssmDouble  melting[NUMVERTICES];
+        IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];;
         IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
 
@@ -2048 +2139 @@
         GetInputListOnVertices(&b[0],BedEnum);
         GetInputListOnVertices(&ba[0],BathymetryEnum);
+        if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
         rho_water=matpar->GetRhoWater();
         rho_ice=matpar->GetRhoIce();
@@ -2087 +2179 @@
         /*If at least one vertex is now floating, the element is now floating*/
         for(i=0;i<NUMVERTICES;i++){
-                if(nodes[i]->IsFloating()){
-                        elementonshelf=true;
-                        break;
+                if(migration_style==SubelementMigrationEnum){
+                        if(nodes[i]->IsGrounded()){
+                                elementonsheet=true;
+                                break;
+                        }
+                        phi[i]=h[i]+ba[i]/density;
+                        elementonshelf=!elementonsheet;
+                }
+                else{
+                        if(nodes[i]->IsFloating()){
+                                elementonshelf=true;
+                                break;
+                        }
                 }
         }
@@ -2101 +2203 @@
         /*Update inputs*/
    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
-
-        /*Update inputs*/    
         this->inputs->AddInput(new TriaP1Input(SurfaceEnum,&s[0]));
         this->inputs->AddInput(new TriaP1Input(BedEnum,&b[0]));
+        if(migration_style==SubelementMigrationEnum) this->inputs->AddInput(new TriaP1Input(GLlevelsetEnum,&phi[0]));
 }
 /*}}}*/
@@ -2995 +3096 @@
         IssmDouble MOUNTAINKEXPONENT = 10;
         IssmDouble slope_magnitude,alpha2;
+        IssmDouble migration_style;
         IssmDouble Jdet;
+        IssmDouble phi=1.0;
         IssmDouble L[2][numdof];
         IssmDouble DL[2][2]  = {{ 0,0 },{0,0}};
@@ -3009 +3112 @@
 
         /*Retrieve all inputs and parameters*/
+        parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-        Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
-        Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
+        Input* surface_input=inputs->GetInput(SurfaceEnum);       _assert_(surface_input);
+        Input* vx_input=inputs->GetInput(VxEnum);                 _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);                 _assert_(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);                 _assert_(vz_input);
         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
         /*build friction object, used later on: */
         friction=new Friction("2d",inputs,matpar,analysis_type);
+
+        /*Recover portion of element that is grounded*/
+        if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);
 
         /* Start  looping on the number of gaussian points: */
@@ -3031 +3138 @@
                 if(slope_magnitude>MAXSLOPE) alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT);
                 else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
+                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
 
                 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);

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

@@ -196 +196 @@
                 void           GetVertexSidList(int* sidlist);
                 void           GetConnectivityList(int* connectivity);
+                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
                 void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
                 void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp ¶

@@ -471 +471 @@
                 case SoftMigrationEnum : return "SoftMigration";
                 case SubelementMigrationEnum : return "SubelementMigration";
+                case GLlevelsetEnum : return "GLlevelset";
                 case StokesSolverEnum : return "StokesSolver";
                 case AdjointEnum : return "Adjoint";

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp ¶

@@ -481 +481 @@
               else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
               else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
+              else if (strcmp(name,"GLlevelset")==0) return GLlevelsetEnum;
               else if (strcmp(name,"StokesSolver")==0) return StokesSolverEnum;
               else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
@@ -506 +507 @@
               else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
               else if (strcmp(name,"Regular")==0) return RegularEnum;
-              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"Separate")==0) return SeparateEnum;
+              if (strcmp(name,"Scaled")==0) return ScaledEnum;
+              else if (strcmp(name,"Separate")==0) return SeparateEnum;
               else if (strcmp(name,"Sset")==0) return SsetEnum;
               else if (strcmp(name,"Verbose")==0) return VerboseEnum;

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

@@ -4549 +4549 @@
         return StringToEnum('SubelementMigration')[0]
 
+def GLlevelsetEnum():
+        """
+        GLLEVELSETENUM - Enum of GLlevelset
+
+           Usage:
+              macro=GLlevelsetEnum()
+        """
+
+        return StringToEnum('GLlevelset')[0]
+
 def StokesSolverEnum():
         """
@@ -4987 +4997 @@
         """
 
-        return 497
-
+        return 498
+

issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m ¶

@@ -9 +9 @@
 %      macro=MaximumNumberOfEnums()
 
-macro=497;
+macro=498;