[cig-commits] r6944 - short/3D/PyLith/trunk/libsrc/materials

[willic3 at geodynamics.org]

Author: willic3
Date: 2007-05-23 12:38:17 -0700 (Wed, 23 May 2007)
New Revision: 6944

Modified:
   short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
Log:
More work on stress computation for Maxwell isotropic material.


Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc	2007-05-23 18:27:50 UTC (rev 6943)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc	2007-05-23 19:38:17 UTC (rev 6944)
@@ -197,8 +197,16 @@
   const double mu = parameters[_MaxwellIsotropic3D::pidMu][0];
   const double lambda = parameters[_MaxwellIsotropic3D::pidLambda][0];
   const double viscosity = parameters[_MaxwellIsotropic3D::pidViscosity][0];
+  const double maxwelltime = parameters[_MaxwellIsotropic3D::pidMaxwellTime][0];
 
   const double lambda2mu = lambda + 2.0 * mu;
+  const double bulkmodulus = lambda + 2.0 * mu/3.0;
+
+  // The following three definitions are dummies until we figure out how to get the
+  // information into the function.
+  const bool elastic = true;
+  const double dt = 1.0;
+  const bool updatestate = true;
   
   const double e11 = totalStrain[0];
   const double e22 = totalStrain[1];
@@ -207,14 +215,45 @@
   const double e23 = totalStrain[4];
   const double e13 = totalStrain[5];
   
-  const double s123 = lambda * (e11 + e22 + e33);
+  const double etraceTpdt = e11 + e22 + e33;
+  const double emeanTpdt = etraceTpdt/3.0;
+  const double s123 = lambda * etraceTpdt;
 
-  (*stress)[0] = s123 + 2.0*mu*e11;
-  (*stress)[1] = s123 + 2.0*mu*e22;
-  (*stress)[2] = s123 + 2.0*mu*e33;
-  (*stress)[3] = 2.0 * mu * e12;
-  (*stress)[4] = 2.0 * mu * e23;
-  (*stress)[5] = 2.0 * mu * e13;
+  if (elastic) {
+    (*stress)[0] = s123 + 2.0*mu*e11;
+    (*stress)[1] = s123 + 2.0*mu*e22;
+    (*stress)[2] = s123 + 2.0*mu*e33;
+    (*stress)[3] = 2.0 * mu * e12;
+    (*stress)[4] = 2.0 * mu * e23;
+    (*stress)[5] = 2.0 * mu * e13;
+  } else {
+    const double diag[6] = {1.0, 1.0, 1.0, 0.0, 0.0, 0.0};
+    const double smeanTpdt = bulkmodulus * etraceTpdt;
+    // The code below should probably be in a separate function since it
+    // is used more than once.  I should also probably cover the possibility
+    // that Maxwell time is zero (although this should never happen).
+    const double timeFrac = 1.0e-5;
+    const double numTerms = 5.0;
+    double dq = 0.0;
+    if(maxwelltime < timeFrac*dt) {
+      double fSign = 1.0;
+      double factorial = 1.0;
+      double fraction = 1.0;
+      dq = 1.0;
+      double term = 2.0;
+      while (term <= numTerms) {
+	factorial = factorial*term;
+	fSign = -1.0 * fSign;
+	fraction = fraction * (dt/maxwelltime);
+	dq = dq + fSign*fraction/factorial;
+	term++;
+      }
+    } else {
+      dq = maxwelltime*(1.0-exp(-dt/maxwelltime))/dt;
+    }
+
+
+
 } // _calcStress
 
 // ----------------------------------------------------------------------