Changeset 6969

Timestamp:

01/06/11 11:13:40 (14 years ago)

Author:

seroussi

Message:

update bed and surface on elements

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (3 diffs)
• Tria.cpp (modified) (2 diffs)

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -4737 +4737 @@
         for(i=0;i<numdof2d;i++){
                 values[i]         =solution[doflist[i]];
+                if(isnan(values[i])) _error_("NaN found in solution vector");
+                /*Constrain thickness to be at least 1m*/
+                if(values[i]<1) values[i]=1;
                 values[i+numdof2d]=values[i];
-                if(isnan(values[i])) _error_("NaN found in solution vector");
         }
 
@@ -4749 +4751 @@
         surface_input->GetValuesPtr(&surface_ptr,&dummy);
 
-        /*If shelf: hydrostatic equilibrium*/
-        if (this->IsOnShelf()){
-
-                /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
-                this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
-
-                /*recover material parameters: */
-                rho_ice=matpar->GetRhoIce();
-                rho_water=matpar->GetRhoWater();
-
-                if(hydroadjustment==AbsoluteEnum){
-                        for(i=0;i<numdof;i++) {
-                                surface[i]=values[i]*(1-rho_ice/rho_water);
-                                bed[i]=-values[i]*rho_ice/rho_water;
+        /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
+        this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
+
+        /*recover material parameters: */
+        rho_ice=matpar->GetRhoIce();
+        rho_water=matpar->GetRhoWater();
+
+        for(i=0;i<numdof;i++) {
+                /*If shelf: hydrostatic equilibrium*/
+                if (this->nodes[i]->IsOnSheet()){
+                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
+                        bed[i]=bed_ptr[i]; //bed does not change
+                }
+                else{ //so it is an ice shelf
+                        if(hydroadjustment==AbsoluteEnum){
+                                        surface[i]=values[i]*(1-rho_ice/rho_water);
+                                        bed[i]=values[i]*(-rho_ice/rho_water);
                         }
-                }
-                else if(hydroadjustment==IncrementalEnum){
-                        for(i=0;i<numdof;i++) {
-                                surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH 
-                                bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
+                        else if(hydroadjustment==IncrementalEnum){
+                                        surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH 
+                                        bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
                         }
-                }
-                else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
-        }
-
-        /*If sheet: surface = bed + thickness*/
-        else{
-                /*Now Compute surface only*/
-                for(i=0;i<numdof;i++) {
-                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
+                        else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
                 }
         }
@@ -4788 +4783 @@
                 penta->inputs->AddInput(new PentaVertexInput(ThicknessEnum,values));
                 penta->inputs->AddInput(new PentaVertexInput(SurfaceEnum,surface));
-                if (penta->IsOnShelf()) penta->inputs->AddInput(new PentaVertexInput(BedEnum,bed));
+                penta->inputs->AddInput(new PentaVertexInput(BedEnum,bed));
 
                 /*Stop if we have reached the surface*/

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -3958 +3958 @@
                 values[i]=solution[doflist[i]];
                 if(isnan(values[i])) _error_("NaN found in solution vector");
+                /*Constrain thickness to be at least 1m*/
+                if(values[i]<1) values[i]=1;
         }
 
@@ -3968 +3970 @@
         surface_input->GetValuesPtr(&surface_ptr,&dummy);
 
-        /*If shelf: hydrostatic equilibrium*/
-        if (this->IsOnShelf()){
-
-                /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
-                this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
-
-                /*recover material parameters: */
-                rho_ice=matpar->GetRhoIce();
-                rho_water=matpar->GetRhoWater();
-
-                if(hydroadjustment==AbsoluteEnum){
-                        for(i=0;i<numdof;i++) {
+        /*Fing HydrostaticAdjustment to figure out how to update the geometry:*/
+        this->parameters->FindParam(&hydroadjustment,HydrostaticAdjustmentEnum);
+
+        /*recover material parameters: */
+        rho_ice=matpar->GetRhoIce();
+        rho_water=matpar->GetRhoWater();
+
+        for(i=0;i<numdof;i++) {
+                /*If shelf: hydrostatic equilibrium*/
+                if (this->nodes[i]->IsOnSheet()){
+                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
+                }
+                else{ //this is an ice shelf
+
+                        if(hydroadjustment==AbsoluteEnum){
                                 surface[i]=values[i]*(1-rho_ice/rho_water);
-                                bed[i]=-values[i]*rho_ice/rho_water;
                         }
-                }
-                else if(hydroadjustment==IncrementalEnum){
-                        for(i=0;i<numdof;i++) {
+                        else if(hydroadjustment==IncrementalEnum){
                                 surface[i]=surface_ptr[i]+(1.0-rho_ice/rho_water)*(values[i]-thickness_ptr[i]); //surface = oldsurface + (1-di) * dH 
                                 bed[i]=bed_ptr[i]-rho_ice/rho_water*(values[i]-thickness_ptr[i]); //bed = oldbed + di * dH
                         }
-                }
-                else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
-        }
-
-        /*If sheet: surface = bed + thickness*/
-        else{
-                /*Now Compute surface only*/
-                for(i=0;i<numdof;i++) {
-                        surface[i]=bed_ptr[i]+values[i]; //surface=oldbed+newthickness
+                        else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToString(hydroadjustment));
                 }
         }