/[escript]/trunk/escript/src/LocalOps.h
ViewVC logotype

Diff of /trunk/escript/src/LocalOps.h

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 583 by gross, Wed Mar 8 08:15:34 2006 UTC revision 585 by gross, Thu Mar 9 23:47:42 2006 UTC
# Line 80  void eigenvalues3(const double A00, cons Line 80  void eigenvalues3(const double A00, cons
80        const register double A02_2=A02*A02;        const register double A02_2=A02*A02;
81        const register double A12_2=A12*A12;        const register double A12_2=A12*A12;
82        const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;        const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
83        const register double q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);        if (p<=0.) {
84        const register double sq_p=sqrt(p/3.);           *ev2=trA;
85        register double z=-q/(2*pow(sq_p,3));           *ev1=trA;
86        if (z<-1.) {           *ev0=trA;
87           z=-1.;  
88        } else if (z>1.) {        } else {
89           z=1.;           const register double q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);
90             const register double sq_p=sqrt(p/3.);
91             register double z=-q/(2*pow(sq_p,3));
92             if (z<-1.) {
93                z=-1.;
94             } else if (z>1.) {
95                z=1.;
96             }
97             const register double alpha_3=acos(z)/3.;
98             *ev2=trA+2.*sq_p*cos(alpha_3);
99             *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
100             *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
101        }        }
       const register double alpha_3=acos(z)/3.;  
       *ev2=trA+2.*sq_p*cos(alpha_3);  
       *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);  
       *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);  
102  }  }
103  /**  /**
104     \brief     \brief
Legend:
Removed from v.583  
changed lines
  Added in v.585
  ViewVC Help
Powered by ViewVC 1.1.26