[cig-commits] r7060 - in
short/3D/PyLith/trunk/unittests/libtests/materials: . data
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon Jun 4 13:44:55 PDT 2007
Author: willic3
Date: 2007-06-04 13:44:55 -0700 (Mon, 04 Jun 2007)
New Revision: 7060
Modified:
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py
Log:
More work on Maxwell unit tests (not yet complete).
Need separate versions of the python for elastic and time-dependent.
Need to check usage of numpy for stress.
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2007-06-04 20:22:56 UTC (rev 7059)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2007-06-04 20:44:55 UTC (rev 7060)
@@ -14,7 +14,7 @@
#include "TestMaxwellIsotropic3D.hh" // Implementation of class methods
-#include "data/MaxwellIsotropic3DData.hh" // USES MaxwellIsotropic3DData
+#include "data/MaxwellIsotropic3DElasticData.hh" // USES MaxwellIsotropic3DElasticData
#include "pylith/materials/MaxwellIsotropic3D.hh" // USES MaxwellIsotropic3D
@@ -27,7 +27,9 @@
pylith::materials::TestMaxwellIsotropic3D::testDBValues(void)
{ // testDBValues
MaxwellIsotropic3D material;
- MaxwellIsotropic3DData data;
+ MaxwellIsotropic3DElasticData data;
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
_testDBValues(&material, data);
} // testDBValues
@@ -37,7 +39,9 @@
pylith::materials::TestMaxwellIsotropic3D::testParameters(void)
{ // testParameters
MaxwellIsotropic3D material;
- MaxwellIsotropic3DData data;
+ MaxwellIsotropic3DElasticData data;
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
_testParameters(&material, data);
} // testParameters
@@ -47,7 +51,9 @@
pylith::materials::TestMaxwellIsotropic3D::testDBToParameters(void)
{ // testDBToParameters
MaxwellIsotropic3D material;
- MaxwellIsotropic3DData data;
+ MaxwellIsotropic3DElasticData data;
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
_testDBToParameters(&material, data);
} // testDBToParameters
@@ -57,42 +63,93 @@
pylith::materials::TestMaxwellIsotropic3D::testCalcDensity(void)
{ // testCalcDensity
MaxwellIsotropic3D material;
- MaxwellIsotropic3DData data;
+ MaxwellIsotropic3DElasticData data;
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
_testCalcDensity(&material, data);
} // testCalcDensity
// ----------------------------------------------------------------------
-// Test calcStress()
+// Test calcStressElastic()
void
-pylith::materials::TestMaxwellIsotropic3D::testCalcStress(void)
-{ // testCalcStress
+pylith::materials::TestMaxwellIsotropic3D::testCalcStressElastic(void)
+{ // testCalcStressElastic
MaxwellIsotropic3D material;
- MaxwellIsotropic3DData data;
+ MaxwellIsotropic3DElasticData data;
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
_testCalcStress(&material, data);
-} // testCalcStress
+} // testCalcStressElastic
// ----------------------------------------------------------------------
-// Test calcElasticConsts()
+// Test calcElasticConstsElastic()
void
-pylith::materials::TestMaxwellIsotropic3D::testCalcElasticConsts(void)
-{ // testElasticConsts
+pylith::materials::TestMaxwellIsotropic3D::testCalcElasticConstsElastic(void)
+{ // testElasticConstsElastic
MaxwellIsotropic3D material;
- MaxwellIsotropic3DData data;
- _testCalcElasticConsts(&material, data);
-} // testElasticConsts
+ MaxwellIsotropic3DElasticData data;
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
+ _testCalcElasticConstsElastic(&material, data);
+} // testElasticConstsElastic
// ----------------------------------------------------------------------
-// Test updateState()
+// Test updateStateElastic()
void
-pylith::materials::TestMaxwellIsotropic3D::testUpdateState(void)
-{ // testUpdateState
+pylith::materials::TestMaxwellIsotropic3D::testUpdateStateElastic(void)
+{ // testUpdateStateElastic
MaxwellIsotropic3D material;
std::vector<double_array> totalStrain;
- material.updateState(totalStrain);
-} // testUpdateState
+ bool elasFlag = true;
+ material.useElasticBehavior(elasFlag);
+ material.updateStateElastic(totalStrain);
+} // testUpdateStateElastic
// ----------------------------------------------------------------------
+// Test calcStressTimeDep()
+void
+pylith::materials::TestMaxwellIsotropic3D::testCalcStressTimeDep(void)
+{ // testCalcStressTimeDep
+ MaxwellIsotropic3D material;
+ MaxwellIsotropic3DTimeDepData data;
+ bool elasFlag = false;
+ material.useElasticBehavior(elasFlag);
+ double dt = 2.0e5;
+ material.timeStep(dt);
+ _testCalcStressTimeDep(&material, data);
+} // testCalcStressTimeDep
+
+// ----------------------------------------------------------------------
+// Test calcElasticConstsTimeDep()
+void
+pylith::materials::TestMaxwellIsotropic3D::testCalcElasticConstsTimeDep(void)
+{ // testElasticConstsTimeDep
+ MaxwellIsotropic3D material;
+ MaxwellIsotropic3DTimeDepData data;
+ bool elasFlag = false;
+ material.useElasticBehavior(elasFlag);
+ double dt = 2.0e5;
+ material.timeStep(dt);
+ _testCalcElasticConstsTimeDep(&material, data);
+} // testElasticConstsTimeDep
+
+// ----------------------------------------------------------------------
+// Test updateStateTimeDep()
+void
+pylith::materials::TestMaxwellIsotropic3D::testUpdateStateTimeDep(void)
+{ // testUpdateStateTimeDep
+ MaxwellIsotropic3D material;
+
+ std::vector<double_array> totalStrain;
+ bool elasFlag = false;
+ material.useElasticBehavior(elasFlag);
+ double dt = 2.0e5;
+ material.timeStep(dt);
+ material.updateStateTimeDep(totalStrain);
+} // testUpdateStateTimeDep
+
+// ----------------------------------------------------------------------
// Test timeStep()
void
pylith::materials::TestMaxwellIsotropic3D::testTimeStep(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py 2007-06-04 20:22:56 UTC (rev 7059)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/MaxwellIsotropic3D.py 2007-06-04 20:44:55 UTC (rev 7060)
@@ -43,31 +43,44 @@
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
- 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]
+
+ # 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)
@@ -84,37 +97,63 @@
dtype=numpy.float64)
(self.elasticConsts[0,:], self.stress[0,:]) = \
- self._calcStress(strainA, densityA, muA, lambdaA)
+ self._calcStress(strainA, muA, lambdaA, maxwellTimeA, strainTA, visStrainA)
(self.elasticConsts[1,:], self.stress[1,:]) = \
- self._calcStress(strainB, densityB, muB, lambdaB)
+ self._calcStress(strainB, muB, lambdaB, maxwellTimeB, strainTB, visStrainB)
return
- def _calcStress(self, strainV, densityV, muV, lambdaV):
+ def _calcStress(self, strainV, muV, lambdaV, maxwellTimeV, strainTV, visStrainV):
"""
Compute stress and derivative of elasticity matrix.
+ This assumes behavior is always viscoelastic.
"""
- C1111 = lambdaV + 2.0*muV
- C1122 = lambdaV
- C1133 = lambdaV
+ 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 = lambdaV + 2.0*muV
- C2233 = lambdaV
+ C2222 = C1111
+ C2233 = C1122
C2212 = 0.0
C2223 = 0.0
C2213 = 0.0
- C3333 = lambdaV + 2.0*muV
+ C3333 = C1111
C3312 = 0.0
C3323 = 0.0
C3313 = 0.0
- C1212 = 2.0*muV
+ C1212 = 3.0*visFac
C1223 = 0.0
C1213 = 0.0
- C2323 = 2.0*muV
+ C2323 = C1212
C2313 = 0.0
- C1313 = 2.0*muV
+ C1313 = C1212
elasticConsts = numpy.array([C1111, C1122, C1133, C1112, C1123, C1113,
C2222, C2233, C2212, C2223, C2213,
C3333, C3312, C3323, C3313,
@@ -122,15 +161,21 @@
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)
+ 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))
More information about the cig-commits
mailing list