Changeset 13065

Timestamp:

08/16/12 11:09:34 (13 years ago)

Author:

Mathieu Morlighem

Message:

CHG: simpler way of computing root of cubic equation

File:

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

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

@@ -8183 +8183 @@
 
         int        i;
-        IssmDouble a,c,d,z,s,viscosity;
+        IssmDouble z,s,viscosity,p,q,delta;
         IssmDouble tau_perp,tau_par,eps_b,A;
         IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/
@@ -8215 +8215 @@
         A=matice->GetA();
 
-        /*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);
+        /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
+        p     = tau_perp *tau_perp;
+        q     = - eps_b/A;
+        delta = q *q + p*p*p*4./27.;
+        _assert_(delta>0);
+        tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
 
         /*Viscosity*/
         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));
+        _assert_(viscosity > 0.);
 
         /*Assign output pointer*/