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

Timestamp:

01/05/15 11:21:22 (10 years ago)

Author:

bdef

Message:

minor clean-up in Hydro Code

File:

• issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (modified) (22 diffs)

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

@@ -24 +24 @@
         iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 
+        /*If not return*/
         if(!isefficientlayer) return;
+
+        /*If yes, initialize a flip flop counter*/
         iomodel->FetchData(&eplflip_lock,HydrologydcEplflipLockEnum);
- 
         parameters->AddObject(new IntParam(HydrologydcEplflipLockEnum,eplflip_lock));
         
-        /*Nothing for now*/
 }/*}}}*/
 void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -35 +36 @@
         bool   isefficientlayer;
         int    hydrology_model;
-
 
         /*Now, do we really want DC?*/
@@ -61 +61 @@
         iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum);
         iomodel->FetchDataToInput(elements,HydrologydcEplMaxThicknessEnum);
-        iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
         iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum);
-        if(iomodel->domaintype!=Domain2DhorizontalEnum){
+                if(iomodel->domaintype!=Domain2DhorizontalEnum){
                 iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
                 iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
@@ -85 +84 @@
         iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
@@ -98 +98 @@
 
         IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
-
-}/*}}}*/
+}/*}}}*/
+
 void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
         /*Nothing for now*/
@@ -131 +131 @@
 }/*}}}*/
 
- ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
+ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
         /*Default, return NULL*/
         return NULL;
 }/*}}}*/
 
- ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
+ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 _error_("Not implemented");
 }/*}}}*/
 
- ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
+ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
         /*Intermediaries*/
@@ -244 +244 @@
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         return Ke;
-
-}/*}}}*/
+}/*}}}*/
+
 ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/
 
@@ -348 +348 @@
         return pe;
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         element->GetSolutionFromInputsOneDof(solution,EplHeadEnum);
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
         _error_("Not implemented yet");
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
@@ -375 +378 @@
         /*Fetch number of nodes for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
-
 
         /*Fetch dof list and allocate solution vector*/
@@ -401 +403 @@
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 } /*}}}*/
+
 void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
         /*Default, do nothing*/
@@ -415 +418 @@
         return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));                 
 }/*}}}*/
+
 IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/
         IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
@@ -535 +539 @@
                 _error_("not Implemented Yet");
                 }
-                        
+                
                 int         numnodes      = element->GetNumberOfNodes();
                 IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
@@ -602 +606 @@
                 xDelete<IssmDouble>(bed);
         }
-}
-/*}}}*/
+}/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
@@ -631 +635 @@
         xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+
 void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
 
@@ -655 +660 @@
         /*Intermediaries*/
 
-        int         penalty_lock;
+        int         eplflip_lock;
         int         numnodes      =basalelement->GetNumberOfNodes();
         IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
@@ -669 +674 @@
         active_element_input->GetInputValue(&active_element);
 
-        basalelement->parameters->FindParam(&penalty_lock,HydrologydcEplflipLockEnum); 
+        basalelement->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
 
         basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); 
@@ -691 +696 @@
                 else if(old_active[i]>0.){
                         vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
-                        /*If epl thickness gets under 10-3 initial thickness, close the layer*/
+                        /*If epl thickness gets under colapse thickness, close the layer*/
                         if(epl_thickness[i]<colapse_thick){
                                 eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
-                                if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<penalty_lock |penalty_lock==0){
+                                /*Avoid flipfloping between open and closed states*/
+                                if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
                                         vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
                                         epl_thickness[i]=init_thick;
@@ -705 +711 @@
                 /*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
                 GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
-
                 if(eplhead[i]>=h_max && active_element){
                         for(j=0;j<numnodes;j++){
@@ -729 +734 @@
 }
 /*}}}*/
+
 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
         /*Constants*/
@@ -750 +756 @@
         IssmDouble  flag     = 0.;
         IssmDouble* active   = xNew<IssmDouble>(numnodes);
-                
+
+        /*Pass the activity mask from elements to nodes*/
         basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
         bool active_element;
@@ -758 +765 @@
         for(int i=0;i<numnodes;i++) flag+=active[i];
 
+        /*If any node is active all the node in the element are active*/
         if(flag>0.){
                 for(int i=0;i<numnodes;i++){
@@ -763 +771 @@
                 }
         }
+        /*If the element is active all its nodes are active*/
         else if(active_element){
                 for(int i=0;i<numnodes;i++){