[cig-commits] r16079 - short/3D/PyLith/branches/pylith-friction/tests/2d/quad4
brad at geodynamics.org
brad at geodynamics.org
Sun Dec 6 14:36:14 PST 2009
Author: brad
Date: 2009-12-06 14:36:13 -0800 (Sun, 06 Dec 2009)
New Revision: 16079
Added:
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction_soln.py
Modified:
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/Makefile.am
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformRigidBody.py
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction.cfg
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/testpylith.py
Log:
Added axial traction full-scale tests for large deformations.
Modified: short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/Makefile.am 2009-12-06 21:25:29 UTC (rev 16078)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/Makefile.am 2009-12-06 22:36:13 UTC (rev 16079)
@@ -28,6 +28,8 @@
TestLgDeformRigidBody.py \
rigidbody_soln.py \
rigidbody_gendb.py \
+ TestLgDeformTraction.py \
+ lgdeformtraction_soln.py \
testpylith.py
dist_noinst_DATA = \
@@ -59,7 +61,10 @@
rigidbody_disp.spatialdb \
lgdeformrigidbody_t0000000.vtk \
lgdeformrigidbody-elastic_info.vtk \
- lgdeformrigidbody-elastic_t0000000.vtk
+ lgdeformrigidbody-elastic_t0000000.vtk \
+ lgdeformtraction_t0000000.vtk \
+ lgdeformtraction-elastic_info.vtk \
+ lgdeformtraction-elastic_t0000000.vtk
TESTS_ENVIRONMENT = $(PYTHON)
Modified: short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformRigidBody.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformRigidBody.py 2009-12-06 21:25:29 UTC (rev 16078)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformRigidBody.py 2009-12-06 22:36:13 UTC (rev 16079)
@@ -49,7 +49,7 @@
class TestRigidBody(TestQuad4):
"""
- Test suite for testing pylith with 2-D axial extension.
+ Test suite for testing pylith with 2-D rigid body motion.
"""
def setUp(self):
Added: short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py 2009-12-06 22:36:13 UTC (rev 16079)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/TestRigidBody.py
+##
+## @brief Test suite for testing pylith with uniform tractions for
+## large deformations in 2-D.
+
+import numpy
+from TestQuad4 import TestQuad4
+from lgdeformtraction_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 TractionApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="lgdeformtraction")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Run PyLith
+ app = TractionApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestTraction(TestQuad4):
+ """
+ Test suite for testing pylith with 2-D axial tractions.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestQuad4.setUp(self)
+ run_pylith()
+ self.outputRoot = "lgdeformtraction"
+ 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
Modified: short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction.cfg 2009-12-06 21:25:29 UTC (rev 16078)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction.cfg 2009-12-06 22:36:13 UTC (rev 16079)
@@ -1,13 +1,13 @@
# -*- Python -*-
-[pylithapp]
+[lgdeformtraction]
-[pylithapp.launcher] # WARNING: THIS IS NOT PORTABLE
+[lgdeformtraction.launcher] # WARNING: THIS IS NOT PORTABLE
command = mpirun -np ${nodes}
# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
-[pylithapp.journal.info]
+[lgdeformtraction.journal.info]
#timedependent = 1
#implicit = 1
#petsc = 1
@@ -20,35 +20,35 @@
# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
-[pylithapp.mesh_generator]
+[lgdeformtraction.mesh_generator]
#debug = 1
reader = pylith.meshio.MeshIOCubit
-[pylithapp.mesh_generator.reader]
+[lgdeformtraction.mesh_generator.reader]
filename = mesh.exo
coordsys.space_dim = 2
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
-[pylithapp.timedependent]
+[lgdeformtraction.timedependent]
dimension = 2
bc = [x_neg,x_pos,y_neg,y_pos]
formulation = pylith.problems.ImplicitLgDeform
formulation.solver = pylith.problems.SolverNonlinear
-[pylithapp.timedependent.formulation.time_step]
+[lgdeformtraction.timedependent.formulation.time_step]
total_time = 0.0*s
# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------
-[pylithapp.timedependent]
+[lgdeformtraction.timedependent]
materials = [elastic]
materials.elastic = pylith.materials.ElasticPlaneStrain
-[pylithapp.timedependent.materials.elastic]
+[lgdeformtraction.timedependent.materials.elastic]
label = Elastic material
id = 1
db_properties.iohandler.filename = matprops.spatialdb
@@ -58,19 +58,19 @@
# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
-[pylithapp.timedependent.bc]
+[lgdeformtraction.timedependent.bc]
x_pos = pylith.bc.Neumann
y_pos = pylith.bc.Neumann
-[pylithapp.timedependent.bc.x_neg]
+[lgdeformtraction.timedependent.bc.x_neg]
bc_dof = [0]
label = 21
-[pylithapp.timedependent.bc.y_neg]
+[lgdeformtraction.timedependent.bc.y_neg]
bc_dof = [1]
label = 23
-[pylithapp.timedependent.bc.x_pos]
+[lgdeformtraction.timedependent.bc.x_pos]
label = 20
db_initial = spatialdata.spatialdb.UniformDB
db_initial.label = Neumann BC +x edge
@@ -79,7 +79,7 @@
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
-[pylithapp.timedependent.bc.y_pos]
+[lgdeformtraction.timedependent.bc.y_pos]
label = 22
db_initial = spatialdata.spatialdb.UniformDB
db_initial.label = Neumann BC +y edge
@@ -91,7 +91,7 @@
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
-[pylithapp.petsc]
+[lgdeformtraction.petsc]
ksp_rtol = 1.0e-8
pc_type = asm
# Change the preconditioner settings (must turn off
@@ -114,9 +114,9 @@
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
-[pylithapp.problem.formulation.output.output.writer]
+[lgdeformtraction.problem.formulation.output.output.writer]
filename = lgdeformtraction.vtk
-[pylithapp.timedependent.materials.elastic.output]
+[lgdeformtraction.timedependent.materials.elastic.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
writer.filename = lgdeformtraction-elastic.vtk
Added: short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction_soln.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction_soln.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/lgdeformtraction_soln.py 2009-12-06 22:36:13 UTC (rev 16079)
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/2d/quad4/lgdeformtraction_soln.py
+##
+## @brief Analytical solution to axial tractions with large
+## deformation formulation.
+
+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)
+sxx0 = -100.0e+6
+syy0 = -10.0e+6
+ux = -12.6626
+uy = 3.37568
+
+sxx = -1.00159e+8
+sxy = 0.0
+syy = -9.99578e+6
+szz = p_lambda/(2*p_lambda+2*p_mu)*(sxx+syy)
+
+# Uniform strain field
+exx = -0.00158157
+eyy = 0.000422049
+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 sxx,syy
+#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] = ux*(4000.0 + locs[:,0]) / 8000.0
+ disp[:,1] = uy*(4000.0 + locs[:,1]) / 8000.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/branches/pylith-friction/tests/2d/quad4/testpylith.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/testpylith.py 2009-12-06 21:25:29 UTC (rev 16078)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/testpylith.py 2009-12-06 22:36:13 UTC (rev 16079)
@@ -35,6 +35,9 @@
from TestLgDeformRigidBody import TestRigidBody
suite.addTest(unittest.makeSuite(TestRigidBody))
+ from TestLgDeformTraction import TestTraction
+ suite.addTest(unittest.makeSuite(TestTraction))
+
return suite
More information about the CIG-COMMITS
mailing list