[cig-commits] r7299 - in short/3D/PyLith/trunk: libsrc/materials
unittests/libtests/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Jun 19 07:31:16 PDT 2007
Author: willic3
Date: 2007-06-19 07:31:16 -0700 (Tue, 19 Jun 2007)
New Revision: 7299
Modified:
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
Log:
Worked on unit test for viscoelastic updateState, but there is still a
problem -- the expected value is the negative of the actual value.
Some debugging will be needed.
Removed some unused variables from MaxwellIsotropic3D.cc.
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-06-19 07:06:17 UTC (rev 7298)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-06-19 14:31:16 UTC (rev 7299)
@@ -267,6 +267,13 @@
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
const double meanStressTpdt = bulkmodulus * traceStrainTpdt;
+ // See what's going on in state variables.
+ std::cout << " pidStrainT, pidVisStrain : " << std::endl;
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
+ std::cout << " " << parameters[_MaxwellIsotropic3D::pidStrainT][iComp]
+ << " " << parameters[_MaxwellIsotropic3D::pidVisStrain][iComp]
+ << std::endl;
+
const double meanStrainT = (parameters[_MaxwellIsotropic3D::pidStrainT][0] +
parameters[_MaxwellIsotropic3D::pidStrainT][1] +
parameters[_MaxwellIsotropic3D::pidStrainT][2])/3.0;
@@ -307,6 +314,13 @@
// Later I will want to put in initial stresses.
(*stress)[iComp] =diag[iComp]*meanStressTpdt+devStressTpdt;
} // for
+ std::cout << " totalStrain: " << std::endl;
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
+ std::cout << " " << totalStrain[iComp];
+ std::cout << std::endl << " stress: " << std::endl;
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
+ std::cout << " " << (*stress)[iComp];
+ std::cout << std::endl;
} // _calcStress
// ----------------------------------------------------------------------
@@ -433,21 +447,14 @@
assert(6 == (*parameters)[_MaxwellIsotropic3D::pidStrainT].size());
assert(6 == (*parameters)[_MaxwellIsotropic3D::pidVisStrain].size());
- const double density = (*parameters)[_MaxwellIsotropic3D::pidDensity][0];
- const double mu = (*parameters)[_MaxwellIsotropic3D::pidMu][0];
- const double lambda = (*parameters)[_MaxwellIsotropic3D::pidLambda][0];
const double maxwelltime = (*parameters)[_MaxwellIsotropic3D::pidMaxwellTime][0];
- const double lambda2mu = lambda + 2.0 * mu;
- const double bulkmodulus = lambda + 2.0 * mu/3.0;
-
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
const double traceStrainTpdt = e11 + e22 + e33;
const double meanStrainTpdt = traceStrainTpdt/3.0;
- const double s123 = lambda * traceStrainTpdt;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
@@ -456,6 +463,12 @@
(*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp] =
totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
} // for
+ std::cout << " pidStrainT, pidVisStrain : " << std::endl;
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
+ std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
+ << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
+ << std::endl;
+ _needNewJacobian = true;
} // _calcStressElastic
// ----------------------------------------------------------------------
@@ -475,21 +488,13 @@
assert(6 == (*parameters)[_MaxwellIsotropic3D::pidStrainT].size());
assert(6 == (*parameters)[_MaxwellIsotropic3D::pidVisStrain].size());
- const double density = (*parameters)[_MaxwellIsotropic3D::pidDensity][0];
- const double mu = (*parameters)[_MaxwellIsotropic3D::pidMu][0];
- const double lambda = (*parameters)[_MaxwellIsotropic3D::pidLambda][0];
const double maxwelltime = (*parameters)[_MaxwellIsotropic3D::pidMaxwellTime][0];
- const double lambda2mu = lambda + 2.0 * mu;
- const double bulkmodulus = lambda + 2.0 * mu/3.0;
-
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
- const double traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
- const double s123 = lambda * traceStrainTpdt;
+ const double meanStrainTpdt = (e11 + e22 + e33)/3.0;
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2007-06-19 07:06:17 UTC (rev 7298)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaxwellIsotropic3D.cc 2007-06-19 14:31:16 UTC (rev 7299)
@@ -226,34 +226,67 @@
MaxwellIsotropic3D material;
MaxwellIsotropic3DTimeDepData data;
- // CHARLES: Setup parameters and totalStrain so this works.
-#if 0
const int numParams = data.numParameters;
+ material.useElasticBehavior(false);
+ const double dt = 2.0e5;
+ material.timeStep(dt);
+ const double viscosity = 1.0e18;
+ const double mu = 3.0e10;
+ const double maxwelltime = viscosity/mu;
+
+ const int tensorSize = 6;
+ double_array totalStrainTpdt(tensorSize);
+ double_array totalStrainT(tensorSize);
+ double_array visStrainT(tensorSize);
+ for (int i=0; i < tensorSize; ++i) {
+ totalStrainTpdt[i] = i;
+ totalStrainT[i] = totalStrainTpdt[i]/2.0;
+ visStrainT[i] = totalStrainTpdt[i]/4.0;
+ } // for
+
+ const double meanStrainTpdt = (totalStrainTpdt[0] + totalStrainTpdt[1] + totalStrainTpdt[2])/3.0;
+ const double meanStrainT = (totalStrainT[0] + totalStrainT[1] + totalStrainT[2])/3.0;
+
+ const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+
std::vector<double_array> parameters(numParams);
- const int paramsSize = 1;
+ std::vector<double_array> paramdata(numParams);
+ const int paramsSize []= { 1, 1, 1, 1, 6, 6};
for (int i=0; i < numParams; ++i) {
- parameters[i].resize(numParams);
- for (int j=0; j < paramsSize; ++j)
+ parameters[i].resize(paramsSize[i]);
+ paramdata[i].resize(paramsSize[i]);
+ for (int j=0; j < paramsSize[i]; ++j)
parameters[i][j] = i+j;
} // for
-
- const int tensorSize = 9;
- double_array totalStrain(tensorSize);
- for (int i=0; i < tensorSize; ++i)
- totalStrain[i] = i;
+
+ parameters[3][0] = maxwelltime;
+ paramdata[3][0] = maxwelltime;
+
+ const double dq = maxwelltime*(1.0-exp(-dt/maxwelltime))/dt;
+ const double expFac = exp(-dt/maxwelltime);
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+
+ for (int i=0; i < tensorSize; ++i) {
+ devStrainTpdt = totalStrainTpdt[i] - diag[i]*meanStrainTpdt;
+ devStrainT = totalStrainT[i] - diag[i]*meanStrainT;
+ parameters[4][i] = totalStrainT[i];
+ parameters[5][i] = visStrainT[i];
+ paramdata[4][i] = totalStrainTpdt[i];
+ paramdata[5][i] = expFac * visStrainT[i] + dq * (devStrainTpdt - devStrainT);
+ } //for
- material._updateState(¶meters, totalStrain);
+ material._updateState(¶meters, totalStrainTpdt);
const double tolerance = 1.0e-06;
- for (int i=0; i < numParams; ++i)
- for (int j=0; j < paramsSize; ++j)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i+j), parameters[i][j], tolerance);
+ // Test vector parameters and Maxwell time.
+ for (int i=3; i < numParams; ++i)
+ for (int j=0; j < paramsSize[i]; ++j)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(paramdata[i][j], parameters[i][j], tolerance);
for (int i=0; i < tensorSize; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrain[i], tolerance);
-#endif
- throw std::logic_error("Unit test not implemented.");
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(double(i), totalStrainTpdt[i], tolerance);
} // testUpdateStateTimeDep
More information about the cig-commits
mailing list