source: issm/oecreview/Archive/16554-17801/ISSM-16779-16780.diff@ 17802

Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16779)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16780)
@@ -1059,11 +1059,23 @@
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
-       int         i,meshtype;
-       IssmDouble  rho_ice,g;
+       int         i;
        int*        doflist=NULL;
        IssmDouble* xyz_list=NULL;
 
+       /*Deal with pressure first*/
+       int numvertices = element->GetNumberOfVertices();
+       IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
+       IssmDouble* surface   = xNew<IssmDouble>(numvertices);
+       IssmDouble  rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+       IssmDouble  g         = element->GetMaterialParameter(ConstantsGEnum);
+       element->GetVerticesCoordinates(&xyz_list);
+       element->GetInputListOnVertices(surface,SurfaceEnum);
+       for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+       element->AddInput(PressureEnum,pressure,P1Enum);
+       xDelete<IssmDouble>(pressure);
+       xDelete<IssmDouble>(surface);
+
        /*Fetch number of nodes and dof for this finite element*/
        int numnodes = element->GetNumberOfNodes();
        int numdof   = numnodes*2;
@@ -1075,15 +1087,12 @@
        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
-       IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
-       IssmDouble* surface   = xNew<IssmDouble>(numnodes);
 
        /*Use the dof list to index into the solution vector: */
        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
        /*Transform solution in Cartesian Space*/
        element->TransformSolutionCoord(&values[0],XYEnum);
-       element->FindParam(&meshtype,MeshTypeEnum);
 
        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
        for(i=0;i<numnodes;i++){
@@ -1099,14 +1108,6 @@
        element->GetInputListOnNodes(&vz[0],VzEnum,0.);
        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
-       /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
-        *so the pressure is just the pressure at the bedrock: */
-       rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
-       g       = element->GetMaterialParameter(ConstantsGEnum);
-       element->GetVerticesCoordinates(&xyz_list);
-       element->GetInputListOnNodes(&surface[0],SurfaceEnum);
-       for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
-
        /*Now, we have to move the previous Vx and Vy inputs  to old 
         * status, otherwise, we'll wipe them off: */
        element->InputChangeName(VxEnum,VxPicardEnum);
@@ -1117,11 +1118,8 @@
        element->AddInput(VxEnum,vx,P1Enum);
        element->AddInput(VyEnum,vy,P1Enum);
        element->AddBasalInput(VelEnum,vel,P1Enum);
-       element->AddBasalInput(PressureEnum,pressure,P1Enum);
 
        /*Free ressources:*/
-       xDelete<IssmDouble>(surface);
-       xDelete<IssmDouble>(pressure);
        xDelete<IssmDouble>(vel);
        xDelete<IssmDouble>(vz);
        xDelete<IssmDouble>(vy);