[cig-commits] r20119 - in short/3D/PyLith/branches/v1.7-trunk/unittests: libtests/materials libtests/materials/data pytests/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon May 14 21:50:33 PDT 2012
Author: willic3
Date: 2012-05-14 21:50:33 -0700 (Mon, 14 May 2012)
New Revision: 20119
Added:
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.hh
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.hh
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.hh
short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/TestPowerLawPlaneStrain.py
Modified:
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/generate.sh
short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/testmaterials.py
Log:
Added unit tests for PowerLawPlaneStrain, but they don't pass yet.
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/Makefile.am 2012-05-15 04:49:18 UTC (rev 20118)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/Makefile.am 2012-05-15 04:50:33 UTC (rev 20119)
@@ -41,6 +41,7 @@
TestGenMaxwellPlaneStrain.cc \
TestGenMaxwellQpQsIsotropic3D.cc \
TestPowerLaw3D.cc \
+ TestPowerLawPlaneStrain.cc \
TestDruckerPrager3D.cc \
TestDruckerPragerPlaneStrain.cc \
TestEffectiveStress.cc \
@@ -62,6 +63,7 @@
TestMaxwellIsotropic3D.hh \
TestMaxwellPlaneStrain.hh \
TestPowerLaw3D.hh \
+ TestPowerLawPlaneStrain.hh \
TestDruckerPrager3D.hh \
TestDruckerPragerPlaneStrain.hh \
TestEffectiveStress.hh
@@ -87,6 +89,8 @@
data/GenMaxwellQpQsIsotropic3DTimeDepData.cc \
data/PowerLaw3DElasticData.cc \
data/PowerLaw3DTimeDepData.cc \
+ data/PowerLawPlaneStrainElasticData.cc \
+ data/PowerLawPlaneStrainTimeDepData.cc \
data/DruckerPrager3DElasticData.cc \
data/DruckerPrager3DTimeDepData.cc \
data/DruckerPragerPlaneStrainElasticData.cc \
@@ -113,6 +117,8 @@
data/MaxwellPlaneStrainTimeDepData.hh \
data/PowerLaw3DElasticData.hh \
data/PowerLaw3DTimeDepData.hh \
+ data/PowerLawPlaneStrainElasticData.hh \
+ data/PowerLawPlaneStrainTimeDepData.hh \
data/DruckerPrager3DElasticData.hh \
data/DruckerPrager3DTimeDepData.hh \
data/DruckerPragerPlaneStrainElasticData.hh \
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.cc (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.cc 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,224 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestPowerLawPlaneStrain.hh" // Implementation of class methods
+
+#include "data/PowerLawPlaneStrainElasticData.hh" // USES PowerLawPlaneStrainElasticData
+#include "data/PowerLawPlaneStrainTimeDepData.hh" // USES PowerLawPlaneStrainTimeDepData
+
+#include "pylith/materials/PowerLawPlaneStrain.hh" // USES PowerLawPlaneStrain
+
+#include <cstring> // USES memcpy()
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::materials::TestPowerLawPlaneStrain );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::materials::TestPowerLawPlaneStrain::setUp(void)
+{ // setUp
+ _material = new PowerLawPlaneStrain();
+ _matElastic = new PowerLawPlaneStrain();
+ _data = new PowerLawPlaneStrainElasticData();
+ _dataElastic = new PowerLawPlaneStrainElasticData();
+ setupNormalizer();
+} // setUp
+
+// ----------------------------------------------------------------------
+// Test timeStep()
+void
+pylith::materials::TestPowerLawPlaneStrain::testTimeStep(void)
+{ // testTimeStep
+ PowerLawPlaneStrain material;
+
+ CPPUNIT_ASSERT_EQUAL(false, material._needNewJacobian);
+
+ const PylithScalar dt1 = 1.0;
+ material.timeStep(dt1);
+ CPPUNIT_ASSERT_EQUAL(dt1, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(true, material.needNewJacobian());
+
+ const PylithScalar dt2 = 2.0;
+ material.timeStep(dt2);
+ CPPUNIT_ASSERT_EQUAL(dt2, material.Material::timeStep());
+ CPPUNIT_ASSERT_EQUAL(true, material.needNewJacobian());
+} // testTimeStep
+
+// ----------------------------------------------------------------------
+// Test useLinearBehavior()
+void
+pylith::materials::TestPowerLawPlaneStrain::testUseLinearBehavior(void)
+{ // testUseLinearBehavior
+ PowerLawPlaneStrain material;
+
+ material.useLinearBehavior(true);
+ // Some compilers/operating systems (cygwin) don't allow comparing
+ // pointers. Use first test to determine if we can compare pointers.
+ if (&pylith::materials::PowerLawPlaneStrain::_calcStressElastic ==
+ material._calcStressFn) {
+ CPPUNIT_ASSERT(&pylith::materials::PowerLawPlaneStrain::_calcStressElastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(&pylith::materials::PowerLawPlaneStrain::_calcElasticConstsElastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(&pylith::materials::PowerLawPlaneStrain::_updateStateVarsElastic ==
+ material._updateStateVarsFn);
+
+ material.useLinearBehavior(false);
+ CPPUNIT_ASSERT(&pylith::materials::PowerLawPlaneStrain::_calcStressViscoelastic ==
+ material._calcStressFn);
+ CPPUNIT_ASSERT(&pylith::materials::PowerLawPlaneStrain::_calcElasticConstsViscoelastic ==
+ material._calcElasticConstsFn);
+ CPPUNIT_ASSERT(&pylith::materials::PowerLawPlaneStrain::_updateStateVarsViscoelastic ==
+ material._updateStateVarsFn);
+ } // if
+} // testUseLinearBehavior
+
+// ----------------------------------------------------------------------
+// Test usesHasStateVars()
+void
+pylith::materials::TestPowerLawPlaneStrain::testHasStateVars(void)
+{ // testHasStateVars
+ PowerLawPlaneStrain material;
+ CPPUNIT_ASSERT_EQUAL(true, material.hasStateVars());
+} // testHasStateVars
+
+// ----------------------------------------------------------------------
+// Test _calcStressElastic()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_calcStressElastic(void)
+{ // test_calcStressElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useLinearBehavior(true);
+
+ test_calcStress();
+} // test_calcStressElastic
+
+// ----------------------------------------------------------------------
+// Test calcElasticConstsElastic()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_calcElasticConstsElastic(void)
+{ // test_calcElasticConstsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useLinearBehavior(true);
+
+ test_calcElasticConsts();
+} // test_calcElasticConstsElastic
+
+// ----------------------------------------------------------------------
+// Test _updateStateVarsElastic()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_updateStateVarsElastic(void)
+{ // test_updateStateVarsElastic
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useLinearBehavior(true);
+
+ test_updateStateVars();
+} // test_updateStateVarsElastic
+
+// ----------------------------------------------------------------------
+// Test _calcStressTimeDep()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_calcStressTimeDep(void)
+{ // test_calcStressTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useLinearBehavior(false);
+
+ delete _dataElastic; _dataElastic = new PowerLawPlaneStrainTimeDepData();
+
+ PylithScalar dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcStress();
+} // test_calcStressTimeDep
+
+// ----------------------------------------------------------------------
+// Test _calcElasticConstsTimeDep()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_calcElasticConstsTimeDep(void)
+{ // test_calcElasticConstsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useLinearBehavior(false);
+
+ delete _dataElastic; _dataElastic = new PowerLawPlaneStrainTimeDepData();
+
+ PylithScalar dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_calcElasticConsts();
+} // test_calcElasticConstsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _updateStateVarsTimeDep()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_updateStateVarsTimeDep(void)
+{ // test_updateStateVarsTimeDep
+ CPPUNIT_ASSERT(0 != _matElastic);
+ _matElastic->useLinearBehavior(false);
+
+ delete _dataElastic; _dataElastic = new PowerLawPlaneStrainTimeDepData();
+
+ PylithScalar dt = 2.0e+5;
+ _matElastic->timeStep(dt);
+ test_updateStateVars();
+
+} // test_updateStateVarsTimeDep
+
+// ----------------------------------------------------------------------
+// Test _stableTimeStepImplicit()
+void
+pylith::materials::TestPowerLawPlaneStrain::test_stableTimeStepImplicit(void)
+{ // test_stableTimeStepImplicit
+ CPPUNIT_ASSERT(0 != _matElastic);
+
+ delete _dataElastic; _dataElastic = new PowerLawPlaneStrainTimeDepData();
+
+ TestElasticMaterial::test_stableTimeStepImplicit();
+} // test_stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
+// Test hasProperty().
+void
+pylith::materials::TestPowerLawPlaneStrain::testHasProperty(void)
+{ // testHasProperty
+ PowerLawPlaneStrain material;
+
+ CPPUNIT_ASSERT(material.hasProperty("mu"));
+ CPPUNIT_ASSERT(material.hasProperty("lambda"));
+ CPPUNIT_ASSERT(material.hasProperty("density"));
+ CPPUNIT_ASSERT(material.hasProperty("reference_strain_rate"));
+ CPPUNIT_ASSERT(material.hasProperty("reference_stress"));
+ CPPUNIT_ASSERT(material.hasProperty("power_law_exponent"));
+ CPPUNIT_ASSERT(!material.hasProperty("aaa"));
+} // testHasProperty
+
+// ----------------------------------------------------------------------
+// Test hasStateVar().
+void
+pylith::materials::TestPowerLawPlaneStrain::testHasStateVar(void)
+{ // testHasStateVar
+ PowerLawPlaneStrain material;
+
+ CPPUNIT_ASSERT(material.hasStateVar("stress_zz_initial"));
+ CPPUNIT_ASSERT(material.hasStateVar("viscous_strain"));
+ CPPUNIT_ASSERT(material.hasStateVar("stress4"));
+ CPPUNIT_ASSERT(!material.hasStateVar("total_strain"));
+} // testHasStateVar
+
+
+// End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.hh (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/TestPowerLawPlaneStrain.hh 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,123 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/materials/TestPowerLawPlaneStrain.hh
+ *
+ * @brief C++ TestPowerLawPlaneStrain object
+ *
+ * C++ unit testing for PowerLawPlaneStrain.
+ */
+
+#if !defined(pylith_materials_testpowerlawplanestrain_hh)
+#define pylith_materials_testpowerlawplanestrain_hh
+
+#include "TestElasticMaterial.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class TestPowerLawPlaneStrain;
+ class PowerLawPlaneStrainElasticData;
+ class PowerLawPlaneStrainTimeDepData;
+ } // materials
+} // pylith
+
+/// C++ unit testing for PowerLawPlaneStrain
+class pylith::materials::TestPowerLawPlaneStrain : public TestElasticMaterial
+{ // class TestPowerLawPlaneStrain
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestPowerLawPlaneStrain );
+
+ CPPUNIT_TEST( testDimension );
+ CPPUNIT_TEST( testTensorSize );
+ CPPUNIT_TEST( testDBToProperties );
+ CPPUNIT_TEST( testNonDimProperties );
+ CPPUNIT_TEST( testDimProperties );
+ CPPUNIT_TEST( testDBToStateVars );
+ CPPUNIT_TEST( testNonDimStateVars );
+ CPPUNIT_TEST( testDimStateVars );
+ CPPUNIT_TEST( test_calcDensity );
+ CPPUNIT_TEST( test_stableTimeStepImplicit );
+
+ // Need to test Power Law viscoelastic specific behavior.
+ CPPUNIT_TEST( testTimeStep );
+ CPPUNIT_TEST( testUseLinearBehavior );
+ CPPUNIT_TEST( testHasStateVars );
+
+ CPPUNIT_TEST( test_calcStressElastic );
+ CPPUNIT_TEST( test_calcStressTimeDep );
+ CPPUNIT_TEST( test_calcElasticConstsElastic );
+ CPPUNIT_TEST( test_calcElasticConstsTimeDep );
+ CPPUNIT_TEST( test_updateStateVarsElastic );
+ CPPUNIT_TEST( test_updateStateVarsTimeDep );
+
+ CPPUNIT_TEST( testHasProperty );
+ CPPUNIT_TEST( testHasStateVar );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Setup testing data.
+ void setUp(void);
+
+ /// Test timeStep()
+ void testTimeStep(void);
+
+ /// Test useLinearBehavior()
+ void testUseLinearBehavior(void);
+
+ /// Test hasStateVars()
+ void testHasStateVars(void);
+
+ /// Test _calcStressElastic()
+ void test_calcStressElastic(void);
+
+ /// Test _calcElasticConstsElastic()
+ void test_calcElasticConstsElastic(void);
+
+ /// Test _updateStateVarsElastic()
+ void test_updateStateVarsElastic(void);
+
+ /// Test _calcStressTimeDep()
+ void test_calcStressTimeDep(void);
+
+ /// Test _calcElasticConstsTimeDep()
+ void test_calcElasticConstsTimeDep(void);
+
+ /// Test _updateStatevarsTimeDep()
+ void test_updateStateVarsTimeDep(void);
+
+ /// Test _stableTimeStepImplicit()
+ void test_stableTimeStepImplicit(void);
+
+ /// Test hasProperty()
+ void testHasProperty(void);
+
+ /// Test hasStateVar()
+ void testHasStateVar(void);
+
+}; // class TestPowerLawPlaneStrain
+
+#endif // pylith_materials_testpowerlawplanestrain_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,288 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ PowerLawPlaneStrain object with elastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+
+# ----------------------------------------------------------------------
+dimension = 2
+numElasticConsts = 9
+tensorSize = 3
+
+# PowerLawPlaneStrainElastic class
+class PowerLawPlaneStrainElastic(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ PowerLawPlaneStrain object with elastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="powerlawplanestrainelastic"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ # import pdb
+ # pdb.set_trace()
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "reference-strain-rate", "reference-stress",
+ "power-law-exponent"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+
+ self.dbStateVarValues = ["stress-zz-initial",
+ "viscous-strain-xx",
+ "viscous-strain-yy",
+ "viscous-strain-zz",
+ "viscous-strain-xy",
+ "stress4-xx",
+ "stress4-yy",
+ "stress4-zz",
+ "stress4-xy"
+ ]
+ self.numStateVarValues = numpy.array([1, 4, 4], dtype=numpy.int32)
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ # Derive new values in based on previous value for power-law coefficient
+ # and viscosity coefficient.
+ powerLawCoeffA = 1.0/3.0e18
+ refStrainRateA = 1.0e-6
+ powerLawExponentA = 1.0
+ strainA = [1.1e-4, 1.2e-4, 1.4e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.4e4]
+ initialStrainA = [3.1e-4, 3.2e-4, 3.4e-4]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ stressZZInitialA = numpy.array([1.0754e4], dtype=numpy.float64)
+
+ viscosityCoeffA = (1.0/((3.0**0.5)**(powerLawExponentA + 1.0) \
+ * powerLawCoeffA))**(1.0/powerLawExponentA)
+ refStressA = viscosityCoeffA * \
+ (2.0 * refStrainRateA) ** (1.0/powerLawExponentA)
+ # refStressA = (refStrainRateA/powerLawCoeffA)**(1.0/powerLawExponentA)
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ powerLawCoeffB = 1.0/9.0e30
+ refStrainRateB = 1.0e-6
+ powerLawExponentB = 3.0
+ strainB = [4.1e-4, 4.2e-4, 4.4e-4]
+ initialStressB = [5.1e4, 5.2e4, 5.4e4]
+ initialStrainB = [6.1e-4, 6.2e-4, 6.4e-4]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ viscosityCoeffB = (1.0/((3.0**0.5)**(powerLawExponentB + 1.0) \
+ * powerLawCoeffB))**(1.0/powerLawExponentB)
+ refStressB = viscosityCoeffB * \
+ (2.0 * refStrainRateB) ** (1.0/powerLawExponentB)
+ # refStressB = (refStrainRateB/powerLawCoeffB)**(1.0/powerLawExponentB)
+ stressZZInitialB = numpy.array([2.575e4], dtype=numpy.float64)
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.timeScale = 1.0
+ self.densityScale = 1.0e+3
+ self.strainRateScale = 1.0/self.timeScale
+
+ self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
+ refStrainRateA, refStressA, \
+ powerLawExponentA],
+ [densityB, vsB, vpB, \
+ refStrainRateB, refStressB, \
+ powerLawExponentB] ],
+ dtype=numpy.float64)
+ self.properties = numpy.array([ [densityA, muA, lambdaA, \
+ refStrainRateA, refStressA, \
+ powerLawExponentA],
+ [densityB, muB, lambdaB, \
+ refStrainRateB, refStressB, \
+ powerLawExponentB] ],
+ dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ self.dbStateVars = numpy.zeros( (numLocs, 1 + 4 + 4),
+ dtype=numpy.float64)
+ self.dbStateVars[0,0] = stressZZInitialA
+ self.dbStateVars[1,0] = stressZZInitialB
+
+ self.stateVars = numpy.zeros( (numLocs, 1 + 4 + 4),
+ dtype=numpy.float64)
+ self.stateVars[0,0] = stressZZInitialA
+ self.stateVars[1,0] = stressZZInitialB
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ time0 = self.timeScale
+ strainRate0 = self.strainRateScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0, \
+ refStrainRateA/strainRate0, refStressA/mu0, \
+ powerLawExponentA],
+ [densityB/density0, muB/mu0, lambdaB/mu0, \
+ refStrainRateB/strainRate0, refStressB/mu0, \
+ powerLawExponentB] ],
+ dtype=numpy.float64)
+
+ stressZZInitialANondim = stressZZInitialA/mu0
+ stressZZInitialBNondim = stressZZInitialB/mu0
+
+ self.stateVarsNondim = self.stateVars # no scaling
+
+ self.stateVarsNondim[0, 0] = stressZZInitialANondim
+ self.stateVarsNondim[1, 0] = stressZZInitialBNondim
+
+ self.initialStress = numpy.array([initialStressA,
+ initialStressB],
+ dtype=numpy.float64)
+ self.initialStrain = numpy.array([initialStrainA,
+ initialStrainB],
+ dtype=numpy.float64)
+
+ self.density = numpy.array([densityA,
+ densityB],
+ dtype=numpy.float64)
+
+ self.strain = numpy.array([strainA,
+ strainB],
+ dtype=numpy.float64)
+
+ self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts), \
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA, \
+ initialStressA, initialStrainA)
+ (self.elasticConsts[1,:], self.stress[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB, \
+ initialStressB, initialStrainB)
+
+ viscousStrainUpdated = numpy.zeros((numLocs, 4), dtype=numpy.float64)
+
+ stressUpdated = numpy.zeros((numLocs, 4), dtype=numpy.float64)
+ stressUpdated[0,0] = self.stress[0,0]
+ stressUpdated[1,0] = self.stress[1,0]
+ stressUpdated[0,1] = self.stress[0,1]
+ stressUpdated[1,1] = self.stress[1,1]
+ stressUpdated[0,3] = self.stress[0,2]
+ stressUpdated[1,3] = self.stress[1,2]
+ stressUpdated[0,2] = lambdaA * (strainA[0] + strainA[1]) + stressZZInitialA
+ stressUpdated[1,2] = lambdaB * (strainB[0] + strainB[1]) + stressZZInitialB
+
+ self.stateVarsUpdated = numpy.zeros( (numLocs, 1 + 4 + 4),
+ dtype=numpy.float64)
+ self.stateVarsUpdated[0,0] = stressZZInitialA
+ self.stateVarsUpdated[1,0] = stressZZInitialB
+ self.stateVarsUpdated[0,5:9] = stressUpdated[0,:]
+ self.stateVarsUpdated[1,5:9] = stressUpdated[1,:]
+
+ maxwellTimeA = self._getMaxwellTime(muA, refStrainRateA, refStressA,
+ powerLawExponentA,
+ stressUpdated[0,:])
+ maxwellTimeB = self._getMaxwellTime(muB, refStrainRateB, refStressB,
+ powerLawExponentB,
+ stressUpdated[1,:])
+
+
+ self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
+ return
+
+
+ def _scalarProduct(self, tensor1, tensor2):
+ """
+ Compute the scalar product of two tensors stored in vector form.
+ """
+ scalarProduct = tensor1[0] * tensor2[0] + \
+ tensor1[1] * tensor2[1] + \
+ tensor1[2] * tensor2[2] + \
+ 2.0 * tensor1[3] * tensor2[3]
+ return scalarProduct
+
+
+ def _getMaxwellTime(self, mu, refStrainRate, refStress, powerLawExponent,
+ stress4):
+ """
+ Compute Maxwell time from stress4, viscosity coefficient, shear modulus, and
+ power-law exponent.
+ """
+ meanStress = (stress4[0] + stress4[1] + stress4[2])/3.0
+ devStress = numpy.array(stress4, dtype=numpy.float64)
+
+ devStress[0] = stress4[0] - meanStress
+ devStress[1] = stress4[1] - meanStress
+ devStress[2] = stress4[2] - meanStress
+
+ devStressProd = self._scalarProduct(devStress, devStress)
+ effStress = (0.5 * devStressProd)**0.5
+ maxwellTime = 1.0
+ if (effStress != 0.0):
+ maxwellTime = 0.5 * (refStress/effStress)**(powerLawExponent - 1.0) * \
+ (refStress/mu)/refStrainRate
+
+ return maxwellTime
+
+ def _calcStress(self, strainV, muV, lambdaV, initialStressV, initialStrainV):
+ """
+ Compute stress and derivative of elasticity matrix.
+ """
+ C1111 = lambdaV + 2.0*muV
+ C1122 = lambdaV
+ C1112 = 0.0
+ C2211 = lambdaV
+ C2222 = lambdaV + 2.0*muV
+ C2212 = 0.0
+ C1211 = 0.0
+ C1222 = 0.0
+ C1212 = 2.0*muV
+ elasticConsts = numpy.array([C1111, C1122, C1112,
+ C2211, C2222, C2212,
+ C1211, C1222, C1212], dtype=numpy.float64)
+
+ strain = numpy.reshape(strainV, (tensorSize,1))
+ initialStress = numpy.reshape(initialStressV, (tensorSize,1))
+ initialStrain = numpy.reshape(initialStrainV, (tensorSize,1))
+ elastic = numpy.array([ [C1111, C1122, C1112],
+ [C2211, C2222, C2212],
+ [C1211, C1222, C1212] ], dtype=numpy.float64)
+ stress = numpy.dot(elastic, strain-initialStrain) + initialStress
+ return (elasticConsts, numpy.ravel(stress))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = PowerLawPlaneStrainElastic()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,315 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application powerlawplanestrainelastic.
+
+#include "PowerLawPlaneStrainElasticData.hh"
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_dimension = 2;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numLocs = 2;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numProperties = 6;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numStateVars = 3;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numDBProperties = 6;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numDBStateVars = 9;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numPropsQuadPt = 6;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numVarsQuadPt = 9;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_lengthScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_timeScale = 1.00000000e+00;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_pressureScale = 2.25000000e+10;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_densityScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_dtStableImplicit = 4.09893495e+06;
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::PowerLawPlaneStrainElasticData::_numStateVarValues[] = {
+1,
+4,
+4,
+};
+
+const char* pylith::materials::PowerLawPlaneStrainElasticData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"reference-strain-rate",
+"reference-stress",
+"power-law-exponent",
+};
+
+const char* pylith::materials::PowerLawPlaneStrainElasticData::_dbStateVarValues[] = {
+"stress-zz-initial",
+"viscous-strain-xx",
+"viscous-strain-yy",
+"viscous-strain-zz",
+"viscous-strain-xy",
+"stress4-xx",
+"stress4-yy",
+"stress4-zz",
+"stress4-xy",
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 1.00000000e-06,
+ 2.00000000e+12,
+ 1.00000000e+00,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 1.00000000e-06,
+ 1.25992105e+08,
+ 3.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_dbStateVars[] = {
+ 1.07540000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.57500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 1.00000000e-06,
+ 2.00000000e+12,
+ 1.00000000e+00,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 1.00000000e-06,
+ 1.25992105e+08,
+ 3.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_stateVars[] = {
+ 4.77955556e-07,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 1.14444444e-06,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 1.00000000e-06,
+ 8.88888889e+01,
+ 1.00000000e+00,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 1.00000000e-06,
+ 5.59964911e-03,
+ 3.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_stateVarsNondim[] = {
+ 4.77955556e-07,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 1.14444444e-06,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_stress[] = {
+ -1.79790000e+07,
+ -1.79780000e+07,
+ -8.97600000e+06,
+ -2.25300000e+06,
+ -2.25200000e+06,
+ -1.09800000e+06,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_elasticConsts[] = {
+ 6.75000000e+10,
+ 2.25000000e+10,
+ 0.00000000e+00,
+ 2.25000000e+10,
+ 6.75000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 8.64000000e+09,
+ 2.88000000e+09,
+ 0.00000000e+00,
+ 2.88000000e+09,
+ 8.64000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.40000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.40000000e+04,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_initialStrain[] = {
+ 3.10000000e-04,
+ 3.20000000e-04,
+ 3.40000000e-04,
+ 6.10000000e-04,
+ 6.20000000e-04,
+ 6.40000000e-04,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_stateVarsUpdated[] = {
+ 1.07540000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ -1.79790000e+07,
+ -1.79780000e+07,
+ 5.18575400e+06,
+ -8.97600000e+06,
+ 2.57500000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ -2.25300000e+06,
+ -2.25200000e+06,
+ 2.41615000e+06,
+ -1.09800000e+06,
+};
+
+pylith::materials::PowerLawPlaneStrainElasticData::PowerLawPlaneStrainElasticData(void)
+{ // constructor
+ dimension = _dimension;
+ numLocs = _numLocs;
+ numProperties = _numProperties;
+ numStateVars = _numStateVars;
+ numDBProperties = _numDBProperties;
+ numDBStateVars = _numDBStateVars;
+ numPropsQuadPt = _numPropsQuadPt;
+ numVarsQuadPt = _numVarsQuadPt;
+ lengthScale = _lengthScale;
+ timeScale = _timeScale;
+ pressureScale = _pressureScale;
+ densityScale = _densityScale;
+ dtStableImplicit = _dtStableImplicit;
+ numPropertyValues = const_cast<int*>(_numPropertyValues);
+ numStateVarValues = const_cast<int*>(_numStateVarValues);
+ dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+ dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+ dbProperties = const_cast<PylithScalar*>(_dbProperties);
+ dbStateVars = const_cast<PylithScalar*>(_dbStateVars);
+ properties = const_cast<PylithScalar*>(_properties);
+ stateVars = const_cast<PylithScalar*>(_stateVars);
+ propertiesNondim = const_cast<PylithScalar*>(_propertiesNondim);
+ stateVarsNondim = const_cast<PylithScalar*>(_stateVarsNondim);
+ density = const_cast<PylithScalar*>(_density);
+ strain = const_cast<PylithScalar*>(_strain);
+ stress = const_cast<PylithScalar*>(_stress);
+ elasticConsts = const_cast<PylithScalar*>(_elasticConsts);
+ initialStress = const_cast<PylithScalar*>(_initialStress);
+ initialStrain = const_cast<PylithScalar*>(_initialStrain);
+ stateVarsUpdated = const_cast<PylithScalar*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::PowerLawPlaneStrainElasticData::~PowerLawPlaneStrainElasticData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.hh (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.hh 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,110 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application powerlawplanestrainelastic.
+
+#if !defined(pylith_materials_powerlawplanestrainelasticdata_hh)
+#define pylith_materials_powerlawplanestrainelasticdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class PowerLawPlaneStrainElasticData;
+ } // pylith
+} // materials
+
+class pylith::materials::PowerLawPlaneStrainElasticData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ PowerLawPlaneStrainElasticData(void);
+
+ /// Destructor
+ ~PowerLawPlaneStrainElasticData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numLocs;
+
+ static const int _numProperties;
+
+ static const int _numStateVars;
+
+ static const int _numDBProperties;
+
+ static const int _numDBStateVars;
+
+ static const int _numPropsQuadPt;
+
+ static const int _numVarsQuadPt;
+
+ static const PylithScalar _lengthScale;
+
+ static const PylithScalar _timeScale;
+
+ static const PylithScalar _pressureScale;
+
+ static const PylithScalar _densityScale;
+
+ static const PylithScalar _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const PylithScalar _dbProperties[];
+
+ static const PylithScalar _dbStateVars[];
+
+ static const PylithScalar _properties[];
+
+ static const PylithScalar _stateVars[];
+
+ static const PylithScalar _propertiesNondim[];
+
+ static const PylithScalar _stateVarsNondim[];
+
+ static const PylithScalar _density[];
+
+ static const PylithScalar _strain[];
+
+ static const PylithScalar _stress[];
+
+ static const PylithScalar _elasticConsts[];
+
+ static const PylithScalar _initialStress[];
+
+ static const PylithScalar _initialStrain[];
+
+ static const PylithScalar _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_powerlawplanestrainelasticdata_hh
+
+// End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,519 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ PowerLawPlaneStrain object with viscoelastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+
+# ----------------------------------------------------------------------
+dimension = 2
+numElasticConsts = 9
+tensorSize = 3
+
+# PowerLawPlaneStrainTimeDep class
+class PowerLawPlaneStrainTimeDep(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ PowerLawPlaneStrain object using viscoelastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="powerlawplanestraintimedep"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ # import pdb
+ # pdb.set_trace()
+
+ numLocs = 2
+
+ self.dimension = dimension
+ self.numLocs = numLocs
+
+ self.dbPropertyValues = ["density", "vs", "vp",
+ "reference-strain-rate", "reference-stress",
+ "power-law-exponent"]
+ self.propertyValues = ["density", "mu", "lambda",
+ "reference-strain-rate", "reference-stress",
+ "power-law-exponent"]
+ self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1], dtype=numpy.int32)
+
+ self.dbStateVarValues = ["stress-zz-initial",
+ "viscous-strain-xx",
+ "viscous-strain-yy",
+ "viscous-strain-zz",
+ "viscous-strain-xy",
+ "stress4-xx",
+ "stress4-yy",
+ "stress4-zz",
+ "stress4-xy"
+ ]
+ self.numStateVarValues = numpy.array([1, 4, 4], dtype=numpy.int32)
+
+ self.alpha = 0.5
+ self.dt = 2.0e5
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ powerLawCoeffA = 1.0/3.0e18
+ refStrainRateA = 1.0e-6
+ powerLawExponentA = 1.0
+ strainA = [1.1e-4, 1.2e-4, 1.4e-4]
+ initialStressA = [2.1e4, 2.2e4, 2.4e4]
+ initialStrainA = [3.6e-5, 3.5e-5, 3.3e-5]
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ viscosityCoeffA = (1.0/((3.0**0.5)**(powerLawExponentA + 1.0) \
+ * powerLawCoeffA))**(1.0/powerLawExponentA)
+ refStressA = viscosityCoeffA * \
+ (2.0 * refStrainRateA) ** (1.0/powerLawExponentA)
+ stressZZInitialA = 2.3e4
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ powerLawCoeffB = 1.0/9.0e36
+ refStrainRateB = 1.0e-6
+ powerLawExponentB = 3.0
+ strainB = [4.1e-4, 4.2e-4, 4.4e-4]
+ initialStressB = [5.1e4, 5.2e4, 5.4e4]
+ initialStrainB = [6.1e-5, 6.2e-5, 6.6e-5]
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ viscosityCoeffB = (1.0/((3.0**0.5)**(powerLawExponentB + 1.0) \
+ * powerLawCoeffB))**(1.0/powerLawExponentB)
+ refStressB = viscosityCoeffB * \
+ (2.0 * refStrainRateB) ** (1.0/powerLawExponentB)
+ stressZZInitialB = 5.3e4
+
+ self.lengthScale = 1.0e+3
+ self.pressureScale = muA
+ self.timeScale = 1.0
+ self.densityScale = 1.0e+3
+ self.strainRateScale = 1.0/self.timeScale
+
+ self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
+ refStrainRateA, refStressA, \
+ powerLawExponentA],
+ [densityB, vsB, vpB, \
+ refStrainRateB, refStressB, \
+ powerLawExponentB] ],
+ dtype=numpy.float64)
+ self.properties = numpy.array([ [densityA, muA, lambdaA, \
+ refStrainRateA, refStressA, \
+ powerLawExponentA],
+ [densityB, muB, lambdaB, \
+ refStrainRateB, refStressB, \
+ powerLawExponentB] ],
+ dtype=numpy.float64)
+
+ # TEMPORARY, need to determine how to use initial state variables
+ self.dbStateVars = numpy.zeros( (numLocs, tensorSize),
+ dtype=numpy.float64)
+ self.dbStateVars[0, 0] = stressZZInitialA
+ self.dbStateVars[1, 0] = stressZZInitialB
+
+ mu0 = self.pressureScale
+ density0 = self.densityScale
+ time0 = self.timeScale
+ strainRate0 = self.strainRateScale
+ self.propertiesNondim = \
+ numpy.array([ [densityA/density0, muA/mu0, \
+ lambdaA/mu0, \
+ refStrainRateA/strainRate0, \
+ refStressA/mu0, \
+ powerLawExponentA], \
+ [densityB/density0, muB/mu0, \
+ lambdaB/mu0, \
+ refStrainRateB/strainRate0, \
+ refStressB/mu0, \
+ powerLawExponentB] ], \
+ dtype=numpy.float64)
+
+ self.initialStress = numpy.array([initialStressA,
+ initialStressB],
+ dtype=numpy.float64)
+ self.initialStrain = numpy.array([initialStrainA,
+ initialStrainB],
+ dtype=numpy.float64)
+
+ self.density = numpy.array([densityA,
+ densityB],
+ dtype=numpy.float64)
+
+ # Define state variables
+ visStrainA = numpy.array([4.1e-5, 4.2e-5, 4.3e-5, 4.4e-5],
+ dtype=numpy.float64)
+ visStrainB = numpy.array([1.1e-5, 1.2e-5, 1.3e-5, 1.4e-5],
+ dtype=numpy.float64)
+ stress4A = numpy.array([3.1e4, 3.2e4, stressZZInitialA, 3.4e4],
+ dtype=numpy.float64)
+ stress4B = numpy.array([5.1e4, 5.2e4, stressZZInitialB, 5.4e4],
+ dtype=numpy.float64)
+ stressNondimA = stress4A/mu0
+ stressNondimB = stress4B/mu0
+ stressZZInitialANondim = stressZZInitialA/mu0
+ stressZZInitialBNondim = stressZZInitialB/mu0
+
+ stateVarsA = numpy.concatenate(([stressZZInitialA], visStrainA, stress4A))
+ stateVarsB = numpy.concatenate(([stressZZInitialB], visStrainB, stress4B))
+ self.stateVars = numpy.array((stateVarsA, stateVarsB), dtype=numpy.float64)
+
+ stateVarsNondimA = numpy.concatenate(([stressZZInitialANondim], visStrainA,
+ stressNondimA))
+ stateVarsNondimB = numpy.concatenate(([stressZZInitialBNondim], visStrainB,
+ stressNondimB))
+ self.stateVarsNondim = numpy.array((stateVarsNondimA, stateVarsNondimB),
+ dtype=numpy.float64)
+
+ self.strain = numpy.array([strainA, strainB],
+ dtype=numpy.float64)
+ self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+ self.stateVarsUpdated = numpy.zeros( (numLocs, 1 + 4 + 4),
+ dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
+ self._calcStress(strainA,
+ muA, lambdaA, refStrainRateA,
+ refStressA,
+ powerLawExponentA,
+ visStrainA, stress4A,
+ initialStressA, initialStrainA,
+ stateVarsA)
+ (self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
+ self._calcStress(strainB,
+ muB, lambdaB, refStrainRateB,
+ refStressB,
+ powerLawExponentB,
+ visStrainB, stress4B,
+ initialStressB, initialStrainB,
+ stateVarsB)
+
+ # Use state variables to compute Maxwell times (and stable time step size).
+ maxwellTimeA = self._getMaxwellTime(muA, refStrainRateA, refStressA, \
+ powerLawExponentA,
+ self.stateVarsUpdated[0,5:])
+
+ maxwellTimeB = self._getMaxwellTime(muB, refStrainRateB, refStressB, \
+ powerLawExponentB,
+ self.stateVarsUpdated[1,5:])
+
+ self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
+ return
+
+
+ def _bracket(self, effStressInitialGuess, ae, b, c, d, alpha, dt, effStressT,
+ powerLawExponentV, refStrainRateV, refStressV):
+ """
+ Function to bracket the effective stress.
+ """
+ maxIterations = 50
+ bracketFactor = 1.6
+
+ x1 = 0.0
+ x2 = 0.0
+
+ if effStressInitialGuess > 0.0:
+ x1 = effStressInitialGuess - 0.5 * effStressInitialGuess
+ x2 = effStressInitialGuess + 0.5 * effStressInitialGuess
+ else:
+ x1 = 500.0
+ x2 = 1500.0
+
+ funcValue1 = self._effStressFunc(x1, ae, b, c, d, alpha, dt,
+ effStressT, powerLawExponentV,
+ refStrainRateV, refStressV)
+ funcValue2 = self._effStressFunc(x2, ae, b, c, d, alpha, dt,
+ effStressT, powerLawExponentV,
+ refStrainRateV, refStressV)
+
+ iteration = 0
+ bracketed = False
+
+ while iteration < maxIterations:
+ if (funcValue1 * funcValue2) < 0.0:
+ bracketed = True
+ break
+ if abs(funcValue1) < abs(funcValue2):
+ x1 += bracketFactor * (x1 - x2)
+ x1 = max(x1, 0.0)
+
+ funcValue1 = self._effStressFunc(x1, ae, b, c, d, alpha, dt,
+ effStressT, powerLawExponentV,
+ refStrainRateV, refStressV)
+ else:
+ x2 += bracketFactor * (x1 - x2)
+ x2 = max(x2, 0.0)
+
+ funcValue2 = self._effStressFunc(x2, ae, b, c, d, alpha, dt,
+ effStressT, powerLawExponentV,
+ refStrainRateV, refStressV)
+ iteration += 1
+
+ if bracketed == False:
+ raise RuntimeError("Unable to bracket root.")
+
+ return x1, x2
+
+
+ def _getMaxwellTime(self, mu, refStrainRate, refStress, powerLawExponent,
+ stress4):
+ """
+ Compute Maxwell time from stress, reference stress and strain rate, shear
+ modulus, and power-law exponent.
+ """
+ meanStress = (stress4[0] + stress4[1] + stress4[2])/3.0
+ devStress = numpy.array(stress4, dtype=numpy.float64)
+
+ devStress[0] = stress4[0] - meanStress
+ devStress[1] = stress4[1] - meanStress
+ devStress[2] = stress4[2] - meanStress
+ devStress[3] = stress4[3]
+
+ devStressProd = self._scalarProduct(devStress, devStress)
+ effStress = (0.5 * devStressProd)**0.5
+ maxwellTime = 1.0
+ if (effStress != 0.0):
+ maxwellTime = 0.5 * (refStress/effStress)**(powerLawExponent - 1.0) * \
+ (refStress/mu)/refStrainRate
+
+ return maxwellTime
+
+
+ def _scalarProduct(self, tensor1, tensor2):
+ """
+ Compute the scalar product of two tensors stored in vector form.
+ """
+ scalarProduct = tensor1[0] * tensor2[0] + \
+ tensor1[1] * tensor2[1] + \
+ tensor1[2] * tensor2[2] + \
+ 2.0 * tensor1[3] * tensor2[3]
+ return scalarProduct
+
+
+ def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
+ muV, lambdaV, refStrainRateV, refStressV,
+ powerLawExponentV, visStrainT, stressT,
+ initialStress, initialStrain, stateVars):
+ """
+ Function to compute a particular stress component as a function of a
+ strain component.
+ """
+ strainTest = numpy.array(strainTpdt, dtype=numpy.float64)
+ strainTest[strainComp] = strainVal
+ stressTpdt, visStrainTpdt = self._computeStress(strainTest, muV, lambdaV,
+ refStrainRateV,
+ refStressV,
+ powerLawExponentV,
+ visStrainT,
+ stressT,
+ initialStress,
+ initialStrain,
+ stateVars)
+ return stressTpdt[stressComp]
+
+
+ def _effStressFunc(self, effStressTpdt, ae, b, c, d, alpha, dt, effStressT,
+ powerLawExponentV, refStrainRateV, refStressV):
+ """
+ Function to compute effective stress function for a given effective stress.
+ """
+
+ factor1 = 1.0 - alpha
+ effStressTau = factor1 * effStressT + alpha * effStressTpdt
+ gammaTau = refStrainRateV * (effStressTau/refStressV)** \
+ (powerLawExponentV - 1.0) / refStressV
+ a = ae + alpha * dt * gammaTau
+ effStressFunc = a * a * effStressTpdt * effStressTpdt - b + \
+ c * gammaTau - d * d * gammaTau * gammaTau
+
+ return effStressFunc
+
+
+ def _computeStress(self, strainTpdt, muV, lambdaV, refStrainRateV, refStressV,
+ powerLawExponentV, visStrainT, stressT,
+ initialStress, initialStrain, stateVars):
+ """
+ Function to compute stresses and viscous strains using the
+ effective stress function algorithm.
+ """
+ import scipy.optimize
+
+ # Constants
+ mu2 = 2.0 * muV
+ lamPlusMu = lambdaV + muV
+ bulkModulus = lambdaV + 2.0 * muV/3.0
+ ae = 1.0/mu2
+ timeFac = self.dt * (1.0 - self.alpha)
+ diag = numpy.array([1.0, 1.0, 1.0, 0.0], dtype=numpy.float64)
+
+ # Initial stress values
+ initialStress4 = numpy.array([initialStress[0], initialStress[1],
+ stateVars[0], initialStress[2]],
+ dtype=numpy.float64)
+ meanStressInitial = (initialStress[0] + initialStress[1] +
+ stateVars[0])/3.0
+ devStressInitial = initialStress4 - meanStressInitial * diag
+ stressInvar2Initial = 0.5 * self._scalarProduct(devStressInitial,
+ devStressInitial)
+
+ # Initial strain values
+ initialStrain4 = numpy.array([initialStrain[0], initialStrain[1],
+ 0.0, initialStrain[2]], dtype=numpy.float64)
+ meanStrainInitial = (initialStrain[0] + initialStrain[1])/3.0
+
+ # Values for current time step
+ strainTpdt4 = numpy.array([strainTpdt[0], strainTpdt[1],
+ 0.0, strainTpdt[2]], dtype=numpy.float64)
+ meanStrainTpdt = (strainTpdt[0] + strainTpdt[1])/3.0 - meanStrainInitial
+ meanStressTpdt = 3.0 * bulkModulus * meanStrainTpdt
+ strainPPTpdt = strainTpdt4 - meanStrainTpdt * diag - \
+ visStrainT - initialStrain4
+ strainPPInvar2Tpdt = 0.5 * self._scalarProduct(strainPPTpdt, strainPPTpdt)
+
+ # Values for previous time step
+ meanStressT = (stressT[0] + stressT[1] + stressT[2])/3.0
+ devStressT = stressT - diag * meanStressT
+ stressInvar2T = 0.5 * self._scalarProduct(devStressT, devStressT)
+ effStressT = stressInvar2T**0.5
+
+ # Finish defining parameters needed for root-finding algorithm.
+ b = strainPPInvar2Tpdt + \
+ ae * self._scalarProduct(strainPPTpdt, devStressInitial) + \
+ ae * ae * stressInvar2Initial
+ c = (self._scalarProduct(strainPPTpdt, devStressT) + \
+ ae * self._scalarProduct(devStressT, devStressInitial)) * timeFac
+ d = timeFac * effStressT
+
+ # Bracket the root
+ effStressInitialGuess = effStressT
+
+ x1, x2 = self._bracket(effStressInitialGuess, ae, b, c, d, self.alpha,
+ self.dt, effStressT, powerLawExponentV,
+ refStrainRateV, refStressV)
+
+ # Find the root using Brent's method (from scipy)
+ rootTolerance = 1.0e-14
+ effStressTpdt = scipy.optimize.brentq(self._effStressFunc, x1, x2,
+ args=(ae, b, c, d, self.alpha,
+ self.dt, effStressT,
+ powerLawExponentV,
+ refStrainRateV, refStressV),
+ xtol=rootTolerance)
+
+ # Compute stresses from the effective stress.
+ effStressTau = (1.0 - self.alpha) * effStressT + self.alpha * effStressTpdt
+ gammaTau = refStrainRateV * ((effStressTau/refStressV)** \
+ (powerLawExponentV - 1.0)) / refStressV
+ factor1 = 1.0/(ae + self.alpha * self.dt * gammaTau)
+ factor2 = timeFac * gammaTau
+ devStressTpdt = 0.0
+ stressTpdt = numpy.zeros( (4), dtype=numpy.float64)
+ visStrainTpdt = numpy.zeros( (4), dtype=numpy.float64)
+
+ for iComp in range(4):
+ devStressTpdt = factor1 * (strainPPTpdt[iComp] - \
+ factor2 * devStressT[iComp] + \
+ ae * devStressInitial[iComp])
+ stressTpdt[iComp] = devStressTpdt + diag[iComp] * \
+ (meanStressTpdt + meanStressInitial)
+ devStressTau = (1.0 - self.alpha) * devStressT[iComp] + \
+ self.alpha * devStressTpdt
+ deltaVisStrain = self.dt * gammaTau * devStressTau
+ visStrainTpdt[iComp] = visStrainT[iComp] + deltaVisStrain
+
+ return stressTpdt, visStrainTpdt
+
+
+ def _calcStress(self, strainV, muV, lambdaV, refStrainRateV, refStressV,
+ powerLawExponentV,visStrainV, stressV,
+ initialStressV, initialStrainV, stateVars):
+ """
+ Compute stress, updated state variables and derivative of elasticity matrix.
+ This assumes behavior is always viscoelastic.
+ """
+ import scipy.misc
+
+ # Define some numpy arrays
+ strainTpdt = numpy.array(strainV, dtype=numpy.float64)
+ visStrainT = numpy.array(visStrainV, dtype=numpy.float64)
+ stressT = numpy.array(stressV, dtype=numpy.float64)
+ initialStress = numpy.array(initialStressV, dtype=numpy.float64)
+ initialStrain = numpy.array(initialStrainV, dtype=numpy.float64)
+
+ stressTpdt, visStrainTpdt = self._computeStress(strainTpdt, muV, lambdaV,
+ refStrainRateV, refStressV,
+ powerLawExponentV,
+ visStrainT, stressT,
+ initialStress,
+ initialStrain,
+ stateVars)
+
+ stateVarsUpdated = numpy.concatenate(([stateVars[0]], visStrainTpdt,
+ stressTpdt))
+
+ # Compute components of tangent constitutive matrix using numerical
+ # derivatives.
+ derivDx = 1.0e-12
+ derivOrder = 3
+ elasticConstsList = []
+ components = [0, 1, 3]
+ stress3 = stressTpdt[components]
+
+ for stressComp in components:
+ for strainComp in range(tensorSize):
+ dStressDStrain = scipy.misc.derivative(self._calcStressComponent,
+ strainTpdt[strainComp],
+ dx=derivDx,
+ args=(strainComp,
+ stressComp,
+ strainTpdt, muV, lambdaV,
+ refStrainRateV, refStressV,
+ powerLawExponentV,
+ visStrainT,
+ stressT, initialStress,
+ initialStrain,
+ stateVars),
+ order=derivOrder)
+ elasticConstsList.append(dStressDStrain)
+
+ elasticConsts = numpy.array(elasticConstsList, dtype=numpy.float64)
+
+ return (elasticConsts, numpy.ravel(stress3),
+ numpy.ravel(stateVarsUpdated))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = PowerLawPlaneStrainTimeDep()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,303 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application powerlawplanestraintimedep.
+
+#include "PowerLawPlaneStrainTimeDepData.hh"
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_dimension = 2;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numLocs = 2;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numProperties = 6;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numStateVars = 3;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numDBProperties = 6;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numDBStateVars = 9;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numPropsQuadPt = 6;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numVarsQuadPt = 9;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_lengthScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_timeScale = 1.00000000e+00;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_pressureScale = 2.25000000e+10;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_densityScale = 1.00000000e+03;
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_dtStableImplicit = 4.44444444e+06;
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numPropertyValues[] = {
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::PowerLawPlaneStrainTimeDepData::_numStateVarValues[] = {
+1,
+4,
+4,
+};
+
+const char* pylith::materials::PowerLawPlaneStrainTimeDepData::_dbPropertyValues[] = {
+"density",
+"vs",
+"vp",
+"reference-strain-rate",
+"reference-stress",
+"power-law-exponent",
+};
+
+const char* pylith::materials::PowerLawPlaneStrainTimeDepData::_dbStateVarValues[] = {
+"stress-zz-initial",
+"viscous-strain-xx",
+"viscous-strain-yy",
+"viscous-strain-zz",
+"viscous-strain-xy",
+"stress4-xx",
+"stress4-yy",
+"stress4-zz",
+"stress4-xy",
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_dbProperties[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 1.00000000e-06,
+ 2.00000000e+12,
+ 1.00000000e+00,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 1.00000000e-06,
+ 1.25992105e+10,
+ 3.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_dbStateVars[] = {
+ 2.30000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.30000000e+04,
+ 0.00000000e+00,
+ 0.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_properties[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 1.00000000e-06,
+ 2.00000000e+12,
+ 1.00000000e+00,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 1.00000000e-06,
+ 1.25992105e+10,
+ 3.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_stateVars[] = {
+ 2.30000000e+04,
+ 4.10000000e-05,
+ 4.20000000e-05,
+ 4.30000000e-05,
+ 4.40000000e-05,
+ 3.10000000e+04,
+ 3.20000000e+04,
+ 2.30000000e+04,
+ 3.40000000e+04,
+ 5.30000000e+04,
+ 1.10000000e-05,
+ 1.20000000e-05,
+ 1.30000000e-05,
+ 1.40000000e-05,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.30000000e+04,
+ 5.40000000e+04,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_propertiesNondim[] = {
+ 2.50000000e+00,
+ 1.00000000e+00,
+ 1.00000000e+00,
+ 1.00000000e-06,
+ 8.88888889e+01,
+ 1.00000000e+00,
+ 2.00000000e+00,
+ 1.28000000e-01,
+ 1.28000000e-01,
+ 1.00000000e-06,
+ 5.59964911e-01,
+ 3.00000000e+00,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_stateVarsNondim[] = {
+ 1.02222222e-06,
+ 4.10000000e-05,
+ 4.20000000e-05,
+ 4.30000000e-05,
+ 4.40000000e-05,
+ 1.37777778e-06,
+ 1.42222222e-06,
+ 1.02222222e-06,
+ 1.51111111e-06,
+ 2.35555556e-06,
+ 1.10000000e-05,
+ 1.20000000e-05,
+ 1.30000000e-05,
+ 1.40000000e-05,
+ 2.26666667e-06,
+ 2.31111111e-06,
+ 2.35555556e-06,
+ 2.40000000e-06,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_strain[] = {
+ 1.10000000e-04,
+ 1.20000000e-04,
+ 1.40000000e-04,
+ 4.10000000e-04,
+ 4.20000000e-04,
+ 4.40000000e-04,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_stress[] = {
+ 5.08551746e+06,
+ 5.53550274e+06,
+ 2.85250536e+06,
+ 4.03404000e+06,
+ 4.08112000e+06,
+ 2.12760000e+06,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_elasticConsts[] = {
+ 6.74326522e+10,
+ 2.25336747e+10,
+ 0.00000000e+00,
+ 2.25336747e+10,
+ 6.74326518e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.48989770e+10,
+ 8.63999990e+09,
+ 2.88000004e+09,
+ 0.00000000e+00,
+ 2.88000004e+09,
+ 8.64000013e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.75999985e+09,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_initialStress[] = {
+ 2.10000000e+04,
+ 2.20000000e+04,
+ 2.40000000e+04,
+ 5.10000000e+04,
+ 5.20000000e+04,
+ 5.40000000e+04,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_initialStrain[] = {
+ 3.60000000e-05,
+ 3.50000000e-05,
+ 3.30000000e-05,
+ 6.10000000e-05,
+ 6.20000000e-05,
+ 6.60000000e-05,
+};
+
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_stateVarsUpdated[] = {
+ 2.30000000e+04,
+ 4.09551675e-05,
+ 4.19777168e-05,
+ 4.27842521e-05,
+ 4.41443253e-05,
+ 5.08551746e+06,
+ 5.53550274e+06,
+ 1.67520866e+06,
+ 2.85250536e+06,
+ 5.30000000e+04,
+ 1.10000000e-05,
+ 1.20000000e-05,
+ 1.29999999e-05,
+ 1.40000002e-05,
+ 4.03404000e+06,
+ 4.08112000e+06,
+ 2.01428000e+06,
+ 2.12760000e+06,
+};
+
+pylith::materials::PowerLawPlaneStrainTimeDepData::PowerLawPlaneStrainTimeDepData(void)
+{ // constructor
+ dimension = _dimension;
+ numLocs = _numLocs;
+ numProperties = _numProperties;
+ numStateVars = _numStateVars;
+ numDBProperties = _numDBProperties;
+ numDBStateVars = _numDBStateVars;
+ numPropsQuadPt = _numPropsQuadPt;
+ numVarsQuadPt = _numVarsQuadPt;
+ lengthScale = _lengthScale;
+ timeScale = _timeScale;
+ pressureScale = _pressureScale;
+ densityScale = _densityScale;
+ dtStableImplicit = _dtStableImplicit;
+ numPropertyValues = const_cast<int*>(_numPropertyValues);
+ numStateVarValues = const_cast<int*>(_numStateVarValues);
+ dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+ dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+ dbProperties = const_cast<PylithScalar*>(_dbProperties);
+ dbStateVars = const_cast<PylithScalar*>(_dbStateVars);
+ properties = const_cast<PylithScalar*>(_properties);
+ stateVars = const_cast<PylithScalar*>(_stateVars);
+ propertiesNondim = const_cast<PylithScalar*>(_propertiesNondim);
+ stateVarsNondim = const_cast<PylithScalar*>(_stateVarsNondim);
+ density = const_cast<PylithScalar*>(_density);
+ strain = const_cast<PylithScalar*>(_strain);
+ stress = const_cast<PylithScalar*>(_stress);
+ elasticConsts = const_cast<PylithScalar*>(_elasticConsts);
+ initialStress = const_cast<PylithScalar*>(_initialStress);
+ initialStrain = const_cast<PylithScalar*>(_initialStrain);
+ stateVarsUpdated = const_cast<PylithScalar*>(_stateVarsUpdated);
+} // constructor
+
+pylith::materials::PowerLawPlaneStrainTimeDepData::~PowerLawPlaneStrainTimeDepData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.hh (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.hh 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,110 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application powerlawplanestraintimedep.
+
+#if !defined(pylith_materials_powerlawplanestraintimedepdata_hh)
+#define pylith_materials_powerlawplanestraintimedepdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class PowerLawPlaneStrainTimeDepData;
+ } // pylith
+} // materials
+
+class pylith::materials::PowerLawPlaneStrainTimeDepData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ PowerLawPlaneStrainTimeDepData(void);
+
+ /// Destructor
+ ~PowerLawPlaneStrainTimeDepData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numLocs;
+
+ static const int _numProperties;
+
+ static const int _numStateVars;
+
+ static const int _numDBProperties;
+
+ static const int _numDBStateVars;
+
+ static const int _numPropsQuadPt;
+
+ static const int _numVarsQuadPt;
+
+ static const PylithScalar _lengthScale;
+
+ static const PylithScalar _timeScale;
+
+ static const PylithScalar _pressureScale;
+
+ static const PylithScalar _densityScale;
+
+ static const PylithScalar _dtStableImplicit;
+
+ static const int _numPropertyValues[];
+
+ static const int _numStateVarValues[];
+
+ static const char* _dbPropertyValues[];
+
+ static const char* _dbStateVarValues[];
+
+ static const PylithScalar _dbProperties[];
+
+ static const PylithScalar _dbStateVars[];
+
+ static const PylithScalar _properties[];
+
+ static const PylithScalar _stateVars[];
+
+ static const PylithScalar _propertiesNondim[];
+
+ static const PylithScalar _stateVarsNondim[];
+
+ static const PylithScalar _density[];
+
+ static const PylithScalar _strain[];
+
+ static const PylithScalar _stress[];
+
+ static const PylithScalar _elasticConsts[];
+
+ static const PylithScalar _initialStress[];
+
+ static const PylithScalar _initialStrain[];
+
+ static const PylithScalar _stateVarsUpdated[];
+
+};
+
+#endif // pylith_materials_powerlawplanestraintimedepdata_hh
+
+// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/generate.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/generate.sh 2012-05-15 04:49:18 UTC (rev 20118)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/materials/data/generate.sh 2012-05-15 04:50:33 UTC (rev 20119)
@@ -83,6 +83,11 @@
--data.object=DruckerPragerPlaneStrainElasticData \
--data.parent=ElasticMaterialData
+ python PowerLawPlaneStrainElastic.py \
+ --data.namespace=pylith,materials \
+ --data.object=PowerLawPlaneStrainElasticData \
+ --data.parent=ElasticMaterialData
+
# 1-D ----------------------------------------------------------------
python ElasticStrain1D.py \
@@ -130,6 +135,11 @@
--data.object=PowerLaw3DTimeDepData \
--data.parent=ElasticMaterialData
+ python PowerLawPlaneStrainTimeDep.py \
+ --data.namespace=pylith,materials \
+ --data.object=PowerLawPlaneStrainTimeDepData \
+ --data.parent=ElasticMaterialData
+
python DruckerPrager3DTimeDep.py \
--data.namespace=pylith,materials \
--data.object=DruckerPrager3DTimeDepData \
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/Makefile.am 2012-05-15 04:49:18 UTC (rev 20118)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/Makefile.am 2012-05-15 04:50:33 UTC (rev 20119)
@@ -41,6 +41,7 @@
TestGenMaxwellPlaneStrain.py \
TestGenMaxwellQpQsIsotropic3D.py \
TestPowerLaw3D.py \
+ TestPowerLawPlaneStrain.py \
TestDruckerPrager3D.py \
TestDruckerPragerPlaneStrain.py
Added: short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/TestPowerLawPlaneStrain.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/TestPowerLawPlaneStrain.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/TestPowerLawPlaneStrain.py 2012-05-15 04:50:33 UTC (rev 20119)
@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/materials/TestPowerLawPlaneStrain.py
+
+## @brief Unit testing of PowerLawPlaneStrain object.
+
+import unittest
+
+from pylith.materials.PowerLawPlaneStrain import PowerLawPlaneStrain
+
+# ----------------------------------------------------------------------
+class TestPowerLawPlaneStrain(unittest.TestCase):
+ """
+ Unit testing of PowerLawPlaneStrain object.
+ """
+
+ def setUp(self):
+ """
+ Setup test subject.
+ """
+ self.material = PowerLawPlaneStrain()
+ return
+
+
+ def test_constructor(self):
+ """
+ Test constructor.
+ """
+ self.assertEqual(2, self.material.dimension())
+ return
+
+
+ def test_useLinearBehavior(self):
+ """
+ Test useLinearBehavior().
+ """
+ self.material.useLinearBehavior(False)
+ return
+
+
+ def testHasStateVars(self):
+ self.failUnless(self.material.hasStateVars())
+ return
+
+
+ def testTensorSize(self):
+ self.assertEqual(3, self.material.tensorSize())
+ return
+
+
+ def testNeedNewJacobian(self):
+ """
+ Test needNewJacobian().
+ """
+ # Default should be False.
+ self.failIf(self.material.needNewJacobian())
+
+ # Should require a new Jacobian even if time step is the same.
+ self.material.timeStep(1.0)
+ self.failUnless(self.material.needNewJacobian())
+ self.material.timeStep(2.0)
+ self.failUnless(self.material.needNewJacobian())
+
+ self.material.timeStep(2.0)
+ self.failUnless(self.material.needNewJacobian())
+ return
+
+
+ def test_factory(self):
+ """
+ Test factory method.
+ """
+ from pylith.materials.PowerLawPlaneStrain import material
+ m = material()
+ return
+
+
+# End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/testmaterials.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/testmaterials.py 2012-05-15 04:49:18 UTC (rev 20118)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/materials/testmaterials.py 2012-05-15 04:50:33 UTC (rev 20119)
@@ -95,6 +95,9 @@
from TestPowerLaw3D import TestPowerLaw3D
suite.addTest(unittest.makeSuite(TestPowerLaw3D))
+ from TestPowerLawPlaneStrain import TestPowerLawPlaneStrain
+ suite.addTest(unittest.makeSuite(TestPowerLawPlaneStrain))
+
from TestDruckerPrager3D import TestDruckerPrager3D
suite.addTest(unittest.makeSuite(TestDruckerPrager3D))
More information about the CIG-COMMITS
mailing list