[cig-commits] r19956 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16

[brad at geodynamics.org]

Author: brad
Date: 2012-04-20 10:41:37 -0700 (Fri, 20 Apr 2012)
New Revision: 19956

Modified:
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py
Log:
Finished script to plot initial stress and final slip.

Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py	2012-04-20 17:32:10 UTC (rev 19955)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py	2012-04-20 17:41:37 UTC (rev 19956)
@@ -12,12 +12,13 @@
 sim = "tpv16"
 cell = "tet4"
 dx = 75
-griddata = True
+griddata = False
 
 # ----------------------------------------------------------------------
 import tables
 import numpy
 import matplotlib.pyplot as pyplot
+import matplotlib.cm as cm
 import sys
 import subprocess
 
@@ -28,7 +29,7 @@
 
 # ----------------------------------------------------------------------
 figure = matplotlibext.Figure()
-figure.open(3.0, 2.5, margins=[[0.4, 0.3, 0.1], [0.35, 0.4, 0.1]], dpi=150)
+figure.open(3.0, 3.0, margins=[[0.4, 0.3, 0.16], [0.35, 0.3, 0.2]], dpi=150)
 
 nrows = 2
 ncols = 1
@@ -43,36 +44,48 @@
         distStrike = vertices[:,1]
         if irow == 0:
             tractionShear = h5.root.vertex_fields.traction[1,:,0].squeeze()
-            data = numpy.transpose((distStrike, dispDip, tractionShear))
+            data = numpy.transpose((distStrike, distDip, -tractionShear/1.0e+6))
         else:
             slip = h5.root.vertex_fields.slip[-1,:,:].squeeze()
             slipMag = (slip[:,0]**2 + slip[:,1]**2)**0.5
             data = numpy.transpose((distStrike, distDip, slipMag))
         h5.close()
-        numpy.savetxt("tmp.dat", data, fmt='%14.63')
-        cmd = "awk '/^j/ {next}; { print $0 }' tmp.dat | triangulate -Gtmp.grd -I75.0/75.0 -R-24.0e+3/+24.0e+3/0.0/19.5e+3 -Jx1.0 > /dev/null "
-                subprocess.call(cmd, shell=True)
-                cmd = "grd2xyz tmp.grd > tmp.txtgrd"
-                subprocess.call(cmd, shell=True)
+        numpy.savetxt("tmp.dat", data, fmt='%14.6f')
+        cmd = "triangulate -Gtmp.grd -I75.0/75.0 -R-24.0e+3/+24.0e+3/0.0/19.5e+3 -Jx1.0 < tmp.dat > /dev/null "
+        subprocess.call(cmd, shell=True)
+        cmd = "grd2xyz tmp.grd > tmp%d.txtgrd" % irow
+        subprocess.call(cmd, shell=True)
 
 
     ax = figure.axes(nrows, ncols, irow+1, icol+1)
-    filename = "tmp.txtgrd"
+    filename = "tmp%d.txtgrd" % irow
     data = numpy.loadtxt(filename)
     data = numpy.reshape(data, (261, 641, 3))
-    contours = ax.contour(data[:,:,0]/1.0e+3, data[:,:,1]/1e+3, data[:,:,2],
-                                  levels=numpy.arange(0.0, 10.01, 0.5),
-                                  colors=lineStyle[0],
-                                  linewidths=1.0,
-                                  linestyles=lineStyle[1])
+    if irow == 0:
+        label = "Initial Shear Traction (MPa)"
+        contours = numpy.arange(0.0, 50.01, 5.0)
+    else:
+        label = "Final Slip (m)"
+        contours = numpy.arange(0.0, 5.01, 0.5)
+    pyplot.set_cmap(cm.hot_r)
+    CS = ax.contourf(data[:,:,0]/1.0e+3, data[:,:,1]/1e+3, data[:,:,2],
+                     levels=contours)
+    pyplot.hold(True)
+    ax.contour(data[:,:,0]/1.0e+3, data[:,:,1]/1e+3, data[:,:,2],
+                levels=contours,
+               colors=lineStyle[0],
+               linewidths=0.5)
 
     ax.set_xlim((-24.0, +24.0))
     ax.set_xlabel("Dist. Along Strike (km)")
     ax.set_ylim((19.5, 0.0))
     ax.set_yticks(numpy.arange(0, 20.01, 5.0))
     ax.set_ylabel("Dist. Down Dip (km)")
+    ax.set_title(label)
     ax.set_aspect('equal')
     
+    pyplot.colorbar(CS, orientation='vertical', fraction=0.05, pad=0.02)
+
     if irow+1 < nrows:
         ax.set_xticklabels([])
         ax.set_xlabel("")