[cig-commits] r14939 - short/2.5D/benchmarks/savageprescott/figures
brad at geodynamics.org
brad at geodynamics.org
Fri May 8 19:31:34 PDT 2009
Author: brad
Date: 2009-05-08 19:31:34 -0700 (Fri, 08 May 2009)
New Revision: 14939
Added:
short/2.5D/benchmarks/savageprescott/figures/plot_profiles.py
short/2.5D/benchmarks/savageprescott/figures/plot_soln.py
short/2.5D/benchmarks/savageprescott/figures/vtk_geometry.py
Log:
Added some plotting files.
Added: short/2.5D/benchmarks/savageprescott/figures/plot_profiles.py
===================================================================
--- short/2.5D/benchmarks/savageprescott/figures/plot_profiles.py (rev 0)
+++ short/2.5D/benchmarks/savageprescott/figures/plot_profiles.py 2009-05-09 02:31:34 UTC (rev 14939)
@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+cycles = [2,9]
+style = "lightbg"
+fileSuffix = "eps"
+plotSize = "poster"
+
+# ======================================================================
+import numpy
+import pylab
+from mypylab.Figure import Figure
+
+# ----------------------------------------------------------------------
+class ProfileSet(object):
+ """
+ Set of profiles.
+ """
+
+ def __init__(self, filename):
+ self.data = numpy.loadtxt(filename, comments="#")
+ ncols = self.data.shape[1]
+ self.data[:,1:ncols+1] /= 400.0 # Coseismic displacement
+ self.data[:,0] /= 40.0 # Thickness of elastic layer
+ return
+
+
+ def plot(self, linestyle):
+ data = self.data
+ ncols = data.shape[1]
+ if ncols == 42:
+ indicesTime = [2, 10, 20, 30, 38]
+ elif ncols == 21:
+ indicesTime = [1, 5, 10, 15, 19]
+
+ colorOrder = ('orange', 'blue', 'red', 'green', 'purple')
+ i = 0
+ for index in indicesTime:
+ h = pylab.plot(data[:,0], data[:,1+index], linestyle, color=colorOrder[i])
+ i += 1
+ return h
+
+
+# ----------------------------------------------------------------------
+class Profiles(Figure):
+ """
+ Figure with velocity time histories for a site.
+ """
+
+ def __init__(self, cycle):
+ if plotSize == "poster":
+ fontsize = 14
+ elif plotSize == "presentation":
+ fontsize = 14
+ elif plotSize == "manual":
+ fontsize = 10
+ else:
+ raise ValueError("Unknown plotSize '%s'." % plotSize)
+ Figure.__init__(self, color=style, fontsize=fontsize)
+ self.cycle = cycle
+ return
+
+
+ def plot(self):
+
+ self.width = 7.25
+ self.height = 5.0
+
+ # Create figure
+ self.open(self.width, self.height, margins=[[0.70, 0, 0.20],
+ [0.55, 0, 0.12]])
+
+ # Plot profiles
+ self.axes(1, 1, 1, 1)
+ analytic = ProfileSet("savpres_displ_c%dref.txt" % self.cycle)
+ ah = analytic.plot(linestyle="-")
+ pylab.hold(True)
+ simulation = ProfileSet("spbm_hex8_graded3_c%dref.txt" % self.cycle)
+ sh = simulation.plot(linestyle="--")
+ pylab.hold(False)
+
+ pylab.xlim(0, 10)
+ pylab.ylim(0, 0.5)
+
+ leg = pylab.legend((ah, sh),
+ ["analytic", "simulation"],
+ loc="upper left")
+ leg.legendPatch.set_facecolor('bg')
+ leg.legendPatch.set_edgecolor('fg')
+ self._annotate()
+ return
+
+
+ def save(self):
+ pylab.figure(self.handle.number)
+ pylab.savefig("savageprescott_profile%d.%s" % (self.cycle, fileSuffix))
+ return
+
+
+ def _annotate(self):
+ """
+ Add title and labels for axes.
+ """
+ pylab.xlabel("Dist. from fault / Elastic thickness")
+ pylab.ylabel("Disp. / Coseismic Disp. ")
+
+ pylab.text(9.9, 0.03, "t=0.05", horizontalalignment="right")
+ pylab.text(9.9, 0.13, "t=0.25", horizontalalignment="right")
+ pylab.text(9.9, 0.255, "t=0.50", horizontalalignment="right")
+ pylab.text(9.9, 0.37, "t=0.75", horizontalalignment="right")
+ pylab.text(9.9, 0.46, "t=0.95", horizontalalignment="right")
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+
+ figures = []
+ for cycle in cycles:
+ figure = Profiles(cycle)
+ figure.plot()
+ figures.append(figure)
+
+ pylab.show()
+ for figure in figures:
+ figure.save()
+
+# End of file
Property changes on: short/2.5D/benchmarks/savageprescott/figures/plot_profiles.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/2.5D/benchmarks/savageprescott/figures/plot_soln.py
===================================================================
--- short/2.5D/benchmarks/savageprescott/figures/plot_soln.py (rev 0)
+++ short/2.5D/benchmarks/savageprescott/figures/plot_soln.py 2009-05-09 02:31:34 UTC (rev 14939)
@@ -0,0 +1,200 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+shape = "hex8"
+scaleFactor = 5.0/40.0
+t = 1900
+pngfile = "savageprescott_soln.png"
+
+style = {'colors': "lightbg",
+ }
+
+from enthought.mayavi.app import Mayavi
+import vtk_geometry
+from enthought.mayavi.sources.vtk_data_source import VTKDataSource
+from enthought.mayavi.filters.warp_vector import WarpVector
+from enthought.mayavi.filters.extract_vector_norm import ExtractVectorNorm
+from enthought.mayavi.modules.surface import Surface
+from enthought.mayavi.modules.outline import Outline
+from enthought.mayavi.modules.axes import Axes
+from enthought.mayavi.modules.surface import Surface
+
+class PlotSoln(Mayavi):
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def run(self):
+
+ self._setup()
+
+ self._renderScene()
+ self._annotateScene()
+ self._setCamera()
+ scene = self.script.engine.current_scene.scene
+ scene.disable_render = False
+ scene.render()
+ scene.save_png(pngfile)
+
+ return
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _setup(self):
+ """
+ Plot axes, fault surface, and materials.
+ """
+
+ dkslate = (47/255.0, 53/255.0, 72/255.0)
+ ltblue = (51/255.0, 187/255.0, 255/255.0)
+ ltred = (1.0, 0.25, 0.25)
+ red = (1.0, 0.40, 0.40)
+ ltyellow = (1.0, 1.0, 0.45)
+
+ if style['colors'] == "darkbg":
+ dkslate = (47/255.0, 53/255.0, 72/255.0)
+ ltblue = (51/255.0, 187/255.0, 255/255.0)
+ ltred = (1.0, 0.25, 0.25)
+ self.colorFg = (1,1,1)
+ self.colorBg = dkslate
+ self.colorSurfTrace = red
+ self.colorCities = ltyellow
+ else:
+ self.colorFg = (0,0,0)
+ self.colorBg = (1,1,1)
+ self.colorSurfTrace = (0,0,1)
+ self.lut = "hot"
+ self.lutReverse = True
+
+ self.windowSize = (1600+16, 675+120)
+ self.aaframes = 4
+
+ return
+
+
+ def _renderScene(self):
+
+ script = self.script
+ w = script.get_active_window()
+ w.size = self.windowSize
+ script.new_scene()
+ scene = script.engine.current_scene.scene
+ scene.disable_render = True
+ scene.anti_aliasing_frames = self.aaframes
+ scene.background = self.colorBg
+ scene.foreground = self.colorFg
+
+
+ script = self.script
+ script.add_source(VTKDataSource(data=self._readData()))
+ script.engine.current_object.name = "Solution"
+
+ warp = WarpVector()
+ warp.filter.scale_factor = scaleFactor
+ script.add_filter(warp)
+
+ norm = ExtractVectorNorm()
+ script.add_filter(norm)
+
+ surf = Surface()
+ script.add_module(surf)
+
+ wire = Surface()
+ script.add_module(wire)
+ wire.actor.actor.property.representation = "wireframe"
+ wire.actor.actor.property.color = (0.0, 0.0, 0.0)
+ wire.actor.mapper.scalar_visibility = False
+
+ colorbar = script.engine.current_object.module_manager.scalar_lut_manager
+ colorbar.scalar_bar.label_format = "%3.1f"
+ colorbar.scalar_bar.label_text_property.font_size = 18
+ colorbar.scalar_bar.label_text_property.shadow = True
+ colorbar.scalar_bar.label_text_property.italic = False
+ colorbar.scalar_bar.title_text_property.italic = False
+ colorbar.scalar_bar.title_text_property.shadow = True
+ colorbar.show_scalar_bar = True
+ colorbar.data_range = (0.0, 18.0)
+ colorbar.number_of_labels = 7
+ w,h = colorbar.scalar_bar.position2
+ colorbar.scalar_bar.position2 = (0.4, 0.11)
+ colorbar.scalar_bar.position = (0.25, 0.02)
+ colorbar.data_name = "Displacement [m]"
+
+ return
+
+ def _annotateScene(self):
+
+ script = self.script
+
+ # Domain (axes and outline)
+ script.add_source(VTKDataSource(data=vtk_geometry.domain()))
+ script.engine.current_object.name = "Domain"
+ outline = Outline()
+ script.add_module(outline)
+ outline.actor.property.opacity = 0.2
+ axes = Axes()
+ axes.axes.x_label = "X"
+ axes.axes.y_label = "Y"
+ axes.axes.z_label = "Z"
+ axes.axes.label_format = "%-0.1f"
+ script.add_module(axes)
+
+ return
+
+
+ def _setCamera(self):
+ script = self.script
+ vtk_geometry.setCamera(script.engine.current_scene.scene.camera)
+ return
+
+
+ def _readData(self):
+ from enthought.tvtk.api import tvtk
+ import tables
+ import numpy
+
+ filename = "../results/spbm_%s_graded3_20km/spbm_%s_graded3_20km_t%04d.h5" % (shape, shape, t)
+ h5 = tables.openFile(filename, 'r')
+
+ cells = h5.root.topology.cells[:]
+ (ncells, ncorners) = cells.shape
+ vertices = h5.root.geometry.vertices[:] / (1e+3 * 40.0)
+ (nvertices, spaceDim) = vertices.shape
+ disp = h5.root.solution.snapshot0.displacements[:]
+ h5.close()
+
+ if shape == "tet4":
+ assert(spaceDim == 3)
+ assert(ncorners == 4)
+ cellType = tvtk.Tetra().cell_type
+ elif shape == "hex8":
+ assert(spaceDim == 3)
+ assert(ncorners == 8)
+ cellType = tvtk.Hexahedron().cell_type
+ else:
+ raise ValueError("Unknown shape '%s'." % shape)
+
+ data = tvtk.UnstructuredGrid()
+ data.points = vertices
+ data.set_cells(cellType, cells)
+ data.point_data.vectors = disp
+ data.point_data.vectors.name = "Displacement [m]"
+ return data
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ app = PlotSoln()
+ app.main()
+
+
+# End of file
+
Property changes on: short/2.5D/benchmarks/savageprescott/figures/plot_soln.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/2.5D/benchmarks/savageprescott/figures/vtk_geometry.py
===================================================================
--- short/2.5D/benchmarks/savageprescott/figures/vtk_geometry.py (rev 0)
+++ short/2.5D/benchmarks/savageprescott/figures/vtk_geometry.py 2009-05-09 02:31:34 UTC (rev 14939)
@@ -0,0 +1,102 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+from enthought.tvtk.api import tvtk
+import numpy
+
+# ----------------------------------------------------------------------
+def domain():
+ domain = tvtk.CubeSource()
+ domain.set_bounds(-1000/40.0, 1000/40.0, -500/40.0, 500/40.0, -400/40.0, 0.0)
+ return domain.get_output()
+
+# ----------------------------------------------------------------------
+def fault(showTaper=False):
+ vertices = numpy.array(
+ [ [12.0, 0.0, 0.0], # 0
+ [12.0, 12.0, 0.0], # 1
+ [12.0, 16.0, 0.0], # 2
+ [12.0, 0.0, -12.0], # 3
+ [12.0, 12.0, -12.0], # 4
+ [12.0, 0.0, -16.0], # 5
+ [12.0, 12.0, -16.0], # 6
+ [12.0, 16.0, -16.0] ], # 7
+ dtype=numpy.float32)
+
+ # Setup VTK vertices
+ if not showTaper:
+ polys = numpy.array(
+ [ [0, 3, 1],
+ [1, 3, 4],
+ [3, 5, 4],
+ [4, 5, 6],
+ [1, 6, 2],
+ [2, 6, 7] ],
+ dtype=numpy.int32)
+
+ # Setup VTK simplices
+ data = [{'name': "fault",
+ 'object': tvtk.PolyData(points=vertices, polys=polys)}]
+ else:
+ polys = numpy.array(
+ [ [0, 3, 1],
+ [1, 3, 4] ],
+ dtype=numpy.int32)
+ polysTaper = numpy.array(
+ [ [3, 5, 4],
+ [4, 5, 6],
+ [1, 6, 2],
+ [2, 6, 7] ],
+ dtype=numpy.int32)
+
+ data= [{'name': "fault",
+ 'object': tvtk.PolyData(points=vertices, polys=polys)},
+ {'name': "taper",
+ 'object': tvtk.PolyData(points=vertices, polys=polysTaper)}]
+ return data
+
+
+# ----------------------------------------------------------------------
+def materials():
+ l = 24.0
+ elastic = tvtk.CubeSource()
+ elastic.set_bounds(0, l, 0, l, -l/2.0, 0.0)
+ viscoelastic = tvtk.CubeSource()
+ viscoelastic.set_bounds(0, l, 0, l, -l, -l/2.0)
+ return [{'name': "elastic",
+ 'object': elastic.get_output()},
+ {'name': "viscoelastic",
+ 'object': viscoelastic.get_output()}]
+
+# ----------------------------------------------------------------------
+def setCamera(camera):
+ dist = 2000.0/40.0
+ angle = 40.0
+ elev = 30.0
+ azimuth = 25.0
+ clipRange = numpy.array( [200/40.0, 5000/40.0] )
+ ptTo = numpy.array( [0.0, 0.0, -500/40.0] )
+ ptFrom = ptTo + numpy.array( [0.0, -dist, 0.0])
+
+ camera.view_angle = angle
+ camera.view_up = (0,0,1)
+ camera.focal_point = ptTo
+ camera.position = ptFrom
+ camera.elevation(elev)
+ camera.azimuth(azimuth)
+ camera.clipping_range = clipRange
+ return
+
+
+# ----------------------------------------------------------------------
+def setWindow(maya):
+ return
More information about the CIG-COMMITS
mailing list