Changeset 22841 for issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

Timestamp:

06/13/18 02:10:04 (7 years ago)

Author:

bdef

Message:

NEW:Minor Hydro bug and fixing coupling with friction

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (17 diffs)

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

@@ -74 +74 @@
         bool   isefficientlayer;
         int    hydrology_model;
-        
+
         /*Fetch data needed: */
         iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
@@ -161 +161 @@
                                 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
                                 loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
-                        }       
+                        }
                 }
         }
@@ -233 +233 @@
                 active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
         }
-        
+
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=basalelement->NewGauss(2);
@@ -250 +250 @@
                 D[0][0]=D_scalar;
                 D[1][1]=D_scalar;
-                GetB(B,basalelement,xyz_list,gauss); 
+                GetB(B,basalelement,xyz_list,gauss);
                 TripleMultiply(B,2,numnodes,1,
                                                                          &D[0][0],2,2,0,
@@ -265 +265 @@
                                                                                  basis,1,numnodes,0,
                                                                                  &Ke->values[0],1);
-                        
+
                         /*Transfer EPL part*/
                         if(isefficientlayer){
@@ -288 +288 @@
         delete gauss;
         if(domaintype!=Domain2DhorizontalEnum){
-                basalelement->DeleteMaterials(); 
+                basalelement->DeleteMaterials();
                 delete basalelement;
         }
@@ -412 +412 @@
         delete gauss;
         if(domaintype!=Domain2DhorizontalEnum){
-                basalelement->DeleteMaterials(); 
+                basalelement->DeleteMaterials();
                 delete basalelement;
         }
@@ -419 +419 @@
 
 void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
          * For node i, Bi can be expressed in the actual coordinate system
-         * by: 
+         * by:
          *       Bi=[ dN/dx ]
          *          [ dN/dy ]
@@ -507 +507 @@
                 IssmDouble rho_ice        = basalelement->GetMaterialParameter(MaterialsRhoIceEnum);
                 IssmDouble g              = basalelement->GetMaterialParameter(ConstantsGEnum);
-                
+
                 basalelement->GetInputListOnVertices(thickness,ThicknessEnum);
                 basalelement->GetInputListOnVertices(base,BaseEnum);
@@ -542 +542 @@
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){
-                basalelement->DeleteMaterials(); 
+                basalelement->DeleteMaterials();
                 delete basalelement;
         }
@@ -554 +554 @@
 IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/
         int unconf_scheme;
-        IssmDouble expfac; 
+        IssmDouble expfac;
         IssmDouble sediment_storing;
         IssmDouble storing,yield;
@@ -610 +610 @@
                 sed_head_input->GetInputValue(&prestep_head,gauss);
                 water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
-                
+
                 //min definition of the if test
                 sediment_transmitivity=FullLayer_transmitivity*min(water_sheet/sediment_thickness,1.);
@@ -625 +625 @@
 
 void  HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
-        
+
         int        hmax_flag;
         IssmDouble h_max;
@@ -632 +632 @@
         /*Get the flag to the limitation method*/
         element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
-        
+
         /*Switch between the different cases*/
         switch(hmax_flag){
@@ -674 +674 @@
         case 1:
                 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
-                transfer=+leakage; 
+                transfer=+leakage;
                 break;
         default:
@@ -698 +698 @@
                 _assert_(epl_head_input);
                 epl_head_input->GetInputValue(&epl_head,gauss);
-                if (element->Id()==42){
-                        _printf_("epl head in sed Pvec transfer is "<<  epl_head <<"\n");
-                }
                 element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
                 transfer=+epl_head*leakage;
@@ -718 +715 @@
                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
                         Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
-                
+
                 if(node_mask_input->Max()>0.){
                         element_active = true;