[cig-commits] r13809 - short/3D/PyLith/trunk/playpen/postproc
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Jan 7 20:31:40 PST 2009
Author: willic3
Date: 2009-01-07 20:31:40 -0800 (Wed, 07 Jan 2009)
New Revision: 13809
Added:
short/3D/PyLith/trunk/playpen/postproc/vtkcff.cfg
short/3D/PyLith/trunk/playpen/postproc/vtkcff.py
Log:
Unfinished python and .cfg files to compute Coulomb failure function from
stresses in VTK files.
Added: short/3D/PyLith/trunk/playpen/postproc/vtkcff.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/postproc/vtkcff.cfg (rev 0)
+++ short/3D/PyLith/trunk/playpen/postproc/vtkcff.cfg 2009-01-08 04:31:40 UTC (rev 13809)
@@ -0,0 +1,11 @@
+# -*- Python -*-
+[vtkcff]
+
+# Top-level info
+vtk_reference_file = axialdisp-statevars_t000.vtk
+vtk_output_root = axialdisp-statevars-cff
+stress_index = 2
+stress_components_order = [0,1,2,3,4,5]
+friction_coeff = 0.6
+skempton_coeff = 0.5
+isotropic_poroelastic = False
Added: short/3D/PyLith/trunk/playpen/postproc/vtkcff.py
===================================================================
--- short/3D/PyLith/trunk/playpen/postproc/vtkcff.py (rev 0)
+++ short/3D/PyLith/trunk/playpen/postproc/vtkcff.py 2009-01-08 04:31:40 UTC (rev 13809)
@@ -0,0 +1,527 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file postproc/vtkcff
+
+## @brief Python application to compute the Coulomb Failure Function difference
+## for either a specified plane orientation or for optimally oriented planes.
+## Differences are computed from a constant initial stress state, a specified
+## state in a time series, or from the previous step in a time series.
+
+import math
+import numpy
+import os
+import re
+import glob
+from pyre.units.time import s
+
+from pyre.applications.Script import Script as Application
+
+class VtkCff(Application):
+ """
+ Python application to compute the Coulomb Failure Function (CFF) difference
+ for either a specified plane orientation or for optimally oriented planes.
+ Differences are computed from a zero initial stress state, a specified state
+ in a time series, or from the previous step in a time series.
+ """
+
+ class Inventory(Application.Inventory):
+ """
+ Python object for managing VtkCff facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing VtkCff facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b stress_ref_mode Whether to reference stresses to a constant state, a selected state, or the previous state.
+ ## @li \b orientation_mode Compute CFF on predefined plane or optimally-oriented planes.
+ ## @li \b vtk_input_root Root filename for VTK input files.
+ ## @li \b vtk_output_root Root filename for VTK output files.
+ ## @li \b vtk_stress_index Index indicating which VTK field array contains stresses.
+ ## @li \b vtk_stress_components_order Indices corresponding to Sxx,Syy,Szz,Sxy,Syz,Sxz.
+ ## @li \b friction_coeff Coefficient of friction.
+ ## @li \b skempton_coeff Skempton's pore pressure coefficient B.
+ ## @li \b initial_state_index Initial state time step number for initial state mode.
+ ## @li \b constant_state_values Initial stress values to use with constant state mode.
+ ## @li \b cff_plane_normal Normal to plane on which to compute CFF
+ ## @li \b isotropic_poroelastic Use isotropic poroelastic model instead of constant apparent friction model.
+
+ import pyre.inventory
+
+ stressRefMode = pyre.inventory.str("stress_ref_mode",
+ default="initial_state",
+ validator=pyre.inventory.choice(["constant_state",
+ "initial_state", "previous_state"]))
+ stressRefMode.meta['tip'] = "Stress state against which to compute differences."
+
+ orientationMode = pyre.inventory.str("orientation_mode",
+ default="optimally_oriented",
+ validator=pyre.inventory.choice(["optimally_oriented",
+ "predefined_plane"]))
+ orientationMode.meta['tip'] = "Compute CFF on predefined plane or optimally-oriented planes."
+
+ vtkInputRoot = pyre.inventory.str("vtk_input_root",
+ default="stress_t0001.vtk")
+ vtkInputRoot.meta['tip'] = "Root filename for VTK input files."
+
+ vtkOutputRoot = pyre.inventory.str("vtk_output_root", default="output.vtk")
+ vtkOutputRoot.meta['tip'] = "Root filename for VTK output files."
+
+ vtkStressIndex = pyre.inventory.int("vtk_stress_index", default=1)
+ vtkStressIndex.meta['tip'] = "Index indicating which VTK field array contains stresses."
+
+ vtkStressComponentsOrder = pyre.inventory.list("vtk_stress_components_order",
+ default=[0, 1, 2, 3, 4, 5])
+ vtkStressComponentsOrder.meta['tip'] = "Indices corresponding to Sxx, Syy, Szz, Sxy, Syz, Sxz."
+
+ frictionCoeff = pyre.inventory.float("friction_coeff", default=0.6)
+ frictionCoeff.meta['tip'] = "Coefficient of friction."
+
+ skemptonCoeff = pyre.inventory.float("skempton_coeff", default=0.5)
+ skemptonCoeff.meta['tip'] = "Skempton's pore pressure coefficient B."
+
+ initialStateIndex = pyre.inventory.int("initial_state_indes", default=0)
+ initialStateIndex.meta['tip'] = "Initial state time step number for initial state mode."
+
+ constantStateValues = pyre.inventory.list("constant_state_values",
+ default=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+ constantStateValues.meta['tip'] = "Initial stress values to use with constant state mode."
+
+ cffPlaneNormal = pyre.inventory.list("cff_plane_normal",
+ default=[1.0, 0.0, 0.0])
+ cffPlaneNormal.meta['tip'] = "Plane normal for predefined CFF plane."
+
+ isotropicPoroelastic = pyre.inventory.bool("isotropic_poroelastic",
+ default=False)
+ isotropicPoroelastic.meta['tip'] = "Use isotropic poroelastic model instead of constant apparent friction model."
+
+
+ class TooFewFilesError(IOError):
+ """
+ Exception raised when not enough VTK input files are found.
+ """
+
+ def __init__(self, value):
+ self.value = value
+
+
+ def __str__(self):
+ return repr(self.value)
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="vtkcff"):
+ Application.__init__(self, name)
+ self.vtkInputList = []
+ self.numVtkInputFiles = 0
+ self.vtkInputTimes = []
+ self.timeStampWidth = 0
+ self.refFileIndex = 0
+
+ self.numVertsPerCell = 0
+ self.numCells = 0
+ self.cellsArray = numpy.array([0])
+ self.verticesArray = numpy.array([0])
+ self.numVerts = 0
+ self.spaceDim = 0
+ self.cellType = ""
+ self.readMesh = False
+
+ self.numStressPoints = 0
+ return
+
+
+ def main(self):
+ import pdb
+ pdb.set_trace()
+ self._getFileInfo()
+ self._cffLoop()
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ Application._configure(self)
+ import pdb
+ pdb.set_trace()
+
+ # Set up info for input files
+ totalInputPath = os.path.normpath(
+ os.path.join(os.getcwd(), self.inventory.vtkInputRoot))
+ self.vtkInputDir = os.path.dirname(totalInputPath)
+ baseInputName = os.path.basename(totalInputPath)
+ baseInputNameLen = len(baseInputName)
+ if baseInputName.endswith(".vtk"):
+ baseInputNameStripped = baseInputName[0:baseInputNameLen-4]
+ else:
+ baseInputNameStripped = baseInputName
+ testFind = re.search('_t[0-9]*$', baseInputNameStripped)
+ if testFind != None:
+ timeInd = baseInputNameStripped.rfind(testFind.group(0))
+ self.vtkInputRoot = baseInputNameStripped[0:timeInd]
+ else:
+ self.vtkInputRoot = baseInputNameStripped
+
+ # Solution mode info
+ self.stressRefMode = self.inventory.stressRefMode
+ self.orientationMode = self.inventory.orientationMode
+ self.initialStateIndex = self.inventory.initialStateIndex
+ self.constantStateValues = self.inventory.constantStateValues
+ self.isotropicPoroelastic = self.inventory.isotropicPoroelastic
+
+ # Index information
+ self.vtkStressIndex = self.inventory.vtkStressIndex
+ self.vtkStressComponentsOrder = self.inventory.vtkStressComponentsOrder
+
+ # Parameters
+ self.frictionCoeff = self.inventory.frictionCoeff
+ self.skemptonCoeff = self.inventory.skemptonCoeff
+ self.cffPlaneNormal = numpy.array([float(self.inventory.cffPlaneNormal[0]),
+ float(self.inventory.cffPlaneNormal[1]),
+ float(self.inventory.cffPlaneNormal[2]),
+ dtype=numpy.float64)
+
+ # Set up info for output files
+ totalOutputPath = os.path.normpath(os.path.join(
+ os.getcwd(), self.inventory.vtkOutputRoot))
+ self.vtkOutputDir = os.path.dirname(totalOutputPath)
+ baseOutputName = os.path.basename(totalOutputPath)
+ baseOutputNameLen = len(baseOutputName)
+ if baseOutputName.endswith(".vtk"):
+ baseOutputNameStripped = baseOutputName[0:baseOutputNameLen-4]
+ else:
+ baseOutputNameStripped = baseOutputName
+ testFind = re.search('_t[0-9]*$', baseOutputNameStripped)
+ if testFind != None:
+ timeInd = baseOutputNameStripped.rfind(testFind.group(0))
+ self.vtkOutputRoot = baseOutputNameStripped[0:timeInd]
+ else:
+ self.vtkOutputRoot = baseOutputNameStripped
+
+ return
+
+
+ def _getFileInfo(self):
+ """
+ Find input files and set up filenames for input and output.
+ """
+
+ # Create list of input files and associated times
+ fileString = self.vtkInputRoot + "_t[0-9]*.vtk"
+ searchString = os.path.join(self.vtkInputDir, fileString)
+ self.vtkInputList = glob.glob(searchString)
+ self.numVtkInputFiles = len(self.vtkInputList)
+ index1 = self.vtkInputList[0].rfind("_t")
+ index2 = self.vtkInputList[0].rfind(".vtk")
+ self.timeStampWidth = index2 - index1 - 2
+ for vtkFile in self.vtkInputList:
+ timeString = vtkFile[index1+2:index2]
+ self.vtkInputTimes.append(float(timeString))
+
+ # Determine time index from which to start processing
+ if self.stressRefMode == "constant_state":
+ self.startIndex = 0
+ elif self.stressRefMode == "previous_state":
+ self.startIndex = 1
+ else:
+ self.startIndex = self.initialStateIndex + 1
+
+ # Raise exception if there aren't enough files to process
+ numFilesToProcess = self.numVtkInputFiles - self.startIndex - 1
+ if numFilesToProcess < 1:
+ try:
+ raise TooFewFilesError(numFilesToProcess)
+ except TooFewFilesError, e:
+ print 'Not enough files found for search string: ', searchString
+ print 'Number of files found: ', e.value
+
+ # Create output directory if it doesn't exist
+ if not os.path.isdir(self.vtkOutputDir):
+ os.mkdir(self.vtkOutputDir)
+
+ return
+
+
+ def _cffLoop(self):
+ """
+ Function to loop over input files, compute Cff relative to reference state,
+ and write the results to output files.
+ """
+
+ # Define function to use, depending on whether we are computing CFF on a
+ # predefined plane or an optimally-oriented plane.
+ if self.orientationMode == "optimally_oriented":
+ cffFunct = self._cffOop
+ else:
+ cffFunct = self._cffPredefined
+
+ # Define reference stresses unless differences are computed between each
+ # time step.
+ if self.stressRefMode == "constant_state":
+ testFile = os.path.join(self.vtkInputDir, self.vtkInputList[0]))
+ testStress = self._readStress(testFile)
+ refStress = numpy.tile(self.constantStateValues,
+ (self.numStressPoints, 1))
+ else if self.stressRefMode == "initial_state":
+ refFile = os.path.normpath(os.path.join(self.vtkInputDir,
+ self.vtkInputList[self.initialStateIndex]))
+ refStress = self._readStress(refFile)
+
+ # Loop over input VTK files.
+ for fileInd in range(self.startIndex, self.numVtkInputFiles):
+ timeStamp = self.vtkInputTimes[fileInd]
+ timeStampOut = int(timeStamp)
+ timeStampOutString = repr(timeStampOut).rjust(self.timeStampWidth, '0')
+ outputFileName = self.vtkOutputRoot + "_t" + timeStampOutString + ".vtk"
+ vtkOutputFile = os.path.join(self.vtkOutputDir, outputFileName)
+ if self.stressRefMode == "previous_state":
+ refFile = os.path.normpath(os.path.join(self.vtkInputDir,
+ self.vtkInputList[fileInd - 1]))
+ refStress = self._readStress(refFile)
+ newFile = os.path.normpath(os.path.join(self.vtkInputDir,
+ self.vtkInputList[fileInd]))
+ newStress = self._readStress(newFile)
+ cffFunct(refStress, newStress, vtkOutputFile)
+
+ return
+
+
+ def _readStress(self, vtkFile):
+ """
+ Function to read stresses from a file and store the info in a numpy array.
+ """
+ from enthought.mayavi.sources.vtk_file_reader import VTKFileReader
+ from enthought.tvtk.api import tvtk
+
+ reader = VTKFileReader()
+ reader.initialize(vtkFile)
+ data = reader.outputs[0]
+
+ # Get vertex and cell info if it hasn't already been done
+ if not self.readMesh:
+ cellVtk = data.get_cells()
+ self.numVertsPerCell = cellVtk._get_max_cell_size()
+ self.numCells = cellVtk.number_of_cells
+ cellArray = cellVtk.to_array()
+ self.cells = tvtk.CellArray()
+ self.cells.set_cells(self.numCells, cellArray)
+ self.vertArray = data._get_points().to_array()
+ self.cellType = data.get_cell_type(0)
+ (self.numVerts, self.spaceDim) = self.vertArray.shape
+ self.readMesh = True
+
+
+ # Get cell fields and extract stresses
+ cellData = data._get_cell_data()
+ numCellDataArrays = cellData._get_number_of_arrays()
+ stress = cellData.get_array(self.stressIndex).to_array()
+ (self.numStressPoints, numCols) = stress.shape
+
+ sxx = stress[:,self.stressComponentsOrder[0]]
+ syy = stress[:,self.stressComponentsOrder[1]]
+ szz = stress[:,self.stressComponentsOrder[2]]
+ sxy = stress[:,self.stressComponentsOrder[3]]
+ syz = stress[:,self.stressComponentsOrder[4]]
+ sxz = stress[:,self.stressComponentsOrder[5]]
+ stressOrdered = numpy.column_stack((sxx, syy, szz, sxy, syz, sxz))
+
+ return stressOrdered
+
+
+ def _princStress(self, stress):
+ """
+ Function to compute 3D principal stresses and sort them.
+ """
+ stressMat = numpy.array([(stress[0], stress[3], stress[5]),
+ (stress[3], stress[1], stress[4]),
+ (stress[5], stress[4], stress[2])],
+ dtype=numpy.float64)
+ (princStress, princAxes) = numpy.linalg.eigh(stressMat)
+ idx = princStress.argsort()
+ princStressOrdered = princStress[idx]
+ princAxesOrdered = princAxes[:,idx]
+ return princStressOrdered, princAxesOrdered
+
+
+ def _CffOop(self, refStress, newStress, vtkOutputFile):
+ """
+ Function to compute CFF for optimally-oriented planes and output results to
+ a file.
+ """
+ from enthought.tvtk.api import tvtk
+
+ beta = math.atan(self.frictionCoeff)/2.0
+ sinBeta = math.sin(beta)
+ cosBeta = math.cos(beta)
+ sin2Beta = math.sin(2.0*beta)
+ cos2Beta = math.cos(2.0*beta)
+ sinCosBeta = sinBeta * cosBeta
+ sinBetaSq = sinBeta * sinBeta
+ cosBetaSq = cosBeta * cosBeta
+ cff = numpy.empty((self.numStressPoints), dtype=numpy.float64)
+ failDir1 = numpy.empty((self.numStressPoints, self.spaceDim),
+ dtype=numpy.float64)
+ failDir2 = numpy.empty((self.numStressPoints, self.spaceDim),
+ dtype=numpy.float64)
+ effFrictionCoeff = self.frictionCoeff * (1.0 - self.skemptonCoeff)
+ presFac = -self.skemptonCoeff/3.0
+ vec1 = numpy.array([cosBeta, 0.0, sinBeta], dtype=numpy.float64)
+ vec2 = numpy.array([cosBeta, 0.0, -sinBeta], dtype=numpy.float64)
+
+ # Loop over stress points, compute total stress and the associated
+ # principal stresses, as well as the stress difference, and then
+ # compute CFF.
+ for point in xrange(self.numStressPoints):
+ totStress = newStress[point, :]
+ deltaStress = newStress[point, :] - refStress[point, :]
+ (totPrincStress, totPrincAxes) = self._princStress(totStress)
+ # Rotate stress changes into principal axis coordinate system for
+ # total stress.
+ deltaStressMat = numpy.array(
+ [(deltaStress[0], deltaStress[3], deltaStress[5]),
+ (deltaStress[3], deltaStress[1], deltaStress[4]),
+ (deltaStress[5], deltaStress[4], deltaStress[2])],
+ dtype=numpy.float64)
+ prod1 = numpy.dot(totPrincAxes, deltaStressMat)
+ deltaStressRot = numpy.dot(prod1, totPrincAxes.transpose())
+ s33 = deltaStressRot[0,0] * sinBetaSq - \
+ 2.0 * deltaStressRot[0,2] * sinCosBeta + \
+ deltaStressRot[2,2] * cosBetaSq
+ s13 = 0.5 * (deltaStressRot[2,2] - deltaStressRot[0,0]) * sin2Beta + \
+ deltaStressRot[0,2] * cos2Beta
+ deltaNormStress = deltaStress[0] + deltaStress[1] + deltaStress[2]
+ # Compute CFF depending on which model is used
+ if self.isotropicPoroelastic:
+ cff[point] = s13 + self.frictionCoeff * \
+ (s33 + presFac * deltaNormStress)
+ else:
+ cff[point] = s13 + effFrictionCoeff * s33
+
+ # Get failure planes by rotating vector in principal axis coordinate
+ # system into global coordinate system.
+ # NOTE: make sure I'm applying the rotation the right way.
+ failDir1[point,:] = numpy.dot(totPrincAxes.transpose(), vec1)
+ failDir2[point,:] = numpy.dot(totPrincAxes.transpose(), vec2)
+
+
+ # Set up mesh info for VTK file
+ mesh = tvtk.UnstructuredGrid(points=self.vertArray)
+ mesh.set_cells(self.cellType, self.cells)
+
+ # Add output fields and write VTK file
+ # For now, output cff as a scalar, failDir1 as a vector, and failDir2 as
+ # a general array.
+ cffName = "cff"
+ failDir1Name = "failure_plane_1"
+ failDir2Name = "failure_plane_2"
+ mesh.cell_data.scalars = cff
+ mesh.cell_data.scalars.name = cffName
+ mesh.cell_data.vectors = failDir1
+ mesh.cell_data.vectors.name = failDir1Name
+ mesh.cell_data.add_array(failDir2)
+ w = tvtk.UnstructuredGridWriter(file_name=vtkOutputFile, input=mesh)
+ w.write()
+
+ return
+
+
+ def _CffPredefined(self, refStress, newStress, vtkOutputFile):
+ """
+ Function to compute CFF for predefined planes and output results to a file.
+ """
+ from enthought.tvtk.api import tvtk
+
+ beta = math.atan(self.frictionCoeff)/2.0
+ sinBeta = math.sin(beta)
+ cosBeta = math.cos(beta)
+ sin2Beta = math.sin(2.0*beta)
+ cos2Beta = math.cos(2.0*beta)
+ sinCosBeta = sinBeta * cosBeta
+ sinBetaSq = sinBeta * sinBeta
+ cosBetaSq = cosBeta * cosBeta
+ cff = numpy.empty((self.numStressPoints), dtype=numpy.float64)
+ failDir1 = numpy.empty((self.numStressPoints, self.spaceDim),
+ dtype=numpy.float64)
+ failDir2 = numpy.empty((self.numStressPoints, self.spaceDim),
+ dtype=numpy.float64)
+ effFrictionCoeff = self.frictionCoeff * (1.0 - self.skemptonCoeff)
+ presFac = -self.skemptonCoeff/3.0
+ vec1 = numpy.array([cosBeta, 0.0, sinBeta], dtype=numpy.float64)
+ vec2 = numpy.array([cosBeta, 0.0, -sinBeta], dtype=numpy.float64)
+
+ # Loop over stress points, compute total stress and the associated
+ # principal stresses, as well as the stress difference, and then
+ # compute CFF.
+ for point in xrange(self.numStressPoints):
+ totStress = newStress[point, :]
+ deltaStress = newStress[point, :] - refStress[point, :]
+ (totPrincStress, totPrincAxes) = self._princStress(totStress)
+ # Rotate stress changes into principal axis coordinate system for
+ # total stress.
+ deltaStressMat = numpy.array(
+ [(deltaStress[0], deltaStress[3], deltaStress[5]),
+ (deltaStress[3], deltaStress[1], deltaStress[4]),
+ (deltaStress[5], deltaStress[4], deltaStress[2])],
+ dtype=numpy.float64)
+ prod1 = numpy.dot(totPrincAxes, deltaStressMat)
+ deltaStressRot = numpy.dot(prod1, totPrincAxes.transpose())
+ s33 = deltaStressRot[0,0] * sinBetaSq - \
+ 2.0 * deltaStressRot[0,2] * sinCosBeta + \
+ deltaStressRot[2,2] * cosBetaSq
+ s13 = 0.5 * (deltaStressRot[2,2] - deltaStressRot[0,0]) * sin2Beta + \
+ deltaStressRot[0,2] * cos2Beta
+ deltaNormStress = deltaStress[0] + deltaStress[1] + deltaStress[2]
+ # Compute CFF depending on which model is used
+ if self.isotropicPoroelastic:
+ cff[point] = s13 + self.frictionCoeff * \
+ (s33 + presFac * deltaNormStress)
+ else:
+ cff[point] = s13 + effFrictionCoeff * s33
+
+ # Get failure planes by rotating vector in principal axis coordinate
+ # system into global coordinate system.
+ # NOTE: make sure I'm applying the rotation the right way.
+ failDir1[point,:] = numpy.dot(totPrincAxes.transpose(), vec1)
+ failDir2[point,:] = numpy.dot(totPrincAxes.transpose(), vec2)
+
+
+ # Set up mesh info for VTK file
+ mesh = tvtk.UnstructuredGrid(points=self.vertArray)
+ mesh.set_cells(self.cellType, self.cells)
+
+ # Add output fields and write VTK file
+ # For now, output cff as a scalar, failDir1 as a vector, and failDir2 as
+ # a general array.
+ cffName = "cff"
+ failDir1Name = "failure_plane_1"
+ failDir2Name = "failure_plane_2"
+ mesh.cell_data.scalars = cff
+ mesh.cell_data.scalars.name = cffName
+ mesh.cell_data.vectors = failDir1
+ mesh.cell_data.vectors.name = failDir1Name
+ mesh.cell_data.add_array(failDir2)
+ w = tvtk.UnstructuredGridWriter(file_name=vtkOutputFile, input=mesh)
+ w.write()
+
+ return
+# ----------------------------------------------------------------------
+if __name__ == '__main__':
+ app = VtkCff()
+ app.run()
+
+# End of file
More information about the CIG-COMMITS
mailing list