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

Timestamp:

12/07/16 02:45:35 (8 years ago)

Author:

bdef

Message:

BUG. Fixing the Epl initialisation which was erased for non active nodes

File:

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

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

@@ -151 +151 @@
 
         int*     eplzigzag_counter=NULL;
-        Element* element=NULL;
         femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
         for(int i=0;i<femmodel->nodes->Size();i++){
@@ -381 +380 @@
                 pe->values[iv]+=residual/connectivity;
         }
-
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
@@ -399 +397 @@
 
 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-
         /*Intermediaries*/
         bool     active_element;
@@ -417 +414 @@
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
-
         /*Intermediary*/
         int* doflist = NULL;
@@ -425 +421 @@
 
         /*Fetch dof list and allocate solution vector*/
-        IssmDouble* sedhead     = xNew<IssmDouble>(numnodes);
         IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
-
         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-        basalelement->GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
 
         Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
@@ -435 +428 @@
 
         /*Use the dof list to index into the solution vector: */
-        for(int i=0;i<numnodes;i++){
-                eplHeads[i]=solution[doflist[i]];
-                if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
-                if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
+        /*If the EPL is not active we revert to the initialisation vallue*/
+        /*       For now we keep OldValue but it is probably not the best*/
+        if(active_element){
+                for(int i=0;i<numnodes;i++){
+                        eplHeads[i]=solution[doflist[i]];
+                        if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
+                        if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
+                }
+        }
+        else{
+                basalelement->GetInputListOnVertices(&eplHeads[0],EplHeadOldEnum);
+                for(int i=0;i<numnodes;i++){
+                        if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
+                        if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
+                }
         }
         /*Add input to the element: */
         element->AddBasalInput(EplHeadEnum,eplHeads,P1Enum);
-
         /*Free ressources:*/
         xDelete<IssmDouble>(eplHeads);
-        xDelete<IssmDouble>(sedhead);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};