Changeset 15182

Timestamp:

05/30/13 18:08:42 (12 years ago)

Author:

bdef

Message:

CHG: treatment of the residual input is changed to avoid division multiplication by the basis function

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (6 diffs)

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

@@ -6200 +6200 @@
 
         /*Intermediaries */
+        IssmDouble connectivity;
         IssmDouble Jdet;
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble dt,scalar,water_head;
-        IssmDouble residual_load,transfer;
+        IssmDouble transfer,residual;
         IssmDouble epl_storing;
         IssmDouble basis[numdof];
@@ -6233 +6234 @@
 
                 /*Loading term*/
-                residual_input->GetInputValue(&residual_load,gauss);
                 transfer_input->GetInputValue(&transfer,gauss);
-                scalar = Jdet*gauss->weight*(residual_load+transfer);
+                scalar = Jdet*gauss->weight*transfer;
                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
                 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6247 +6247 @@
         }
 
+        /*      Add residual if necessary*/
+        delete gauss;
+        gauss=new GaussTria();
+        for(int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                connectivity = IssmDouble(nodes[iv]->GetConnectivity());
+                residual_input->GetInputValue(&residual,gauss);
+                pe->values[iv]+=residual/connectivity;
+        }
         /*Clean up and return*/
         delete gauss;
@@ -6387 +6397 @@
         IssmDouble  values[numdof];
         IssmDouble  residual[numdof];
-        IssmDouble  intbasis[numdof];   
         IssmDouble  penalty_factor, dt;
         IssmDouble  kmax, kappa, h_max;
@@ -6402 +6411 @@
         /*If converged keep the residual in mind*/
         this->inputs->GetInputValue(&converged,ConvergedEnum);
-        GetInputListOnVertices(&intbasis[0],BasisIntegralEnum);
 
         if(converged){
@@ -6413 +6421 @@
                         this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
                         if(values[i]>h_max){
-                                residual[i]=kappa*(values[i]-h_max)/intbasis[i];
-                                if(reCast<bool,IssmDouble>(dt))residual[i]=residual[i]/dt;
+                                residual[i]=kappa*(values[i]-h_max);
+                                if(reCast<bool,IssmDouble>(dt))residual[i]=residual[i];
                         }
                         else
                                 residual[i]=0.0;
                 }
-                if(this->id==1)printf("residual, %g \n", residual[1]);
-                if(this->id==1)printf("hSed, %g \n", values[1]);
         }
         /*Add input to the element: */