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

Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 26289)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 26290)
@@ -193,13 +193,11 @@
                                        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,7 +3063,7 @@
                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,8 +3088,7 @@
 
        /*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,17 +3121,15 @@
                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;