Changeset 21341 for issm/trunk/src/m/mech/thomasparams.py

Timestamp:

Nov 4, 2016, 1:48:43 PM (10 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 21337

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/mech/thomasparams.py (modified) (6 diffs)

issm/trunk/src/m/mech/thomasparams.py

@@ -1 @@
-import numpy as npy
+import numpy as  np
 from averaging import averaging
 
@@ -75 +75 @@
         
         # checks: any of e1 or e2 equal to zero?
-        pos=npy.nonzero(e1==0)
-        if npy.any(pos==1):
+        pos=np.nonzero(e1==0)
+        if np.any(pos==1):
                 print 'WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1'
                 e1[pos]=1.e-13
-        pos=npy.nonzero(e2==0)
-        if npy.any(pos==1):
+        pos=np.nonzero(e2==0)
+        if np.any(pos==1):
                 print 'WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1'
                 e2[pos]=1.e-13
@@ -89 +89 @@
         
         if coordsys=='principal':
-                b=npy.zeros((md.mesh.numberofvertices,))
+                b=np.zeros((md.mesh.numberofvertices,))
                 ex=e1
                 a=e2/e1
-                pos=npy.nonzero(npy.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
+                pos=np.nonzero(np.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
                 a[pos]=e1[pos]/e2[pos]
                 ex[pos]=e2[pos]
-                pos2=npy.nonzero(e1<0 & e2<0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal compression
+                pos2=np.nonzero(e1<0 & e2<0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal compression
                 a[pos2]=e1[pos2]/e2[pos2]
                 ex[pos2]=e2[pos2]
-                pos3=npy.nonzero(e1>0 & e2>0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal tension
+                pos3=np.nonzero(e1>0 & e2>0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal tension
                 a[pos3]=e1[pos3]/e2[pos3]
                 ex[pos3]=e2[pos3]
-                ind=npy.nonzero(e1<0 & e2<0)
+                ind=np.nonzero(e1<0 & e2<0)
                 a[ind]=-a[ind] # where both strain rates are compressive, enforce negative alpha
-                sigxx=(npy.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
+                sigxx=(np.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
         elif coordsys=='xy':
                 ex=exx
@@ -110 +110 @@
         elif coordsys=='longitudinal':
                 # using longitudinal strain rates defined by observed velocity vector
-                velangle=npy.arctan(md.initialization.vy/md.initialization.vx)
-                pos=npy.nonzero(md.initialization.vx==0)
-                velangle[pos]=npy.pi/2
-                ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
+                velangle=np.arctan(md.initialization.vy/md.initialization.vx)
+                pos=np.nonzero(md.initialization.vx==0)
+                velangle[pos]=np.pi/2
+                ex=0.5*(exx+eyy)+0.5*(exx-eyy)*np.cos(2.*velangle)+exy*np.sin(2.*velangle)
                 ey=exx+eyy-ex # trace of strain rate tensor is invariant
-                exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
+                exy=-0.5*(exx-eyy)*np.sin(2.*velangle)+exy*np.cos(2.*velangle)
                 a=ey/ex
                 b=exy/ex
@@ -124 +124 @@
         # a < -1 in areas of strong lateral compression or longitudinal compression and
         # theta flips sign at a = -2
-        pos=npy.nonzero(npy.abs((npy.abs(a)-2.))<1.e-3)
+        pos=np.nonzero(np.abs((np.abs(a)-2.))<1.e-3)
         if len(pos)>0:
                 print 'Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2'
@@ -131 +131 @@
         if eq=='Weertman1D':
                 theta=1./8
-                a=npy.zeros((md.mesh.numberofvertices,))
+                a=np.zeros((md.mesh.numberofvertices,))
         elif eq=='Weertman2D':
                 theta=1./9
-                a=npy.ones((md.mesh.numberofvertices,))
+                a=np.ones((md.mesh.numberofvertices,))
         elif eq=='Thomas':
-                theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(npy.abs(2.+a)**n)
+                theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(np.abs(2.+a)**n)
         else:
                 raise ValueError('argument passed to "eq" not valid')