Changeset 20650

Timestamp:

05/25/16 20:40:09 (9 years ago)

Author:

Mathieu Morlighem

Message:

BUG: missing powers and wrong sign in eps_prime

File:

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

issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp

@@ -331 +331 @@
         }
         else{
+                dvx[2] = 0.;
                 vy = 0.;
                 dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
         }
         vz = 0.;
-        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
+        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
 
         /*Get enhancement factors and ko*/
@@ -371 +372 @@
         }
         else{
+                dvx[1] = 0.;
+                dvx[2] = 0.;
                 vy = 0.;
                 dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
         }
         vz = 0.;
-        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
+        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
 
         /*Get enhancement factors and ko*/
@@ -382 +385 @@
         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
         ec_input->GetInputValue(&Ec,gauss);
-        es_input->GetInputValue(&Es,gauss);
-        ko_input->GetInputValue(&ko,gauss);
-
-        /*Compute viscosity*/
-        *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}/*}}}*/
-void  Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
-void  Matbafl::ResetHooks(){/*{{{*/
-
-        this->element=NULL;
-
-        /*Get Element type*/
-        this->helement->reset();
-
-}
-/*}}}*/
-IssmDouble Matbafl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
-
-        /*Intermediaries*/
-        IssmDouble viscosity;
-        IssmDouble vorticity[3],vorticity_norm;
-        IssmDouble nrsp[3],nrsp_norm;
-        IssmDouble eps[3][3],epso;
-        IssmDouble epsprime[3],epsprime_norm;
-        IssmDouble E,lambdas;
-
-        /*Create vorticity vector*/
-        _assert_(dvx && dvy && dvz);
-        vorticity[0] =  dvz[1] - dvy[2];
-        vorticity[1] =  dvx[2] - dvz[0];
-        vorticity[2] =  dvy[0] - dvx[1];
-
-        /*Normalize*/
-        vorticity_norm = sqrt(vorticity[0]*vorticity[0] + vorticity[1]*vorticity[1] + vorticity[2]*vorticity[2]);
-        if(vorticity_norm==0){
-                vorticity[0] = 0.;
-                vorticity[1] = 0.;
-                vorticity[2] = 1.;
-        }
         else{
                 vorticity[0] =vorticity[0]/vorticity_norm;
@@ -472 +434 @@
                 for(int j=0;j<3;j++){
                         for(int k=0;k<3;k++){
-                                epsprime[j] += nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
+                                epsprime[j] += -nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
                         }
                 }
@@ -480 +442 @@
                 for(int j=0;j<3;j++){
                         for(int k=0;k<3;k++){
-                                epsprime[j] += nrsp[i]*eps[i][k]*vorticity[k]*vorticity[j];
+                                epsprime[j] += -nrsp[i]*eps[i][k]*vorticity[k]*vorticity[j];
                         }
                 }
@@ -490 +452 @@
         /*Compute lambda_s*/
         _assert_(epso>0.);
-        lambdas = sqrt(2*epsprime_norm/(3*epso));
+        lambdas = sqrt(2*epsprime_norm*epsprime_norm/(3*epso*epso));
 
         /*Get total enhancement factor E(lambdas)*/