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

Timestamp:

02/10/14 14:38:22 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: simplifying getviscosity calls, now only provide effective strain rate

File:

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

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

@@ -283 +283 @@
 
                 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-                this->material->GetViscosity3dFS(&viscosity,&epsilon[0]);
+                this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
                 this->GetPhi(&phi,&epsilon[0],viscosity);
 
@@ -303 +303 @@
         IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
+        IssmDouble eps_eff;
+        IssmDouble eps0=1.e-27;
 
         if(dim==3){
-                /*3D*/
+                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
                 this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-                material->GetViscosity3dFS(&viscosity, &epsilon3d[0]);
+                eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
         }
         else{
-                /*2D*/
+                /* eps_eff^2 = exx^2 + eyy^2 + 2*exy^2 */
                 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-                material->GetViscosity2dvertical(&viscosity,&epsilon2d[0]);
-        }
+                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
+        }
+
+        /*Get viscosity*/
+        material->GetViscosity(&viscosity,eps_eff);
 
         /*Assign output pointer*/
@@ -384 +389 @@
         IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
         IssmDouble epsilon2d[2]; /* epsilon=[exx,exy];           */
+        IssmDouble eps_eff;
 
         if(dim==3){
+                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
                 this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
-                material->GetViscosity3d(&viscosity, &epsilon3d[0]);
+                eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
         }
         else{
+                /* eps_eff^2 = exx^2 + exy^2 */
                 this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-                material->GetViscosity2dverticalHO(&viscosity, &epsilon2d[0]);
-        }
+                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1]);
+        }
+
+        /*Get viscosity*/
+        material->GetViscosity(&viscosity,eps_eff);
 
         /*Assign output pointer*/
@@ -402 +413 @@
         IssmDouble viscosity;
         IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];    */
-
+        IssmDouble eps_eff;
+
+        /*Get effective strain rate
+         * eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/
         this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-        material->GetViscosity2d(&viscosity, &epsilon[0]);
+        eps_eff = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + epsilon[0]*epsilon[1]);
+
+        /*Get viscosity*/
+        material->GetViscosityBar(&viscosity,eps_eff);
 
         /*Assign output pointer*/