Changeset 22308

Timestamp:

12/22/17 10:10:23 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: do not call element->GetMaterialParameter(MaterialsRheologyBEnum) anymore since this is NOT a constant

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCEfficientAnalysis.cpp (modified) (4 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (3 diffs)
• classes/Elements/Tria.cpp (modified) (6 diffs)

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

@@ -522 +522 @@
         int         iseplthickcomp;
         int         domaintype;
-        IssmDouble  dt,A,B;
+        IssmDouble  dt,A;
         IssmDouble  EPLgrad2;
         IssmDouble  EPL_N;
@@ -534 +534 @@
                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
                 
-                switch(domaintype){
-                case Domain2DhorizontalEnum:
-                        if(!element->IsOnBase()) return;        
-                        B = element->GetMaterialParameter(MaterialsRheologyBbarEnum);
-                        break;
-                case Domain3DEnum:
-                        B = element->GetMaterialParameter(MaterialsRheologyBEnum);
-                        break;
-                default:
-                _error_("not Implemented Yet");
-                }
+                /*skip element if 3d and not on base*/
+                if(domaintype==Domain3DEnum && !element->IsOnBase()) continue;
                 
                 int         numnodes      = element->GetNumberOfNodes();
                 IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
+                IssmDouble* B             = xNew<IssmDouble>(numnodes);
                 IssmDouble* eplhead       = xNew<IssmDouble>(numnodes);
                 IssmDouble* epl_slopeX    = xNew<IssmDouble>(numnodes);
@@ -569 +561 @@
                 IssmDouble init_thick       =   element->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
                 IssmDouble max_thick        =   element->GetMaterialParameter(HydrologydcEplMaxThicknessEnum);
-                
-                A=pow(B,-n);
+
+                switch(domaintype){
+                        case Domain2DhorizontalEnum: element->GetInputListOnVertices(&B[0],MaterialsRheologyBbarEnum); break;
+                        case Domain3DEnum:           element->GetInputListOnVertices(&B[0],MaterialsRheologyBEnum); break;
+                        default: _error_("not Implemented Yet");
+                }
                 
                 element->GetInputListOnVertices(&eplhead[0],EplHeadEnum);
@@ -587 +583 @@
                 else{
                         for(int i=0;i<numnodes;i++){
+                                A=pow(B[i],-n);
                                 /*Compute first the effective pressure in the EPL*/
                                 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));

issm/trunk-jpl/src/c/classes/Elements/Element.cpp ¶

@@ -1132 +1132 @@
 
         _assert_(this->matpar);
-        switch(enum_in){ // FIXME: change this to material and replace NULL by gauss
+        switch(enum_in){ // FIXME: change this to material
                 case MaterialsRheologyNEnum:
                         return this->material->GetN();
-                case MaterialsRheologyBEnum:
-                        return this->material->GetB(NULL);
-                case MaterialsRheologyBbarEnum:
-                        return this->material->GetBbar(NULL);
                 default:
                         return this->matpar->GetMaterialParameter(enum_in);

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶

@@ -192 +192 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-        IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
+        IssmDouble  B,sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
         IssmDouble  epse_2,groundedice,bed;
 
@@ -208 +208 @@
         Input* gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
         Input* bs_input = inputs->GetInput(BaseEnum);                    _assert_(bs_input);
+        Input* B_input  = inputs->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
         Input* smax_fl_input = inputs->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
         Input* smax_gr_input = inputs->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
-        IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
         IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
 
@@ -219 +219 @@
 
                 /*Get velocity components and thickness*/
+                B_input->GetInputValue(&B,gauss);
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶

@@ -217 +217 @@
         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
         IssmDouble  sigma_vm[NUMVERTICES];
-        IssmDouble  sigma_max,sigma_max_floating,sigma_max_grounded;
+        IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded;
         IssmDouble  epse_2,groundedice,bed;
 
@@ -226 +226 @@
         Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
         Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
+        Input* B_input  = inputs->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
         Input* gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
         Input* bs_input = inputs->GetInput(BaseEnum);                    _assert_(bs_input);
         Input* smax_fl_input = inputs->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
         Input* smax_gr_input = inputs->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
-        IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
         IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
 
@@ -239 +239 @@
 
                 /*Get velocity components and thickness*/
+                B_input->GetInputValue(&B,gauss);
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
@@ -301 +302 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  vx,vy,vel;
-        IssmDouble  critical_fraction,water_height;
+        IssmDouble  rheology_B,critical_fraction,water_height;
         IssmDouble  bathymetry,Ho,thickness,float_depth,groundedice;
         IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
@@ -318 +319 @@
         IssmDouble rho_freshwater = this->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
         IssmDouble constant_g     = this->GetMaterialParameter(ConstantsGEnum);
-        IssmDouble rheology_B     = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
         IssmDouble rheology_n     = this->GetMaterialParameter(MaterialsRheologyNEnum);
-        
+
+        Input*   B_input       = inputs->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
         Input*   H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
         Input*   bed           = inputs->GetInput(BedEnum); _assert_(bed);
@@ -336 +337 @@
                 gauss->GaussVertex(iv);
         
+                B_input->GetInputValue(&rheology_B,gauss);
                 H_input->GetInputValue(&thickness,gauss);
                 bed->GetInputValue(&bathymetry,gauss);