Changeset 17093

Timestamp:

01/10/14 13:28:59 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: simplifying strain rate computation

File:

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

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

@@ -413 +413 @@
          */
 
-        IssmDouble epsilonvx[6];
-        IssmDouble epsilonvy[6];
-        IssmDouble epsilonvz[6];
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+        IssmDouble dvy[3];
+        IssmDouble dvz[3];
 
         /*Check that both inputs have been found*/
@@ -423 +424 @@
 
         /*Get strain rate assuming that epsilon has been allocated*/
-        vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
-        vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
-        vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
-
-        /*Sum all contributions*/
-        for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
+        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = dvy[1];
+        epsilon[2] = dvz[2];
+        epsilon[3] = 0.5*(dvx[1] + dvy[0]);
+        epsilon[4] = 0.5*(dvx[2] + dvz[0]);
+        epsilon[5] = 0.5*(dvy[2] + dvz[1]);
+
 }/*}}}*/
 void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
@@ -441 +446 @@
          */
 
-        int i;
-        IssmDouble epsilonvx[5];
-        IssmDouble epsilonvy[5];
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+        IssmDouble dvy[3];
 
         /*Check that both inputs have been found*/
@@ -451 +456 @@
 
         /*Get strain rate assuming that epsilon has been allocated*/
-        vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
-        vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
-
-        /*Sum all contributions*/
-        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = dvy[1];
+        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
+        epsilon[3] = 0.5*dvx[2];
+        epsilon[4] = 0.5*dvy[2];
+
 }/*}}}*/
 void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
-        int i;
-        IssmDouble epsilonvx[3];
-        IssmDouble epsilonvy[3];
+        /*Intermediaries*/
+        IssmDouble dvx[3];
+        IssmDouble dvy[3];
 
         /*Check that both inputs have been found*/
@@ -469 +477 @@
 
         /*Get strain rate assuming that epsilon has been allocated*/
-        vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
-        vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
-
-        /*Sum all contributions*/
-        for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+        epsilon[0] = dvx[0];
+        epsilon[1] = dvy[1];
+        epsilon[2] = 0.5*(dvx[1] + dvy[0]);
+
 }/*}}}*/
 void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/