/[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 576 by gross, Fri Mar 3 08:28:42 2006 UTC revision 580 by gross, Wed Mar 8 05:45:51 2006 UTC
# Line 15  Line 15 
15    
16  #if !defined escript_LocalOps_H  #if !defined escript_LocalOps_H
17  #define escript_LocalOps_H  #define escript_LocalOps_H
18    #include <math.h>
19  namespace escript {  namespace escript {
20    
21    
# Line 29  namespace escript { Line 29  namespace escript {
29  inline  inline
30  void eigenvalues1(const double A00,double* ev0) {  void eigenvalues1(const double A00,double* ev0) {
31    
32     *ev0=1.;     *ev0=A00;
33    
34  }  }
35  /**  /**
# Line 46  inline Line 46  inline
46  void eigenvalues2(const double A00,const double A01  void eigenvalues2(const double A00,const double A01
47                                   ,const double A11,                                   ,const double A11,
48                   double* ev0, double* ev1) {                   double* ev0, double* ev1) {
49          const register double trA=(A00+A11)/2.;
50     *ev0=1.;        const register double A_00=A00-trA;
51     *ev1=2.;        const register double A_11=A11-trA;
52          const register double s=sqrt(A01*A01-A_00*A_11);
53          *ev0=trA-s;
54          *ev1=trA+s;
55  }  }
56  /**  /**
57     \brief     \brief
# Line 70  void eigenvalues3(const double A00, cons Line 73  void eigenvalues3(const double A00, cons
73                                                       const double A22,                                                       const double A22,
74                   double* ev0, double* ev1,double* ev2) {                   double* ev0, double* ev1,double* ev2) {
75    
76     *ev0=1.;        const register double trA=(A00+A11+A22)/3.;
77     *ev1=2.;        const register double A_00=A00-trA;
78     *ev2=3.;        const register double A_11=A11-trA;
79          const register double A_22=A22-trA;
80          const register double A01_2=A01*A01;
81          const register double A02_2=A02*A02;
82          const register double A12_2=A12*A12;
83          const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
84          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);
85          const register double sq_p=sqrt(p/3.);
86          register double z=-q/(2*pow(sq_p,3));
87          if (z<-1.) {
88             z=-1.;
89          } else if (z>1.) {
90             z=1.;
91          }
92          const register double alpha_3=acos(z)/3.;
93          *ev2=trA+2.*sq_p*cos(alpha_3);
94          *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
95          *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
96  }  }
   
97  } // end of namespace  } // end of namespace
98  #endif  #endif
Legend:
Removed from v.576  
changed lines
  Added in v.580
  ViewVC Help
Powered by ViewVC 1.1.26