[cig-commits] r16082 - in short/3D/PyLith/branches/pylith-friction: . libsrc/feassemble pylith/tests pylith/utils tests tests/2d/quad4 tests/3dnew/hex8
brad at geodynamics.org
brad at geodynamics.org
Mon Dec 7 12:52:02 PST 2009
Author: brad
Date: 2009-12-07 12:52:00 -0800 (Mon, 07 Dec 2009)
New Revision: 16082
Added:
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/Makefile.am
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestHex8.py
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformRigidBody.py
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformTraction.py
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformrigidbody.cfg
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction.cfg
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction_soln.py
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/matprops.spatialdb
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_gendb.py
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_soln.py
short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/testpylith.py
Modified:
short/3D/PyLith/branches/pylith-friction/configure.ac
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian2d_lgdeform.wxm
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian3d_lgdeform.wxm
short/3D/PyLith/branches/pylith-friction/pylith/tests/Solution.py
short/3D/PyLith/branches/pylith-friction/pylith/utils/VTKDataReader.py
short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py
short/3D/PyLith/branches/pylith-friction/tests/Makefile.am
short/3D/PyLith/branches/pylith-friction/tests/README
Log:
Added full-scale tests for large deformations in 3-D.
Modified: short/3D/PyLith/branches/pylith-friction/configure.ac
===================================================================
--- short/3D/PyLith/branches/pylith-friction/configure.ac 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/configure.ac 2009-12-07 20:52:00 UTC (rev 16082)
@@ -284,6 +284,7 @@
tests/2d/Makefile
tests/2d/tri3/Makefile
tests/2d/quad4/Makefile
+ tests/3dnew/hex8/Makefile
tests/petsc/Makefile
doc/Makefile
doc/developer/Makefile
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2009-12-07 20:52:00 UTC (rev 16082)
@@ -386,12 +386,19 @@
iBasis < numBasis;
++iBasis) {
const int iB = iBasis*spaceDim;
- const double N1 = basisDeriv[iQ+iB ];
- const double N2 = basisDeriv[iQ+iB+1];
+ const double Nip = basisDeriv[iQ+iB ];
+ const double Niq = basisDeriv[iQ+iB+1];
+
+ // Generated using Maxima (see jacobian2d_lgdeform.wxm)
_cellVector[iB ] -=
- wt*((1.0+l11)*N1*s11 + N2*l12*s22 + ((1.0+l11)*N2 + l12*N1)*s12);
+ wt*(l12*Niq*s22 +
+ ((l11+1)*Niq+l12*Nip)*s12 +
+ (l11+1)*Nip*s11);
+
_cellVector[iB+1] -=
- wt*(l21*N1*s11 + (1.0+l22)*N2*s22 + ((1.0+l22)*N1 + l21*N2)*s12);
+ wt*((l22+1)*Niq*s22 +
+ (l21*Niq+(l22+1)*Nip)*s12 +
+ l21*Nip*s11);
} // for
} // for
PetscLogFlops(numQuadPts*(1+numBasis*(numBasis*8+14)));
@@ -452,25 +459,34 @@
iBasis < numBasis;
++iBasis) {
const int iB = iBasis*spaceDim;
- const double N1 = basisDeriv[iQ+iB ];
- const double N2 = basisDeriv[iQ+iB+1];
- const double N3 = basisDeriv[iQ+iB+2];
+ const double Nip = basisDeriv[iQ+iB ];
+ const double Niq = basisDeriv[iQ+iB+1];
+ const double Nir = basisDeriv[iQ+iB+2];
+ // Generated using Maxima (see jacobian3d_lgdeform.wxm)
_cellVector[iB ] -=
- wt*((1.0+l11)*N1*s11 + l12*N2*s22 + l13*N3*s33 +
- ((1.0+l11)*N2 + l12*N1)*s12 +
- (l12*N3 + l13*N2)*s23 +
- ((1.0+l11)*N3 + l13*N1)*s13);
+ wt*(l13*Nir*s33 +
+ (l12*Nir+l13*Niq)*s23 +
+ l12*Niq*s22 +
+ ((l11+1)*Nir + l13*Nip)*s13 +
+ ((l11+1)*Niq+l12*Nip)*s12 +
+ (l11+1)*Nip*s11);
+
_cellVector[iB+1] -=
- wt*(l21*N1*s11 + (1.0+l22)*N2*s22 + l23*N3*s33 +
- ((1.0+l22)*N1 + l21*N2)*s12 +
- ((1.0+l22)*N3 + l23*N2)*s23 +
- (l21*N3 + l23*N1)*s13);
+ wt*(l23*Nir*s33 +
+ ((l22+1)*Nir+l23*Niq)*s23 +
+ (l22+1)*Niq*s22 +
+ (l21*Nir+l23*Nip)*s13 +
+ (l21*Niq+(l22+1)*Nip)*s12 +
+ l21*Nip*s11);
+
_cellVector[iB+2] -=
- wt*(l31*N1*s11 + l32*N2*s22 + (1.0+l33)*N3*s33 +
- (l31*N2 + l32*N1)*s12 +
- ((1.0+l33)*N2 + l32*N3)*s23 +
- ((1.0+l33)*N1 + l31*N3)*s13);
+ wt*((l33+1)*Nir*s33 +
+ (l32*Nir+(l33+1)*Niq)*s23 +
+ l32*Niq*s22 +
+ (l31*Nir+(l33+1)*Nip)*s13 +
+ (l31*Niq+l32*Nip)*s12 +
+ l31*Nip*s11);
} // for
} // for
PetscLogFlops(numQuadPts*(1+numBasis*(numBasis*18+3*27)));
@@ -668,6 +684,7 @@
assert(3 == cellDim);
assert(quadWts.size() == numQuadPts);
+ assert(6 == tensorSize);
const int numConsts = 21;
// Compute Jacobian for consistent tangent matrix
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian2d_lgdeform.wxm
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian2d_lgdeform.wxm 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian2d_lgdeform.wxm 2009-12-07 20:52:00 UTC (rev 16082)
@@ -41,5 +41,13 @@
y: transpose(Bi) . C . Bj;
/* [wxMaxima: input end ] */
+/* [wxMaxima: input start ] */
+Svec: matrix([s11, s22, s12]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+r: transpose(Bi) . transpose(Svec);
+/* [wxMaxima: input end ] */
+
/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian3d_lgdeform.wxm
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian3d_lgdeform.wxm 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/jacobian3d_lgdeform.wxm 2009-12-07 20:52:00 UTC (rev 16082)
@@ -1,5 +1,5 @@
/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
-/* [ Created with wxMaxima version 0.8.3a ] */
+/* [ Created with wxMaxima version 0.8.2 ] */
/* [wxMaxima: input start ] */
Bnli: matrix([Nip,0,0],
@@ -72,5 +72,13 @@
y: transpose(Bi) . C . Bj;
/* [wxMaxima: input end ] */
+/* [wxMaxima: input start ] */
+Svec: matrix([s11,s22,s33,s12,s23,s13]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+r: transpose(Bi) . transpose(Svec);
+/* [wxMaxima: input end ] */
+
/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$
Modified: short/3D/PyLith/branches/pylith-friction/pylith/tests/Solution.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/tests/Solution.py 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/pylith/tests/Solution.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -22,11 +22,6 @@
"""
data = testcase.reader.read(filename)
- # Check cells
- (ncells, ncorners) = data['cells'].shape
- testcase.assertEqual(mesh['ncells'], ncells)
- testcase.assertEqual(mesh['ncorners'], ncorners)
-
# Check vertices
(nvertices, spaceDim) = data['vertices'].shape
testcase.assertEqual(mesh['nvertices'], nvertices)
Modified: short/3D/PyLith/branches/pylith-friction/pylith/utils/VTKDataReader.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/utils/VTKDataReader.py 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/pylith/utils/VTKDataReader.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -74,6 +74,8 @@
raise ValueError("Expecting cells to all be the same type.")
if cellId == 5: # tri3
ncorners = 3
+ elif cellId == 12: # tri3
+ ncorners = 8
elif cellId == 3: # line2
ncorners = 2
elif cellId == 9: # quad4
Modified: short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/tests/2d/quad4/TestLgDeformTraction.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -10,7 +10,7 @@
# ----------------------------------------------------------------------
#
-## @file tests/2d/quad4/TestRigidBody.py
+## @file tests/2d/quad4/TestLgDeformTraction.py
##
## @brief Test suite for testing pylith with uniform tractions for
## large deformations in 2-D.
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/Makefile.am (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/Makefile.am 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,64 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+TESTS = testpylith.py
+
+check_SCRIPTS = testpylith.py
+
+dist_noinst_PYTHON = \
+ TestHex8.py \
+ TestLgDeformRigidBody.py \
+ rigidbody_soln.py \
+ rigidbody_gendb.py \
+ TestLgDeformTraction.py \
+ lgdeformtraction_soln.py \
+ testpylith.py
+
+dist_noinst_DATA = \
+ geometry.jou \
+ mesh.jou \
+ mesh.exo \
+ matprops.spatialdb \
+ lgdeformrigidbody.cfg \
+ lgdeformtraction.cfg
+
+noinst_TMP = \
+ rigidbody_disp.spatialdb \
+ lgdeformrigidbody_t0000000.vtk \
+ lgdeformrigidbody-elastic_info.vtk \
+ lgdeformrigidbody-elastic_t0000000.vtk \
+ lgdeformrigidbody-viscoelastic_info.vtk \
+ lgdeformrigidbody-viscoelastic_t0000000.vtk \
+ lgdeformtraction_t0000000.vtk \
+ lgdeformtraction-elastic_info.vtk \
+ lgdeformtraction-elastic_t0000000.vtk \
+ lgdeformtraction-viscoelastic_info.vtk \
+ lgdeformtraction-viscoelastic_t0000000.vtk
+
+
+TESTS_ENVIRONMENT = $(PYTHON)
+
+
+# 'export' the input files by performing a mock install
+export_datadir = $(top_builddir)/tests/3dnew/hex8
+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
+
+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
+
+
+BUILT_SOURCES = export-data
+clean-local: clean-data
+CLEANFILES = *.pyc
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestHex8.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestHex8.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestHex8.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,108 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3dnew/hex8/TestHex8.py
+##
+## @brief Generic tests for problems using 3-D mesh.
+
+import unittest
+import numpy
+
+class TestHex8(unittest.TestCase):
+ """
+ Generic tests for problems using 3-D mesh.
+ """
+
+ def setUp(self):
+ """
+ Setup for tests.
+ """
+ self.mesh = {'ncells_elastic': 128,
+ 'ncells_viscoelastic': 256,
+ 'ncorners': 8,
+ 'nvertices': 567,
+ 'spaceDim': 3,
+ 'tensorSize': 6}
+ return
+
+
+ def test_elastic_info(self):
+ """
+ Check elastic info.
+ """
+ if self.reader is None:
+ return
+
+ from pylith.tests.PhysicalProperties import check_properties
+ from rigidbody_soln import p_mu,p_lambda,p_density
+
+ self.mesh['ncells'] = self.mesh['ncells_elastic']
+ ncells= self.mesh['ncells']
+ filename = "%s-elastic_info.vtk" % self.outputRoot
+ 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}
+ check_properties(self, filename, self.mesh, properties)
+
+ self.mesh['ncells'] = self.mesh['ncells_viscoelastic']
+ ncells= self.mesh['ncells']
+ filename = "%s-viscoelastic_info.vtk" % self.outputRoot
+ 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}
+ 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
+
+ from pylith.tests.StateVariables import check_state_variables
+ stateVars = ["total_strain", "stress"]
+
+ filename = "%s-elastic_t0000000.vtk" % self.outputRoot
+ self.mesh['ncells'] = self.mesh['ncells_elastic']
+ check_state_variables(self, filename, self.mesh, stateVars)
+
+ filename = "%s-viscoelastic_t0000000.vtk" % self.outputRoot
+ self.mesh['ncells'] = self.mesh['ncells_viscoelastic']
+ check_state_variables(self, filename, self.mesh, stateVars)
+
+ return
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformRigidBody.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformRigidBody.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformRigidBody.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,93 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3dnew/hex8/TestRigidBody.py
+##
+## @brief Test suite for testing pylith with rigid body motion for
+## large deformations in 3-D.
+
+import numpy
+from TestHex8 import TestHex8
+from rigidbody_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 RigidBodyApp(PyLithApp):
+ def __init__(self):
+ PyLithApp.__init__(self, name="lgdeformrigidbody")
+ return
+
+
+# Helper function to run PyLith
+def run_pylith():
+ """
+ Run pylith.
+ """
+ if not "done" in dir(run_pylith):
+ # Generate spatial databases
+ from rigidbody_gendb import GenerateDB
+ db = GenerateDB()
+ db.run()
+
+ # Run PyLith
+ app = RigidBodyApp()
+ app.run()
+ run_pylith.done = True
+ return
+
+
+class TestRigidBody(TestHex8):
+ """
+ Test suite for testing pylith with 2-D rigid body motion.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestHex8.setUp(self)
+ run_pylith()
+ self.outputRoot = "lgdeformrigidbody"
+ 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/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformTraction.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformTraction.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/TestLgDeformTraction.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3dnew/hex8/TestLgDeformTraction.py
+##
+## @brief Test suite for testing pylith with uniform tractions for
+## large deformations in 3-D.
+
+import numpy
+from TestHex8 import TestHex8
+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(TestHex8):
+ """
+ Test suite for testing pylith with 3-D axial tractions.
+ """
+
+ def setUp(self):
+ """
+ Setup for test.
+ """
+ TestHex8.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
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformrigidbody.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformrigidbody.cfg (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformrigidbody.cfg 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,119 @@
+# -*- Python -*-
+[lgdeformrigidbody]
+
+[lgdeformrigidbody.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solverlinear = 1
+#meshiocubit = 1
+#implicitelasticity = 1
+#quadrature3d = 1
+#fiatlagrange = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+[lgdeformrigidbody.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 3
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.timedependent]
+dimension = 3
+bc = [x_neg,x_pos]
+
+formulation = pylith.problems.ImplicitLgDeform
+formulation.solver = pylith.problems.SolverNonlinear
+
+[lgdeformrigidbody.timedependent.formulation.time_step]
+total_time = 0.0*s
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.timedependent]
+materials = [elastic,viscoelastic]
+materials.elastic = pylith.materials.ElasticIsotropic3D
+materials.viscoelastic = pylith.materials.ElasticIsotropic3D
+
+[lgdeformrigidbody.timedependent.materials.elastic]
+label = Elastic material
+id = 1
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 3
+
+[lgdeformrigidbody.timedependent.materials.viscoelastic]
+label = Elastic material 2
+id = 2
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 3
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.timedependent.bc.x_neg]
+bc_dof = [0,1,2]
+label = 21
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC -x face
+db_initial.iohandler.filename = rigidbody_disp.spatialdb
+
+[lgdeformrigidbody.timedependent.bc.x_pos]
+bc_dof = [0,1,2]
+label = 20
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Dirichlet BC +x face
+db_initial.iohandler.filename = rigidbody_disp.spatialdb
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.petsc]
+ksp_rtol = 1.0e-12
+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 = 200
+ksp_gmres_restart = 100
+snes_max_it = 100
+#ksp_monitor = true
+#ksp_view = true
+#snes_monitor = true
+#snes_view = true
+#ksp_converged_reason = true
+#snes_converged_reason = true
+#log_summary = true
+#start_in_debugger = true
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[lgdeformrigidbody.problem.formulation.output.output.writer]
+filename = lgdeformrigidbody.vtk
+
+[lgdeformrigidbody.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = lgdeformrigidbody-elastic.vtk
+
+[lgdeformrigidbody.timedependent.materials.viscoelastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = lgdeformrigidbody-viscoelastic.vtk
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction.cfg (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction.cfg 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,138 @@
+# -*- Python -*-
+[lgdeformtraction]
+
+[lgdeformtraction.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpirun -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+[lgdeformtraction.journal.info]
+#timedependent = 1
+#implicit = 1
+#petsc = 1
+#solverlinear = 1
+#meshiocubit = 1
+#implicitelasticity = 1
+#quadrature3d = 1
+#fiatlagrange = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[lgdeformtraction.mesh_generator]
+#debug = 1
+reader = pylith.meshio.MeshIOCubit
+
+[lgdeformtraction.mesh_generator.reader]
+filename = mesh.exo
+coordsys.space_dim = 3
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[lgdeformtraction.timedependent]
+dimension = 3
+bc = [x_neg,x_pos,y_neg,y_pos,z_neg]
+
+formulation = pylith.problems.ImplicitLgDeform
+formulation.solver = pylith.problems.SolverNonlinear
+
+[lgdeformtraction.timedependent.formulation.time_step]
+total_time = 0.0*s
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[lgdeformtraction.timedependent]
+materials = [elastic,viscoelastic]
+materials.elastic = pylith.materials.ElasticIsotropic3D
+materials.viscoelastic = pylith.materials.ElasticIsotropic3D
+
+[lgdeformtraction.timedependent.materials.elastic]
+label = Elastic material
+id = 1
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 3
+
+[lgdeformtraction.timedependent.materials.viscoelastic]
+label = Elastic material 2
+id = 2
+db_properties.iohandler.filename = matprops.spatialdb
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 3
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+[lgdeformtraction.timedependent.bc]
+x_pos = pylith.bc.Neumann
+y_pos = pylith.bc.Neumann
+
+[lgdeformtraction.timedependent.bc.x_neg]
+bc_dof = [0]
+label = 21
+
+[lgdeformtraction.timedependent.bc.y_neg]
+bc_dof = [1]
+label = 23
+
+[lgdeformtraction.timedependent.bc.z_neg]
+bc_dof = [2]
+label = 25
+
+[lgdeformtraction.timedependent.bc.x_pos]
+label = 20
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC +x face
+db_initial.values = [traction-shear-horiz,traction-shear-vert,traction-normal]
+db_initial.data = [0.0*MPa,0.0*MPa,-100*MPa]
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+[lgdeformtraction.timedependent.bc.y_pos]
+label = 22
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Neumann BC +y face
+db_initial.values = [traction-shear-horiz,traction-shear-vert,traction-normal]
+db_initial.data = [0.0*MPa,0.0*MPa,-10*MPa]
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[lgdeformtraction.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 = 200
+ksp_gmres_restart = 50
+#ksp_monitor = true
+#ksp_view = true
+#snes_monitor = true
+#snes_view = true
+#ksp_converged_reason = true
+#snes_converged_reason = true
+#log_summary = true
+#start_in_debugger = true
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[lgdeformtraction.problem.formulation.output.output.writer]
+filename = lgdeformtraction.vtk
+
+[lgdeformtraction.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = lgdeformtraction-elastic.vtk
+
+[lgdeformtraction.timedependent.materials.viscoelastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = lgdeformtraction-viscoelastic.vtk
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction_soln.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction_soln.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/lgdeformtraction_soln.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,102 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3dnew/hex8/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
+sxx0 = -100.0e+6
+syy0 = -10.0e+6
+ux = -13.7398
+uy = 2.30274
+uz = 3.05947
+
+sxx = -1.00172e+8
+syy = -9.99712e+6
+szz = 0.0
+sxy = 0.0
+syz = 0.0
+sxz = 0.0
+
+# Uniform strain field
+exx = -0.001716
+eyy = 0.000287884
+ezz = 0.000510042
+exy = 0.0
+eyz = 0.0
+exz = 0.0
+
+# ----------------------------------------------------------------------
+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
+ disp[:,2] = uz*(6000.0 + locs[:,2]) / 6000.0
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (npts, dim) = locs.shape
+ strain = numpy.zeros( (npts, 6), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = ezz
+ strain[:,3] = exy
+ strain[:,4] = eyz
+ strain[:,5] = exz
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (npts, dim) = locs.shape
+ stress = numpy.zeros( (npts, 6), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = szz
+ stress[:,3] = sxy
+ stress[:,4] = syz
+ stress[:,5] = sxz
+ return stress
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/matprops.spatialdb
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/matprops.spatialdb (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/matprops.spatialdb 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,14 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 3
+ value-names = density vs vp
+ value-units = kg/m**3 m/s m/s
+ num-locs = 1
+ data-dim = 0
+ space-dim = 3
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 3
+ }
+}
+0.0 0.0 0.0 2500.0 3000.0 5291.502622129181
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_gendb.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_gendb.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_gendb.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3dnew/hex8/rigidbody_gendb.py
+##
+## @brief Python script to generate spatial database with displacement
+## boundary conditions for the rigid body motion test.
+
+import numpy
+
+class GenerateDB(object):
+ """
+ Python object to generate spatial database with displacement
+ boundary conditions for the rigid body motion 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)
+ z = numpy.arange(-6000.0, 0000.1, 1000.0)
+ nptsX = x.shape[0]
+ nptsY = y.shape[0]
+ nptsZ = z.shape[0]
+
+ xx = x * numpy.ones( (nptsY*nptsZ, 1), dtype=numpy.float64)
+ yy0 = y * numpy.ones( (nptsZ, 1), dtype=numpy.float64)
+ yy = yy0.ravel() * numpy.ones( (nptsX, 1), dtype=numpy.float64)
+ zz = z * numpy.ones( (nptsX*nptsY, 1), dtype=numpy.float64)
+ xyz = numpy.zeros( (nptsX*nptsY*nptsZ, 3), dtype=numpy.float64)
+ xyz[:,0] = numpy.ravel(xx)
+ xyz[:,1] = numpy.ravel(numpy.transpose(yy))
+ xyz[:,2] = numpy.ravel(numpy.transpose(zz))
+
+ from rigidbody_soln import AnalyticalSoln
+ soln = AnalyticalSoln()
+ disp = soln.displacement(xyz)
+
+ from spatialdata.geocoords.CSCart import CSCart
+ cs = CSCart()
+ cs.inventory.spaceDim = 3
+ cs._configure()
+ data = {'locs': xyz,
+ '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])},
+ {'name': "displacement-z",
+ 'units': "m",
+ 'data': numpy.ravel(disp[:,2])}]}
+
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ io = SimpleIOAscii()
+ io.inventory.filename = "rigidbody_disp.spatialdb"
+ io._configure()
+ io.write(data)
+ return
+
+# ======================================================================
+if __name__ == "__main__":
+ app = GenerateDB()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_soln.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_soln.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/rigidbody_soln.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tests/3dnew/hex8/rigidbody_soln.py
+##
+## @brief Analytical solution to rigid body motion 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
+syy = 0.0
+szz = 0.0
+sxy = 0.0
+syz = 0.0
+sxz = 0.0
+
+# Uniform strain field
+exx = 0.0
+eyy = 0.0
+ezz = 0.0
+exy = 0.0
+eyz = 0.0
+exz = 0.0
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+ """
+ Analytical solution to axial/shear displacement problem.
+ """
+
+ def __init__(self):
+ return
+
+
+ def displacement(self, locs):
+ """
+ Compute displacement field at locations.
+ """
+ u0 = 400.0
+ v0 = -300.0
+ w0 = 500.0
+ from math import pi
+ theta = 60.0/180.0*pi
+
+ (npts, dim) = locs.shape
+ disp = numpy.zeros( (npts, 3), dtype=numpy.float64)
+ disp[:,0] = u0 -locs[:,0] + \
+ u0 + numpy.cos(theta)*locs[:,0] + numpy.sin(theta)*locs[:,2]
+ disp[:,1] = v0
+ disp[:,2] = -locs[:,2] + \
+ w0 - numpy.sin(theta)*locs[:,0] + numpy.cos(theta)*locs[:,2]
+ return disp
+
+
+ def strain(self, locs):
+ """
+ Compute strain field at locations.
+ """
+ (npts, dim) = locs.shape
+ strain = numpy.zeros( (npts, 6), dtype=numpy.float64)
+ strain[:,0] = exx
+ strain[:,1] = eyy
+ strain[:,2] = ezz
+ strain[:,3] = exy
+ strain[:,4] = eyz
+ strain[:,5] = exz
+ return strain
+
+
+ def stress(self, locs):
+ """
+ Compute stress field at locations.
+ """
+ (npts, dim) = locs.shape
+ stress = numpy.zeros( (npts, 6), dtype=numpy.float64)
+ stress[:,0] = sxx
+ stress[:,1] = syy
+ stress[:,2] = szz
+ stress[:,3] = sxy
+ stress[:,4] = syz
+ stress[:,5] = sxz
+ return stress
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/testpylith.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/testpylith.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/testpylith.py 2009-12-07 20:52:00 UTC (rev 16082)
@@ -0,0 +1,45 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+__requires__ = "PyLith"
+
+import unittest
+
+def suite():
+ """
+ Create test suite.
+ """
+ suite = unittest.TestSuite()
+
+ from TestLgDeformRigidBody import TestRigidBody
+ suite.addTest(unittest.makeSuite(TestRigidBody))
+
+ from TestLgDeformTraction import TestTraction
+ suite.addTest(unittest.makeSuite(TestTraction))
+
+ return suite
+
+
+def main():
+ """
+ Run test suite.
+ """
+ unittest.TextTestRunner(verbosity=2).run(suite())
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == '__main__':
+ main()
+
+
+# End of file
Property changes on: short/3D/PyLith/branches/pylith-friction/tests/3dnew/hex8/testpylith.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: short/3D/PyLith/branches/pylith-friction/tests/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/Makefile.am 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/tests/Makefile.am 2009-12-07 20:52:00 UTC (rev 16082)
@@ -13,6 +13,7 @@
SUBDIRS = \
1d \
2d \
+ 3dnew \
petsc
Modified: short/3D/PyLith/branches/pylith-friction/tests/README
===================================================================
--- short/3D/PyLith/branches/pylith-friction/tests/README 2009-12-07 17:16:25 UTC (rev 16081)
+++ short/3D/PyLith/branches/pylith-friction/tests/README 2009-12-07 20:52:00 UTC (rev 16082)
@@ -72,13 +72,13 @@
* dislocationdyn --TODO--
[single fault, dynamic]
- * lgdeformrigidbody --TODO--
+ * lgdeformrigidbody --TODO-- (need higher precision output)
- Dirichlet BC w/initial displacement, UniformDB
(rigid body translation and rotation)
- Large deformation
- * lgdeformtraction --TODO--
+ * lgdeformtraction --TODO-- (need higher precision output)
- Dirichlet BC w/initial displacement, UniformDB
- Neuman BC w/initial value, UniformDB
@@ -115,3 +115,15 @@
Gravity with faults, Neumann BC.
+ * lgdeformrigidbody --TODO-- (need higher precision output)
+
+ - Dirichlet BC w/initial displacement, UniformDB
+ (rigid body translation and rotation)
+ - Large deformation
+
+ * lgdeformtraction --TODO-- (need higher precision output)
+
+ - Dirichlet BC w/initial displacement, UniformDB
+ - Neuman BC w/initial value, UniformDB
+ - Large deformation
+
More information about the CIG-COMMITS
mailing list