[cig-commits] r15229 - in short/3D/PyLith/trunk: libsrc/materials unittests/libtests/materials unittests/libtests/materials/data

brad at geodynamics.org brad at geodynamics.org
Sat Jun 13 14:11:46 PDT 2009


Author: brad
Date: 2009-06-13 14:11:46 -0700 (Sat, 13 Jun 2009)
New Revision: 15229

Modified:
   short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py
   short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc
   short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh
Log:
Updated unit tests and fixed some small bugs in update to GenMaxwellIsotropic3D. Still need to update state variables test.

Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc	2009-06-13 05:52:06 UTC (rev 15228)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc	2009-06-13 21:11:46 UTC (rev 15229)
@@ -194,7 +194,7 @@
   _updateStateVarsFn(0)  
 { // constructor
   useElasticBehavior(true);
-  _viscousStrain.resize(_tensorSize);
+  _viscousStrain.resize(_GenMaxwellIsotropic3D::numMaxwellModels*_tensorSize);
 } // constructor
 
 // ----------------------------------------------------------------------
@@ -512,7 +512,7 @@
   const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
   
   double visFrac = 0.0;
-  for (int imodel = 0; imodel < numMaxwellModels; ++imodel) 
+  for (int imodel=0; imodel < numMaxwellModels; ++imodel) 
     visFrac += muRatio[imodel];
   assert(visFrac <= 1.0);
   const double elasFrac = 1.0 - visFrac;
@@ -529,8 +529,6 @@
   } else {
     int index = 0;
     for (int iComp=0; iComp < tensorSize; ++iComp)
-      _viscousStrain[index++] = stateVars[s_totalStrain+iComp];
-    for (int iComp=0; iComp < tensorSize; ++iComp)
       _viscousStrain[index++] = stateVars[s_viscousStrain1+iComp];
     for (int iComp=0; iComp < tensorSize; ++iComp)
       _viscousStrain[index++] = stateVars[s_viscousStrain2+iComp];
@@ -545,14 +543,11 @@
     devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
     devStressTpdt = elasFrac*devStrainTpdt;
     for (int model=0; model < numMaxwellModels; ++model) {
-      if (muRatio[model] != 0.0) {
-	int pindex = iComp + model * tensorSize;
-	devStressTpdt += muRatio[model] * _viscousStrain[pindex];
-      } // if
+      devStressTpdt += muRatio[model] * _viscousStrain[model*tensorSize+iComp];
     } // for
 
     devStressTpdt = mu2*devStressTpdt;
-    stress[iComp] =diag[iComp] * meanStressTpdt + devStressTpdt +
+    stress[iComp] = diag[iComp] * meanStressTpdt + devStressTpdt +
 	    initialStress[iComp];
   } // for
 
@@ -922,9 +917,8 @@
 
     // Maxwell model 1
     int imodel = 0;
-    int vindex = imodel*tensorSize + iComp;
     if (0.0 != muRatio[imodel]) {
-      _viscousStrain[vindex] = exp(-_dt/maxwellTime[imodel])*
+      _viscousStrain[imodel*tensorSize+iComp] = exp(-_dt/maxwellTime[imodel])*
 	stateVars[s_viscousStrain1 + iComp] + dq[imodel] * deltaStrain;
       PetscLogFlops(6);
     } // if
@@ -932,7 +926,7 @@
     // Maxwell model 2
     imodel = 1;
     if (0.0 != muRatio[imodel]) {
-      _viscousStrain[vindex] = exp(-_dt/maxwellTime[imodel])*
+      _viscousStrain[imodel*tensorSize+iComp] = exp(-_dt/maxwellTime[imodel])*
 	stateVars[s_viscousStrain2 + iComp] + dq[imodel] * deltaStrain;
       PetscLogFlops(6);
     } // if
@@ -940,7 +934,7 @@
     // Maxwell model 3
     imodel = 2;
     if (0.0 != muRatio[imodel]) {
-      _viscousStrain[vindex] = exp(-_dt/maxwellTime[imodel])*
+      _viscousStrain[imodel*tensorSize+iComp] = exp(-_dt/maxwellTime[imodel])*
 	stateVars[s_viscousStrain3 + iComp] + dq[imodel] * deltaStrain;
       PetscLogFlops(6);
     } // if

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am	2009-06-13 05:52:06 UTC (rev 15228)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/Makefile.am	2009-06-13 21:11:46 UTC (rev 15229)
@@ -61,10 +61,10 @@
 	data/MaxwellIsotropic3DElasticData.cc \
 	data/MaxwellIsotropic3DTimeDepData.cc \
 	data/GenMaxwellIsotropic3DElasticData.cc \
+	data/GenMaxwellIsotropic3DTimeDepData.cc \
 	data/PowerLaw3DElasticData.cc \
 	data/PowerLaw3DTimeDepData.cc
 
-#	data/GenMaxwellIsotropic3DTimeDepData.cc 
 
 noinst_HEADERS += \
 	data/MaterialData.hh \

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc	2009-06-13 05:52:06 UTC (rev 15228)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestGenMaxwellIsotropic3D.cc	2009-06-13 21:11:46 UTC (rev 15229)
@@ -185,7 +185,6 @@
 void
 pylith::materials::TestGenMaxwellIsotropic3D::test_calcStressTimeDep(void)
 { // testCalcStressTimeDep
-#if 0
   CPPUNIT_ASSERT(0 != _matElastic);
   _matElastic->useElasticBehavior(false);
 
@@ -194,7 +193,6 @@
   double dt = 2.0e+5;
   _matElastic->timeStep(dt);
   test_calcStress();
-#endif
 } // testCalcStressTimeDep
 
 // ----------------------------------------------------------------------
@@ -202,7 +200,6 @@
 void
 pylith::materials::TestGenMaxwellIsotropic3D::test_calcElasticConstsTimeDep(void)
 { // testElasticConstsTimeDep
-#if 0
   CPPUNIT_ASSERT(0 != _matElastic);
   _matElastic->useElasticBehavior(false);
 
@@ -211,7 +208,6 @@
   double dt = 2.0e+5;
   _matElastic->timeStep(dt);
   test_calcElasticConsts();
-#endif
 } // testElasticConstsTimeDep
 
 // ----------------------------------------------------------------------
@@ -410,13 +406,11 @@
 void
 pylith::materials::TestGenMaxwellIsotropic3D::test_stableTimeStepImplicit(void)
 { // test_stableTimeStepImplicit
-#if 0
   CPPUNIT_ASSERT(0 != _matElastic);
 
   delete _dataElastic; _dataElastic = new GenMaxwellIsotropic3DTimeDepData();
 
   TestElasticMaterial::test_stableTimeStepImplicit();
-#endif
 } // test_stableTimeStepImplicit
 
 

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py	2009-06-13 05:52:06 UTC (rev 15228)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDep.py	2009-06-13 21:11:46 UTC (rev 15229)
@@ -20,6 +20,11 @@
 import numpy
 
 # ----------------------------------------------------------------------
+dimension = 3
+numElasticConsts = 21
+tensorSize = 6
+numMaxwellModels = 3
+
 # GenMaxwellIsotropic3DTimeDep class
 class GenMaxwellIsotropic3DTimeDep(ElasticMaterialApp):
   """
@@ -36,126 +41,162 @@
     """
     ElasticMaterialApp.__init__(self, name)
 
-    # import pdb
-    # pdb.set_trace()
+    numLocs = 2
+    self.dt = 2.0e5
+    self.dimension = dimension
+    self.numLocs = numLocs
 
-    self.dimension = 3
-    self.tensorSize = 6
+    self.dbPropertyValues = ["density", "vs", "vp",
+                             "shear_ratio_1", "shear_ratio_2", "shear_ratio_3",
+                             "viscosity_1", "viscosity_2", "viscosity_3"]
+    self.numPropertyValues = numpy.array([1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=numpy.int32)
 
-    self.numMaxwellModels = 3
-    self.numDBValues = 3+2*self.numMaxwellModels
-    self.numInitialStateValues = self.tensorSize
-    self.dbValues = ["density", "vs", "vp" ,
-                     "shear_ratio_1", "shear_ratio_2", "shear_ratio_3",
-                     "viscosity_1", "viscosity_2", "viscosity_3"]
-    self.initialStateDBValues =  ["stress_xx", "stress_yy", "stress_zz",
-                                  "stress_xy", "stress_yz", "stress_xy"]
-    self.numParameters = 7
-    self.numParamValues = [1, 1, 1,
-                           self.numMaxwellModels, self.numMaxwellModels,
-                           self.tensorSize,
-                           self.tensorSize*self.numMaxwellModels]
-    self.parameterNames = ["density", "mu", "lambda",
-                           "shearRatio", "maxwellTime",
-                           "strainT", "visStrain"]
+    self.dbStateVarValues = ["total-strain-xx",
+                             "total-strain-yy",
+                             "total-strain-zz",
+                             "total-strain-xy",
+                             "total-strain-yz",
+                             "total-strain-xz",
+                             "viscous-strain-1-xx",
+                             "viscous-strain-1-yy",
+                             "viscous-strain-1-zz",
+                             "viscous-strain-1-xy",
+                             "viscous-strain-1-yz",
+                             "viscous-strain-1-xz",
+                             "viscous-strain-2-xx",
+                             "viscous-strain-2-yy",
+                             "viscous-strain-2-zz",
+                             "viscous-strain-2-xy",
+                             "viscous-strain-2-yz",
+                             "viscous-strain-2-xz",
+                             "viscous-strain-3-xx",
+                             "viscous-strain-3-yy",
+                             "viscous-strain-3-zz",
+                             "viscous-strain-3-xy",
+                             "viscous-strain-3-yz",
+                             "viscous-strain-3-xz",
+                             ]
+    self.numStateVarValues = numpy.array([tensorSize]*(1+numMaxwellModels),
+                                         dtype=numpy.int32)
 
-    self.dt = 2.0e5
-
     densityA = 2500.0
     vsA = 3000.0
     vpA = vsA*3**0.5
     shearRatioA = [0.5, 0.1, 0.2]
     viscosityA = [1.0e18, 1.0e17, 1.0e19]
-    visStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
     strainA = [1.1e-4, 2.2e-4, 3.3e-4, 4.4e-4, 5.5e-4, 6.6e-4]
-    meanStrainA = (strainA[0] + strainA[1] + strainA[2])/3.0
-    elasDataA = [densityA, vsA, vpA]
-    matDbA = elasDataA + shearRatioA + viscosityA
+    initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
+    #initialStrainA = [3.1e-4, 3.2e-4, 3.3e-4, 3.4e-4, 3.5e-4, 3.6e-4]
     muA = vsA*vsA*densityA
     lambdaA = vpA*vpA*densityA - 2.0*muA
-    initialStateA = [1.2e4, 2.3e4, 3.4e4, 4.5e4, 5.6e4, 6.7e4]
-    
+
     densityB = 2000.0
     vsB = 1200.0
     vpB = vsB*3**0.5
     shearRatioB = [0.2, 0.2, 0.2]
     viscosityB = [1.0e18, 1.0e19, 1.0e20]
-    visStrainB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
     strainB = [1.2e-4, 2.3e-4, 3.4e-4, 4.5e-4, 5.6e-4, 6.7e-4]
-    meanStrainB = (strainB[0] + strainB[1] + strainB[2])/3.0
-    elasDataB = [densityB, vsB, vpB]
-    matDbB = elasDataB + shearRatioB + viscosityB
+    initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
+    #initialStrainB = [6.1e-4, 6.2e-4, 6.3e-4, 6.4e-4, 6.5e-4, 6.6e-4]
     muB = vsB*vsB*densityB
     lambdaB = vpB*vpB*densityB - 2.0*muB
-    initialStateB = [2.1e4, 3.2e4, 4.3e4, 5.4e4, 6.5e4, 7.6e4]
 
-    matDb = matDbA + matDbB
-    self.dbData = numpy.array(matDb, dtype=numpy.float64)
-
-    diag = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
-
+    # TEMPORARY, need to determine how to use initial strain
+    initialStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    initialStrainB = [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[:]
+    diag = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0], dtype=numpy.float64)
+    strainTA = numpy.array(strainA)
+    strainTB = numpy.array(strainB)
+    meanStrainA = (strainA[0] + strainA[1] + strainA[2])/3.0
+    meanStrainB = (strainB[0] + strainB[1] + strainB[2])/3.0
 
     maxwellTimeA = [0.0, 0.0, 0.0]
     maxwellTimeB = [0.0, 0.0, 0.0]
-    for model in range(self.numMaxwellModels):
-      if shearRatioA[model] != 0.0:
-        maxwellTimeA[model] = viscosityA[model]/(muA*shearRatioA[model])
-        for i in range(self.tensorSize):
-          visStrainA[i+self.tensorSize*model] = strainA[i] - diag[i] * meanStrainA
-      if shearRatioB[model] != 0.0:
-        maxwellTimeB[model] = viscosityB[model]/(muB*shearRatioB[model])
-        for i in range(self.tensorSize):
-          visStrainB[i+self.tensorSize*model] = strainB[i] - diag[i] * meanStrainB
-        
-    vecParamsA = numpy.hstack((shearRatioA, maxwellTimeA, strainTA, visStrainA))
-    vecParamsB = numpy.hstack((shearRatioB, maxwellTimeB, strainTB, visStrainB))
-    vecParams = numpy.vstack((vecParamsA, vecParamsB))
-    scalarParams = numpy.array([ [densityA, muA, lambdaA],
-                                 [densityB, muB, lambdaB] ],
+    visStrainA = numpy.zeros( (numMaxwellModels, tensorSize), dtype=numpy.float64)
+    visStrainB = numpy.zeros( (numMaxwellModels, tensorSize), dtype=numpy.float64)
+    for imodel in xrange(numMaxwellModels):
+      if shearRatioA[imodel] != 0.0:
+        maxwellTimeA[imodel] = viscosityA[imodel]/(muA*shearRatioA[imodel])
+        visStrainA[imodel,:] = strainA[:] - diag[:] * meanStrainA
+      if shearRatioB[imodel] != 0.0:
+        maxwellTimeB[imodel] = viscosityB[imodel]/(muB*shearRatioB[imodel])
+        visStrainB[imodel,:] = strainB[:] - diag[:] * meanStrainB
+
+    self.lengthScale = 1.0e+3
+    self.pressureScale = muA
+    self.timeScale = 1.0
+    self.densityScale = 1.0e+3
+
+    propA = [densityA, vsA, vpA] + shearRatioA + viscosityA
+    propB = [densityB, vsB, vpB] + shearRatioB + viscosityB
+    self.dbProperties = numpy.array([propA, propB], dtype=numpy.float64)
+    propA = [densityA, muA, lambdaA] + shearRatioA + maxwellTimeA
+    propB = [densityB, muB, lambdaB] + shearRatioB + maxwellTimeB
+    self.properties = numpy.array([propA, propB], dtype=numpy.float64)
+
+    # TEMPORARY, need to determine how to use initial state variables
+    self.dbStateVars = numpy.zeros( (numLocs, tensorSize+numMaxwellModels*tensorSize),
+                                    dtype=numpy.float64)
+    self.stateVars = numpy.zeros( (numLocs, tensorSize+numMaxwellModels*tensorSize),
+                                  dtype=numpy.float64)
+    self.stateVars[0,0:tensorSize] = strainTA
+    self.stateVars[0,tensorSize:(1+numMaxwellModels)*tensorSize] = visStrainA.ravel()
+    self.stateVars[1,0:tensorSize] = strainTB
+    self.stateVars[1,tensorSize:(1+numMaxwellModels)*tensorSize] = visStrainB.ravel()
+
+    mu0 = self.pressureScale
+    density0 = self.densityScale
+    time0 = self.timeScale
+    self.propertiesNondim = \
+        numpy.array([ [densityA/density0, muA/mu0, lambdaA/mu0,
+                       shearRatioA[0], shearRatioA[1], shearRatioA[2],
+                       maxwellTimeA[0]/time0, maxwellTimeA[1]/time0, maxwellTimeA[2]/time0],
+                      [densityB/density0, muB/mu0, lambdaB/mu0,
+                       shearRatioB[0], shearRatioB[1], shearRatioB[2],
+                       maxwellTimeB[0]/time0, maxwellTimeB[1]/time0, maxwellTimeB[2]/time0] ],
+                    dtype=numpy.float64)
+
+    self.stateVarsNondim = self.stateVars # no scaling
+
+    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.parameterData = numpy.hstack((scalarParams, vecParams))
+
+    self.strain = numpy.array([strainA,
+                               strainB],
+                               dtype=numpy.float64)
     
-    self.numLocs = 2
-    numElasticConsts = 21
-    self.density = numpy.array([densityA, densityB], dtype=numpy.float64)
+    self.stress = numpy.zeros( (numLocs, tensorSize), dtype=numpy.float64)
+    self.elasticConsts = numpy.zeros( (numLocs, numElasticConsts), \
+                                        dtype=numpy.float64)
 
-    self.strain = numpy.array([strainA, strainB], dtype=numpy.float64)
-    self.initialStateDBData = numpy.array([initialStateA, initialStateB],
-                                          dtype=numpy.float64)
-    self.initialState = numpy.array([initialStateA, initialStateB],
-                                    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,
-                                               shearRatioA, maxwellTimeA,
-                                               strainTA, visStrainA,
-                                               initialStateA)
+        self._calcStress(strainA, muA, lambdaA, shearRatioA, maxwellTimeA,
+                         strainTA, visStrainA,
+                         initialStressA, initialStrainA)
     (self.elasticConsts[1,:], self.stress[1,:]) = \
-                              self._calcStress(strainB, muB, lambdaB,
-                                               shearRatioB, maxwellTimeB,
-                                               strainTB, visStrainB,
-                                               initialStateB)
-    self.dtStableImplicit = 0.1*min(min(maxwellTimeA),
-                                    min(maxwellTimeB))
+        self._calcStress(strainB, muB, lambdaB, shearRatioB, maxwellTimeB,
+                         strainTB, visStrainB,
+                         initialStressB, initialStrainB)
+    self.dtStableImplicit = 0.1*min(min(maxwellTimeA), min(maxwellTimeB))
+
     return
 
 
-  def _calcStress(self, strainV, muV, lambdaV,
-                  shearRatioV, maxwellTimeV, strainTV, visStrainV,
-                  initialStateV):
+  def _calcStress(self, strainV, muV, lambdaV, shearRatioV, maxwellTimeV,
+                  strainTV, visStrainV, initialStressV, initialStrainV):
     """
     Compute stress and derivative of elasticity matrix.
     This assumes behavior is always viscoelastic.
@@ -177,25 +218,25 @@
     visFrac = 0.0
     visFac = 0.0
     dq = [0.0, 0.0, 0.0]
-    for model in range(self.numMaxwellModels):
-      visFrac += shearRatioV[model]
-      if shearRatioV[model] != 0.0:
-        if self.dt < timeFrac*maxwellTimeV[model]:
+    for imodel in range(numMaxwellModels):
+      visFrac += shearRatioV[imodel]
+      if shearRatioV[imodel] != 0.0:
+        if self.dt < timeFrac*maxwellTimeV[imodel]:
           fSign = 1.0
           factorial = 1.0
           fraction = 1.0
-          dq[model] = 1.0
+          dq[imodel] = 1.0
           for iTerm in range(2, numTerms + 1):
             factorial *= iTerm
             fSign *= -1.0
-            fraction *= self.dt/maxwellTimeV[model]
-            dq[model] += fSign*fraction/factorial
-        elif maxwellTimeV[model] < timeFrac*self.dt:
-          dq[model] = maxwellTimeV[model]/dt
+            fraction *= self.dt/maxwellTimeV[imodel]
+            dq[imodel] += fSign*fraction/factorial
+        elif maxwellTimeV[imodel] < timeFrac*self.dt:
+          dq[imodel] = maxwellTimeV[imodel]/dt
         else:
-          dq[model] = maxwellTimeV[model] * \
-          (1.0-math.exp(-self.dt/maxwellTimeV[model]))/self.dt
-        visFac += shearRatioV[model] * dq[model]
+          dq[imodel] = maxwellTimeV[imodel] * \
+          (1.0-math.exp(-self.dt/maxwellTimeV[imodel]))/self.dt
+        visFac += shearRatioV[imodel] * dq[imodel]
 
     elasFrac = 1.0 - visFrac
     shearFac = muV*(elasFrac + visFac)/3.0
@@ -236,20 +277,20 @@
     deltaStrain = 0.0
     devStressTpdt = 0.0
     visStrain = 0.0
-    for iComp in range(self.tensorSize):
+    for iComp in xrange(tensorSize):
       devStrainTpdt = strainV[iComp] - diag[iComp]*meanStrainTpdt
       devStrainT = strainTV[iComp] - diag[iComp]*meanStrainT
       deltaStrain = devStrainTpdt - devStrainT
       devStressTpdt = elasFrac*devStrainTpdt
-      for model in range(self.numMaxwellModels):
-        if shearRatioV[model] != 0.0:
-          visStrain = math.exp(-self.dt/maxwellTimeV[model]) * \
-                      visStrainV[iComp + self.tensorSize * model] + \
-                      dq[model] * deltaStrain
-          devStressTpdt += shearRatioV[model] * visStrain
+      for imodel in range(numMaxwellModels):
+        if shearRatioV[imodel] != 0.0:
+          visStrain = math.exp(-self.dt/maxwellTimeV[imodel]) * \
+                      visStrainV[imodel,iComp] + \
+                      dq[imodel] * deltaStrain
+          devStressTpdt += shearRatioV[imodel] * visStrain
       devStressTpdt = elasFac * devStressTpdt
       stressV[iComp] = diag[iComp]*meanStressTpdt + devStressTpdt + \
-                       initialStateV[iComp]
+                       initialStressV[iComp]
       
     stress = numpy.reshape(stressV, (6,1))
     return (elasticConsts, numpy.ravel(stress))

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc	2009-06-13 05:52:06 UTC (rev 15228)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.cc	2009-06-13 21:11:46 UTC (rev 15229)
@@ -17,29 +17,50 @@
 
 const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dimension = 3;
 
-const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numDBValues = 9;
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numLocs = 2;
 
-const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numInitialStateValues = 6;
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numProperties = 9;
 
-const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numParameters = 7;
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numStateVars = 4;
 
-const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numParamsQuadPt = 33;
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numDBProperties = 9;
 
-const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numLocs = 2;
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numDBStateVars = 24;
 
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numPropsQuadPt = 9;
+
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numVarsQuadPt = 24;
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_lengthScale =   1.00000000e+03;
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_timeScale =   1.00000000e+00;
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_pressureScale =   2.25000000e+10;
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_densityScale =   1.00000000e+03;
+
 const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dtStableImplicit =   4.44444444e+06;
 
-const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numParamValues[] = {
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numPropertyValues[] = {
 1,
 1,
 1,
-3,
-3,
+1,
+1,
+1,
+1,
+1,
+1,
+};
+
+const int pylith::materials::GenMaxwellIsotropic3DTimeDepData::_numStateVarValues[] = {
 6,
-18,
+6,
+6,
+6,
 };
 
-const char* pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbValues[] = {
+const char* pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbPropertyValues[] = {
 "density",
 "vs",
 "vp",
@@ -51,16 +72,34 @@
 "viscosity_3",
 };
 
-const char* pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStateDBValues[] = {
-"stress_xx",
-"stress_yy",
-"stress_zz",
-"stress_xy",
-"stress_yz",
-"stress_xy",
+const char* pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbStateVarValues[] = {
+"total-strain-xx",
+"total-strain-yy",
+"total-strain-zz",
+"total-strain-xy",
+"total-strain-yz",
+"total-strain-xz",
+"viscous-strain-1-xx",
+"viscous-strain-1-yy",
+"viscous-strain-1-zz",
+"viscous-strain-1-xy",
+"viscous-strain-1-yz",
+"viscous-strain-1-xz",
+"viscous-strain-2-xx",
+"viscous-strain-2-yy",
+"viscous-strain-2-zz",
+"viscous-strain-2-xy",
+"viscous-strain-2-yz",
+"viscous-strain-2-xz",
+"viscous-strain-3-xx",
+"viscous-strain-3-yy",
+"viscous-strain-3-zz",
+"viscous-strain-3-xy",
+"viscous-strain-3-yz",
+"viscous-strain-3-xz",
 };
 
-const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbData[] = {
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbProperties[] = {
   2.50000000e+03,
   3.00000000e+03,
   5.19615242e+03,
@@ -81,22 +120,58 @@
   1.00000000e+20,
 };
 
-const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStateDBData[] = {
-  1.20000000e+04,
-  2.30000000e+04,
-  3.40000000e+04,
-  4.50000000e+04,
-  5.60000000e+04,
-  6.70000000e+04,
-  2.10000000e+04,
-  3.20000000e+04,
-  4.30000000e+04,
-  5.40000000e+04,
-  6.50000000e+04,
-  7.60000000e+04,
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_dbStateVars[] = {
+  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,
+  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,
+  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,
+  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::GenMaxwellIsotropic3DTimeDepData::_parameterData[] = {
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_properties[] = {
   2.50000000e+03,
   2.25000000e+10,
   2.25000000e+10,
@@ -106,6 +181,18 @@
   8.88888889e+07,
   4.44444444e+07,
   2.22222222e+09,
+  2.00000000e+03,
+  2.88000000e+09,
+  2.88000000e+09,
+  2.00000000e-01,
+  2.00000000e-01,
+  2.00000000e-01,
+  1.73611111e+09,
+  1.73611111e+10,
+  1.73611111e+11,
+};
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_stateVars[] = {
   1.10000000e-04,
   2.20000000e-04,
   3.30000000e-04,
@@ -130,15 +217,78 @@
   4.40000000e-04,
   5.50000000e-04,
   6.60000000e-04,
-  2.00000000e+03,
-  2.88000000e+09,
-  2.88000000e+09,
+  1.20000000e-04,
+  2.30000000e-04,
+  3.40000000e-04,
+  4.50000000e-04,
+  5.60000000e-04,
+  6.70000000e-04,
+ -1.10000000e-04,
+ -2.71050543e-20,
+  1.10000000e-04,
+  4.50000000e-04,
+  5.60000000e-04,
+  6.70000000e-04,
+ -1.10000000e-04,
+ -2.71050543e-20,
+  1.10000000e-04,
+  4.50000000e-04,
+  5.60000000e-04,
+  6.70000000e-04,
+ -1.10000000e-04,
+ -2.71050543e-20,
+  1.10000000e-04,
+  4.50000000e-04,
+  5.60000000e-04,
+  6.70000000e-04,
+};
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_propertiesNondim[] = {
+  2.50000000e+00,
+  1.00000000e+00,
+  1.00000000e+00,
+  5.00000000e-01,
+  1.00000000e-01,
   2.00000000e-01,
+  8.88888889e+07,
+  4.44444444e+07,
+  2.22222222e+09,
+  2.00000000e+00,
+  1.28000000e-01,
+  1.28000000e-01,
   2.00000000e-01,
   2.00000000e-01,
+  2.00000000e-01,
   1.73611111e+09,
   1.73611111e+10,
   1.73611111e+11,
+};
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_stateVarsNondim[] = {
+  1.10000000e-04,
+  2.20000000e-04,
+  3.30000000e-04,
+  4.40000000e-04,
+  5.50000000e-04,
+  6.60000000e-04,
+ -1.10000000e-04,
+  0.00000000e+00,
+  1.10000000e-04,
+  4.40000000e-04,
+  5.50000000e-04,
+  6.60000000e-04,
+ -1.10000000e-04,
+  0.00000000e+00,
+  1.10000000e-04,
+  4.40000000e-04,
+  5.50000000e-04,
+  6.60000000e-04,
+ -1.10000000e-04,
+  0.00000000e+00,
+  1.10000000e-04,
+  4.40000000e-04,
+  5.50000000e-04,
+  6.60000000e-04,
   1.20000000e-04,
   2.30000000e-04,
   3.40000000e-04,
@@ -165,21 +315,6 @@
   6.70000000e-04,
 };
 
-const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialState[] = {
-  1.20000000e+04,
-  2.30000000e+04,
-  3.40000000e+04,
-  4.50000000e+04,
-  5.60000000e+04,
-  6.70000000e+04,
-  2.10000000e+04,
-  3.20000000e+04,
-  4.30000000e+04,
-  5.40000000e+04,
-  6.50000000e+04,
-  7.60000000e+04,
-};
-
 const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_density[] = {
   2.50000000e+03,
   2.00000000e+03,
@@ -201,18 +336,18 @@
 };
 
 const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_stress[] = {
-  1.98198741e+07,
-  2.47730000e+07,
-  2.97261259e+07,
-  1.98135037e+07,
-  2.47666296e+07,
-  2.97197555e+07,
-  2.69941620e+06,
-  3.34400000e+06,
-  3.98858380e+06,
+  1.98288741e+07,
+  2.47720000e+07,
+  2.97151259e+07,
+  1.97925037e+07,
+  2.47356296e+07,
+  2.96787555e+07,
+  2.72941620e+06,
+  3.36400000e+06,
+  3.99858380e+06,
   2.64593371e+06,
-  3.29051751e+06,
-  3.93510131e+06,
+  3.28051751e+06,
+  3.91510131e+06,
 };
 
 const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_elasticConsts[] = {
@@ -260,26 +395,70 @@
   5.75992635e+09,
 };
 
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStress[] = {
+  2.10000000e+04,
+  2.20000000e+04,
+  2.30000000e+04,
+  2.40000000e+04,
+  2.50000000e+04,
+  2.60000000e+04,
+  5.10000000e+04,
+  5.20000000e+04,
+  5.30000000e+04,
+  5.40000000e+04,
+  5.50000000e+04,
+  5.60000000e+04,
+};
+
+const double pylith::materials::GenMaxwellIsotropic3DTimeDepData::_initialStrain[] = {
+  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::GenMaxwellIsotropic3DTimeDepData::_stateVarsUpdated = 0;
+
 pylith::materials::GenMaxwellIsotropic3DTimeDepData::GenMaxwellIsotropic3DTimeDepData(void)
 { // constructor
   dimension = _dimension;
-  numDBValues = _numDBValues;
-  numInitialStateValues = _numInitialStateValues;
-  numParameters = _numParameters;
-  numParamsQuadPt = _numParamsQuadPt;
   numLocs = _numLocs;
+  numProperties = _numProperties;
+  numStateVars = _numStateVars;
+  numDBProperties = _numDBProperties;
+  numDBStateVars = _numDBStateVars;
+  numPropsQuadPt = _numPropsQuadPt;
+  numVarsQuadPt = _numVarsQuadPt;
+  lengthScale = _lengthScale;
+  timeScale = _timeScale;
+  pressureScale = _pressureScale;
+  densityScale = _densityScale;
   dtStableImplicit = _dtStableImplicit;
-  numParamValues = const_cast<int*>(_numParamValues);
-  dbValues = const_cast<char**>(_dbValues);
-  initialStateDBValues = const_cast<char**>(_initialStateDBValues);
-  dbData = const_cast<double*>(_dbData);
-  initialStateDBData = const_cast<double*>(_initialStateDBData);
-  parameterData = const_cast<double*>(_parameterData);
-  initialState = const_cast<double*>(_initialState);
+  numPropertyValues = const_cast<int*>(_numPropertyValues);
+  numStateVarValues = const_cast<int*>(_numStateVarValues);
+  dbPropertyValues = const_cast<char**>(_dbPropertyValues);
+  dbStateVarValues = const_cast<char**>(_dbStateVarValues);
+  dbProperties = const_cast<double*>(_dbProperties);
+  dbStateVars = const_cast<double*>(_dbStateVars);
+  properties = const_cast<double*>(_properties);
+  stateVars = const_cast<double*>(_stateVars);
+  propertiesNondim = const_cast<double*>(_propertiesNondim);
+  stateVarsNondim = const_cast<double*>(_stateVarsNondim);
   density = const_cast<double*>(_density);
   strain = const_cast<double*>(_strain);
   stress = const_cast<double*>(_stress);
   elasticConsts = const_cast<double*>(_elasticConsts);
+  initialStress = const_cast<double*>(_initialStress);
+  initialStrain = const_cast<double*>(_initialStrain);
+  stateVarsUpdated = const_cast<double*>(_stateVarsUpdated);
 } // constructor
 
 pylith::materials::GenMaxwellIsotropic3DTimeDepData::~GenMaxwellIsotropic3DTimeDepData(void)

Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh	2009-06-13 05:52:06 UTC (rev 15228)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/GenMaxwellIsotropic3DTimeDepData.hh	2009-06-13 21:11:46 UTC (rev 15229)
@@ -39,32 +39,50 @@
 
   static const int _dimension;
 
-  static const int _numDBValues;
+  static const int _numLocs;
 
-  static const int _numInitialStateValues;
+  static const int _numProperties;
 
-  static const int _numParameters;
+  static const int _numStateVars;
 
-  static const int _numParamsQuadPt;
+  static const int _numDBProperties;
 
-  static const int _numLocs;
+  static const int _numDBStateVars;
 
+  static const int _numPropsQuadPt;
+
+  static const int _numVarsQuadPt;
+
+  static const double _lengthScale;
+
+  static const double _timeScale;
+
+  static const double _pressureScale;
+
+  static const double _densityScale;
+
   static const double _dtStableImplicit;
 
-  static const int _numParamValues[];
+  static const int _numPropertyValues[];
 
-  static const char* _dbValues[];
+  static const int _numStateVarValues[];
 
-  static const char* _initialStateDBValues[];
+  static const char* _dbPropertyValues[];
 
-  static const double _dbData[];
+  static const char* _dbStateVarValues[];
 
-  static const double _initialStateDBData[];
+  static const double _dbProperties[];
 
-  static const double _parameterData[];
+  static const double _dbStateVars[];
 
-  static const double _initialState[];
+  static const double _properties[];
 
+  static const double _stateVars[];
+
+  static const double _propertiesNondim[];
+
+  static const double _stateVarsNondim[];
+
   static const double _density[];
 
   static const double _strain[];
@@ -73,6 +91,12 @@
 
   static const double _elasticConsts[];
 
+  static const double _initialStress[];
+
+  static const double _initialStrain[];
+
+  static const double* _stateVarsUpdated;
+
 };
 
 #endif // pylith_materials_genmaxwellisotropic3dtimedepdata_hh



More information about the CIG-COMMITS mailing list