[cig-commits] r7065 - in
short/3D/PyLith/trunk/unittests/libtests/materials: . data
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon Jun 4 20:47:24 PDT 2007
Author: willic3
Date: 2007-06-04 20:47:23 -0700 (Mon, 04 Jun 2007)
New Revision: 7065
Added:
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
Removed:
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py
Modified:
short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.hh
short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh
Log:
More work on unit tests for Maxwell material.
Everything compiles and python generates data, but testing is needed.
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am 2007-06-05 03:47:23 UTC (rev 7065)
@@ -28,6 +28,8 @@
TestElasticPlaneStress.cc \
TestElasticStrain1D.cc \
TestElasticStress1D.cc \
+ TestMaxwellIsotropic3D.cc \
+ TestElasticPlaneStrain.cc \
test_materials.cc
noinst_HEADERS = \
@@ -41,7 +43,9 @@
data/ElasticPlaneStrainData.cc \
data/ElasticPlaneStressData.cc \
data/ElasticStrain1DData.cc \
- data/ElasticStress1DData.cc
+ data/ElasticStress1DData.cc \
+ data/MaxwellIsotropic3DElasticData.cc \
+ data/MaxwellIsotropic3DTimeDepData.cc
noinst_HEADERS += \
data/MaterialData.hh \
@@ -51,6 +55,8 @@
data/ElasticPlaneStressData.hh \
data/ElasticStrain1DData.hh \
data/ElasticStress1DData.hh \
+ data/MaxwellIsotropic3DElasticData.hh \
+ data/MaxwellIsotropic3DTimeDepData.hh \
data/header.hh
testmaterials_LDFLAGS = $(PETSC_LIB) $(PYTHON_BLDLIBRARY)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.hh 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.hh 2007-06-05 03:47:23 UTC (rev 7065)
@@ -28,7 +28,8 @@
namespace materials {
class MaxwellIsotropic3D;
class TestMaxwellIsotropic3D;
- class MaxwellIsotropic3DData;
+ class MaxwellIsotropic3DElasticData;
+ class MaxwellIsotropic3DTimeDepData;
} // materials
} // pylith
@@ -42,9 +43,12 @@
CPPUNIT_TEST( testDBValues );
CPPUNIT_TEST( testParameters );
CPPUNIT_TEST( testCalcDensity );
- CPPUNIT_TEST( testCalcStress );
- CPPUNIT_TEST( testCalcElasticConsts );
- CPPUNIT_TEST( testUpdateState );
+ CPPUNIT_TEST( testCalcStressElastic );
+ CPPUNIT_TEST( testCalcStressTimeDep );
+ CPPUNIT_TEST( testCalcElasticConstsElastic );
+ CPPUNIT_TEST( testCalcElasticConstsTimeDep );
+ CPPUNIT_TEST( testUpdateStateElastic );
+ CPPUNIT_TEST( testUpdateStateTimeDep );
CPPUNIT_TEST( testTimeStep );
CPPUNIT_TEST_SUITE_END();
@@ -63,15 +67,24 @@
/// Test calcDensity()
void testCalcDensity(void);
- /// Test calcStress()
- void testCalcStress(void);
+ /// Test calcStressElastic()
+ void testCalcStressElastic(void);
- /// Test calcElasticConsts()
- void testCalcElasticConsts(void);
+ /// Test calcStressTimeDep()
+ void testCalcStressTimeDep(void);
- /// Test updateState()
- void testUpdateState(void);
+ /// Test calcElasticConstsElastic()
+ void testCalcElasticConstsElastic(void);
+ /// Test calcElasticConstsTimeDep()
+ void testCalcElasticConstsTimeDep(void);
+
+ /// Test updateStateElastic()
+ void testUpdateStateElastic(void);
+
+ /// Test updateStateTimeDep()
+ void testUpdateStateTimeDep(void);
+
/// Test timeStep()
void testTimeStep(void);
Deleted: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py 2007-06-05 03:47:23 UTC (rev 7065)
@@ -1,189 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file unittests/libtests/materials/data/MaxwellIsotropic3D.py
-
-## @brief Python application for generating C++ data files for testing
-## C++ MaxwellIsotropic3D object.
-
-from ElasticMaterialApp import ElasticMaterialApp
-
-import numpy
-
-# ----------------------------------------------------------------------
-# MaxwellIsotropic3D class
-class MaxwellIsotropic3D(ElasticMaterialApp):
- """
- Python application for generating C++ data files for testing C++
- MaxwellIsotropic3D object.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="maxwellisotropic3d"):
- """
- Constructor.
- """
- ElasticMaterialApp.__init__(self, name)
-
- self.dimension = 3
-
- self.numDBValues = 4
- self.dbValues = ["density", "vs", "vp" , "viscosity"]
- self.numParameters = 6
- self.numParamValues = [1, 1, 1, 1, 6, 6]
- self.parameterNames = ["density", "mu", "lambda", "maxwellTime", "strainT", "visStrain"]
-
- self.dt = 2.0e5
-
- densityA = 2500.0
- vsA = 3000.0
- vpA = vsA*3**0.5
- viscosityA = 1.0e18
- strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
- meanStrainA = (strainA[1] + strainA[2] + strainA[3])/3.0
-
- densityB = 2000.0
- vsB = 1200.0
- vpB = vsB*3**0.5
- viscosityB = 1.0e19
- strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
- meanStrainB = (strainB[1] + strainB[2] + strainB[3])/3.0
-
- diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
-
- self.dbData = numpy.array([ [densityA, vsA, vpA, viscosityA],
- [densityB, vsB, vpB, viscosityB] ],
- dtype=numpy.float64)
- muA = vsA*vsA*densityA
- lambdaA = vpA*vpA*densityA - 2.0*muA
- maxwellTimeA = viscosityA/muA
-
- muB = vsB*vsB*densityB
- lambdaB = vpB*vpB*densityB - 2.0*muB
- maxwellTimeB = viscosityB/muB
-
- # Simplest approach for now is to assume this is the first step after the elastic solution.
- # In that case, both the total strain from the last step (strainT) and the total viscous
- # strain (visStrain) are defined by the assigned elastic strain.
- strainTA = strainA.copy()
- strainTB = strainB.copy()
- for i in range(6):
- visStrainA[i] = strainA[i] - diag[i] * meanStrainA[i]
- visStrainB[i] = strainB[i] - diag[i] * meanStrainB[i]
-
- self.parameterData = numpy.array([ [densityA, muA, lambdaA, maxwellTimeA, strainTA, visStrainA],
- [densityB, muB, lambdaB, maxwellTimeB, strainTB, visStrainB] ],
- dtype=numpy.float64)
-
- self.numLocs = 2
- numElasticConsts = 21
- self.density = numpy.array([densityA, densityB],
- dtype=numpy.float64)
-
- self.strain = numpy.array([strainA, strainB],
- dtype=numpy.float64)
- self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
- self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
- dtype=numpy.float64)
-
- (self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, muA, lambdaA, maxwellTimeA, strainTA, visStrainA)
- (self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, muB, lambdaB, maxwellTimeB, strainTB, visStrainB)
- return
-
-
- def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, strainTV, visStrainV):
- """
- Compute stress and derivative of elasticity matrix.
- This assumes behavior is always viscoelastic.
- """
- bulkModulus = lambdaV + 2.0 * muV/3.0
-
- diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
- traceStrainT = strainTV[0] + strainTV[1] + strainTV[2]
- traceStrainTpdt = strainV[0] + strainV[1] + strainV[2]
- meanStrainT = traceStrainT/3.0
- meanStrainTpdt = traceStrainTpdt/3.0
- timeFrac = 1.0e-5
- numTerms = 5
- dq = 0.0
- if maxwellTimeV < timeFrac*self.dt:
- fSign = 1.0
- factorial = 1.0
- fraction = 1.0
- dq = 1.0
- for iTerm in range(2, numTerms + 1):
- factorial *= iTerm
- fSign *= -1.0
- fraction *= self.dt/maxwellTimeV
- dq += fSign*fraction/factorial
- else:
- dq = maxwellTimeV*(1.0-exp(-self.dt/maxwellTimeV))/self.dt
-
- visFac = muV*dq/3.0
-
- C1111 = bulkModulus + 4.0*visFac
- C1122 = bulkModulus - 2.0*visFac
- C1133 = C1122
- C1112 = 0.0
- C1123 = 0.0
- C1113 = 0.0
- C2222 = C1111
- C2233 = C1122
- C2212 = 0.0
- C2223 = 0.0
- C2213 = 0.0
- C3333 = C1111
- C3312 = 0.0
- C3323 = 0.0
- C3313 = 0.0
- C1212 = 3.0*visFac
- C1223 = 0.0
- C1213 = 0.0
- C2323 = C1212
- C2313 = 0.0
- C1313 = C1212
- elasticConsts = numpy.array([C1111, C1122, C1133, C1112, C1123, C1113,
- C2222, C2233, C2212, C2223, C2213,
- C3333, C3312, C3323, C3313,
- C1212, C1223, C1213,
- C2323, C2313,
- C1313], dtype=numpy.float64)
-
- stress = numpy.zeros( 6, dtype=numpy.float64)
-
- expFac = exp(-self.dt/maxwellTimeV)
- elasFac = 2.0*muV
- devStrainTpdt = 0.0
- devStrainT = 0.0
- devStressTpdt = 0.0
- visStrain = 0.0
- for iComp in range(6):
- devStrainTpdt = strainV[iComp] - diag[iComp]*meanStrainTpdt
- devStrainT = strainTV[iComp] - diag[iComp]*meanStrainT
- visStrain = expFac*visStrainV[iComp] + dq*(devStrainTpdt - devStrainT)
- devStressTpdt = elasFac*visStrain
- stress[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt
-
- return (elasticConsts, numpy.ravel(stress))
-
-
-# MAIN /////////////////////////////////////////////////////////////////
-if __name__ == "__main__":
-
- app = MaxwellIsotropic3D()
- app.run()
-
-
-# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py 2007-06-05 03:47:23 UTC (rev 7065)
@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/MaxwellIsotropic3DElastic.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ MaxwellIsotropic3D object with elastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+
+# ----------------------------------------------------------------------
+# MaxwellIsotropic3DElastic class
+class MaxwellIsotropic3DElastic(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ MaxwellIsotropic3D object with elastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="maxwellisotropic3delastic"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ self.dimension = 3
+
+ self.numDBValues = 4
+ self.dbValues = ["density", "vs", "vp", "viscosity"]
+ self.numParameters = 6
+ self.numParamValues = [1, 1, 1, 1, 6, 6]
+ self.parameterNames = ["density", "mu", "lambda", "maxwellTime", "strainT", "visStrain"]
+
+ self.dt = 2.0e5
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ viscosityA = 1.0e18
+ strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ viscosityB = 1.0e18
+ strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
+
+ self.dbData = numpy.array([ [densityA, vsA, vpA, viscosityA],
+ [densityB, vsB, vpB, viscosityB] ],
+ dtype=numpy.float64)
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ maxwellTimeA = viscosityA/muA
+ strainTA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ visStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ maxwellTimeB = viscosityB/muB
+ strainTB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+ visStrainB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+ vecParamsA = numpy.hstack((strainTA, visStrainA))
+ vecParamsB = numpy.hstack((strainTB, visStrainB))
+ vecParams = numpy.vstack((vecParamsA, vecParamsB))
+ scalarParams = numpy.array([ [densityA, muA, lambdaA, maxwellTimeA],
+ [densityB, muB, lambdaB, maxwellTimeB] ],
+ dtype=numpy.float64)
+ self.parameterData = numpy.hstack((scalarParams, vecParams))
+
+ self.numLocs = 2
+ numElasticConsts = 21
+ self.density = numpy.array([densityA, densityB],
+ dtype=numpy.float64)
+
+ self.strain = numpy.array([strainA, strainB],
+ dtype=numpy.float64)
+ self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA)
+ (self.elasticConsts[1,:], self.stress[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB)
+ return
+
+
+ def _calcStress(self, strainV, muV, lambdaV):
+ """
+ Compute stress and derivative of elasticity matrix.
+ """
+ C1111 = lambdaV + 2.0*muV
+ C1122 = lambdaV
+ C1133 = lambdaV
+ C1112 = 0.0
+ C1123 = 0.0
+ C1113 = 0.0
+ C2222 = lambdaV + 2.0*muV
+ C2233 = lambdaV
+ C2212 = 0.0
+ C2223 = 0.0
+ C2213 = 0.0
+ C3333 = lambdaV + 2.0*muV
+ C3312 = 0.0
+ C3323 = 0.0
+ C3313 = 0.0
+ C1212 = 2.0*muV
+ C1223 = 0.0
+ C1213 = 0.0
+ C2323 = 2.0*muV
+ C2313 = 0.0
+ C1313 = 2.0*muV
+ elasticConsts = numpy.array([C1111, C1122, C1133, C1112, C1123, C1113,
+ C2222, C2233, C2212, C2223, C2213,
+ C3333, C3312, C3323, C3313,
+ C1212, C1223, C1213,
+ C2323, C2313,
+ C1313], dtype=numpy.float64)
+
+ strain = numpy.reshape(strainV, (6,1))
+ elastic = numpy.array([ [C1111, C1122, C1133, C1112, C1123, C1113],
+ [C1122, C2222, C2233, C2212, C2223, C2213],
+ [C1133, C2233, C3333, C3312, C3323, C3313],
+ [C1112, C2212, C3312, C1212, C1223, C1213],
+ [C1123, C2223, C3323, C1223, C2323, C2313],
+ [C1113, C2213, C3313, C1213, C2313, C1313] ],
+ dtype=numpy.float64)
+ stress = numpy.dot(elastic, strain)
+ return (elasticConsts, numpy.ravel(stress))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = MaxwellIsotropic3DElastic()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.cc 2007-06-05 03:47:23 UTC (rev 7065)
@@ -0,0 +1,198 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application maxwellisotropic3delastic.
+
+#include "MaxwellIsotropic3DElasticData.hh"
+
+const int pylith::materials::MaxwellIsotropic3DElasticData::_dimension = 3;
+
+const int pylith::materials::MaxwellIsotropic3DElasticData::_numDBValues = 4;
+
+const int pylith::materials::MaxwellIsotropic3DElasticData::_numParameters = 6;
+
+const int pylith::materials::MaxwellIsotropic3DElasticData::_numLocs = 2;
+
+const int pylith::materials::MaxwellIsotropic3DElasticData::_numParamValues[] = {
+1,
+1,
+1,
+1,
+6,
+6,
+};
+
+const char* pylith::materials::MaxwellIsotropic3DElasticData::_dbValues[] = {
+"density",
+"vs",
+"vp",
+"viscosity",
+};
+
+const char* pylith::materials::MaxwellIsotropic3DElasticData::_parameterNames[] = {
+"density",
+"mu",
+"lambda",
+"maxwellTime",
+"strainT",
+"visStrain",
+};
+
+const double pylith::materials::MaxwellIsotropic3DElasticData::_dbData[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 1.00000000e+18,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 1.00000000e+18,
+};
+
+const double pylith::materials::MaxwellIsotropic3DElasticData::_parameterData[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 4.44444444e+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,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 3.47222222e+08,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 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 double pylith::materials::MaxwellIsotropic3DElasticData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const double pylith::materials::MaxwellIsotropic3DElasticData::_strain[] = {
+ 1.10000000e-04,
+ 2.20000000e-04,
+ 3.30000000e-04,
+ 4.40000000e-04,
+ 5.50000000e-04,
+ 6.60000000e-04,
+ 1.20000000e-04,
+ 2.30000000e-04,
+ 3.40000000e-04,
+ 4.50000000e-04,
+ 5.60000000e-04,
+ 6.70000000e-04,
+};
+
+const double pylith::materials::MaxwellIsotropic3DElasticData::_stress[] = {
+ 1.98000000e+07,
+ 2.47500000e+07,
+ 2.97000000e+07,
+ 1.98000000e+07,
+ 2.47500000e+07,
+ 2.97000000e+07,
+ 2.67840000e+06,
+ 3.31200000e+06,
+ 3.94560000e+06,
+ 2.59200000e+06,
+ 3.22560000e+06,
+ 3.85920000e+06,
+};
+
+const double pylith::materials::MaxwellIsotropic3DElasticData::_elasticConsts[] = {
+ 6.75000000e+10,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 6.75000000e+10,
+ 2.25000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 6.75000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 0.00000000e+00,
+ 4.50000000e+10,
+ 8.64000000e+09,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 8.64000000e+09,
+ 2.88000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 8.64000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 5.76000000e+09,
+ 0.00000000e+00,
+ 5.76000000e+09,
+};
+
+pylith::materials::MaxwellIsotropic3DElasticData::MaxwellIsotropic3DElasticData(void)
+{ // constructor
+ dimension = _dimension;
+ numDBValues = _numDBValues;
+ numParameters = _numParameters;
+ numLocs = _numLocs;
+ numParamValues = const_cast<int*>(_numParamValues);
+ dbValues = const_cast<char**>(_dbValues);
+ parameterNames = const_cast<char**>(_parameterNames);
+ dbData = const_cast<double*>(_dbData);
+ parameterData = const_cast<double*>(_parameterData);
+ density = const_cast<double*>(_density);
+ strain = const_cast<double*>(_strain);
+ stress = const_cast<double*>(_stress);
+ elasticConsts = const_cast<double*>(_elasticConsts);
+} // constructor
+
+pylith::materials::MaxwellIsotropic3DElasticData::~MaxwellIsotropic3DElasticData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DElasticData.hh 2007-06-05 03:47:23 UTC (rev 7065)
@@ -0,0 +1,70 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application maxwellisotropic3delastic.
+
+#if !defined(pylith_materials_maxwellisotropic3delasticdata_hh)
+#define pylith_materials_maxwellisotropic3delasticdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class MaxwellIsotropic3DElasticData;
+ } // pylith
+} // materials
+
+class pylith::materials::MaxwellIsotropic3DElasticData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ MaxwellIsotropic3DElasticData(void);
+
+ /// Destructor
+ ~MaxwellIsotropic3DElasticData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numDBValues;
+
+ static const int _numParameters;
+
+ static const int _numLocs;
+
+ static const int _numParamValues[];
+
+ static const char* _dbValues[];
+
+ static const char* _parameterNames[];
+
+ static const double _dbData[];
+
+ static const double _parameterData[];
+
+ static const double _density[];
+
+ static const double _strain[];
+
+ static const double _stress[];
+
+ static const double _elasticConsts[];
+
+};
+
+#endif // pylith_materials_maxwellisotropic3delasticdata_hh
+
+// End of file
Copied: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py (from rev 7064, short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py 2007-06-05 03:47:23 UTC (rev 7065)
@@ -0,0 +1,199 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/materials/data/MaxwellIsotropic3DTimeDep.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ MaxwellIsotropic3D object with viscoelastic behavior.
+
+from ElasticMaterialApp import ElasticMaterialApp
+
+import numpy
+
+# ----------------------------------------------------------------------
+# MaxwellIsotropic3DTimeDep class
+class MaxwellIsotropic3DTimeDep(ElasticMaterialApp):
+ """
+ Python application for generating C++ data files for testing C++
+ MaxwellIsotropic3D object using viscoelastic behavior.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="maxwellisotropic3dtimedep"):
+ """
+ Constructor.
+ """
+ ElasticMaterialApp.__init__(self, name)
+
+ self.dimension = 3
+
+ self.numDBValues = 4
+ self.dbValues = ["density", "vs", "vp" , "viscosity"]
+ self.numParameters = 6
+ self.numParamValues = [1, 1, 1, 1, 6, 6]
+ self.parameterNames = ["density", "mu", "lambda", "maxwellTime", "strainT", "visStrain"]
+
+ self.dt = 2.0e5
+
+ densityA = 2500.0
+ vsA = 3000.0
+ vpA = vsA*3**0.5
+ viscosityA = 1.0e18
+ strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
+ meanStrainA = (strainA[1] + strainA[2] + strainA[3])/3.0
+
+ densityB = 2000.0
+ vsB = 1200.0
+ vpB = vsB*3**0.5
+ viscosityB = 1.0e19
+ strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
+ meanStrainB = (strainB[1] + strainB[2] + strainB[3])/3.0
+
+ diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
+
+ self.dbData = numpy.array([ [densityA, vsA, vpA, viscosityA],
+ [densityB, vsB, vpB, viscosityB] ],
+ dtype=numpy.float64)
+ muA = vsA*vsA*densityA
+ lambdaA = vpA*vpA*densityA - 2.0*muA
+ maxwellTimeA = viscosityA/muA
+ visStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+ muB = vsB*vsB*densityB
+ lambdaB = vpB*vpB*densityB - 2.0*muB
+ maxwellTimeB = viscosityB/muB
+ visStrainB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+ # Simplest approach for now is to assume this is the first step after the elastic solution.
+ # In that case, both the total strain from the last step (strainT) and the total viscous
+ # strain (visStrain) are defined by the assigned elastic strain.
+ strainTA = strainA[:]
+ strainTB = strainB[:]
+ for i in range(6):
+ visStrainA[i] = strainA[i] - diag[i] * meanStrainA
+ visStrainB[i] = strainB[i] - diag[i] * meanStrainB
+
+ vecParamsA = numpy.hstack((strainTA, visStrainA))
+ vecParamsB = numpy.hstack((strainTB, visStrainB))
+ vecParams = numpy.vstack((vecParamsA, vecParamsB))
+ scalarParams = numpy.array([ [densityA, muA, lambdaA, maxwellTimeA],
+ [densityB, muB, lambdaB, maxwellTimeB] ],
+ dtype=numpy.float64)
+ self.parameterData = numpy.hstack((scalarParams, vecParams))
+
+ self.numLocs = 2
+ numElasticConsts = 21
+ self.density = numpy.array([densityA, densityB],
+ dtype=numpy.float64)
+
+ self.strain = numpy.array([strainA, strainB],
+ dtype=numpy.float64)
+ self.stress = numpy.zeros( (self.numLocs, 6), dtype=numpy.float64)
+ self.elasticConsts = numpy.zeros( (self.numLocs, numElasticConsts),
+ dtype=numpy.float64)
+
+ (self.elasticConsts[0,:], self.stress[0,:]) = \
+ self._calcStress(strainA, muA, lambdaA, maxwellTimeA, strainTA, visStrainA)
+ (self.elasticConsts[1,:], self.stress[1,:]) = \
+ self._calcStress(strainB, muB, lambdaB, maxwellTimeB, strainTB, visStrainB)
+ return
+
+
+ def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, strainTV, visStrainV):
+ """
+ Compute stress and derivative of elasticity matrix.
+ This assumes behavior is always viscoelastic.
+ """
+ import math
+
+ bulkModulus = lambdaV + 2.0 * muV/3.0
+
+ diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
+ traceStrainT = strainTV[0] + strainTV[1] + strainTV[2]
+ traceStrainTpdt = strainV[0] + strainV[1] + strainV[2]
+ meanStrainT = traceStrainT/3.0
+ meanStrainTpdt = traceStrainTpdt/3.0
+ meanStressTpdt = bulkModulus * traceStrainTpdt
+ timeFrac = 1.0e-5
+ numTerms = 5
+ dq = 0.0
+ if maxwellTimeV < timeFrac*self.dt:
+ fSign = 1.0
+ factorial = 1.0
+ fraction = 1.0
+ dq = 1.0
+ for iTerm in range(2, numTerms + 1):
+ factorial *= iTerm
+ fSign *= -1.0
+ fraction *= self.dt/maxwellTimeV
+ dq += fSign*fraction/factorial
+ else:
+ dq = maxwellTimeV*(1.0-math.exp(-self.dt/maxwellTimeV))/self.dt
+
+ visFac = muV*dq/3.0
+
+ C1111 = bulkModulus + 4.0*visFac
+ C1122 = bulkModulus - 2.0*visFac
+ C1133 = C1122
+ C1112 = 0.0
+ C1123 = 0.0
+ C1113 = 0.0
+ C2222 = C1111
+ C2233 = C1122
+ C2212 = 0.0
+ C2223 = 0.0
+ C2213 = 0.0
+ C3333 = C1111
+ C3312 = 0.0
+ C3323 = 0.0
+ C3313 = 0.0
+ C1212 = 3.0*visFac
+ C1223 = 0.0
+ C1213 = 0.0
+ C2323 = C1212
+ C2313 = 0.0
+ C1313 = C1212
+ elasticConsts = numpy.array([C1111, C1122, C1133, C1112, C1123, C1113,
+ C2222, C2233, C2212, C2223, C2213,
+ C3333, C3312, C3323, C3313,
+ C1212, C1223, C1213,
+ C2323, C2313,
+ C1313], dtype=numpy.float64)
+
+ stressV = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+ expFac = math.exp(-self.dt/maxwellTimeV)
+ elasFac = 2.0*muV
+ devStrainTpdt = 0.0
+ devStrainT = 0.0
+ devStressTpdt = 0.0
+ visStrain = 0.0
+ for iComp in range(6):
+ devStrainTpdt = strainV[iComp] - diag[iComp]*meanStrainTpdt
+ devStrainT = strainTV[iComp] - diag[iComp]*meanStrainT
+ visStrain = expFac*visStrainV[iComp] + dq*(devStrainTpdt - devStrainT)
+ devStressTpdt = elasFac*visStrain
+ stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt
+
+ stress = numpy.reshape(stressV, (6,1))
+ return (elasticConsts, numpy.ravel(stress))
+
+
+# MAIN /////////////////////////////////////////////////////////////////
+if __name__ == "__main__":
+
+ app = MaxwellIsotropic3DTimeDep()
+ app.run()
+
+
+# End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.cc 2007-06-05 03:47:23 UTC (rev 7065)
@@ -0,0 +1,198 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application maxwellisotropic3dtimedep.
+
+#include "MaxwellIsotropic3DTimeDepData.hh"
+
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_dimension = 3;
+
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numDBValues = 4;
+
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParameters = 6;
+
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numLocs = 2;
+
+const int pylith::materials::MaxwellIsotropic3DTimeDepData::_numParamValues[] = {
+1,
+1,
+1,
+1,
+6,
+6,
+};
+
+const char* pylith::materials::MaxwellIsotropic3DTimeDepData::_dbValues[] = {
+"density",
+"vs",
+"vp",
+"viscosity",
+};
+
+const char* pylith::materials::MaxwellIsotropic3DTimeDepData::_parameterNames[] = {
+"density",
+"mu",
+"lambda",
+"maxwellTime",
+"strainT",
+"visStrain",
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_dbData[] = {
+ 2.50000000e+03,
+ 3.00000000e+03,
+ 5.19615242e+03,
+ 1.00000000e+18,
+ 2.00000000e+03,
+ 1.20000000e+03,
+ 2.07846097e+03,
+ 1.00000000e+19,
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_parameterData[] = {
+ 2.50000000e+03,
+ 2.25000000e+10,
+ 2.25000000e+10,
+ 4.44444444e+07,
+ 1.10000000e-04,
+ 2.20000000e-04,
+ 3.30000000e-04,
+ 4.40000000e-04,
+ 5.50000000e-04,
+ 6.60000000e-04,
+ -2.20000000e-04,
+ -1.10000000e-04,
+ 0.00000000e+00,
+ 4.40000000e-04,
+ 5.50000000e-04,
+ 6.60000000e-04,
+ 2.00000000e+03,
+ 2.88000000e+09,
+ 2.88000000e+09,
+ 3.47222222e+09,
+ 1.20000000e-04,
+ 2.30000000e-04,
+ 3.40000000e-04,
+ 4.50000000e-04,
+ 5.60000000e-04,
+ 6.70000000e-04,
+ -2.20000000e-04,
+ -1.10000000e-04,
+ 0.00000000e+00,
+ 4.50000000e-04,
+ 5.60000000e-04,
+ 6.70000000e-04,
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_density[] = {
+ 2.50000000e+03,
+ 2.00000000e+03,
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_strain[] = {
+ 1.10000000e-04,
+ 2.20000000e-04,
+ 3.30000000e-04,
+ 4.40000000e-04,
+ 5.50000000e-04,
+ 6.60000000e-04,
+ 1.20000000e-04,
+ 2.30000000e-04,
+ 3.40000000e-04,
+ 4.50000000e-04,
+ 5.60000000e-04,
+ 6.70000000e-04,
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_stress[] = {
+ 1.48944499e+07,
+ 1.98222250e+07,
+ 2.47500000e+07,
+ 1.97111002e+07,
+ 2.46388752e+07,
+ 2.95666503e+07,
+ 2.04487299e+06,
+ 2.67843649e+06,
+ 3.31200000e+06,
+ 2.59185071e+06,
+ 3.22541421e+06,
+ 3.85897772e+06,
+};
+
+const double pylith::materials::MaxwellIsotropic3DTimeDepData::_elasticConsts[] = {
+ 6.74326011e+10,
+ 2.25336994e+10,
+ 2.25336994e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 6.74326011e+10,
+ 2.25336994e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 6.74326011e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.24494509e+10,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.24494509e+10,
+ 0.00000000e+00,
+ 2.24494509e+10,
+ 8.63988941e+09,
+ 2.88005529e+09,
+ 2.88005529e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 8.63988941e+09,
+ 2.88005529e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 8.63988941e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.87991706e+09,
+ 0.00000000e+00,
+ 0.00000000e+00,
+ 2.87991706e+09,
+ 0.00000000e+00,
+ 2.87991706e+09,
+};
+
+pylith::materials::MaxwellIsotropic3DTimeDepData::MaxwellIsotropic3DTimeDepData(void)
+{ // constructor
+ dimension = _dimension;
+ numDBValues = _numDBValues;
+ numParameters = _numParameters;
+ numLocs = _numLocs;
+ numParamValues = const_cast<int*>(_numParamValues);
+ dbValues = const_cast<char**>(_dbValues);
+ parameterNames = const_cast<char**>(_parameterNames);
+ dbData = const_cast<double*>(_dbData);
+ parameterData = const_cast<double*>(_parameterData);
+ density = const_cast<double*>(_density);
+ strain = const_cast<double*>(_strain);
+ stress = const_cast<double*>(_stress);
+ elasticConsts = const_cast<double*>(_elasticConsts);
+} // constructor
+
+pylith::materials::MaxwellIsotropic3DTimeDepData::~MaxwellIsotropic3DTimeDepData(void)
+{}
+
+
+// End of file
Added: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3DTimeDepData.hh 2007-06-05 03:47:23 UTC (rev 7065)
@@ -0,0 +1,70 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// DO NOT EDIT THIS FILE
+// This file was generated from python application maxwellisotropic3dtimedep.
+
+#if !defined(pylith_materials_maxwellisotropic3dtimedepdata_hh)
+#define pylith_materials_maxwellisotropic3dtimedepdata_hh
+
+#include "ElasticMaterialData.hh"
+
+namespace pylith {
+ namespace materials {
+ class MaxwellIsotropic3DTimeDepData;
+ } // pylith
+} // materials
+
+class pylith::materials::MaxwellIsotropic3DTimeDepData : public ElasticMaterialData
+{
+
+public:
+
+ /// Constructor
+ MaxwellIsotropic3DTimeDepData(void);
+
+ /// Destructor
+ ~MaxwellIsotropic3DTimeDepData(void);
+
+private:
+
+ static const int _dimension;
+
+ static const int _numDBValues;
+
+ static const int _numParameters;
+
+ static const int _numLocs;
+
+ static const int _numParamValues[];
+
+ static const char* _dbValues[];
+
+ static const char* _parameterNames[];
+
+ static const double _dbData[];
+
+ static const double _parameterData[];
+
+ static const double _density[];
+
+ static const double _strain[];
+
+ static const double _stress[];
+
+ static const double _elasticConsts[];
+
+};
+
+#endif // pylith_materials_maxwellisotropic3dtimedepdata_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh 2007-06-05 00:26:06 UTC (rev 7064)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/generate.sh 2007-06-05 03:47:23 UTC (rev 7065)
@@ -25,6 +25,11 @@
--data.object=ElasticIsotropic3DData \
--data.parent=ElasticMaterialData
+ python MaxwellIsotropic3DElastic.py \
+ --data.namespace=pylith,materials \
+ --data.object=MaxwellIsotropic3DElasticData \
+ --data.parent=ElasticMaterialData
+
# 2-D ----------------------------------------------------------------
python ElasticPlaneStrain.py \
@@ -53,8 +58,12 @@
# //////////////////////////////////////////////////////////////////////
if [ $1 == "viscoelastic" ] || [ $1 == "all" ]; then
- echo "" >& /dev/null
+ python MaxwellIsotropic3DTimeDep.py \
+ --data.namespace=pylith,materials \
+ --data.object=MaxwellIsotropic3DTimeDepData \
+ --data.parent=ElasticMaterialData
+
fi
More information about the cig-commits
mailing list