source: issm/oecreview/Archive/25834-26739/ISSM-26289-26290.diff

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -193 +193 @@
                                         IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,1);
                                 }
                         }
-                        else{//itapopo FIXME only base velocities are constrained here
+                        else{//MLHO itapopo FIXME base and shear velocities are constrained
                                 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
-                                //IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
-                                if(iomodel->domaintype!=Domain2DverticalEnum){
-                                        IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
-                                        //IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
-                                }
+                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
+                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
+                                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
                         }
                 }
 
@@ -3065 +3063 @@
                 thickness_input->GetInputValue(&thickness,gauss);
                 surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
                 n_input->GetInputValue(&n,gauss);
-
+                
                 for(int i=0;i<numnodes;i++){// per node: vx (basal vx), vshx, vy (basal vy), vshy
                         pe->values[i*4+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1
                         pe->values[i*4+1]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F2
@@ -3090 +3088 @@
 
         /*Intermediaries*/
         IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
-        IssmDouble  water_pressure_sh,ice_pressure_sh,pressure_sh,n,s,b;
-        IssmDouble  surface_under_water,base_under_water,pressure;
+        IssmDouble  water_pressure_sh,ice_pressure_sh,pressure_sh,pressure,n,s,b;
         IssmDouble *xyz_list = NULL;
         IssmDouble *xyz_list_front = NULL;
         IssmDouble  normal[2];
@@ -3124 +3121 @@
                 sealevel_input->GetInputValue(&sealevel,gauss);
                 element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
                 element->NodalFunctions(basis,gauss);
-
-                surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level
-                base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
+                
+                b = base-sealevel; // ice base shifted by the sea level
+                s = thickness+b;   // ice surface shifted by the sea level
                 /*Vertically integrated pressure - SSA type*/
-                water_pressure = (1.0/2.0)*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
-                ice_pressure   = (1.0/2.0)*gravity*rho_ice*thickness*thickness;
+                water_pressure = -(1.0/2.0)*gravity*rho_water*(b*b);
+                ice_pressure   =  (1.0/2.0)*gravity*rho_ice*thickness*thickness;
                 pressure = ice_pressure + water_pressure;
                 /*Vertically integrated pressure - HO type*/
-                b=min(0.,base-sealevel); // 0 if the bottom of the glacier is above water level
-                s=thickness+b; // ice surface regardless of whether the top of the glacier is above water level or not
-                water_pressure_sh = gravity*rho_water*( -b*b/2 + (pow(-b,n+2)/pow(thickness,n+1))*(-s/(n+2)-b/(n+3)) );
+                water_pressure_sh = gravity*rho_water*( -b*b/2 - (s*thickness/(n+2))*(1-pow(s/thickness,n+2)) + (thickness*thickness/(n+3))*(1-pow(s/thickness,n+3)) );
                 ice_pressure_sh   = gravity*rho_ice*thickness*thickness*(n+1)/(2*(n+3));
                 pressure_sh = ice_pressure_sh + water_pressure_sh;