[cig-commits] r15201 - in short/3D/PyLith/trunk: . pylith/tests pylith/utils tests/2d/quad4
brad at geodynamics.org
brad at geodynamics.org
Thu Jun 11 16:34:31 PDT 2009
Author: brad
Date: 2009-06-11 16:34:30 -0700 (Thu, 11 Jun 2009)
New Revision: 15201
Added:
short/3D/PyLith/trunk/tests/2d/quad4/dislocation_soln.py
Removed:
short/3D/PyLith/trunk/tests/2d/quad4/TestAxialPlaneStrain.py
short/3D/PyLith/trunk/tests/2d/quad4/TestShearPlaneStrain.py
short/3D/PyLith/trunk/tests/2d/quad4/axialplanestrain.cfg
short/3D/PyLith/trunk/tests/2d/quad4/axialshear_soln.py
short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp.cfg
short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp_gendb.py
short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipLL.spatialdb
short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipRL.spatialdb
short/3D/PyLith/trunk/tests/2d/quad4/dislocation_sliptime.spatialdb
short/3D/PyLith/trunk/tests/2d/quad4/shearplanestrain.cfg
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/pylith/tests/Fault.py
short/3D/PyLith/trunk/pylith/tests/StateVariables.py
short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py
short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am
short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py
short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg
short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py
Log:
More 2d/quad4 full-scale tests.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/TODO 2009-06-11 23:34:30 UTC (rev 15201)
@@ -2,6 +2,9 @@
CURRENT ISSUES/PRIORITIES
======================================================================
+BUGS
+ malloc during assembly when running in parallel
+
----------------------------------------------------------------------
RELEASE 1.4
----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/pylith/tests/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/tests/Fault.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/pylith/tests/Fault.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -35,6 +35,10 @@
# Check fault information
tolerance = 1.0e-5
+ from spatialdata.units.NondimElasticQuasistatic import NondimElasticQuasistatic
+ normalizer = NondimElasticQuasistatic()
+ normalizer._configure()
+
for name in fieldNames:
valuesE = testcase.calcFaultField(name, data['vertices'])
values = data['vertex_fields'][name]
@@ -47,9 +51,13 @@
testcase.assertEqual(nverticesE, nvertices)
testcase.assertEqual(dimE, dim)
+ scale = 1.0
+ if name == "traction_change":
+ scale *= normalizer.pressureScale().value
+
for i in xrange(dim):
ratio = numpy.abs(1.0 - values[:,i]/valuesE[:,i])
- diff = numpy.abs(values[:,i] - valuesE[:,i])
+ diff = numpy.abs(values[:,i] - valuesE[:,i]) / scale
mask = valuesE[:,i] != 0.0
okay = mask*(ratio < tolerance) + ~mask*(diff < tolerance)
if numpy.sum(okay) != nvertices:
Modified: short/3D/PyLith/trunk/pylith/tests/StateVariables.py
===================================================================
--- short/3D/PyLith/trunk/pylith/tests/StateVariables.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/pylith/tests/StateVariables.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -32,6 +32,10 @@
testcase.assertEqual(mesh['nvertices'], nvertices)
testcase.assertEqual(mesh['spaceDim'], spaceDim)
+ from spatialdata.units.NondimElasticQuasistatic import NondimElasticQuasistatic
+ normalizer = NondimElasticQuasistatic()
+ normalizer._configure()
+
# Check state variables
tolerance = 1.0e-5
@@ -47,9 +51,13 @@
testcase.assertEqual(ncellsE, ncells)
testcase.assertEqual(dimE, dim)
+ scale = 1.0
+ if name == "stress":
+ scale *= normalizer.pressureScale().value
+
for i in xrange(dim):
ratio = numpy.abs(1.0 - values[:,i]/valuesE[:,i])
- diff = numpy.abs(values[:,i] - valuesE[:,i])
+ diff = numpy.abs(values[:,i] - valuesE[:,i]) / scale
mask = valuesE[:,i] != 0.0
okay = mask*(ratio < tolerance) + ~mask*(diff < tolerance)
if numpy.sum(okay) != ncells:
Modified: short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py
===================================================================
--- short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/pylith/utils/VTKDataReader.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -72,11 +72,13 @@
cellId = cellTypes[0]
if numpy.sum(cellTypes-cellId) != 0:
raise ValueError("Expecting cells to all be the same type.")
- if cellId == 5:
+ if cellId == 5: # line3
ncorners = 3
- elif cellId == 3:
+ elif cellId == 3: # line2
ncorners = 2
- elif cellId == 255:
+ elif cellId == 9: # quad4
+ ncorners = 4
+ elif cellId == 255: # unknown?
ncorners = 1
else:
raise ValueError("Unknown VTK cell type '%d'." % cellId)
Modified: short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/Makefile.am 2009-06-11 23:34:30 UTC (rev 15201)
@@ -21,13 +21,16 @@
axialdisp_gendb.py \
TestShearDisp.py \
sheardisp_soln.py \
- sheardisp_gendb.py
+ sheardisp_gendb.py \
+ TestDislocation.py \
+ dislocation_soln.py
dist_noinst_DATA = \
mesh.exo \
matprops.spatialdb \
axialdisp.cfg \
- sheardisp.cfg
+ sheardisp.cfg \
+ dislocation.cfg
noinst_TMP = \
axial_disp.spatialdb \
@@ -37,7 +40,12 @@
shear_disp.spatialdb \
sheardisp_t0000000.vtk \
sheardisp-elastic_info.vtk \
- sheardisp-elastic_t0000000.vtk
+ sheardisp-elastic_t0000000.vtk \
+ dislocation_t0000000.vtk \
+ dislocation-elastic_info.vtk \
+ dislocation-elastic_t0000000.vtk \
+ dislocation-fault_info.vtk \
+ dislocation-fault_t0000000.vtk
TESTS_ENVIRONMENT = $(PYTHON)
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/TestAxialPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestAxialPlaneStrain.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestAxialPlaneStrain.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,179 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file tests/2d/quad4/TestAxialPlaneStrain.py
-##
-## @brief Test suite for testing pylith with axial compression in
-## y-direction for 2-D box.
-
-import unittest
-import numpy
-from pylith.utils.VTKDataReader import has_vtk
-from pylith.utils.VTKDataReader import VTKDataReader
-
-# Local version of PyLithApp
-from pylith.PyLithApp import PyLithApp
-class AxialPlaneStrainApp(PyLithApp):
- def __init__(self):
- PyLithApp.__init__(self, name="axialplanestrain")
- return
-
-
-# Helper function to run PyLith
-def run_pylith():
- """
- Run pylith.
- """
- if not "done" in dir(run_pylith):
- app = AxialPlaneStrainApp()
- app.run()
- run_pylith.done = True
- return
-
-
-class TestAxialPlaneStrain(unittest.TestCase):
- """
- Test suite for testing pylith with axial extension 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")
-
- # Check cells
- ncellsE = 30
- ncornersE = 4
- (ncells, ncorners) = data['cells'].shape
- self.assertEqual(ncellsE, ncells)
- self.assertEqual(ncornersE, ncorners)
-
- # Check vertices
- nverticesE = 42
- spaceDimE = 3
- (nvertices, spaceDim) = data['vertices'].shape
- self.assertEqual(nverticesE, nvertices)
- self.assertEqual(spaceDimE, spaceDim)
-
- # Check physical properties
- tolerance = 1.0e-5
- vsE = 3000.0
- vpE = 5291.502622129181
- densityE = 2500.0
-
- # Lame's constant mu (shear modulus)
- muE = densityE*vsE**2
- diff = numpy.abs(1.0 - data['cell_fields']['mu']/muE)
- okay = diff < tolerance
- if numpy.sum(okay) != ncells:
- print "Lame's constant mu: ",data['cell_fields']['mu']
- self.assertEqual(ncells, numpy.sum(okay))
-
- # Lame's constant lambda
- lambdaE = densityE*vpE**2 - 2*muE
- diff = numpy.abs(1.0 - data['cell_fields']['lambda']/lambdaE)
- okay = diff < tolerance
- if numpy.sum(okay) != ncells:
- print "Lame's constant lambda: ",data['cell_fields']['lambda']
- self.assertEqual(ncells, numpy.sum(okay))
-
- # Density
- diff = numpy.abs(1.0 - data['cell_fields']['density']/densityE)
- okay = diff < tolerance
- if numpy.sum(okay) != ncells:
- print "Density: ",data['cell_fields']['density']
- self.assertEqual(ncells, numpy.sum(okay))
- return
-
-
- def test_soln(self):
- """
- Check solution (displacement) field.
- """
- if self.reader is None:
- return
-
- data = self.reader.read("axialplanestrain_t0000000.vtk")
-
- # Check cells
- ncellsE = 30
- ncornersE = 4
- (ncells, ncorners) = data['cells'].shape
- self.assertEqual(ncellsE, ncells)
- self.assertEqual(ncornersE, ncorners)
-
- # Check vertices
- nverticesE = 42
- spaceDimE = 3
- vertices = data['vertices']
- (nvertices, spaceDim) = vertices.shape
- self.assertEqual(nverticesE, nvertices)
- self.assertEqual(spaceDimE, spaceDim)
-
- # Check displacement solution
- tolerance = 1.0e-5
- dispE = numpy.zeros( (nvertices, spaceDim), dtype=numpy.float64)
- dispE[:,1] = -0.004 * vertices[:,1]
-
- disp = data['vertex_fields']['displacements']
-
- # Check x displacements
- diff = numpy.abs(disp[:,0] - dispE[:,0])
- okay = diff < tolerance
- if numpy.sum(okay) != nvertices:
- print "Displacement field: ",disp
- self.assertEqual(nvertices, numpy.sum(okay))
-
- # Check y displacements
- mask = dispE[:,1] > 0.0
- diff = mask * numpy.abs(1.0 - disp[:,1] / dispE[:,1]) + \
- ~mask * numpy.abs(disp[:,1] - dispE[:,1])
- okay = diff < tolerance
- if numpy.sum(okay) != nvertices:
- print "Displacement field: ",disp
- self.assertEqual(nvertices, numpy.sum(okay))
-
- # Check z displacements
- diff = numpy.abs(disp[:,2] - dispE[:,2])
- okay = diff < tolerance
- if numpy.sum(okay) != nvertices:
- print "Displacement field: ",disp
- self.assertEqual(nvertices, numpy.sum(okay))
-
- return
-
-
- def test_elastic_statevars(self):
- """
- Check elastic state variables.
- """
- return
-
-
-# End of file
Modified: short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestDislocation.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -12,14 +12,18 @@
## @file tests/2d/quad4/TestDislocation.py
##
-## @brief Test suite for testing pylith with shear dislocation for 2-D
-## box.
+## @brief Test suite for testing pylith with 2-D axial extension.
-import unittest
+import numpy
+from TestQuad4 import TestQuad4
+from dislocation_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 application
-from pylith.PyLithApp import PyLithApp
-class DislocationApp(PyLithApp):
+# Local version of PyLithApp
+from pylith.apps.PyLithApp import PyLithApp
+class AxialApp(PyLithApp):
def __init__(self):
PyLithApp.__init__(self, name="dislocation")
return
@@ -31,30 +35,127 @@
Run pylith.
"""
if not "done" in dir(run_pylith):
- app = DislocationApp()
+ # Run PyLith
+ app = AxialApp()
app.run()
run_pylith.done = True
return
-class TestDislocation(unittest.TestCase):
+class TestDislocation(TestQuad4):
"""
- Test suite for testing pylith with shear dislocation for 2-D box.
+ Test suite for testing pylith with 2-D axial extension.
"""
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 = "dislocation"
+ if has_vtk():
+ self.reader = VTKDataReader()
+ self.soln = AnalyticalSoln()
+ else:
+ self.reader = None
+
return
- def test_disp(self):
+ def test_fault_info(self):
"""
- Check displacement field.
+ Check fault information.
"""
+ if self.reader is None:
+ return
+
+ filename = "%s-fault_info.vtk" % self.outputRoot
+ fields = ["normal_dir", "final_slip", "slip_time"]
+ 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 = ["cumulative_slip", "traction_change"]
+ 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.
+ """
+
+ normalDir = (-1.0, 0.0)
+ finalSlip = -2.0
+ slipTime = 0.0
+
+ nvertices = self.faultMesh['nvertices']
+
+ if name == "normal_dir":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = normalDir[0]
+ field[:,1] = normalDir[1]
+
+ elif name == "final_slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = finalSlip
+
+ elif name == "slip_time":
+ field = slipTime*numpy.zeros( (nvertices, 1), dtype=numpy.float64)
+
+ elif name == "cumulative_slip":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = finalSlip
+
+ elif name == "traction_change":
+ field = numpy.zeros( (nvertices, 3), dtype=numpy.float64)
+ field[:,0] = 0.0
+
+ else:
+ raise ValueError("Unknown fault field '%s'." % name)
+
+ return field
+
+
# End of file
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/TestShearPlaneStrain.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/TestShearPlaneStrain.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/TestShearPlaneStrain.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,84 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file tests/2d/quad4/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.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
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/axialplanestrain.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialplanestrain.cfg 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialplanestrain.cfg 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,105 +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
-importer = pylith.meshio.MeshIOCubit
-
-[axialplanestrain.mesh_generator.importer]
-filename = box_quad4_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.iohandler.filename = matprops.spatialdb
-quadrature = pylith.feassemble.quadrature.Quadrature2D
-quadrature.cell = pylith.feassemble.FIATLagrange
-quadrature.cell.dimension = 2
-quadrature.cell.quad_order = 2
-
-# ----------------------------------------------------------------------
-# 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/quad4/axialshear_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialshear_soln.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialshear_soln.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,82 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file tests/2d/quad4/axialsheardisp_soln.py
-##
-## @brief Analytical solution to axial/shear 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
-sxx = 1.0e+7
-sxy = 8.0e+6
-syy = 0.0
-
-# Uniform strain field
-exx = 1.0/(2*p_mu) * (sxx - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy))
-eyy = 1.0/(2*p_mu) * (syy - p_lambda/(3*p_lambda+2*p_mu) * (sxx+syy))
-exy = 1.0/(2*p_mu) * (sxy)
-
-# ----------------------------------------------------------------------
-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, 2), 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
-
-
- def stress(self, locs):
- """
- Compute stress field at locations.
- """
- (npts, dim) = locs.shape
- stress = numpy.zeros( (npts, 3), dtype=numpy.float64)
- stress[:,0] = self.sxx
- stress[:,1] = self.syy
- stress[:,2] = self.sxy
- return
-
-
-# End of file
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp.cfg 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp.cfg 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,113 +0,0 @@
-# -*- Python -*-
-[pylithapp]
-
-[pylithapp.launcher] # WARNING: THIS IS NOT PORTABLE
-command = mpirun -np ${nodes}
-
-# ----------------------------------------------------------------------
-# journal
-# ----------------------------------------------------------------------
-[pylithapp.journal.info]
-#timedependent = 1
-#implicit = 1
-#petsc = 1
-#solverlinear = 1
-#meshiocubit = 1
-#implicitelasticity = 1
-#quadrature2d = 1
-#fiatlagrange = 1
-
-# ----------------------------------------------------------------------
-# mesh_generator
-# ----------------------------------------------------------------------
-[pylithapp.mesh_generator]
-#debug = 1
-reader = pylith.meshio.MeshIOCubit
-
-[pylithapp.mesh_generator.reader]
-filename = mesh.exo
-coordsys.space_dim = 2
-
-# ----------------------------------------------------------------------
-# problem
-# ----------------------------------------------------------------------
-[pylithapp.timedependent]
-dimension = 2
-bc = [x_neg,x_pos,y_neg,y_pos]
-
-[pylithapp.timedependent.formulation.time_step]
-total_time = 0.0*s
-
-# ----------------------------------------------------------------------
-# materials
-# ----------------------------------------------------------------------
-[pylithapp.timedependent]
-materials = [elastic]
-materials.elastic = pylith.materials.ElasticPlaneStrain
-
-[pylithapp.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
-# ----------------------------------------------------------------------
-[pylithapp.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 = axialshear_disp.spatialdb
-
-[pylithapp.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 = axialshear_disp.spatialdb
-
-[pylithapp.timedependent.bc.y_pos]
-fixed_dof = [1]
-label = 22
-db_initial = spatialdata.spatialdb.SimpleDB
-db_initial.label = Dirichlet BC +y edge
-db_initial.iohandler.filename = axialshear_disp.spatialdb
-
-[pylithapp.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 = axialshear_disp.spatialdb
-
-# ----------------------------------------------------------------------
-# PETSc
-# ----------------------------------------------------------------------
-[pylithapp.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
-# ----------------------------------------------------------------------
-[pylithapp.problem.formulation.output.output.writer]
-filename = axialsheardisp.vtk
-
-[pylithapp.timedependent.materials.elastic.output]
-cell_filter = pylith.meshio.CellFilterAvgMesh
-writer.filename = axialsheardisp-elastic.vtk
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp_gendb.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp_gendb.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/axialsheardisp_gendb.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,76 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file tests/2d/quad4/gendb_axialsheardisp.py
-##
-## @brief Python script to generate spatial database with displacement
-## boundary conditions for the axial/shear displacement test.
-
-import numpy
-
-class GenerateDB(object):
- """
- Python object to generate spatial database with displacement
- boundary conditions for the axial/shear displacement test.
- """
-
- def __init__(self):
- 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 axialshear_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 = "axialshear_disp.spatialdb"
- io._configure()
- io.write(data)
- return
-
-# ======================================================================
-if __name__ == "__main__":
- app = GenerateDB()
- app.run()
-
-
-# End of file
Modified: short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation.cfg 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,6 +1,9 @@
# -*- Python -*-
[dislocation]
+[dislocation.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
@@ -12,17 +15,17 @@
#meshiocubit = 1
#implicitelasticity = 1
#quadrature2d = 1
-#fiatsimplex = 1
+#fiatlagrange = 1
# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
[dislocation.mesh_generator]
#debug = 1
-importer = pylith.meshio.MeshIOCubit
+reader = pylith.meshio.MeshIOCubit
-[dislocation.mesh_generator.importer]
-filename = box_quad4_100m.exo
+[dislocation.mesh_generator.reader]
+filename = mesh.exo
coordsys.space_dim = 2
# ----------------------------------------------------------------------
@@ -46,30 +49,28 @@
[dislocation.timedependent.materials.elastic]
label = Elastic material
id = 1
-db.iohandler.filename = matprops.spatialdb
-quadrature = pylith.feassemble.quadrature.Quadrature2D
+db_properties.iohandler.filename = matprops.spatialdb
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 2
-quadrature.cell.quad_order = 2
# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
[dislocation.timedependent.bc.x_pos]
-fixed_dof = [0, 1]
+fixed_dof = [0,1]
label = 20
-db = spatialdata.spatialdb.UniformDB
-db.label = Dirichlet BC +x edge
-db.values = [dof-0,dof-1]
-db.data = [0.0,+0.5]
+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,-1.0*m]
[dislocation.timedependent.bc.x_neg]
-fixed_dof = [0, 1]
+fixed_dof = [0,1]
label = 21
-db = spatialdata.spatialdb.UniformDB
-db.label = Dirichlet BC -x edge
-db.values = [dof-0,dof-1]
-db.data = [0.0,-0.5]
+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,+1.0*m]
# ----------------------------------------------------------------------
# faults
@@ -80,23 +81,32 @@
[dislocation.timedependent.interfaces.fault]
id = 2
label = 10
-quadrature = pylith.feassemble.quadrature.Quadrature1Din2D
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
-quadrature.cell.quad_order = 2
-mat_db.iohandler.filename = matprops.spatialdb
[dislocation.timedependent.interfaces.fault.eq_srcs.rupture.slip_function]
-slip.iohandler.filename = dislocation_slipLL.spatialdb
-slip_time.iohandler.filename = dislocation_sliptime.spatialdb
+slip = spatialdata.spatialdb.UniformDB
+slip.label = Final slip
+slip.values = [left-lateral-slip,fault-opening]
+slip.data = [-2.0*m,0.0*m]
+slip_time = spatialdata.spatialdb.UniformDB
+slip_time.label = Slip start time
+slip_time.values = [slip-time]
+slip_time.data = [0.0*s]
+
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
[dislocation.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
@@ -104,17 +114,16 @@
#log_summary = true
# start_in_debugger = true
+
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
-# Give basename for VTK domain output of solution over domain.
[dislocation.problem.formulation.output.output.writer]
filename = dislocation.vtk
-# Give basename for VTK output of state variables.
[dislocation.timedependent.materials.elastic.output]
-cell_filter = pylith.meshio.CellFilterAvg
-writer.filename = dislocation-statevars-elastic.vtk
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = dislocation-elastic.vtk
[dislocation.timedependent.interfaces.fault.output.writer]
filename = dislocation-fault.vtk
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipLL.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipLL.spatialdb 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipLL.spatialdb 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,14 +0,0 @@
-#SPATIAL.ascii 1
-SimpleDB {
- num-values = 2
- value-names = left-lateral-slip fault-opening
- value-units = m m
- num-locs = 1
- data-dim = 0
- space-dim = 2
- cs-data = cartesian {
- to-meters = 1.0
- space-dim = 2
- }
-}
-0.0 0.0 1.0 0.0
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipRL.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipRL.spatialdb 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation_slipRL.spatialdb 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,14 +0,0 @@
-#SPATIAL.ascii 1
-SimpleDB {
- num-values = 2
- value-names = left-lateral-slip fault-opening
- value-units = m m
- num-locs = 1
- data-dim = 0
- space-dim = 2
- cs-data = cartesian {
- to-meters = 1.0
- space-dim = 2
- }
-}
-0.0 0.0 -1.0 0.0
Deleted: short/3D/PyLith/trunk/tests/2d/quad4/dislocation_sliptime.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation_sliptime.spatialdb 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation_sliptime.spatialdb 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,14 +0,0 @@
-#SPATIAL.ascii 1
-SimpleDB {
- num-values = 1
- value-names = slip-time
- value-units = s
- num-locs = 1
- data-dim = 0
- space-dim = 2
- cs-data = cartesian {
- to-meters = 1.0
- space-dim = 2
- }
-}
-0.0 0.0 -1.0
Added: short/3D/PyLith/trunk/tests/2d/quad4/dislocation_soln.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/dislocation_soln.py (rev 0)
+++ short/3D/PyLith/trunk/tests/2d/quad4/dislocation_soln.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -0,0 +1,93 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/dislocation_soln.py
+##
+## @brief Analytical solution to dislocation (rigid motion).
+
+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
+#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, nlocsO):
+ """
+ Compute displacement field at locations.
+ """
+ (nlocs, dim) = locs.shape
+
+ disp = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ maskP = locs[:,0] >= 0.0
+ maskP[nlocsO:nlocs] = False
+ maskN = numpy.bitwise_and(locs[:,0] <= 0.0, ~maskP)
+ disp[:,1] = \
+ maskN*(+1.0) + \
+ maskP*(-1.0)
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (nlocs, dim) = locs.shape
+ strain = numpy.zeros( (nlocs, 3), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = exy
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (nlocs, dim) = locs.shape
+ stress = numpy.zeros( (nlocs, 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/quad4/shearplanestrain.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/shearplanestrain.cfg 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/shearplanestrain.cfg 2009-06-11 23:34:30 UTC (rev 15201)
@@ -1,96 +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
-importer = pylith.meshio.MeshIOCubit
-
-[shearplanestrain.mesh_generator.importer]
-filename = box_quad4_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.iohandler.filename = matprops.spatialdb
-quadrature = pylith.feassemble.quadrature.Quadrature2D
-quadrature.cell = pylith.feassemble.FIATLagrange
-quadrature.cell.dimension = 2
-quadrature.cell.quad_order = 2
-
-# ----------------------------------------------------------------------
-# 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/quad4/testpylith.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py 2009-06-11 21:03:39 UTC (rev 15200)
+++ short/3D/PyLith/trunk/tests/2d/quad4/testpylith.py 2009-06-11 23:34:30 UTC (rev 15201)
@@ -26,12 +26,9 @@
from TestShearDisp import TestShearDisp
suite.addTest(unittest.makeSuite(TestShearDisp))
- #from TestAxialShearTract import TestAxislShearTract
- #suite.addTest(unittest.makeSuite(TestAxialShearTract))
+ from TestDislocation import TestDislocation
+ suite.addTest(unittest.makeSuite(TestDislocation))
- #from TestDislocation import TestDislocation
- #suite.addTest(unittest.makeSuite(TestDislocation))
-
#from TestDislocation2 import TestDislocation2
#suite.addTest(unittest.makeSuite(TestDislocation2))
More information about the CIG-COMMITS
mailing list