/[escript]/trunk/escript/py_src/util.py
ViewVC logotype

Diff of /trunk/escript/py_src/util.py

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

revision 584 by gross, Wed Mar 8 05:45:51 2006 UTC revision 585 by gross, Thu Mar 9 23:47:42 2006 UTC
# Line 3566  def eigenvalues(arg): Line 3566  def eigenvalues(arg):
3566        if s[0]==1:        if s[0]==1:
3567            return arg[0]            return arg[0]
3568        elif s[0]==2:        elif s[0]==2:
3569            A11=arg[0,0]            arg1=symmetric(arg)
3570            A12=arg[0,1]            A11=arg1[0,0]
3571            A22=arg[1,1]            A12=arg1[0,1]
3572              A22=arg1[1,1]
3573            trA=(A11+A22)/2.            trA=(A11+A22)/2.
3574            A11-=trA            A11-=trA
3575            A22-=trA            A22-=trA
3576            s=sqrt(A12**2-A11*A22)            s=sqrt(A12**2-A11*A22)
3577            return trA+s*numarray.array([-1.,1.],type=numarray.Float64)            return trA+s*numarray.array([-1.,1.],type=numarray.Float64)
3578        elif s[0]==3:        elif s[0]==3:
3579            A11=arg[0,0]            arg1=symmetric(arg)
3580            A12=arg[0,1]            A11=arg1[0,0]
3581            A22=arg[1,1]            A12=arg1[0,1]
3582            A13=arg[0,2]            A22=arg1[1,1]
3583            A23=arg[1,2]            A13=arg1[0,2]
3584            A33=arg[2,2]            A23=arg1[1,2]
3585              A33=arg1[2,2]
3586            trA=(A11+A22+A33)/3.            trA=(A11+A22+A33)/3.
3587            A11-=trA            A11-=trA
3588            A22-=trA            A22-=trA
# Line 3591  def eigenvalues(arg): Line 3593  def eigenvalues(arg):
3593            p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.            p=A13_2+A23_2+A12_2+(A11**2+A22**2+A33**2)/2.
3594            q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13            q=A13_2*A22+A23_2*A11+A12_2*A33-A11*A22*A33-2*A12*A23*A13
3595            sq_p=sqrt(p/3.)            sq_p=sqrt(p/3.)
3596            alpha_3=acos(clip(-q*sq_p**(-3.)/2.,-1.,1.))/3.            alpha_3=acos(clip(-q*(sq_p+whereZero(p,0.)*1.e-15)**(-3.)/2.,-1.,1.))/3.  # whereZero is protection against divison by zero
3597            sq_p*=2.            sq_p*=2.
3598            f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \            f=cos(alpha_3)               *numarray.array([0.,0.,1.],type=numarray.Float64) \
3599             -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \             -cos(alpha_3+numarray.pi/3.)*numarray.array([0.,1.,0.],type=numarray.Float64) \
# Line 3929  def power(arg0,arg1): Line 3931  def power(arg0,arg1):
3931         """         """
3932         args=matchShape(arg0,arg1)         args=matchShape(arg0,arg1)
3933         if testForZero(args[0]):         if testForZero(args[0]):
3934            return numarray.zeros(args[0],numarray.Float64)            return numarray.zeros(pokeShape(args[0]),numarray.Float64)
3935         elif testForZero(args[1]):         elif testForZero(args[1]):
3936            return numarray.ones(args[0],numarray.Float64)            return numarray.ones(pokeShape(args[1]),numarray.Float64)
3937         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):         elif isinstance(args[0],Symbol) or isinstance(args[1],Symbol):
3938            return Power_Symbol(args[0],args[1])            return Power_Symbol(args[0],args[1])
3939         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):         elif isinstance(args[0],numarray.NumArray) and not isinstance(args[1],numarray.NumArray):
Legend:
Removed from v.584  
changed lines
  Added in v.585
  ViewVC Help
Powered by ViewVC 1.1.26