[cig-commits] r16514 - short/3D/PyLith/trunk/tests/2d/quad4
surendra at geodynamics.org
surendra at geodynamics.org
Fri Apr 9 03:41:09 PDT 2010
Author: surendra
Date: 2010-04-09 03:41:08 -0700 (Fri, 09 Apr 2010)
New Revision: 16514
Added:
short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionOpening.py
short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearSliding.py
short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearStick.py
short/3D/PyLith/trunk/tests/2d/quad4/friction_opening.cfg
short/3D/PyLith/trunk/tests/2d/quad4/friction_opening_soln.py
short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding.cfg
short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding_soln.py
short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick.cfg
short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick_soln.py
Modified:
short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am
short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionCompression.py
short/3D/PyLith/trunk/tests/2d/quad4/friction_compression.cfg
short/3D/PyLith/trunk/tests/2d/quad4/friction_compression_soln.py
short/3D/PyLith/trunk/tests/2d/quad4/mesh.jou
short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py
Log:
More work on tests for .vtk files of friction full-scale 2D quad4 examples
Modified: short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am 2010-04-09 01:16:31 UTC (rev 16513)
+++ short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am 2010-04-09 10:41:08 UTC (rev 16514)
@@ -32,6 +32,9 @@
lgdeformtraction_soln.py \
TestFrictionCompression.py \
friction_compression_soln.py \
+ friction_shear_stick_soln.py \
+ friction_shear_sliding_soln.py \
+ friction_opening_soln.py \
testpylith.py
dist_noinst_DATA = \
@@ -45,7 +48,10 @@
dislocation2.cfg \
lgdeformrigidbody.cfg \
lgdeformtraction.cfg \
- friction_compression.cfg
+ friction_compression.cfg \
+ friction_shear_stick.cfg \
+ friction_shear_sliding.cfg \
+ friction_opening.cfg
noinst_TMP = \
axial_disp.spatialdb \
@@ -70,7 +76,16 @@
lgdeformtraction-elastic_t0000000.vtk \
friction_compression_t0000000.vtk \
friction_compression-elastic_info.vtk \
- friction_compression-elastic_t0000000.vtk
+ friction_compression-elastic_t0000000.vtk \
+ friction_shear_stick_t0000000.vtk \
+ friction_shear_stick-elastic_info.vtk \
+ friction_shear_stick-elastic_t0000000.vtk \
+ friction_shear_sliding_t0000000.vtk \
+ friction_shear_sliding-elastic_info.vtk \
+ friction_shear_sliding-elastic_t0000000.vtk \
+ friction_opening_t0000000.vtk \
+ friction_opening-elastic_info.vtk \
+ friction_opening-elastic_t0000000.vtk
TESTS_ENVIRONMENT = $(PYTHON)
Modified: short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionCompression.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionCompression.py 2010-04-09 01:16:31 UTC (rev 16513)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionCompression.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -19,6 +19,7 @@
from friction_compression_soln import AnalyticalSoln
from pylith.utils.VTKDataReader import has_vtk
from pylith.utils.VTKDataReader import VTKDataReader
+from pylith.tests.Fault import check_vertex_fields
# Local version of PyLithApp
from pylith.apps.PyLithApp import PyLithApp
@@ -51,6 +52,13 @@
Setup for test.
"""
TestQuad4.setUp(self)
+ self.mesh['nvertices'] = 81+9
+ self.nverticesO = 81
+ self.faultMesh = {'nvertices': 9,
+ 'spaceDim': 3,
+ 'ncells': 8,
+ 'ncorners': 2}
+
run_pylith()
self.outputRoot = "friction_compression"
if has_vtk():
@@ -61,11 +69,39 @@
return
+ def test_fault_info(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_info.vtk" % self.outputRoot
+ fields = ["strike_dir", "normal_dir", "initial_traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def test_fault_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_t0000000.vtk" % self.outputRoot
+ fields = ["slip", "traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
def calcDisplacements(self, vertices):
"""
Calculate displacement field given coordinates of vertices.
"""
- return self.soln.displacement(vertices)
+ return self.soln.displacement(vertices, self.nverticesO)
def calcStateVar(self, name, vertices, cells):
@@ -84,4 +120,44 @@
return stateVar
+ def calcFaultField(self, name, vertices):
+ """
+ Calculate fault info.
+ """
+
+ strikeDir = (0.0, -1.0)
+ normalDir = (-1.0, 0.0)
+ initialTraction = (0.0, -1.0e+10)
+
+ nvertices = self.faultMesh['nvertices']
+
+ if name == "strike_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = strikeDir[0]
+ field[:,1] = strikeDir[1]
+
+ elif name == "normal_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir[0]
+ field[:,1] = normalDir[1]
+
+ elif name == "initial_traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = initialTraction[0]
+ field[:,1] = initialTraction[1]
+
+ elif name == "slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+
+ elif name == "traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = 0.0
+ field[:,1] = -2.2e+10
+
+ else:
+ raise ValueError("Unknown fault field '%s'." % name)
+
+ return field
+
+
# End of file
Added: short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionOpening.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionOpening.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionOpening.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,164 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestFrictionOpening.py
+##
+## @brief Test suite for testing pylith with 2-D opening with friction.
+
+import numpy
+from TestQuad4 import TestQuad4
+from friction_opening_soln import AnalyticalSoln
+from pylith.utils.VTKDataReader import has_vtk
+from pylith.utils.VTKDataReader import VTKDataReader
+from pylith.tests.Fault import check_vertex_fields
+
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class OpeningApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="friction_opening")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Run PyLith
+ app = OpeningApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestFrictionOpening(TestQuad4):
+ """
+ Test suite for testing pylith with 2-D opening with friction.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestQuad4.setUp(self)
+ self.mesh['nvertices'] = 81+9
+ self.nverticesO = 81
+ self.faultMesh = {'nvertices': 9,
+ 'spaceDim': 3,
+ 'ncells': 8,
+ 'ncorners': 2}
+
+ run_pylith()
+ self.outputRoot = "friction_opening"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+ return
+
+
+ def test_fault_info(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_info.vtk" % self.outputRoot
+ fields = ["strike_dir", "normal_dir", "initial_traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def test_fault_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_t0000000.vtk" % self.outputRoot
+ fields = ["slip", "traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def calcDisplacements(self, vertices):
+ """
+ Calculate displacement field given coordinates of vertices.
+ """
+ return self.soln.displacement(vertices, self.nverticesO)
+
+
+ def calcStateVar(self, name, vertices, cells):
+ """
+ Calculate state variable.
+ """
+ ncells = self.mesh['ncells']
+ pts = numpy.zeros( (ncells, 3), dtype=numpy.float64)
+ if name == "total_strain":
+ stateVar = self.soln.strain(pts)
+ elif name == "stress":
+ stateVar = self.soln.stress(pts)
+ else:
+ raise ValueError("Unknown state variable '%s'." % name)
+
+ return stateVar
+
+
+ def calcFaultField(self, name, vertices):
+ """
+ Calculate fault info.
+ """
+
+ strikeDir = (0.0, -1.0)
+ normalDir = (-1.0, 0.0)
+ initialTraction = (0.0, -1.0e+10)
+ slip = (0.0, 4000.0)
+
+ nvertices = self.faultMesh['nvertices']
+
+ if name == "strike_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = strikeDir[0]
+ field[:,1] = strikeDir[1]
+
+ elif name == "normal_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir[0]
+ field[:,1] = normalDir[1]
+
+ elif name == "initial_traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = initialTraction[0]
+ field[:,1] = initialTraction[1]
+
+ elif name == "slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = slip[0]
+ field[:,1] = slip[1]
+
+ elif name == "traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+
+ else:
+ raise ValueError("Unknown fault field '%s'." % name)
+
+ return field
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearSliding.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearSliding.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearSliding.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,163 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestFrictionShearSliding.py
+##
+## @brief Test suite for testing pylith with 2-D shear sliding with friction.
+
+import numpy
+from TestQuad4 import TestQuad4
+from friction_shear_sliding_soln import AnalyticalSoln
+from pylith.utils.VTKDataReader import has_vtk
+from pylith.utils.VTKDataReader import VTKDataReader
+from pylith.tests.Fault import check_vertex_fields
+
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class ShearSlidingApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="friction_shear_sliding")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Run PyLith
+ app = ShearSlidingApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestFrictionShearSliding(TestQuad4):
+ """
+ Test suite for testing pylith with 2-D shear sliding with friction.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestQuad4.setUp(self)
+ self.mesh['nvertices'] = 81+9
+ self.nverticesO = 81
+ self.faultMesh = {'nvertices': 9,
+ 'spaceDim': 3,
+ 'ncells': 8,
+ 'ncorners': 2}
+
+ run_pylith()
+ self.outputRoot = "friction_shear_sliding"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+ return
+
+
+ def test_fault_info(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_info.vtk" % self.outputRoot
+ fields = ["strike_dir", "normal_dir", "initial_traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def test_fault_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_t0000000.vtk" % self.outputRoot
+ fields = ["slip", "traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def calcDisplacements(self, vertices):
+ """
+ Calculate displacement field given coordinates of vertices.
+ """
+ return self.soln.displacement(vertices, self.nverticesO)
+
+
+ def calcStateVar(self, name, vertices, cells):
+ """
+ Calculate state variable.
+ """
+ ncells = self.mesh['ncells']
+ pts = numpy.zeros( (ncells, 3), dtype=numpy.float64)
+ if name == "total_strain":
+ stateVar = self.soln.strain(pts)
+ elif name == "stress":
+ stateVar = self.soln.stress(pts)
+ else:
+ raise ValueError("Unknown state variable '%s'." % name)
+
+ return stateVar
+
+
+ def calcFaultField(self, name, vertices):
+ """
+ Calculate fault info.
+ """
+
+ strikeDir = (0.0, -1.0)
+ normalDir = (-1.0, 0.0)
+ initialTraction = (0.0, -1.0e+10)
+
+ nvertices = self.faultMesh['nvertices']
+
+ if name == "strike_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = strikeDir[0]
+ field[:,1] = strikeDir[1]
+
+ elif name == "normal_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir[0]
+ field[:,1] = normalDir[1]
+
+ elif name == "initial_traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = initialTraction[0]
+ field[:,1] = initialTraction[1]
+
+ elif name == "slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+
+ elif name == "traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = 0.0
+ field[:,1] = -2.2e+10
+
+ else:
+ raise ValueError("Unknown fault field '%s'." % name)
+
+ return field
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearStick.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearStick.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestFrictionShearStick.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,163 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestFrictionShearStick.py
+##
+## @brief Test suite for testing pylith with 2-D shear stick with friction.
+
+import numpy
+from TestQuad4 import TestQuad4
+from friction_shear_stick_soln import AnalyticalSoln
+from pylith.utils.VTKDataReader import has_vtk
+from pylith.utils.VTKDataReader import VTKDataReader
+from pylith.tests.Fault import check_vertex_fields
+
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class ShearStickApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="friction_shear_stick")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Run PyLith
+ app = ShearStickApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestFrictionShearStick(TestQuad4):
+ """
+ Test suite for testing pylith with 2-D shear stick with friction.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestQuad4.setUp(self)
+ self.mesh['nvertices'] = 81+9
+ self.nverticesO = 81
+ self.faultMesh = {'nvertices': 9,
+ 'spaceDim': 3,
+ 'ncells': 8,
+ 'ncorners': 2}
+
+ run_pylith()
+ self.outputRoot = "friction_shear_stick"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+ return
+
+
+ def test_fault_info(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_info.vtk" % self.outputRoot
+ fields = ["strike_dir", "normal_dir", "initial_traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def test_fault_data(self):
+ """
+ Check fault information.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_t0000000.vtk" % self.outputRoot
+ fields = ["slip", "traction"]
+ check_vertex_fields(self, filename, self.faultMesh, fields)
+
+ return
+
+
+ def calcDisplacements(self, vertices):
+ """
+ Calculate displacement field given coordinates of vertices.
+ """
+ return self.soln.displacement(vertices, self.nverticesO)
+
+
+ def calcStateVar(self, name, vertices, cells):
+ """
+ Calculate state variable.
+ """
+ ncells = self.mesh['ncells']
+ pts = numpy.zeros( (ncells, 3), dtype=numpy.float64)
+ if name == "total_strain":
+ stateVar = self.soln.strain(pts)
+ elif name == "stress":
+ stateVar = self.soln.stress(pts)
+ else:
+ raise ValueError("Unknown state variable '%s'." % name)
+
+ return stateVar
+
+
+ def calcFaultField(self, name, vertices):
+ """
+ Calculate fault info.
+ """
+
+ strikeDir = (0.0, -1.0)
+ normalDir = (-1.0, 0.0)
+ initialTraction = (0.0, -1.0e+10)
+
+ nvertices = self.faultMesh['nvertices']
+
+ if name == "strike_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = strikeDir[0]
+ field[:,1] = strikeDir[1]
+
+ elif name == "normal_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir[0]
+ field[:,1] = normalDir[1]
+
+ elif name == "initial_traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = initialTraction[0]
+ field[:,1] = initialTraction[1]
+
+ elif name == "slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+
+ elif name == "traction":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = 0.0
+ field[:,1] = -2.2e+10
+
+ else:
+ raise ValueError("Unknown fault field '%s'." % name)
+
+ return field
+
+
+# End of file
Modified: short/3D/PyLith/trunk/tests/2d/quad4/friction_compression.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_compression.cfg 2010-04-09 01:16:31 UTC (rev 16513)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_compression.cfg 2010-04-09 10:41:08 UTC (rev 16514)
@@ -173,4 +173,4 @@
# Give basename for VTK output of state variables.
[friction_compression.timedependent.materials.material.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
-writer.filename = friction_compression-statevars.vtk
+writer.filename = friction_compression-elastic.vtk
Modified: short/3D/PyLith/trunk/tests/2d/quad4/friction_compression_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_compression_soln.py 2010-04-09 01:16:31 UTC (rev 16513)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_compression_soln.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -37,7 +37,7 @@
exy = 1.0/(2*p_mu) * (sxy)
-#print exx,eyy,exy,ezz
+#print exx,eyy,exy,ezz,szz
#print -exx*p_lambda/(p_lambda+2*p_mu)
# ----------------------------------------------------------------------
@@ -50,15 +50,15 @@
return
- def displacement(self, locs):
+ def displacement(self, locs, nlocsO):
"""
Compute displacement field at locations.
"""
(nlocs, dim) = locs.shape
disp = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
- disp[:,0] = exx*locs[:,0] + exy*locs[:,1]
- disp[:,1] = eyy*locs[:,1] + exy*locs[:,0]
+ disp[:,0] = exx*(locs[:,0]+max(abs(locs[:,0])))
+ disp[:,1] = eyy*(locs[:,1]+max(abs(locs[:,1])))
return disp
Added: short/3D/PyLith/trunk/tests/2d/quad4/friction_opening.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_opening.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_opening.cfg 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,203 @@
+# -*- Python -*-
+[friction_opening]
+
+[friction_opening.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+# The settings below turn on journal info for the specified components.
+# If you want less output to stdout, you can turn these off.
+[friction_opening.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solvernonlinear = 1
+#meshioascii = 1
+#homogeneous = 1
+#elasticityimplicit = 1
+#fiatlagrange = 1
+#quadrature1d = 1
+#faultcohesivedyn = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+# The settings below control the mesh generation (importing mesh info).
+# Turn on debugging output for mesh generation.
+[friction_opening.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+# This component specification means we are using PyLith ASCII format,
+# and we then specify the filename and number of space dimensions for
+# the mesh.
+[friction_opening.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 2
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+# Specify the problem settings.
+# This is a time-dependent problem, so we select this as our problem type.
+# We select a total time of 0 sec, and a time step size of 1 sec, so we
+# are performing a single time step.
+# The spatial dimension for this problem is 2.
+# For an implicit formulation (using implicit.cfg), we will perform 1
+# implicit time step from t = -1.0 to t = 0.0 (elastic solution step).
+[friction_opening.timedependent]
+dimension = 2
+normalizer.length_scale = 1.0*m
+formulation = pylith.problems.Implicit
+formulation.solver = pylith.problems.SolverNonlinear
+
+# Set bc to an array with 2 boundary conditions: 'x_neg', 'x_pos'
+bc = [x_neg,x_pos]
+bc.x_pos=pylith.bc.DirichletBoundary
+
+# Set interfaces to an array with 1 fault: 'fault'.
+interfaces = [fault]
+
+
+[friction_opening.timedependent.formulation.time_step]
+total_time = 0.0*s
+dt = 1.0*s
+
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+# Specify the material information for the problem.
+# The material type is isotropic elastic formulated for plane strain.
+[friction_opening.timedependent.materials]
+material = pylith.materials.ElasticPlaneStrain
+
+[friction_opening.timedependent.materials.material]
+
+# We give a label of 'elastic material' to this material.
+label = elastic material
+
+# The cells associated with this material are given a material ID of 0
+# in the mesh file.
+id = 1
+
+# The properties for this material are given in the spatial database file
+# 'matelastic2D.spatialdb'.
+db_properties.iohandler.filename = matprops.spatialdb
+
+# Set cell type to quadrilateral (2-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+# Provide information on the boundary conditions.
+
+# Boundary conditions to be applied to the negative x-side of the mesh.
+[friction_opening.timedependent.bc.x_neg]
+
+# We are fixing the 0 (x) and 1 (y) degree of freedom.
+bc_dof = [0,1]
+
+# The nodes associated with this boundary condition have the name
+# 'x_neg' in the mesh file.
+label = 21
+
+# Boundary conditions to be applied to the positive x-side of the mesh.
+[friction_opening.timedependent.bc.x_pos]
+
+# We are fixing the 1 (y) and prescribing 0 (x) degree of freedom.
+bc_dof = [0,1]
+
+# The nodes associated with this boundary condition have the name
+# 'x_pos' in the mesh file.
+label = 20
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Dirichlet BC +x edge
+db_initial.values = [displacement-x,displacement-y]
+db_initial.data = [4000.0*m,0.0*m]
+
+# ----------------------------------------------------------------------
+# faults
+# ----------------------------------------------------------------------
+# Provide information on the fault (interface).
+[friction_opening.timedependent.interfaces]
+
+fault = pylith.faults.FaultCohesiveDyn
+
+# Define fault properties.
+[friction_opening.timedependent.interfaces.fault]
+
+# The nodes associated with this fault have the name 'fault' in the mesh file.
+label = 10
+
+# NOTE: It is possible to assign an ID number to a fault (e.g.,
+# 'id = 10'). Care must be taken when doing this, however, because the
+# assigned ID will become the material ID for the cohesive element.
+# This ID must not conflict with any of the material ID numbers for
+# volume elements. The default ID for a fault is 100. If you have a
+# fault in your mesh you must:
+# 1. If you create your own fault ID, make sure it does not conflict
+# with any of you material ID's.
+# 2. If you use the default fault ID, make sure that none of your
+# material ID's are equal to 100.
+
+# The quadrature for a 2D fault is 1D with a linear shape.
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+db_initial_tractions = spatialdata.spatialdb.UniformDB
+db_initial_tractions.label = "Initial fault tractions"
+db_initial_tractions.values = [traction-shear,traction-normal]
+db_initial_tractions.data = [0.0*Pa, -10000.0*MPa]
+
+friction.db_properties = spatialdata.spatialdb.UniformDB
+friction.db_properties.label = Static friction
+friction.db_properties.values = [friction-coefficient,cohesion]
+friction.db_properties.data = [0.6,0.0*Pa]
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[friction_opening.petsc]
+pc_type = asm
+
+# Change the preconditioner settings.
+sub_pc_factor_shift_type = nonzero
+
+ksp_rtol = 1.0e-8
+ksp_max_it = 100
+ksp_gmres_restart = 50
+snes_max_it = 200
+
+#ksp_monitor = true
+#ksp_view = true
+#ksp_converged_reason = true
+
+#snes_monitor = true
+#snes_view = true
+#snes_converged_reason = true
+
+#log_summary = true
+#start_in_debugger = true
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+# Give basename for VTK domain output of solution over domain.
+[friction_opening.problem.formulation.output.output.writer]
+filename = friction_opening.vtk
+
+# Give basename for VTK fault output.
+[friction_opening.timedependent.interfaces.fault.output]
+writer.filename = friction_opening-fault.vtk
+vertex_info_fields = [strike_dir,normal_dir,initial_traction]
+
+# Give basename for VTK output of state variables.
+[friction_opening.timedependent.materials.material.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = friction_opening-elastic.vtk
Added: short/3D/PyLith/trunk/tests/2d/quad4/friction_opening_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_opening_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_opening_soln.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/friction_compression_soln.py
+##
+## @brief Analytical solution to compression problem with friction.
+
+import numpy
+
+# Physical properties
+p_density = 2500.0
+p_vs = 3000.0
+p_vp = 5291.502622129181
+
+p_mu = p_density*p_vs**2
+p_lambda = p_density*p_vp**2 - 2*p_mu
+
+# Uniform stress field (plane strain)
+sxx = 0.0
+sxy = 0.0
+syy = 0.0
+szz = p_lambda/(2*p_lambda+2*p_mu)*(sxx+syy)
+
+# Uniform strain field
+exx = 1.0/(2*p_mu) * (sxx - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+eyy = 1.0/(2*p_mu) * (syy - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+ezz = 1.0/(2*p_mu) * (szz - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+
+exy = 1.0/(2*p_mu) * (sxy)
+
+#print exx,eyy,exy,ezz,szz
+#print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to friction_compression problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs, nlocsO):
+ """
+ Compute displacement field at locations.
+ """
+ (nlocs, dim) = locs.shape
+
+ disp = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ disp[0:nlocs/2,0] = 4000.0
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (npts, dim) = locs.shape
+ strain = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = exy
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (npts, dim) = locs.shape
+ stress = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = sxy
+ return stress
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding.cfg 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,239 @@
+# -*- Python -*-
+[friction_shear_sliding]
+
+[friction_shear_sliding.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+# The settings below turn on journal info for the specified components.
+# If you want less output to stdout, you can turn these off.
+[friction_shear_sliding.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solvernonlinear = 1
+#meshioascii = 1
+#homogeneous = 1
+#elasticityimplicit = 1
+#fiatlagrange = 1
+#quadrature1d = 1
+#faultcohesivedyn = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+# The settings below control the mesh generation (importing mesh info).
+# Turn on debugging output for mesh generation.
+[friction_shear_sliding.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+# This component specification means we are using PyLith ASCII format,
+# and we then specify the filename and number of space dimensions for
+# the mesh.
+[friction_shear_sliding.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 2
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+# Specify the problem settings.
+# This is a time-dependent problem, so we select this as our problem type.
+# We select a total time of 0 sec, and a time step size of 1 sec, so we
+# are performing a single time step.
+# The spatial dimension for this problem is 2.
+# For an implicit formulation (using implicit.cfg), we will perform 1
+# implicit time step from t = -1.0 to t = 0.0 (elastic solution step).
+[friction_shear_sliding.timedependent]
+dimension = 2
+normalizer.length_scale = 1.0*m
+formulation = pylith.problems.Implicit
+formulation.solver = pylith.problems.SolverNonlinear
+
+# Set bc to an array with 2 boundary conditions: 'x_neg', 'x_pos'
+bc = [x_neg,x_pos,y_pos_fault,y_neg_fault]
+bc.x_pos = pylith.bc.DirichletBoundary
+bc.y_pos_fault = pylith.bc.Neumann
+bc.y_neg_fault = pylith.bc.Neumann
+
+# Set interfaces to an array with 1 fault: 'fault'.
+interfaces = [fault]
+
+
+[friction_shear_sliding.timedependent.formulation.time_step]
+total_time = 0.0*s
+dt = 1.0*s
+
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+# Specify the material information for the problem.
+# The material type is isotropic elastic formulated for plane strain.
+[friction_shear_sliding.timedependent.materials]
+material = pylith.materials.ElasticPlaneStrain
+
+[friction_shear_sliding.timedependent.materials.material]
+
+# We give a label of 'elastic material' to this material.
+label = elastic material
+
+# The cells associated with this material are given a material ID of 0
+# in the mesh file.
+id = 1
+
+# The properties for this material are given in the spatial database file
+# 'matelastic2D.spatialdb'.
+db_properties.iohandler.filename = matprops.spatialdb
+
+# Set cell type to quadrilateral (2-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+# Provide information on the boundary conditions.
+
+# Boundary conditions to be applied to the negative x-side of the mesh.
+[friction_shear_sliding.timedependent.bc.x_neg]
+
+# We are fixing the 0 (x) and 1 (y) degrees of freedom.
+bc_dof = [0,1]
+
+# The nodes associated with this boundary condition have the name
+# 'x_neg' in the mesh file.
+label = 21
+
+# Boundary conditions to be applied to the positive x-side of the mesh.
+[friction_shear_sliding.timedependent.bc.x_pos]
+
+# We are fixing the 0 (x) and prescribing 1 (y) degree of freedom.
+bc_dof = [0,1]
+
+# The nodes associated with this boundary condition have the name
+# 'x_pos' in the mesh file.
+label = 20
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Dirichlet BC +x edge
+db_initial.values = [displacement-x,displacement-y]
+db_initial.data = [0.0*m,4000.0*m]
+
+
+
+# Boundary conditions to be applied to the positive y-side of the mesh.
+[friction_shear_sliding.timedependent.bc.y_pos_fault]
+
+# The nodes associated with this boundary condition have the name
+# 'y_pos_fault' in the mesh file.
+label = 25
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC +y edge
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [-6000*MPa, 0.0*MPa]
+
+# Set cell type to quadrilateral (1-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+# Boundary conditions to be applied to the negative y-side of the mesh.
+[friction_shear_sliding.timedependent.bc.y_neg_fault]
+
+# The nodes associated with this boundary condition have the name
+# 'y_neg_fault' in the mesh file.
+label = 24
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC -y edge
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [-6000*MPa, 0.0*MPa]
+
+# Set cell type to quadrilateral (1-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+# ----------------------------------------------------------------------
+# faults
+# ----------------------------------------------------------------------
+# Provide information on the fault (interface).
+[friction_shear_sliding.timedependent.interfaces]
+
+fault = pylith.faults.FaultCohesiveDyn
+
+# Define fault properties.
+[friction_shear_sliding.timedependent.interfaces.fault]
+
+# The nodes associated with this fault have the name 'fault' in the mesh file.
+label = 10
+
+# NOTE: It is possible to assign an ID number to a fault (e.g.,
+# 'id = 10'). Care must be taken when doing this, however, because the
+# assigned ID will become the material ID for the cohesive element.
+# This ID must not conflict with any of the material ID numbers for
+# volume elements. The default ID for a fault is 100. If you have a
+# fault in your mesh you must:
+# 1. If you create your own fault ID, make sure it does not conflict
+# with any of you material ID's.
+# 2. If you use the default fault ID, make sure that none of your
+# material ID's are equal to 100.
+
+# The quadrature for a 2D fault is 1D with a linear shape.
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+db_initial_tractions = spatialdata.spatialdb.UniformDB
+db_initial_tractions.label = "Initial fault tractions"
+db_initial_tractions.values = [traction-shear,traction-normal]
+db_initial_tractions.data = [0.0*Pa, -10000.0*MPa]
+
+friction.db_properties = spatialdata.spatialdb.UniformDB
+friction.db_properties.label = Static friction
+friction.db_properties.values = [friction-coefficient,cohesion]
+friction.db_properties.data = [0.6,0.0*Pa]
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[friction_shear_sliding.petsc]
+pc_type = asm
+
+# Change the preconditioner settings.
+sub_pc_factor_shift_type = nonzero
+
+ksp_rtol = 1.0e-8
+ksp_max_it = 100
+ksp_gmres_restart = 50
+snes_max_it = 200
+
+#ksp_monitor = true
+#ksp_view = true
+#ksp_converged_reason = true
+
+#snes_monitor = true
+#snes_view = true
+#snes_converged_reason = true
+
+#log_summary = true
+#start_in_debugger = true
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+# Give basename for VTK domain output of solution over domain.
+[friction_shear_sliding.problem.formulation.output.output.writer]
+filename = friction_shear_sliding.vtk
+
+# Give basename for VTK fault output.
+[friction_shear_sliding.timedependent.interfaces.fault.output]
+writer.filename = friction_shear_sliding-fault.vtk
+vertex_info_fields = [strike_dir,normal_dir,initial_traction]
+
+# Give basename for VTK output of state variables.
+[friction_shear_sliding.timedependent.materials.material.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = friction_shear_sliding-elastic.vtk
Added: short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_sliding_soln.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/friction_compression_soln.py
+##
+## @brief Analytical solution to compression problem with friction.
+
+import numpy
+
+# Physical properties
+p_density = 2500.0
+p_vs = 3000.0
+p_vp = 5291.502622129181
+
+p_mu = p_density*p_vs**2
+p_lambda = p_density*p_vp**2 - 2*p_mu
+
+# Uniform stress field (plane strain)
+sxx = -1.2e+10
+sxy = 0.0
+syy = 0.0
+szz = p_lambda/(2*p_lambda+2*p_mu)*(sxx+syy)
+
+# Uniform strain field
+exx = 1.0/(2*p_mu) * (sxx - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+eyy = 1.0/(2*p_mu) * (syy - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+ezz = 1.0/(2*p_mu) * (szz - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+
+exy = 1.0/(2*p_mu) * (sxy)
+
+print exx,eyy,exy,ezz,szz
+print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to friction_compression problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs, nlocsO):
+ """
+ Compute displacement field at locations.
+ """
+ (nlocs, dim) = locs.shape
+
+ disp = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ disp[:,0] = exx*(locs[:,0]+max(abs(locs[:,0])))
+ disp[:,1] = eyy*(locs[:,1]+max(abs(locs[:,1])))
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (npts, dim) = locs.shape
+ strain = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = exy
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (npts, dim) = locs.shape
+ stress = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = sxy
+ return stress
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick.cfg 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,253 @@
+# -*- Python -*-
+[friction_shear_stick]
+
+[friction_shear_stick.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+# The settings below turn on journal info for the specified components.
+# If you want less output to stdout, you can turn these off.
+[friction_shear_stick.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solvernonlinear = 1
+#meshioascii = 1
+#homogeneous = 1
+#elasticityimplicit = 1
+#fiatlagrange = 1
+#quadrature1d = 1
+#faultcohesivedyn = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+# The settings below control the mesh generation (importing mesh info).
+# Turn on debugging output for mesh generation.
+[friction_shear_stick.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+# This component specification means we are using PyLith ASCII format,
+# and we then specify the filename and number of space dimensions for
+# the mesh.
+[friction_shear_stick.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 2
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+# Specify the problem settings.
+# This is a time-dependent problem, so we select this as our problem type.
+# We select a total time of 0 sec, and a time step size of 1 sec, so we
+# are performing a single time step.
+# The spatial dimension for this problem is 2.
+# For an implicit formulation (using implicit.cfg), we will perform 1
+# implicit time step from t = -1.0 to t = 0.0 (elastic solution step).
+[friction_shear_stick.timedependent]
+dimension = 2
+normalizer.length_scale = 1.0*m
+formulation = pylith.problems.Implicit
+formulation.solver = pylith.problems.SolverNonlinear
+
+# Set bc to an array with 3 boundary conditions: 'x_neg', 'x_pos_disp', 'x_pos_tract'
+bc = [x_neg,x_pos_disp,x_pos_tract,y_pos_fault,y_neg_fault]
+bc.x_pos_tract = pylith.bc.Neumann
+bc.y_pos_fault = pylith.bc.Neumann
+bc.y_neg_fault = pylith.bc.Neumann
+
+# Set interfaces to an array with 1 fault: 'fault'.
+interfaces = [fault]
+
+[friction_shear_stick.timedependent.formulation.time_step]
+total_time = 0.0*s
+dt = 1.0*s
+
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+# Specify the material information for the problem.
+# The material type is isotropic elastic formulated for plane strain.
+[friction_shear_stick.timedependent.materials]
+material = pylith.materials.ElasticPlaneStrain
+
+[friction_shear_stick.timedependent.materials.material]
+
+# We give a label of 'elastic material' to this material.
+label = elastic material
+
+# The cells associated with this material are given a material ID of 0
+# in the mesh file.
+id = 1
+
+# The properties for this material are given in the spatial database file
+# 'matelastic2D.spatialdb'.
+db_properties.iohandler.filename = matprops.spatialdb
+
+# Set cell type to quadrilateral (2-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+# Provide information on the boundary conditions.
+
+# Boundary conditions to be applied to the negative x-side of the mesh.
+[friction_shear_stick.timedependent.bc.x_neg]
+
+# We are fixing the 0 (x) and 1 (y) degrees of freedom.
+bc_dof = [0,1]
+
+# The nodes associated with this boundary condition have the name
+# 'x_neg' in the mesh file.
+label = 21
+
+# Boundary conditions to be applied to the positive x-side of the mesh.
+[friction_shear_stick.timedependent.bc.x_pos_disp]
+
+# We are fixing the 0 (x) degree of freedom.
+bc_dof = [0]
+
+# The nodes associated with this boundary condition have the name
+# 'x_pos' in the mesh file.
+label = 20
+
+# Boundary conditions to be applied to the positive x-side of the mesh.
+[friction_shear_stick.timedependent.bc.x_pos_tract]
+
+# The nodes associated with this boundary condition have the name
+# 'x_pos' in the mesh file.
+label = 20
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC +x edge
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [1000.0*MPa,0.0*MPa]
+
+# Set cell type to quadrilateral (1-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+
+
+
+# Boundary conditions to be applied to the positive y-side of the mesh.
+[friction_shear_stick.timedependent.bc.y_pos_fault]
+
+# The nodes associated with this boundary condition have the name
+# 'y_pos_fault' in the mesh file.
+label = 25
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC +y edge
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [-1000.0*MPa, 0.0*MPa]
+
+# Set cell type to quadrilateral (1-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+# Boundary conditions to be applied to the negative y-side of the mesh.
+[friction_shear_stick.timedependent.bc.y_neg_fault]
+
+# The nodes associated with this boundary condition have the name
+# 'y_neg_fault' in the mesh file.
+label = 24
+
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC -y edge
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [-1000.0*MPa, 0.0*MPa]
+
+# Set cell type to quadrilateral (1-d Lagrange).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+
+
+# ----------------------------------------------------------------------
+# faults
+# ----------------------------------------------------------------------
+# Provide information on the fault (interface).
+[friction_shear_stick.timedependent.interfaces]
+
+fault = pylith.faults.FaultCohesiveDyn
+
+# Define fault properties.
+[friction_shear_stick.timedependent.interfaces.fault]
+
+# The nodes associated with this fault have the name 'fault' in the mesh file.
+label = 10
+
+# NOTE: It is possible to assign an ID number to a fault (e.g.,
+# 'id = 10'). Care must be taken when doing this, however, because the
+# assigned ID will become the material ID for the cohesive element.
+# This ID must not conflict with any of the material ID numbers for
+# volume elements. The default ID for a fault is 100. If you have a
+# fault in your mesh you must:
+# 1. If you create your own fault ID, make sure it does not conflict
+# with any of you material ID's.
+# 2. If you use the default fault ID, make sure that none of your
+# material ID's are equal to 100.
+
+# The quadrature for a 2D fault is 1D with a linear shape.
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 1
+
+db_initial_tractions = spatialdata.spatialdb.UniformDB
+db_initial_tractions.label = "Initial fault tractions"
+db_initial_tractions.values = [traction-shear,traction-normal]
+db_initial_tractions.data = [0.0*Pa, -10000.0*MPa]
+
+friction.db_properties = spatialdata.spatialdb.UniformDB
+friction.db_properties.label = Static friction
+friction.db_properties.values = [friction-coefficient,cohesion]
+friction.db_properties.data = [0.6,0.0*Pa]
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[friction_shear_stick.petsc]
+pc_type = asm
+
+# Change the preconditioner settings.
+sub_pc_factor_shift_type = nonzero
+
+ksp_rtol = 1.0e-8
+ksp_atol = 1.0e-15
+ksp_max_it = 100
+ksp_gmres_restart = 50
+snes_max_it = 200
+
+#ksp_monitor = true
+#ksp_view = true
+#ksp_converged_reason = true
+
+#snes_monitor = true
+#snes_view = true
+#snes_converged_reason = true
+
+#log_summary = true
+#start_in_debugger = true
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+# Give basename for VTK domain output of solution over domain.
+[friction_shear_stick.problem.formulation.output.output.writer]
+filename = friction_shear_stick.vtk
+
+# Give basename for VTK fault output.
+[friction_shear_stick.timedependent.interfaces.fault.output]
+writer.filename = friction_shear_stick-fault.vtk
+vertex_info_fields = [strike_dir,normal_dir,initial_traction]
+
+# Give basename for VTK output of state variables.
+[friction_shear_stick.timedependent.materials.material.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = friction_shear_stick-elastic.vtk
Added: short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/friction_shear_stick_soln.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/friction_compression_soln.py
+##
+## @brief Analytical solution to compression problem with friction.
+
+import numpy
+
+# Physical properties
+p_density = 2500.0
+p_vs = 3000.0
+p_vp = 5291.502622129181
+
+p_mu = p_density*p_vs**2
+p_lambda = p_density*p_vp**2 - 2*p_mu
+
+# Uniform stress field (plane strain)
+sxx = -1.2e+10
+sxy = 0.0
+syy = 0.0
+szz = p_lambda/(2*p_lambda+2*p_mu)*(sxx+syy)
+
+# Uniform strain field
+exx = 1.0/(2*p_mu) * (sxx - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+eyy = 1.0/(2*p_mu) * (syy - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+ezz = 1.0/(2*p_mu) * (szz - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy+szz))
+
+exy = 1.0/(2*p_mu) * (sxy)
+
+print exx,eyy,exy,ezz,szz
+print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to friction_compression problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs, nlocsO):
+ """
+ Compute displacement field at locations.
+ """
+ (nlocs, dim) = locs.shape
+
+ disp = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ disp[:,0] = exx*(locs[:,0]+max(abs(locs[:,0])))
+ disp[:,1] = eyy*(locs[:,1]+max(abs(locs[:,1])))
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (npts, dim) = locs.shape
+ strain = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = exy
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (npts, dim) = locs.shape
+ stress = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = sxy
+ return stress
+
+
+# End of file
Modified: short/3D/PyLith/trunk/tests/2d/quad4/mesh.jou
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/mesh.jou 2010-04-09 01:16:31 UTC (rev 16513)
+++ short/3D/PyLith/trunk/tests/2d/quad4/mesh.jou 2010-04-09 10:41:08 UTC (rev 16514)
@@ -72,6 +72,14 @@
nodeset 24 name "edge yneg nofault"
# ----------------------------------------------------------------------
+# Create nodeset for +y edge w/o overlap with fault
+# ----------------------------------------------------------------------
+group "edge_ypos_nofault" add node in edge_ypos
+group "edge_ypos_nofault" remove node in fault_x
+nodeset 25 group edge_ypos_nofault
+nodeset 25 name "edge ypos nofault"
+
+# ----------------------------------------------------------------------
# Export exodus file
# ----------------------------------------------------------------------
export mesh "mesh.exo" dimension 2 overwrite
Modified: short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py 2010-04-09 01:16:31 UTC (rev 16513)
+++ short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py 2010-04-09 10:41:08 UTC (rev 16514)
@@ -38,9 +38,18 @@
#from TestLgDeformTraction import TestTraction
#suite.addTest(unittest.makeSuite(TestTraction))
- from TestFrictionCompression import TestFrictionCompression
- suite.addTest(unittest.makeSuite(TestFrictionCompression))
+ #from TestFrictionCompression import TestFrictionCompression
+ #suite.addTest(unittest.makeSuite(TestFrictionCompression))
+
+ #from TestFrictionShearStick import TestFrictionShearStick
+ #suite.addTest(unittest.makeSuite(TestFrictionShearStick))
+ #from TestFrictionShearSliding import TestFrictionShearSliding
+ #suite.addTest(unittest.makeSuite(TestFrictionShearSliding))
+
+ from TestFrictionOpening import TestFrictionOpening
+ suite.addTest(unittest.makeSuite(TestFrictionOpening))
+
return suite
More information about the CIG-COMMITS
mailing list