Changeset 13064 for issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp

Timestamp:

08/16/12 10:27:48 (13 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on L1L2, almost there

File:

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

issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp

@@ -8182 +8182 @@
          * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
 
-        _error_("Not supported yet");
         int        i;
         IssmDouble a,c,d,z,s,viscosity;
@@ -8200 +8199 @@
         surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
         PentaRef::GetInputValue(&z,&z_list[0],gauss);
-        tau_perp = matpar->GetRhoIce() * matpar->GetG() * (s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
+        tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
 
         /* Get eps_b*/
@@ -8208 +8207 @@
         eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
         if(eps_b==0.){
-                *pviscosity=2.5*pow(10.,17.);
+                *pviscosity = 2.5e+17;
                 return;
         }
@@ -8216 +8215 @@
         A=matice->GetA();
 
-        /*Solve for tau_par (http://en.wikipedia.org/wiki/Cubic_function)*/
-        a=A;
-        c=A*tau_perp*tau_perp;
-        d=-eps_b;
-        tau_par = -1./(3.*a) * pow(0.5*(27.*a*a*d + sqrt(27.*a*a*d*27.*a*a*d -4.*pow(-3.*a*c,3.))),1./3.);
-
-        //printf("=======================================\n");
-        //printf("tau_par = %g (%g=0?)\n",tau_par,a*tau_par*tau_par*tau_par+c*tau_par+d);
-        //printf("a=%g c=%g d=%g\n",a,c,d);
-        //IssmDouble coeff[4];
-        //coeff[0]=d;
-        //coeff[1]=c;
-        //coeff[2]=0.;
-        //coeff[3]=a;
-        //int numroots;
-        //IssmDouble roots[3];
-        //cubic(coeff,roots,&numroots);
-        //tau_par=roots[0];
-        //for(i=0;i<numroots;i++) printf(" %g ",roots[i]);
+        /*Solve for tau_perp*/
+        int        numroots;
+        IssmDouble roots[3];
+        a = A;
+        c = A *tau_perp*tau_perp;
+        d = - eps_b;
+        cubic(a,0.,c,d,roots,&numroots);
+        tau_par=roots[0];
         //printf(" (%g =0?)\n",a*roots[0]*roots[0]*roots[0]+c*roots[0]+d);
-        
 
         /*Viscosity*/
-        viscosity = 1./(2.*A)*pow(tau_par*tau_par + tau_perp*tau_perp ,-2.);
+        viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
+        //printf("denom = %g (%g,%g)\n",tau_par*tau_par + tau_perp*tau_perp,tau_par,tau_perp);
         _assert_(!isnan(viscosity));