Changeset 21140

Timestamp:

08/16/16 09:48:33 (9 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed pressure in HO

File:

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

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

@@ -1807 +1807 @@
                         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]);
+                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+1]);
                         dim=1;
                         break;
@@ -2777 +2777 @@
 void           StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
-        int         i,dim;
+        int         i,domaintype,dim;
         int*        doflist=NULL;
         IssmDouble* xyz_list=NULL;
-
-        /*Get mesh dimension*/
-        element->FindParam(&dim,DomainDimensionEnum);
 
         /*Deal with pressure first*/
@@ -2788 +2785 @@
         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->FindParam(&domaintype,DomainTypeEnum);
+        IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble g       =element->GetMaterialParameter(ConstantsGEnum);
+        switch(domaintype){
+                case Domain3DEnum:
+                        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]);
+                        dim=3;
+                        break;
+                case Domain2DverticalEnum:
+                        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+1]);
+                        dim=2;
+                        break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
         element->AddInput(PressureEnum,pressure,P1Enum);
         xDelete<IssmDouble>(pressure);
@@ -2843 +2854 @@
 
         /*Add vx and vy as inputs to the element: */
-        //element->AddInput(VxEnum,  vx,element->GetElementType());
-        //element->AddInput(VyEnum,  vy,element->GetElementType());
-        //element->AddInput(VelEnum,vel,element->GetElementType());
         element->AddInput(VxEnum,vx,element->GetElementType());
         if(dim==3)element->AddInput(VyEnum,vy,element->GetElementType());