Changeset 22308 for issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

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

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (4 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])));