[cig-commits] r15200 - in short/3D/PyLith/trunk/tests: . 1d/line2 1d/line3 2d 2d/quad4
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 11 14:03:40 PDT 2009
Author: brad
Date: 2009-06-11 14:03:39 -0700 (Thu, 11 Jun 2009)
New Revision: 15200
Added:
short/3D/PyLith/trunk/tests/2d/quad4/TestAxialDisp.py
short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py
short/3D/PyLith/trunk/tests/2d/quad4/TestShearDisp.py
short/3D/PyLith/trunk/tests/2d/quad4/axialdisp.cfg
short/3D/PyLith/trunk/tests/2d/quad4/axialdisp_gendb.py
short/3D/PyLith/trunk/tests/2d/quad4/axialdisp_soln.py
short/3D/PyLith/trunk/tests/2d/quad4/sheardisp.cfg
short/3D/PyLith/trunk/tests/2d/quad4/sheardisp_gendb.py
short/3D/PyLith/trunk/tests/2d/quad4/sheardisp_soln.py
Modified:
short/3D/PyLith/trunk/tests/1d/line2/Makefile.am
short/3D/PyLith/trunk/tests/1d/line3/Makefile.am
short/3D/PyLith/trunk/tests/2d/Makefile.am
short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am
short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py
short/3D/PyLith/trunk/tests/Makefile.am
Log:
Finished implementing some 2d/quad4 full-scale tests.
Modified: short/3D/PyLith/trunk/tests/1d/line2/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/Makefile.am 2009-06-11 20:35:36 UTC (rev 15199)
+++ short/3D/PyLith/trunk/tests/1d/line2/Makefile.am 2009-06-11 21:03:39 UTC (rev 15200)
@@ -38,6 +38,7 @@
dislocation-elastic_t0000000.vtk \
dislocation_t0000000.vtk
+
TESTS_ENVIRONMENT = $(PYTHON)
@@ -51,6 +52,7 @@
BUILT_SOURCES = export-data
clean-local: clean-data
+CLEANFILES = *.pyc
# End of file
Modified: short/3D/PyLith/trunk/tests/1d/line3/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line3/Makefile.am 2009-06-11 20:35:36 UTC (rev 15199)
+++ short/3D/PyLith/trunk/tests/1d/line3/Makefile.am 2009-06-11 21:03:39 UTC (rev 15200)
@@ -52,6 +52,7 @@
BUILT_SOURCES = export-data
clean-local: clean-data
+CLEANFILES = *.pyc
# End of file
Modified: short/3D/PyLith/trunk/tests/2d/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/2d/Makefile.am 2009-06-11 20:35:36 UTC (rev 15199)
+++ short/3D/PyLith/trunk/tests/2d/Makefile.am 2009-06-11 21:03:39 UTC (rev 15200)
@@ -11,7 +11,6 @@
#
SUBDIRS = \
- tri3 \
quad4
Modified: short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am 2009-06-11 20:35:36 UTC (rev 15199)
+++ short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am 2009-06-11 21:03:39 UTC (rev 15200)
@@ -14,42 +14,30 @@
check_SCRIPTS = testpylith.py
-noinst_PYTHON = \
- TestAxialPlaneStrain.py \
- TestShearPlaneStrain.py \
- TestDislocation.py \
- TestDislocation2.py
+dist_noinst_PYTHON = \
+ TestQuad4.py \
+ TestAxialDisp.py \
+ axialdisp_soln.py \
+ axialdisp_gendb.py \
+ TestShearDisp.py \
+ sheardisp_soln.py \
+ sheardisp_gendb.py
-noinst_DATA = \
- box_quad4_100m.exo \
+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
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 \
- dislocation_t0000000.vtk \
- dislocation-statevars-elastic_info.vtk \
- dislocation-statevars-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
+ 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
TESTS_ENVIRONMENT = $(PYTHON)
@@ -57,16 +45,15 @@
# 'export' the input files by performing a mock install
export_datadir = $(top_builddir)/tests/2d/quad4
-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/quad4/TestAxialDisp.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestAxialDisp.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestAxialDisp.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,92 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestAxialDisp.py
+##
+## @brief Test suite for testing pylith with 2-D axial extension.
+
+import numpy
+from TestQuad4 import TestQuad4
+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(TestQuad4):
+ """
+ Test suite for testing pylith with 2-D axial extension.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestQuad4.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/quad4/TestQuad4.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestQuad4.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,93 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestQuad4.py
+##
+## @brief Generic tests for problems using 1-D bar mesh.
+
+import unittest
+import numpy
+
+class TestQuad4(unittest.TestCase):
+ """
+ Generic tests for problems using 1-D bar mesh.
+ """
+
+ def setUp(self):
+ """
+ Setup for tests.
+ """
+ self.mesh = {'ncells': 64,
+ 'ncorners': 4,
+ 'nvertices': 81,
+ '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/quad4/TestShearDisp.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestShearDisp.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestShearDisp.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,92 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestShearDisp.py
+##
+## @brief Test suite for testing pylith with 2-D shear.
+
+import numpy
+from TestQuad4 import TestQuad4
+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(TestQuad4):
+ """
+ Test suite for testing pylith with 2-D shear extension.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestQuad4.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
Added: short/3D/PyLith/trunk/tests/2d/quad4/axialdisp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialdisp.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialdisp.cfg 2009-06-11 21:03:39 UTC (rev 15200)
@@ -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
+#fiatlagrange = 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.ElasticPlaneStrain
+
+[axialdisp.timedependent.materials.elastic]
+label = Elastic material
+id = 1
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# ----------------------------------------------------------------------
+# 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/quad4/axialdisp_gendb.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialdisp_gendb.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialdisp_gendb.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/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/quad4/axialdisp_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialdisp_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialdisp_soln.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/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 strain)
+sxx = 1.0e+7
+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
+#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
Added: short/3D/PyLith/trunk/tests/2d/quad4/sheardisp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/sheardisp.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/sheardisp.cfg 2009-06-11 21:03:39 UTC (rev 15200)
@@ -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
+#fiatlagrange = 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.FIATLagrange
+quadrature.cell.dimension = 2
+
+# ----------------------------------------------------------------------
+# 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/quad4/sheardisp_gendb.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/sheardisp_gendb.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/sheardisp_gendb.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/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/quad4/sheardisp_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/sheardisp_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/sheardisp_soln.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/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 strain)
+sxx = 0.0
+sxy = 1.0e+07
+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
+#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
Modified: short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py 2009-06-11 20:35:36 UTC (rev 15199)
+++ short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py 2009-06-11 21:03:39 UTC (rev 15200)
@@ -20,9 +20,12 @@
"""
suite = unittest.TestSuite()
- from TestAxialShearDisp import TestAxialShearDisp
- suite.addTest(unittest.makeSuite(TestAxialShearDisp))
+ from TestAxialDisp import TestAxialDisp
+ suite.addTest(unittest.makeSuite(TestAxialDisp))
+ from TestShearDisp import TestShearDisp
+ suite.addTest(unittest.makeSuite(TestShearDisp))
+
#from TestAxialShearTract import TestAxislShearTract
#suite.addTest(unittest.makeSuite(TestAxialShearTract))
Modified: short/3D/PyLith/trunk/tests/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/Makefile.am 2009-06-11 20:35:36 UTC (rev 15199)
+++ short/3D/PyLith/trunk/tests/Makefile.am 2009-06-11 21:03:39 UTC (rev 15200)
@@ -11,10 +11,8 @@
#
SUBDIRS = \
- petsc \
1d \
- 2d \
- 3d
+ 2d
# End of file
More information about the CIG-COMMITS
mailing list