Changeset 16780

Timestamp:

11/15/13 10:20:45 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed pressure HO issue for higher order element (need to be on vertices and not nodes)

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)

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

@@ -1060 +1060 @@
 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*/
@@ -1076 +1088 @@
         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: */
@@ -1084 +1094 @@
         /*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: */
@@ -1099 +1108 @@
         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 
@@ -1118 +1119 @@
         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);