[cig-commits] r15202 - in short/3D/PyLith/trunk: . pylith/utils tests tests/2d tests/2d/quad4 tests/2d/tri3
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 11 17:24:05 PDT 2009
Author: brad
Date: 2009-06-11 17:24:04 -0700 (Thu, 11 Jun 2009)
New Revision: 15202
Added:
short/3D/PyLith/trunk/tests/2d/tri3/TestAxialDisp.py
short/3D/PyLith/trunk/tests/2d/tri3/TestShearDisp.py
short/3D/PyLith/trunk/tests/2d/tri3/TestTri3.py
short/3D/PyLith/trunk/tests/2d/tri3/axialdisp.cfg
short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_gendb.py
short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_soln.py
short/3D/PyLith/trunk/tests/2d/tri3/mesh.exo
short/3D/PyLith/trunk/tests/2d/tri3/mesh.jou
short/3D/PyLith/trunk/tests/2d/tri3/sheardisp.cfg
short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_gendb.py
short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_soln.py
Removed:
short/3D/PyLith/trunk/tests/2d/tri3/TestShearPlaneStrain.py
short/3D/PyLith/trunk/tests/2d/tri3/axialplanestrain.cfg
short/3D/PyLith/trunk/tests/2d/tri3/box_tri3_100m.exo
short/3D/PyLith/trunk/tests/2d/tri3/mesh_tri3_100m.jou
short/3D/PyLith/trunk/tests/2d/tri3/shearplanestrain.cfg
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py
short/3D/PyLith/trunk/tests/2d/Makefile.am
short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py
short/3D/PyLith/trunk/tests/2d/tri3/Makefile.am
short/3D/PyLith/trunk/tests/2d/tri3/geometry.jou
short/3D/PyLith/trunk/tests/2d/tri3/matprops.spatialdb
short/3D/PyLith/trunk/tests/2d/tri3/testpylith.py
short/3D/PyLith/trunk/tests/README
Log:
Added some 2d/tri3 full-scale tests.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/TODO 2009-06-12 00:24:04 UTC (rev 15202)
@@ -32,6 +32,15 @@
Update Neumann
ArbitratySlipFn (final slip, start time)
full-scale testing
+ 2d/quad4
+ axiialtract
+ sheartract
+ dislocation2
+ 2d/tri3
+ axialdisp
+ sheardisp
+ dislocation
+ dislocation2
cleanup
Charles
Modified: short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py
===================================================================
--- short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -72,7 +72,7 @@
cellId = cellTypes[0]
if numpy.sum(cellTypes-cellId) != 0:
raise ValueError("Expecting cells to all be the same type.")
- if cellId == 5: # line3
+ if cellId == 5: # tri3
ncorners = 3
elif cellId == 3: # line2
ncorners = 2
Modified: short/3D/PyLith/trunk/tests/2d/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/2d/Makefile.am 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/Makefile.am 2009-06-12 00:24:04 UTC (rev 15202)
@@ -11,7 +11,8 @@
#
SUBDIRS = \
- quad4
+ quad4 \
+ tri3
# End of file
Modified: short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -12,14 +12,14 @@
## @file tests/2d/quad4/TestQuad4.py
##
-## @brief Generic tests for problems using 1-D bar mesh.
+## @brief Generic tests for problems using 2-D mesh.
import unittest
import numpy
class TestQuad4(unittest.TestCase):
"""
- Generic tests for problems using 1-D bar mesh.
+ Generic tests for problems using 2-D mesh.
"""
def setUp(self):
Modified: short/3D/PyLith/trunk/tests/2d/tri3/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/Makefile.am 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/Makefile.am 2009-06-12 00:24:04 UTC (rev 15202)
@@ -14,42 +14,40 @@
check_SCRIPTS = testpylith.py
-noinst_PYTHON = \
- TestAxialPlaneStrain.py \
- TestShearPlaneStrain.py \
- TestDislocation.py \
- TestDislocation2.py
+dist_noinst_PYTHON = \
+ TestTri3.py \
+ TestAxialDisp.py \
+ axialdisp_soln.py \
+ axialdisp_gendb.py \
+ TestShearDisp.py \
+ sheardisp_soln.py \
+ sheardisp_gendb.py
-noinst_DATA = \
- box_tri3_100m.exo \
+# TestDislocation.py \
+# dislocation_soln.py
+
+dist_noinst_DATA = \
+ mesh.exo \
matprops.spatialdb \
- axialplanestrain.cfg \
- shearplanestrain.cfg \
- dislocation.cfg \
- dislocation2.cfg \
- dislocation_slipLL.spatialdb \
- dislocation_slipRL.spatialdb \
- dislocation_sliptime.spatialdb
+ axialdisp.cfg \
+ sheardisp.cfg
+# dislocation.cfg
+
noinst_TMP = \
- axialplanestrain_t0000000.vtk \
- axialplanestrain-statevars-elastic_info.vtk \
- axialplanestrain-statevars-elastic_t0000000.vtk \
- shearplanestrain_t0000000.vtk \
- shearplanestrain-statevars-elastic_info.vtk \
- shearplanestrain-statevars-elastic_t0000000.vtk \
+ axial_disp.spatialdb \
+ axialdisp_t0000000.vtk \
+ axialdisp-elastic_info.vtk \
+ axialdisp-elastic_t0000000.vtk \
+ shear_disp.spatialdb \
+ sheardisp_t0000000.vtk \
+ sheardisp-elastic_info.vtk \
+ sheardisp-elastic_t0000000.vtk \
dislocation_t0000000.vtk \
- dislocation-statevars-elastic_info.vtk \
- dislocation-statevars-elastic_t0000000.vtk \
+ dislocation-elastic_info.vtk \
+ dislocation-elastic_t0000000.vtk \
dislocation-fault_info.vtk \
- dislocation-fault_t0000000.vtk \
- dislocation2_t0000000.vtk \
- dislocation2-statevars-elastic_info.vtk \
- dislocation2-statevars-elastic_t0000000.vtk \
- dislocation2-fault-pos_info.vtk \
- dislocation2-fault-pos_t0000000.vtk \
- dislocation2-fault-neg_info.vtk \
- dislocation2-fault-neg_t0000000.vtk
+ dislocation-fault_t0000000.vtk
TESTS_ENVIRONMENT = $(PYTHON)
@@ -57,16 +55,15 @@
# 'export' the input files by performing a mock install
export_datadir = $(top_builddir)/tests/2d/tri3
-export-data: $(noinst_PYTHON) $(noinst_DATA)
- for f in $(noinst_PYTHON) $(noinst_DATA); do $(install_sh_DATA) $(srcdir)/$$f $(export_datadir); done
+export-data: $(dist_noinst_PYTHON) $(dist_noinst_DATA)
+ for f in $(dist_noinst_PYTHON) $(dist_noinst_DATA); do $(install_sh_DATA) $(srcdir)/$$f $(export_datadir); done
-BUILT_SOURCES = export-data
+clean-data:
+ if [ "X$(top_srcdir)" != "X$(top_builddir)" ]; then for f in $(dist_noinst_PYTHON) $(dist_noinst_DATA) $(noinst_TMP); do $(RM) $(RM_FLAGS) $(export_datadir)/$$f; done; fi
-CLEANFILES = \
- $(export_datadir)/$(noinst_PYTHON) \
- $(export_datadir)/$(noinst_DATA) \
- $(export_datadir)/$(noinst_TMP) \
- *.pyc
+BUILT_SOURCES = export-data
+clean-local: clean-data
+CLEANFILES = *.pyc
# End of file
Added: short/3D/PyLith/trunk/tests/2d/tri3/TestAxialDisp.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/TestAxialDisp.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/TestAxialDisp.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,92 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/TestAxialDisp.py
+##
+## @brief Test suite for testing pylith with 2-D axial extension.
+
+import numpy
+from TestTri3 import TestTri3
+from axialdisp_soln import AnalyticalSoln
+from pylith.utils.VTKDataReader import has_vtk
+from pylith.utils.VTKDataReader import VTKDataReader
+
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class AxialApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="axialdisp")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Generate spatial databases
+ from axialdisp_gendb import GenerateDB
+ db = GenerateDB()
+ db.run()
+
+ # Run PyLith
+ app = AxialApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestAxialDisp(TestTri3):
+ """
+ Test suite for testing pylith with 2-D axial extension.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestTri3.setUp(self)
+ run_pylith()
+ self.outputRoot = "axialdisp"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+ return
+
+
+ def calcDisplacements(self, vertices):
+ """
+ Calculate displacement field given coordinates of vertices.
+ """
+ return self.soln.displacement(vertices)
+
+
+ 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
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/tri3/TestShearDisp.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/TestShearDisp.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/TestShearDisp.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,92 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/TestShearDisp.py
+##
+## @brief Test suite for testing pylith with 2-D shear.
+
+import numpy
+from TestTri3 import TestTri3
+from sheardisp_soln import AnalyticalSoln
+from pylith.utils.VTKDataReader import has_vtk
+from pylith.utils.VTKDataReader import VTKDataReader
+
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class ShearApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="sheardisp")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Generate spatial databases
+ from sheardisp_gendb import GenerateDB
+ db = GenerateDB()
+ db.run()
+
+ # Run PyLith
+ app = ShearApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestShearDisp(TestTri3):
+ """
+ Test suite for testing pylith with 2-D shear extension.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestTri3.setUp(self)
+ run_pylith()
+ self.outputRoot = "sheardisp"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+ return
+
+
+ def calcDisplacements(self, vertices):
+ """
+ Calculate displacement field given coordinates of vertices.
+ """
+ return self.soln.displacement(vertices)
+
+
+ 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
+
+
+# End of file
Deleted: short/3D/PyLith/trunk/tests/2d/tri3/TestShearPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/TestShearPlaneStrain.py 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/TestShearPlaneStrain.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -1,84 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file tests/2d/tri3/TestShearPlaneStrain.py
-##
-## @brief Test suite for testing pylith with shear in y-direction for
-## 2-D box.
-
-import unittest
-from pylith.utils.VTKDataReader import has_vtk
-from pylith.utils.VTKDataReader import VTKDataReader
-
-# Local version of PyLithApp
-from pylith.apps.PyLithApp import PyLithApp
-class ShearPlaneStrainApp(PyLithApp):
- def __init__(self):
- PyLithApp.__init__(self, name="shearplanestrain")
- return
-
-
-# Helper function to run PyLith
-def run_pylith():
- """
- Run pylith.
- """
- if not "done" in dir(run_pylith):
- app = ShearPlaneStrainApp()
- app.run()
- run_pylith.done = True
- return
-
-
-class TestShearPlaneStrain(unittest.TestCase):
- """
- Test suite for testing pylith with shear in y-direction for 2-D box.
- """
-
- def setUp(self):
- """
- Setup for test.
- """
- run_pylith()
- if has_vtk():
- self.reader = VTKDataReader()
- else:
- self.reader = None
- return
-
-
- def test_elastic_info(self):
- """
- Check elastic info.
- """
- if self.reader is None:
- return
-
- data = self.reader.read("axialplanestrain-statevars-elastic_info.vtk")
- return
-
-
- def test_soln(self):
- """
- Check solution (displacement) field.
- """
- return
-
-
- def test_elastic_statevars(self):
- """
- Check elastic state variables.
- """
- return
-
-
-# End of file
Added: short/3D/PyLith/trunk/tests/2d/tri3/TestTri3.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/TestTri3.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/TestTri3.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,93 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/TestTri3.py
+##
+## @brief Generic tests for problems using 2-D mesh.
+
+import unittest
+import numpy
+
+class TestTri3(unittest.TestCase):
+ """
+ Generic tests for problems using 2-D mesh.
+ """
+
+ def setUp(self):
+ """
+ Setup for tests.
+ """
+ self.mesh = {'ncells': 124,
+ 'ncorners': 3,
+ 'nvertices': 79,
+ 'spaceDim': 3,
+ 'tensorSize': 3}
+ return
+
+
+ def test_elastic_info(self):
+ """
+ Check elastic info.
+ """
+ if self.reader is None:
+ return
+
+ ncells= self.mesh['ncells']
+
+ filename = "%s-elastic_info.vtk" % self.outputRoot
+ from axialdisp_soln import p_mu,p_lambda,p_density
+
+ propMu = p_mu*numpy.ones( (ncells, 1), dtype=numpy.float64)
+ propLambda = p_lambda*numpy.ones( (ncells, 1), dtype=numpy.float64)
+ propDensity = p_density*numpy.ones( (ncells, 2), dtype=numpy.float64)
+
+ properties = {'mu': propMu,
+ 'lambda': propLambda,
+ 'density': propDensity}
+
+ from pylith.tests.PhysicalProperties import check_properties
+ check_properties(self, filename, self.mesh, properties)
+
+ return
+
+
+ def test_soln(self):
+ """
+ Check solution (displacement) field.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s_t0000000.vtk" % self.outputRoot
+ from pylith.tests.Solution import check_displacements
+ check_displacements(self, filename, self.mesh)
+
+ return
+
+
+ def test_elastic_statevars(self):
+ """
+ Check elastic state variables.
+ """
+ if self.reader is None:
+ return
+
+ filename = "%s-elastic_t0000000.vtk" % self.outputRoot
+
+ from pylith.tests.StateVariables import check_state_variables
+ stateVars = ["total_strain", "stress"]
+ check_state_variables(self, filename, self.mesh, stateVars)
+
+ return
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/tri3/axialdisp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/axialdisp.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/axialdisp.cfg 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,106 @@
+# -*- Python -*-
+[axialdisp]
+
+[axialdisp.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+[axialdisp.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solverlinear = 1
+#meshiocubit = 1
+#implicitelasticity = 1
+#quadrature2d = 1
+#fiatsimplex = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[axialdisp.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+[axialdisp.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 2
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[axialdisp.timedependent]
+dimension = 2
+bc = [x_neg,x_pos,y_neg]
+
+[axialdisp.timedependent.formulation.time_step]
+total_time = 0.0*s
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[axialdisp.timedependent]
+materials = [elastic]
+materials.elastic = pylith.materials.ElasticPlaneStress
+
+[axialdisp.timedependent.materials.elastic]
+label = Elastic material
+id = 1
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATSimplex
+quadrature.cell.shape = triangle
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+[axialdisp.timedependent.bc.x_pos]
+fixed_dof = [0]
+label = 20
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC +x edge
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+[axialdisp.timedependent.bc.x_neg]
+fixed_dof = [0]
+label = 21
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -x edge
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+[axialdisp.timedependent.bc.y_neg]
+fixed_dof = [1]
+label = 23
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -y edge
+db_initial.iohandler.filename = axial_disp.spatialdb
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[axialdisp.petsc]
+ksp_rtol = 1.0e-8
+pc_type = asm
+# Change the preconditioner settings (must turn off
+# shift_positive_definite and turn on shift_nonzero).
+sub_pc_factor_shift_positive_definite = 0
+sub_pc_factor_shift_nonzero =
+
+ksp_max_it = 100
+ksp_gmres_restart = 50
+#ksp_monitor = true
+#ksp_view = true
+#log_summary = true
+# start_in_debugger = true
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[axialdisp.problem.formulation.output.output.writer]
+filename = axialdisp.vtk
+
+[axialdisp.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = axialdisp-elastic.vtk
Added: short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_gendb.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_gendb.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_gendb.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/axialdisp_gendb.py
+##
+## @brief Python script to generate spatial database with displacement
+## boundary conditions for the axial displacement test.
+
+import numpy
+
+class GenerateDB(object):
+ """
+ Python object to generate spatial database with displacement
+ boundary conditions for the axial displacement test.
+ """
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ return
+
+
+ def run(self):
+ """
+ Generate the database.
+ """
+ # Domain
+ x = numpy.arange(-4000.0, 4000.1, 1000.0)
+ y = numpy.arange(-4000.0, 4000.1, 1000.0)
+ npts = x.shape[0]
+
+ xx = x * numpy.ones( (npts, 1), dtype=numpy.float64)
+ yy = y * numpy.ones( (npts, 1), dtype=numpy.float64)
+ xy = numpy.zeros( (npts**2, 2), dtype=numpy.float64)
+ xy[:,0] = numpy.ravel(xx)
+ xy[:,1] = numpy.ravel(numpy.transpose(yy))
+
+ from axialdisp_soln import AnalyticalSoln
+ soln = AnalyticalSoln()
+ disp = soln.displacement(xy)
+
+ from spatialdata.geocoords.CSCart import CSCart
+ cs = CSCart()
+ cs.inventory.spaceDim = 2
+ cs._configure()
+ data = {'locs': xy,
+ 'coordsys': cs,
+ 'data_dim': 2,
+ 'values': [{'name': "displacement-x",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,0])},
+ {'name': "displacement-y",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,1])}]}
+
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ io = SimpleIOAscii()
+ io.inventory.filename = "axial_disp.spatialdb"
+ io._configure()
+ io.write(data)
+ return
+
+# ======================================================================
+if __name__ == "__main__":
+ app = GenerateDB()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/axialdisp_soln.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/axialdisp_soln.py
+##
+## @brief Analytical solution to axial displacement problem.
+
+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 stress)
+sxx = 1.0e+7
+sxy = 0.0
+syy = 0.0
+szz = 0.0
+
+# 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
+#print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to axial/shear displacement problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs):
+ """
+ Compute displacement field at locations.
+ """
+ (npts, dim) = locs.shape
+ disp = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ disp[:,0] = exx*locs[:,0] + exy*locs[:,1]
+ disp[:,1] = eyy*locs[:,1] + exy*locs[:,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
Deleted: short/3D/PyLith/trunk/tests/2d/tri3/axialplanestrain.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/axialplanestrain.cfg 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/axialplanestrain.cfg 2009-06-12 00:24:04 UTC (rev 15202)
@@ -1,102 +0,0 @@
-# -*- Python -*-
-[axialplanestrain]
-
-# ----------------------------------------------------------------------
-# journal
-# ----------------------------------------------------------------------
-[axialplanestrain.journal.info]
-#timedependent = 1
-#implicit = 1
-#petsc = 1
-#solverlinear = 1
-#meshiocubit = 1
-#implicitelasticity = 1
-#quadrature2d = 1
-#fiatsimplex = 1
-
-# ----------------------------------------------------------------------
-# mesh_generator
-# ----------------------------------------------------------------------
-[axialplanestrain.mesh_generator]
-#debug = 1
-reader = pylith.meshio.MeshIOCubit
-
-[axialplanestrain.mesh_generator.reader]
-filename = box_tri3_100m.exo
-coordsys.space_dim = 2
-
-# ----------------------------------------------------------------------
-# problem
-# ----------------------------------------------------------------------
-[axialplanestrain.timedependent]
-dimension = 2
-bc = [x_neg,x_pos,y_neg,y_pos]
-
-[axialplanestrain.timedependent.formulation.time_step]
-total_time = 0.0*s
-
-# ----------------------------------------------------------------------
-# materials
-# ----------------------------------------------------------------------
-[axialplanestrain.timedependent]
-materials = [elastic]
-materials.elastic = pylith.materials.ElasticPlaneStrain
-
-[axialplanestrain.timedependent.materials.elastic]
-label = Elastic material
-id = 1
-db_properties.iohandler.filename = matprops.spatialdb
-quadrature.cell.shape = triangle
-
-# ----------------------------------------------------------------------
-# boundary conditions
-# ----------------------------------------------------------------------
-[axialplanestrain.timedependent.bc.x_pos]
-fixed_dof = [0]
-label = 20
-db.label = Dirichlet BC +x edge
-
-[axialplanestrain.timedependent.bc.x_neg]
-fixed_dof = [0]
-label = 21
-db.label = Dirichlet BC -x edge
-
-[axialplanestrain.timedependent.bc.y_pos]
-fixed_dof = [1]
-label = 22
-db = spatialdata.spatialdb.UniformDB
-db.label = Dirichlet BC +y edge
-db.values = [dof-1]
-db.data = [-1.0]
-
-[axialplanestrain.timedependent.bc.y_neg]
-fixed_dof = [1]
-label = 23
-db = spatialdata.spatialdb.UniformDB
-db.label = Dirichlet BC -y edge
-db.values = [dof-1]
-db.data = [+1.0]
-
-# ----------------------------------------------------------------------
-# PETSc
-# ----------------------------------------------------------------------
-[axialplanestrain.petsc]
-ksp_rtol = 1.0e-8
-pc_type = asm
-ksp_max_it = 100
-ksp_gmres_restart = 50
-#ksp_monitor = true
-#ksp_view = true
-#log_summary = true
-# start_in_debugger = true
-
-
-# ----------------------------------------------------------------------
-# output
-# ----------------------------------------------------------------------
-[axialplanestrain.problem.formulation.output.output.writer]
-filename = axialplanestrain.vtk
-
-[axialplanestrain.timedependent.materials.elastic.output]
-cell_filter = pylith.meshio.CellFilterAvg
-writer.filename = axialplanestrain-statevars-elastic.vtk
Deleted: short/3D/PyLith/trunk/tests/2d/tri3/box_tri3_100m.exo
===================================================================
(Binary files differ)
Modified: short/3D/PyLith/trunk/tests/2d/tri3/geometry.jou
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/geometry.jou 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/geometry.jou 2009-06-12 00:24:04 UTC (rev 15202)
@@ -2,14 +2,14 @@
# Create surface using vertices
# ----------------------------------------------------------------------
-# Block is 800m x 800m
-# -400 m <= x <= 400 km
-# -400 km <= y <= 400 km
+# Block is 8000m x 8000m
+# -4000 m <= x <= 4000 m
+# -4000 m <= y <= 4000 m
reset
-create vertex -400.0 -400.0 0.0
-create vertex -400.0 +400.0 0.0
-create vertex +400.0 +400.0 0.0
-create vertex +400.0 -400.0 0.0
+create vertex -4000.0 -4000.0 0.0
+create vertex -4000.0 +4000.0 0.0
+create vertex +4000.0 +4000.0 0.0
+create vertex +4000.0 -4000.0 0.0
create surface vertex 1 2 3 4
delete vertex all
@@ -18,20 +18,14 @@
# ----------------------------------------------------------------------
# Fault at x=0
-create vertex 0.0 -400.0 0.0
-create vertex 0.0 +400.0 0.0
+create vertex 0.0 -4000.0 0.0
+create vertex 0.0 +4000.0 0.0
split surface 1 across location vertex 9 location vertex 10
-# Fault at y=-100
-create vertex -400.0 -100.0 0.0
-create vertex +400.0 -100.0 0.0
+# Fault at x=-2000.0
+create vertex -2000.0 +4000.0 0.0
+create vertex -2000.0 -4000.0 0.0
split surface 2 across location vertex 13 location vertex 14
split surface 3 across location vertex 13 location vertex 14
-# Fault at y=+100
-create vertex -400.0 +100.0 0.0
-create vertex +400.0 +100.0 0.0
-split surface 4 across location vertex 18 location vertex 19
-split surface 6 across location vertex 18 location vertex 19
-
delete vertex all
Modified: short/3D/PyLith/trunk/tests/2d/tri3/matprops.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/matprops.spatialdb 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/matprops.spatialdb 2009-06-12 00:24:04 UTC (rev 15202)
@@ -2,7 +2,7 @@
SimpleDB {
num-values = 3
value-names = density vs vp
- value-units = kg/m^3 m/s m/s
+ value-units = kg/m**3 m/s m/s
num-locs = 1
data-dim = 0
space-dim = 2
Added: short/3D/PyLith/trunk/tests/2d/tri3/mesh.exo
===================================================================
(Binary files differ)
Property changes on: short/3D/PyLith/trunk/tests/2d/tri3/mesh.exo
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: short/3D/PyLith/trunk/tests/2d/tri3/mesh.jou
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/mesh.jou (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/mesh.jou 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,69 @@
+# ----------------------------------------------------------------------
+# Generate geometry
+# ----------------------------------------------------------------------
+playback 'geometry.jou'
+
+# ----------------------------------------------------------------------
+# Set discretization size
+# ----------------------------------------------------------------------
+surface all size 1000
+
+# ----------------------------------------------------------------------
+# Generate the mesh
+# ----------------------------------------------------------------------
+surface all scheme trimesh
+mesh surface all
+
+# ----------------------------------------------------------------------
+# Create blocks for materials
+# ----------------------------------------------------------------------
+block 1 surface 3 4 5
+block 1 name "elastic"
+
+# ----------------------------------------------------------------------
+# Create nodeset for faults
+# ----------------------------------------------------------------------
+group "fault_x" add node in curve 5
+nodeset 10 group fault_x
+nodeset 10 name "fault_x"
+
+group "fault_x2" add node in curve 10
+nodeset 11 group fault_x2
+nodeset 11 name "fault_x2"
+
+# ----------------------------------------------------------------------
+# Create nodeset for +x edge
+# ----------------------------------------------------------------------
+group "edge_xpos" add node in curve 3
+nodeset 20 group edge_xpos
+nodeset 20 name "edge xpos"
+
+# ----------------------------------------------------------------------
+# Create nodeset for -x edge
+# ----------------------------------------------------------------------
+group "edge_xneg" add node in curve 1
+nodeset 21 group edge_xneg
+nodeset 21 name "edge xneg"
+
+# ----------------------------------------------------------------------
+# Create nodeset for +y edge
+# ----------------------------------------------------------------------
+group "edge_ypos" add node in curve 8
+group "edge_ypos" add node in curve 11
+group "edge_ypos" add node in curve 14
+nodeset 22 group edge_ypos
+nodeset 22 name "edge ypos"
+
+# ----------------------------------------------------------------------
+# Create nodeset for -y edge
+# ----------------------------------------------------------------------
+group "edge_yneg" add node in curve 9
+group "edge_yneg" add node in curve 12
+group "edge_yneg" add node in curve 13
+nodeset 23 group edge_yneg
+nodeset 23 name "edge yneg"
+
+# ----------------------------------------------------------------------
+# Export exodus file
+# ----------------------------------------------------------------------
+export mesh "mesh.exo" dimension 2 overwrite
Deleted: short/3D/PyLith/trunk/tests/2d/tri3/mesh_tri3_100m.jou
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/mesh_tri3_100m.jou 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/mesh_tri3_100m.jou 2009-06-12 00:24:04 UTC (rev 15202)
@@ -1,85 +0,0 @@
-# ----------------------------------------------------------------------
-# Generate geometry
-# ----------------------------------------------------------------------
-playback 'geometry.jou'
-
-# ----------------------------------------------------------------------
-# Set discretization size
-# ----------------------------------------------------------------------
-surface all size 100
-
-# ----------------------------------------------------------------------
-# Generate the mesh
-# ----------------------------------------------------------------------
-surface all scheme tridelaunay
-mesh surface all
-
-# ----------------------------------------------------------------------
-# Improve triangle quality with smoothing
-# ----------------------------------------------------------------------
-surface all smooth scheme condition number beta 1.05 cpu 2
-smooth surface all
-
-# ----------------------------------------------------------------------
-# Create blocks for materials
-# ----------------------------------------------------------------------
-block 1 surface 5 7 9 11 8 10
-block 1 name "elastic"
-
-# ----------------------------------------------------------------------
-# Create nodeset for faults
-# ----------------------------------------------------------------------
-group "fault_x" add node in curve 13
-group "fault_x" add node in curve 21
-group "fault_x" add node in curve 20
-nodeset 10 group fault_x
-nodeset 10 name "fault_x"
-
-group "fault_ypos" add node in curve 18
-group "fault_ypos" add node in curve 23
-nodeset 11 group fault_ypos
-nodeset 11 name "fault_ypos"
-
-group "fault_yneg" add node in curve 10
-group "fault_yneg" add node in curve 15
-nodeset 12 group fault_yneg
-nodeset 12 name "fault_yneg"
-
-# ----------------------------------------------------------------------
-# Create nodeset for +x edge
-# ----------------------------------------------------------------------
-group "edge_xpos" add node in curve 17
-group "edge_xpos" add node in curve 25
-group "edge_xpos" add node in curve 24
-nodeset 20 group edge_xpos
-nodeset 20 name "edge xpos"
-
-# ----------------------------------------------------------------------
-# Create nodeset for -x edge
-# ----------------------------------------------------------------------
-group "edge_xneg" add node in curve 14
-group "edge_xneg" add node in curve 22
-group "edge_xneg" add node in curve 19
-nodeset 21 group edge_xneg
-nodeset 21 name "edge xneg"
-
-# ----------------------------------------------------------------------
-# Create nodeset for +y edge
-# ----------------------------------------------------------------------
-group "edge_ypos" add node in curve 7
-group "edge_ypos" add node in curve 8
-nodeset 22 group edge_ypos
-nodeset 22 name "edge ypos"
-
-# ----------------------------------------------------------------------
-# Create nodeset for -y edge
-# ----------------------------------------------------------------------
-group "edge_yneg" add node in curve 6
-group "edge_yneg" add node in curve 9
-nodeset 23 group edge_yneg
-nodeset 23 name "edge yneg"
-
-# ----------------------------------------------------------------------
-# Export exodus file
-# ----------------------------------------------------------------------
-export mesh "box_tri3_100m.exo" dimension 2 overwrite
Added: short/3D/PyLith/trunk/tests/2d/tri3/sheardisp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/sheardisp.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/sheardisp.cfg 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,113 @@
+# -*- Python -*-
+[sheardisp]
+
+[sheardisp.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+[sheardisp.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solverlinear = 1
+#meshiocubit = 1
+#implicitelasticity = 1
+#quadrature2d = 1
+#fiatsimplex = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[sheardisp.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+[sheardisp.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 2
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[sheardisp.timedependent]
+dimension = 2
+bc = [x_neg,x_pos,y_neg,y_pos]
+
+[sheardisp.timedependent.formulation.time_step]
+total_time = 0.0*s
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[sheardisp.timedependent]
+materials = [elastic]
+materials.elastic = pylith.materials.ElasticPlaneStrain
+
+[sheardisp.timedependent.materials.elastic]
+label = Elastic material
+id = 1
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATSimplex
+quadrature.cell.shape = triangle
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+[sheardisp.timedependent.bc.x_pos]
+fixed_dof = [1]
+label = 20
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC +x edge
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+[sheardisp.timedependent.bc.x_neg]
+fixed_dof = [1]
+label = 21
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -x edge
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+[sheardisp.timedependent.bc.y_pos]
+fixed_dof = [0]
+label = 22
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC +y edge
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+[sheardisp.timedependent.bc.y_neg]
+fixed_dof = [0]
+label = 23
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -y edge
+db_initial.iohandler.filename = shear_disp.spatialdb
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[sheardisp.petsc]
+ksp_rtol = 1.0e-8
+pc_type = asm
+# Change the preconditioner settings (must turn off
+# shift_positive_definite and turn on shift_nonzero).
+sub_pc_factor_shift_positive_definite = 0
+sub_pc_factor_shift_nonzero =
+
+ksp_max_it = 100
+ksp_gmres_restart = 50
+#ksp_monitor = true
+#ksp_view = true
+#log_summary = true
+# start_in_debugger = true
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[sheardisp.problem.formulation.output.output.writer]
+filename = sheardisp.vtk
+
+[sheardisp.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = sheardisp-elastic.vtk
Added: short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_gendb.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_gendb.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_gendb.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/sheardisp_gendb.py
+##
+## @brief Python script to generate spatial database with displacement
+## boundary conditions for the axial displacement test.
+
+import numpy
+
+class GenerateDB(object):
+ """
+ Python object to generate spatial database with displacement
+ boundary conditions for the axial displacement test.
+ """
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ return
+
+
+ def run(self):
+ """
+ Generate the database.
+ """
+ # Domain
+ x = numpy.arange(-4000.0, 4000.1, 1000.0)
+ y = numpy.arange(-4000.0, 4000.1, 1000.0)
+ npts = x.shape[0]
+
+ xx = x * numpy.ones( (npts, 1), dtype=numpy.float64)
+ yy = y * numpy.ones( (npts, 1), dtype=numpy.float64)
+ xy = numpy.zeros( (npts**2, 2), dtype=numpy.float64)
+ xy[:,0] = numpy.ravel(xx)
+ xy[:,1] = numpy.ravel(numpy.transpose(yy))
+
+ from sheardisp_soln import AnalyticalSoln
+ soln = AnalyticalSoln()
+ disp = soln.displacement(xy)
+
+ from spatialdata.geocoords.CSCart import CSCart
+ cs = CSCart()
+ cs.inventory.spaceDim = 2
+ cs._configure()
+ data = {'locs': xy,
+ 'coordsys': cs,
+ 'data_dim': 2,
+ 'values': [{'name': "displacement-x",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,0])},
+ {'name': "displacement-y",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,1])}]}
+
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ io = SimpleIOAscii()
+ io.inventory.filename = "shear_disp.spatialdb"
+ io._configure()
+ io.write(data)
+ return
+
+# ======================================================================
+if __name__ == "__main__":
+ app = GenerateDB()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/tri3/sheardisp_soln.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/tri3/axialdisp_soln.py
+##
+## @brief Analytical solution to axial displacement problem.
+
+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 stress)
+sxx = 0.0
+sxy = 1.0e+07
+syy = 0.0
+szz = 0.0
+
+# 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
+#print -exx*p_lambda/(p_lambda+2*p_mu)
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to axial/shear displacement problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs):
+ """
+ Compute displacement field at locations.
+ """
+ (npts, dim) = locs.shape
+ disp = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ disp[:,0] = exx*locs[:,0] + exy*locs[:,1]
+ disp[:,1] = eyy*locs[:,1] + exy*locs[:,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
Deleted: short/3D/PyLith/trunk/tests/2d/tri3/shearplanestrain.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/shearplanestrain.cfg 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/shearplanestrain.cfg 2009-06-12 00:24:04 UTC (rev 15202)
@@ -1,93 +0,0 @@
-# -*- Python -*-
-[shearplanestrain]
-
-# ----------------------------------------------------------------------
-# journal
-# ----------------------------------------------------------------------
-[shearplanestrain.journal.info]
-#timedependent = 1
-#implicit = 1
-#petsc = 1
-#solverlinear = 1
-#meshiocubit = 1
-#implicitelasticity = 1
-#quadrature2d = 1
-#fiatsimplex = 1
-
-# ----------------------------------------------------------------------
-# mesh_generator
-# ----------------------------------------------------------------------
-[shearplanestrain.mesh_generator]
-#debug = 1
-reader = pylith.meshio.MeshIOCubit
-
-[shearplanestrain.mesh_generator.reader]
-filename = box_tri3_100m.exo
-coordsys.space_dim = 2
-
-# ----------------------------------------------------------------------
-# problem
-# ----------------------------------------------------------------------
-[shearplanestrain.timedependent]
-dimension = 2
-bc = [y_neg,y_pos]
-
-[shearplanestrain.timedependent.formulation.time_step]
-total_time = 0.0*s
-
-# ----------------------------------------------------------------------
-# materials
-# ----------------------------------------------------------------------
-[shearplanestrain.timedependent]
-materials = [elastic]
-materials.elastic = pylith.materials.ElasticPlaneStrain
-
-[shearplanestrain.timedependent.materials.elastic]
-label = Elastic material
-id = 1
-db_properties.iohandler.filename = matprops.spatialdb
-quadrature.cell.shape = triangle
-
-# ----------------------------------------------------------------------
-# boundary conditions
-# ----------------------------------------------------------------------
-[shearplanestrain.timedependent.bc.y_neg]
-fixed_dof = [0, 1]
-label = 23
-db = spatialdata.spatialdb.UniformDB
-db.label = Dirichlet BC -y edge
-db.values = [dof-0,dof-1]
-db.data = [-1.0,0.0]
-
-[shearplanestrain.timedependent.bc.y_pos]
-fixed_dof = [0, 1]
-label = 22
-db = spatialdata.spatialdb.UniformDB
-db.label = Dirichlet BC +y edge
-db.values = [dof-0,dof-1]
-db.data = [+1.0,0.0]
-
-# ----------------------------------------------------------------------
-# PETSc
-# ----------------------------------------------------------------------
-[shearplanestrain.petsc]
-ksp_rtol = 1.0e-8
-pc_type = asm
-ksp_max_it = 100
-ksp_gmres_restart = 50
-#ksp_monitor = true
-#ksp_view = true
-#log_summary = true
-# start_in_debugger = true
-
-# ----------------------------------------------------------------------
-# output
-# ----------------------------------------------------------------------
-# Give basename for VTK domain output of solution over domain.
-[shearplanestrain.problem.formulation.output.output.writer]
-filename = shearplanestrain.vtk
-
-# Give basename for VTK output of state variables.
-[shearplanestrain.timedependent.materials.elastic.output]
-cell_filter = pylith.meshio.CellFilterAvg
-writer.filename = shearplanestrain-statevars-elastic.vtk
Modified: short/3D/PyLith/trunk/tests/2d/tri3/testpylith.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/tri3/testpylith.py 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/2d/tri3/testpylith.py 2009-06-12 00:24:04 UTC (rev 15202)
@@ -1,4 +1,4 @@
-#!/usr/bin/env nemesis
+#!/usr/bin/env python
#
# ======================================================================
#
@@ -20,11 +20,11 @@
"""
suite = unittest.TestSuite()
- from TestAxialPlaneStrain import TestAxialPlaneStrain
- suite.addTest(unittest.makeSuite(TestAxialPlaneStrain))
+ from TestAxialDisp import TestAxialDisp
+ suite.addTest(unittest.makeSuite(TestAxialDisp))
- from TestShearPlaneStrain import TestShearPlaneStrain
- suite.addTest(unittest.makeSuite(TestShearPlaneStrain))
+ from TestShearDisp import TestShearDisp
+ suite.addTest(unittest.makeSuite(TestShearDisp))
#from TestDislocation import TestDislocation
#suite.addTest(unittest.makeSuite(TestDislocation))
Modified: short/3D/PyLith/trunk/tests/README
===================================================================
--- short/3D/PyLith/trunk/tests/README 2009-06-11 23:34:30 UTC (rev 15201)
+++ short/3D/PyLith/trunk/tests/README 2009-06-12 00:24:04 UTC (rev 15202)
@@ -21,18 +21,18 @@
* axialextension
- - Dirichlet w/initial displacement
+ - Dirichlet w/initial displacement, SimpleDB
* dislocation
- - Dirichlet w/initial displacement
+ - Dirichlet w/initial displacement, SimpleDB
- Fault slip (1 fault, 1 event)
line3 (quadratic cells)
* axialextension
- - Dirichlet w/initial displacement
+ - Dirichlet w/initial displacement, SimpleDB
* dislocation
@@ -43,15 +43,48 @@
2-D tests
----------------------------------------------------------------------
-1. axial/shear (DirichletBC)
-2. axial/shear (DirichletBC, parallel)
-3. axial traction (Neumann, DirichletBC)
-4. single fault (static)
-5. single fault (quasi-static w/multiple ruptures)
-6. single fault (dynamic)
-7. two fault (static)
-8. two fault (static, parallel)
+quad4 (linear cells, plane strain)
+ * axialdisp
+
+ - Dirichlet w/initial displacement, SimpleDB
+
+ * sheardisp
+
+ - Dirichlet w/initial displacement, SimpleDB
+
+ * axialtract --TODO--
+ * sheartract --TODO--
+
+ * dislocation
+
+ - Dirichlet w/initial displacement, UniformDB
+ - Fault slip (1 fault, 1 event), UniformDB
+
+ * dislocation2 --TODO--
+
+ - Dirichlet w/initial displacement, UniformDB
+ - Fault slip (2 faults, 1 event), UniformDB
+
+ * dislocationdyn --TODO--
+ [single fault, dynamic]
+
+tri3 (linear cells, plane stress)
+
+ * axialdisp
+
+ - Dirichlet w/initial displacement, SimpleDB
+
+ * sheardisp
+
+ - Dirichlet w/initial displacement, SimpleDB
+
+ * dislocation --TODO--
+ [1 fault, quasi-static w/multiple ruptures, diff slip]
+
+ * dislocation2 --TODO--
+ [2 faults, static, parallel]
+
----------------------------------------------------------------------
3-D tests
----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list