[cig-commits] r18575 - in short/3D/PyLith/trunk/doc/userguide: . alternativeformul governingeqns materials

willic3 at geodynamics.org willic3 at geodynamics.org
Thu Jun 9 22:16:18 PDT 2011


Author: willic3
Date: 2011-06-09 22:16:17 -0700 (Thu, 09 Jun 2011)
New Revision: 18575

Added:
   short/3D/PyLith/trunk/doc/userguide/alternativeformul/
   short/3D/PyLith/trunk/doc/userguide/alternativeformul/alternativeformul.lyx
Modified:
   short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx
   short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx
   short/3D/PyLith/trunk/doc/userguide/userguide.lyx
Log:
Moved unused formulation to an appendix.
Added some background info for linear viscoelastic models.
Cleaned up the descriptions and equations a bit.
Changed equations to describe generalized Maxwell model rather than
regular Maxwell model.
I still need to go through and make sure there aren't any mistakes or
inconsistencies.



Added: short/3D/PyLith/trunk/doc/userguide/alternativeformul/alternativeformul.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/alternativeformul/alternativeformul.lyx	                        (rev 0)
+++ short/3D/PyLith/trunk/doc/userguide/alternativeformul/alternativeformul.lyx	2011-06-10 05:16:17 UTC (rev 18575)
@@ -0,0 +1,320 @@
+#LyX 2.0 created this file. For more info see http://www.lyx.org/
+\lyxformat 413
+\begin_document
+\begin_header
+\textclass book
+\begin_preamble
+
+\end_preamble
+\use_default_options false
+\maintain_unincluded_children false
+\language english
+\language_package default
+\inputencoding latin1
+\fontencoding global
+\font_roman default
+\font_sans default
+\font_typewriter default
+\font_default_family default
+\use_non_tex_fonts false
+\font_sc false
+\font_osf false
+\font_sf_scale 100
+\font_tt_scale 100
+
+\graphics default
+\default_output_format default
+\output_sync 0
+\bibtex_command default
+\index_command default
+\paperfontsize default
+\spacing single
+\use_hyperref false
+\papersize default
+\use_geometry true
+\use_amsmath 1
+\use_esint 0
+\use_mhchem 0
+\use_mathdots 1
+\cite_engine basic
+\use_bibtopic false
+\use_indices false
+\paperorientation portrait
+\suppress_date false
+\use_refstyle 0
+\index Index
+\shortcut idx
+\color #008000
+\end_index
+\leftmargin 1in
+\topmargin 1in
+\rightmargin 1in
+\bottommargin 2in
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\paragraph_indentation default
+\quotes_language english
+\papercolumns 1
+\papersides 2
+\paperpagestyle default
+\tracking_changes false
+\output_changes false
+\html_math_output 0
+\html_css_as_file 0
+\html_be_strict false
+\end_header
+
+\begin_body
+
+\begin_layout Chapter
+\begin_inset CommandInset label
+LatexCommand label
+name "cha:Alternative-Formulations"
+
+\end_inset
+
+Alternative Material Model Formulations
+\end_layout
+
+\begin_layout Section
+\begin_inset CommandInset label
+LatexCommand label
+name "sec:ViscoelasticFormulations"
+
+\end_inset
+
+Viscoelastic Formulations
+\end_layout
+
+\begin_layout Standard
+The viscoelastic formulations presently used in PyLith are described in
+ 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sec:Viscoelastic-Materials"
+
+\end_inset
+
+.
+ In some cases there are alternative formulations that may be used in future
+ versions of PyLith, and those are described here.
+\end_layout
+
+\begin_layout Subsection
+\begin_inset CommandInset label
+LatexCommand label
+name "sub:Effective-Stress-Formulation-Maxwell"
+
+\end_inset
+
+Effective Stress Formulation for a Linear Maxwell Viscoelastic Material
+\end_layout
+
+\begin_layout Standard
+An alternative technique for solving the equations for a Maxwell viscoelastic
+ material is based on the effective stress formulation described in 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sub:Effective-Stress-Formulations-Viscoelastic"
+
+\end_inset
+
+.
+ A linear Maxwell viscoelastic material may be characterized by the same
+ elastic parameters as an isotropic elastic material (
+\begin_inset Formula $E$
+\end_inset
+
+ and 
+\begin_inset Formula $\nu$
+\end_inset
+
+), as well as the viscosity, 
+\begin_inset Formula $\eta$
+\end_inset
+
+.
+ The creep strain increment is
+\begin_inset Formula 
+\begin{gather}
+\underline{\Delta e}^{C}=\frac{\Delta t\phantom{}^{\tau}\underline{S}}{2\eta}\,\,.\label{eq:D1}
+\end{gather}
+
+\end_inset
+
+Therefore,
+\begin_inset Formula 
+\begin{gather}
+\Delta\overline{e}^{C}=\frac{\Delta t\sqrt{^{\tau}J_{2}^{\prime}}}{\sqrt{3\eta}}=\frac{\Delta t\phantom{}^{\tau}\overline{\sigma}}{3\eta}\,,\,\mathrm{and}\,^{\tau}\gamma=\frac{1}{2\eta}\,\,.\label{eq:D2}
+\end{gather}
+
+\end_inset
+
+Substituting Equations 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:46"
+
+\end_inset
+
+, 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:D1"
+
+\end_inset
+
+, and 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:D2"
+
+\end_inset
+
+ into 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:43"
+
+\end_inset
+
+, we obtain
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\underline{S}=\frac{1}{a_{E}}\left\{ ^{t+\Delta t}\underline{e}^{\prime}-\frac{\Delta t}{2\eta}\left[(1-\alpha)^{t}\underline{S}+\alpha\phantom{}^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,\,.\label{eq:D3}
+\end{gather}
+
+\end_inset
+
+Solving for 
+\begin_inset Formula $^{t+\Delta t}\underline{S}$
+\end_inset
+
+,
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\underline{S}=\frac{1}{a_{E}+\frac{\alpha\Delta t}{2\eta}}\left[^{t+\Delta t}\underline{e}^{\prime}-\frac{\Delta t}{2\eta}(1-\alpha)^{t}\underline{S}+\frac{1+\mathrm{\nu}}{E}\underline{S}^{I}\right]\,\,.\label{eq:D4}
+\end{gather}
+
+\end_inset
+
+In this case it is possible to solve directly for the deviatoric stresses,
+ and the effective stress function approach is not needed.
+ To obtain the total stress, we simply use
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\sigma_{ij}=\phantom{}^{t+\Delta t}S_{ij}+\frac{\mathit{1}}{a_{m}}\left(\,^{t+\Delta t}\theta-\theta^{I}\right)\delta_{ij}+P^{I}\delta_{ij}\,\,.\label{eq:D5}
+\end{gather}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+To compute the viscoelastic tangent material matrix relating stress and
+ strain, we need to compute the first term in 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:58"
+
+\end_inset
+
+.
+ From Equation 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:D4"
+
+\end_inset
+
+, we have
+\begin_inset Formula 
+\begin{gather}
+\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=\frac{\delta_{ik}}{a_{E}+\frac{\alpha\Delta t}{2\eta}}\,\,.\label{eq:D12}
+\end{gather}
+
+\end_inset
+
+Using this, along with 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:58"
+
+\end_inset
+
+, 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:59"
+
+\end_inset
+
+, and 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:60"
+
+\end_inset
+
+, the final material matrix relating stress and tensor strain is
+\begin_inset Formula 
+\begin{gather}
+C_{ij}^{VE}=\frac{1}{3a_{m}}\left[\begin{array}{cccccc}
+1 & 1 & 1 & 0 & 0 & 0\\
+1 & 1 & 1 & 0 & 0 & 0\\
+1 & 1 & 1 & 0 & 0 & 0\\
+0 & 0 & 0 & 0 & 0 & 0\\
+0 & 0 & 0 & 0 & 0 & 0\\
+0 & 0 & 0 & 0 & 0 & 0
+\end{array}\right]+\frac{1}{3\left(a_{E}+\frac{\alpha\Delta t}{2\eta}\right)}\left[\begin{array}{cccccc}
+2 & -1 & -1 & 0 & 0 & 0\\
+-1 & 2 & -1 & 0 & 0 & 0\\
+-1 & -1 & 2 & 0 & 0 & 0\\
+0 & 0 & 0 & 3 & 0 & 0\\
+0 & 0 & 0 & 0 & 3 & 0\\
+0 & 0 & 0 & 0 & 0 & 3
+\end{array}\right]\,.\label{eq:D13}
+\end{gather}
+
+\end_inset
+
+Note that the coefficient of the second matrix approaches 
+\begin_inset Formula $E/3(1+\nu)=1/3a_{E}$
+\end_inset
+
+ as 
+\begin_inset Formula $\eta$
+\end_inset
+
+ goes to infinity.
+ To check the results we make sure that the regular elastic constitutive
+ matrix is obtained for selected terms in the case where 
+\begin_inset Formula $\eta$
+\end_inset
+
+ goes to infinity.
+\begin_inset Formula 
+\begin{gather}
+C_{11}^{E}=\frac{E(1-\nu)}{(1+\nu)(1-2\nu)}\,\,\nonumber \\
+C_{12}^{E}=\frac{E\nu}{(1+\nu)(1-2\nu)}\,.\label{eq:D14}\\
+C_{44}^{E}=\frac{E}{1+\nu}\,\,\nonumber 
+\end{gather}
+
+\end_inset
+
+This is consistent with the regular elasticity matrix, and Equation 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:D13"
+
+\end_inset
+
+ should thus be used when forming the stiffness matrix.
+ We do not presently use this formulation, but it may be included in future
+ versions.
+\end_layout
+
+\end_body
+\end_document

Modified: short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx	2011-06-10 01:29:23 UTC (rev 18574)
+++ short/3D/PyLith/trunk/doc/userguide/governingeqns/governingeqns.lyx	2011-06-10 05:16:17 UTC (rev 18575)
@@ -1,21 +1,29 @@
-#LyX 1.6.4 created this file. For more info see http://www.lyx.org/
-\lyxformat 345
+#LyX 2.0 created this file. For more info see http://www.lyx.org/
+\lyxformat 413
 \begin_document
 \begin_header
 \textclass book
 \use_default_options false
+\maintain_unincluded_children false
 \language english
+\language_package default
 \inputencoding auto
+\fontencoding global
 \font_roman default
 \font_sans default
 \font_typewriter default
 \font_default_family default
+\use_non_tex_fonts false
 \font_sc false
 \font_osf false
 \font_sf_scale 100
 \font_tt_scale 100
 
 \graphics default
+\default_output_format default
+\output_sync 0
+\bibtex_command default
+\index_command default
 \paperfontsize default
 \spacing single
 \use_hyperref false
@@ -23,9 +31,18 @@
 \use_geometry true
 \use_amsmath 1
 \use_esint 0
+\use_mhchem 1
+\use_mathdots 1
 \cite_engine basic
 \use_bibtopic false
+\use_indices false
 \paperorientation portrait
+\suppress_date false
+\use_refstyle 0
+\index Index
+\shortcut idx
+\color #008000
+\end_index
 \leftmargin 1in
 \topmargin 1in
 \rightmargin 1in
@@ -33,15 +50,16 @@
 \secnumdepth 3
 \tocdepth 3
 \paragraph_separation indent
-\defskip medskip
+\paragraph_indentation default
 \quotes_language english
 \papercolumns 1
 \papersides 1
 \paperpagestyle default
 \tracking_changes false
 \output_changes false
-\author "" 
-\author "" 
+\html_math_output 0
+\html_css_as_file 0
+\html_be_strict false
 \end_header
 
 \begin_body
@@ -110,7 +128,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="11" columns="3">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="0">
 <column alignment="center" valignment="top" width="0">
 <column alignment="center" valignment="top" width="0">
@@ -521,14 +539,18 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial u_{i}}{\partial t}\, dV=\int_{V}f_{i}\, dV+\int_{S}T_{i}\, dS.\label{eqn:momentum:index}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial u_{i}}{\partial t}\, dV=\int_{V}f_{i}\, dV+\int_{S}T_{i}\, dS.\label{eqn:momentum:index}
+\end{equation}
 
 \end_inset
 
 The traction vector field is related to the stress tensor through
-\begin_inset Formula \begin{equation}
-T_{i}=\sigma_{ij}n_{j},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+T_{i}=\sigma_{ij}n_{j},
+\end{equation}
 
 \end_inset
 
@@ -549,26 +571,34 @@
 \end_inset
 
  yields
-\begin_inset Formula \begin{equation}
-\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial u_{i}}{\partial t}\, dV=\int_{V}f_{i}\, dV+\int_{S}\sigma_{ij}n_{j}\, dS.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial u_{i}}{\partial t}\, dV=\int_{V}f_{i}\, dV+\int_{S}\sigma_{ij}n_{j}\, dS.
+\end{equation}
 
 \end_inset
 
 Applying the divergence theorem,
-\begin_inset Formula \begin{equation}
-\int_{V}a_{i,j}\: dV=\int_{S}a_{j}n_{j}\: dS,\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}a_{i,j}\: dV=\int_{S}a_{j}n_{j}\: dS,
+\end{equation}
 
 \end_inset
 
 to the surface integral results in
-\begin_inset Formula \begin{equation}
-\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial u_{i}}{\partial t}\, dV=\int_{V}f_{i}\, dV+\int_{V}\sigma_{ij,j}\, dV,\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial u_{i}}{\partial t}\, dV=\int_{V}f_{i}\, dV+\int_{V}\sigma_{ij,j}\, dV,
+\end{equation}
 
 \end_inset
 
 which we can rewrite as
-\begin_inset Formula \begin{equation}
-\int_{V}\left(\rho\frac{\partial^{2}u_{i}}{\partial t^{2}}-f_{i}-\sigma_{ij,j}\right)\, dV=0.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}\left(\rho\frac{\partial^{2}u_{i}}{\partial t^{2}}-f_{i}-\sigma_{ij,j}\right)\, dV=0.
+\end{equation}
 
 \end_inset
 
@@ -578,11 +608,13 @@
 
  is arbitrary, the integrand must be zero at every location in the volume,
  so that we end up with
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \rho\frac{\partial^{2}u_{i}}{\partial t^{2}}-f_{i}-\sigma_{ij,j}=0\text{ in }V,\\
 \sigma_{ij}n_{j}=T_{i}\text{ on }S_{T}\text{,}\\
 u_{i}=u_{i}^{o}\text{ on }S_{u}\text{, and}\\
-R_{ki}(u_{i}^{+}-u_{i}^{-})=d_{k}\text{ on }S_{f}.\end{gather}
+R_{ki}(u_{i}^{+}-u_{i}^{-})=d_{k}\text{ on }S_{f}.
+\end{gather}
 
 \end_inset
 
@@ -663,14 +695,18 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial\vec{u}}{\partial t}\, dV=\int_{V}\overrightarrow{f}\, dV+\int_{S}\overrightarrow{T}\, dS.\label{eqn:momentum:vec}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial\vec{u}}{\partial t}\, dV=\int_{V}\overrightarrow{f}\, dV+\int_{S}\overrightarrow{T}\, dS.\label{eqn:momentum:vec}
+\end{equation}
 
 \end_inset
 
 The traction vector field is related to the stress tensor through
-\begin_inset Formula \begin{equation}
-\overrightarrow{T}=\underline{\sigma}\cdot\overrightarrow{n},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\overrightarrow{T}=\underline{\sigma}\cdot\overrightarrow{n},
+\end{equation}
 
 \end_inset
 
@@ -691,26 +727,34 @@
 \end_inset
 
  yields
-\begin_inset Formula \begin{equation}
-\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial\overrightarrow{u}}{\partial t}\, dV=\int_{V}\overrightarrow{f}\, dV+\int_{S}\underline{\sigma}\cdot\overrightarrow{n}\, dS.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial\overrightarrow{u}}{\partial t}\, dV=\int_{V}\overrightarrow{f}\, dV+\int_{S}\underline{\sigma}\cdot\overrightarrow{n}\, dS.
+\end{equation}
 
 \end_inset
 
 Applying the divergence theorem,
-\begin_inset Formula \begin{equation}
-\int_{V}\nabla\cdot\overrightarrow{a}\: dV=\int_{S}\overrightarrow{a}\cdot\overrightarrow{n}\: dS,\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}\nabla\cdot\overrightarrow{a}\: dV=\int_{S}\overrightarrow{a}\cdot\overrightarrow{n}\: dS,
+\end{equation}
 
 \end_inset
 
 to the surface integral results in
-\begin_inset Formula \begin{equation}
-\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial\overrightarrow{u}}{\partial t}\, dV=\int_{V}\overrightarrow{f}\, dV+\int_{V}\nabla\cdot\underline{\sigma}\, dV,\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial}{\partial t}\int_{V}\rho\frac{\partial\overrightarrow{u}}{\partial t}\, dV=\int_{V}\overrightarrow{f}\, dV+\int_{V}\nabla\cdot\underline{\sigma}\, dV,
+\end{equation}
 
 \end_inset
 
 which we can rewrite as
-\begin_inset Formula \begin{equation}
-\int_{V}\left(\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}-\overrightarrow{f}-\nabla\cdot\overrightarrow{\sigma}\right)\, dV=\vec{0}.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}\left(\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}-\overrightarrow{f}-\nabla\cdot\overrightarrow{\sigma}\right)\, dV=\vec{0}.
+\end{equation}
 
 \end_inset
 
@@ -720,11 +764,13 @@
 
  is arbitrary, the integrand must be the zero vector at every location in
  the volume, so that we end up with
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}-\overrightarrow{f}-\nabla\cdot\overrightarrow{\sigma}=\vec{0}\text{ in }V,\\
 \underline{\sigma}\cdot\overrightarrow{n}=\overrightarrow{T}\text{ on }S_{T}\text{,}\\
 \overrightarrow{u}=\overrightarrow{u^{o}}\text{ on }S_{u},\text{ and}\\
-\underbar{R}\cdot(\vec{u^{+}}-\vec{u^{-}})=\vec{d}\text{ on }S_{f}.\end{gather}
+\underbar{R}\cdot(\vec{u^{+}}-\vec{u^{-}})=\vec{d}\text{ on }S_{f}.
+\end{gather}
 
 \end_inset
 
@@ -822,20 +868,24 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \sigma_{ij,j}+f_{i}=\rho\ddot{u_{i}}\text{ in }V,\\
 \sigma_{ij}n_{j}=T_{i}\text{ on }S_{T},\\
 u_{i}=u_{i}^{o}\text{ on }S_{u},\\
 R_{ki}(u_{i}^{+}-u_{i}^{-})=d_{k}\text{ on }S_{f},\text{ and}\\
-\sigma_{ij}=\sigma_{ji}\text{ (symmetric).}\end{gather}
+\sigma_{ij}=\sigma_{ji}\text{ (symmetric).}
+\end{gather}
 
 \end_inset
 
 We construct the weak form by computing the dot product of the wave equation
  and weighting function and setting the integral over the domain to zero:
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \int_{V}\left(\sigma_{ij,j}+f_{i}-\rho\ddot{u}_{i}\right)\phi_{i}\, dV=0\text{, or }\\
-\int_{V}\sigma_{ij,j}\phi_{i}\: dV+\int_{V}f_{i}\phi_{i}\: dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\: dV=0.\end{gather}
+\int_{V}\sigma_{ij,j}\phi_{i}\: dV+\int_{V}f_{i}\phi_{i}\: dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\: dV=0.
+\end{gather}
 
 \end_inset
 
@@ -845,21 +895,27 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{equation}
-\int_{V}(\sigma_{ij}\phi_{i})_{,j}\, dV=\int_{S}(\sigma_{ij}\phi_{i})n_{i}\, dS.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}(\sigma_{ij}\phi_{i})_{,j}\, dV=\int_{S}(\sigma_{ij}\phi_{i})n_{i}\, dS.
+\end{equation}
 
 \end_inset
 
 Expanding the left-hand side yields
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \int_{V}\sigma_{ij,j}\phi_{i}\: dV+\int_{V}\sigma_{ij}\phi_{i,j}\: dV=\int_{S}\sigma_{ij}\phi_{i}n_{i}\: dS,\text{ or}\\
-\int_{V}\sigma_{ij,j}\phi_{i}\: dV=-\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S}\sigma_{ij}\phi_{i}n_{i}\, dS.\end{gather}
+\int_{V}\sigma_{ij,j}\phi_{i}\: dV=-\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S}\sigma_{ij}\phi_{i}n_{i}\, dS.
+\end{gather}
 
 \end_inset
 
 Substituting into the weak form gives
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.
+\end{equation}
 
 \end_inset
 
@@ -888,29 +944,37 @@
 \end_inset
 
 ),
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S_{T}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{S_{u}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0,\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}\phi_{i,j}\, dV+\int_{S_{T}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{S_{u}}\sigma_{ij}\phi_{i}n_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0,
+\end{equation}
 
 \end_inset
 
 and recognize that
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \sigma_{ij}n_{i}=T_{i}\text{ on }S_{T}\text{ and}\\
-\phi_{i}=0\text{ on }S_{u},\end{gather}
+\phi_{i}=0\text{ on }S_{u},
+\end{gather}
 
 \end_inset
 
 so that the equation reduces to
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}\phi_{i,j}\: dV+\int_{S_{T}}T_{i}\phi_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\label{eq:elasticity:integral}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}\phi_{i,j}\: dV+\int_{S_{T}}T_{i}\phi_{i}\, dS+\int_{V}f_{i}\phi_{i}\, dV-\int_{V}\rho\ddot{u}_{i}\phi_{i}\, dV=0.\label{eq:elasticity:integral}
+\end{equation}
 
 \end_inset
 
 We express the trial solution and weighting function as linear combinations
  of basis functions,
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 u_{i}=\sum_{m}a_{i}^{m}N^{m},\\
-\phi_{i}=\sum_{n}c_{i}^{n}N^{n}.\end{gather}
+\phi_{i}=\sum_{n}c_{i}^{n}N^{n}.
+\end{gather}
 
 \end_inset
 
@@ -930,9 +994,11 @@
 .
  Substituting in the expressions for the trial solution and weighting function
  yields
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 -\int_{V}\sigma_{ij}\sum_{n}c_{i}^{n}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}\sum_{n}c_{i}^{n}N^{n}\, dS+\int_{V}f_{i}\sum_{n}c_{i}^{n}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}\sum_{n}c_{i}^{n}N^{n}\ dV=0,\text{ or}\\
-\sum_{n}c_{i}^{n}(-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV)=0.\end{gather}
+\sum_{n}c_{i}^{n}(-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV)=0.
+\end{gather}
 
 \end_inset
 
@@ -946,8 +1012,10 @@
 \end_inset
 
 
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV=\vec{0}.\label{eq:elasticity:integral:discretized}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}N^{m}N^{n}\ dV=\vec{0}.\label{eq:elasticity:integral:discretized}
+\end{equation}
 
 \end_inset
 
@@ -968,9 +1036,11 @@
 \bar no
 \noun off
 \color none
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 u_{i}=u_{i}^{o}\text{ on }S_{u},\text{ and}\\
-R_{ki}(u_{i}^{+}-u_{i}^{-})=d_{k}\text{ on }S_{f},\end{gather}
+R_{ki}(u_{i}^{+}-u_{i}^{-})=d_{k}\text{ on }S_{f},
+\end{gather}
 
 \end_inset
 
@@ -986,12 +1056,14 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \nabla\cdot\underline{\sigma}+\overrightarrow{f}=\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\text{ in }V,\\
 \underline{\sigma}\cdot\overrightarrow{n}=\overrightarrow{T}\text{ on }S_{T},\\
 \overrightarrow{u}=\overrightarrow{u^{o}}\text{ on }S_{u},\\
 \underbar{R}\cdot(\overrightarrow{u^{+}}-\overrightarrow{u^{-}})=\vec{d}\text{ on }S_{f}\\
-\underline{\sigma}=\underline{\sigma}^{T}\text{ (symmetric).}\end{gather}
+\underline{\sigma}=\underline{\sigma}^{T}\text{ (symmetric).}
+\end{gather}
 
 \end_inset
 
@@ -1010,9 +1082,11 @@
 \end_inset
 
  Hence our weak form is
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \int_{V}\left(\nabla\cdot\underline{\sigma}+\overrightarrow{f}-\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\right)\cdot\overrightarrow{\phi}\, dV=0\text{, or }\\
-\int_{V}(\nabla\cdot\underline{\sigma})\cdot\overrightarrow{\phi}\: dV+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\: dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\: dV=0.\end{gather}
+\int_{V}(\nabla\cdot\underline{\sigma})\cdot\overrightarrow{\phi}\: dV+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\: dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\: dV=0.
+\end{gather}
 
 \end_inset
 
@@ -1022,26 +1096,34 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{equation}
-\int_{V}\nabla\cdot(\underline{\sigma}\cdot\overrightarrow{\phi})\, dV=\int_{S}(\underline{\sigma}\cdot\overrightarrow{\phi})\cdot\overrightarrow{n}\, dS.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}\nabla\cdot(\underline{\sigma}\cdot\overrightarrow{\phi})\, dV=\int_{S}(\underline{\sigma}\cdot\overrightarrow{\phi})\cdot\overrightarrow{n}\, dS.
+\end{equation}
 
 \end_inset
 
 Expanding the left-hand side yields
-\begin_inset Formula \begin{equation}
-\int_{V}(\nabla\cdot\underline{\sigma})\cdot\overrightarrow{\phi}\: dV+\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\: dV=\int_{S}(\underline{\sigma}\cdot\overrightarrow{\phi})\cdot\overrightarrow{n}\: dS,\text{ or}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}(\nabla\cdot\underline{\sigma})\cdot\overrightarrow{\phi}\: dV+\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\: dV=\int_{S}(\underline{\sigma}\cdot\overrightarrow{\phi})\cdot\overrightarrow{n}\: dS,\text{ or}
+\end{equation}
 
 \end_inset
 
 
-\begin_inset Formula \begin{equation}
-\int_{V}(\nabla\cdot\underline{\sigma})\cdot\overrightarrow{\phi}\: dV=-\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\, dV+\int_{S}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}(\nabla\cdot\underline{\sigma})\cdot\overrightarrow{\phi}\: dV=-\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\, dV+\int_{S}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS.
+\end{equation}
 
 \end_inset
 
 Substituting into the weak form gives
-\begin_inset Formula \begin{equation}
--\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\, dV+\int_{S}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\, dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\, dV=0.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\, dV+\int_{S}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\, dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\, dV=0.
+\end{equation}
 
 \end_inset
 
@@ -1058,29 +1140,37 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{multline}
--\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\, dV+\int_{S_{T}}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS+\int_{S_{u}}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\, dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\, dV=0,\end{multline}
+\begin_inset Formula 
+\begin{multline}
+-\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\, dV+\int_{S_{T}}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS+\int_{S_{u}}\underline{\sigma}\cdot\overrightarrow{n}\cdot\overrightarrow{\phi}\, dS+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\, dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\, dV=0,
+\end{multline}
 
 \end_inset
 
 and recognize that
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \underline{\sigma}\cdot\overrightarrow{n}=\overrightarrow{T}\text{ on }S_{T}\text{ and}\\
-\overrightarrow{\phi}=0\text{ on }S_{u},\end{gather}
+\overrightarrow{\phi}=0\text{ on }S_{u},
+\end{gather}
 
 \end_inset
 
 so that the equation reduces to
-\begin_inset Formula \begin{equation}
--\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\: dV+\int_{S_{T}}\overrightarrow{T}\cdot\overrightarrow{\phi}\, dS+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\, dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\, dV=0.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\underline{\sigma}:\nabla\overrightarrow{\phi}\: dV+\int_{S_{T}}\overrightarrow{T}\cdot\overrightarrow{\phi}\, dS+\int_{V}\overrightarrow{f}\cdot\overrightarrow{\phi}\, dV-\int_{V}\rho\frac{\partial^{2}\overrightarrow{u}}{\partial t^{2}}\cdot\overrightarrow{\phi}\, dV=0.
+\end{equation}
 
 \end_inset
 
 We express the trial solution and weighting function as linear combinations
  of basis functions,
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \vec{u}=\sum_{m}\overrightarrow{a^{m}}N^{m},\\
-\vec{\phi}=\sum_{n}\overrightarrow{c^{n}}N^{n}.\end{gather}
+\vec{\phi}=\sum_{n}\overrightarrow{c^{n}}N^{n}.
+\end{gather}
 
 \end_inset
 
@@ -1100,9 +1190,11 @@
 .
  Substituting in the expressions for the trial solution and weighting function
  yields
-\begin_inset Formula \begin{multline}
+\begin_inset Formula 
+\begin{multline}
 -\int_{V}\underline{\sigma}:\sum_{n}\overrightarrow{c^{n}}\nabla N_{,}^{n}\, dV+\int_{S_{T}}\vec{T}\cdot\sum_{n}\overrightarrow{c^{n}}N^{n}\, dS+\int_{V}\vec{f}\cdot\sum_{n}\overrightarrow{c^{n}}N^{n}\, dV\\
--\int_{V}\rho\sum_{m}\frac{\partial^{2}\overrightarrow{a^{m}}}{\partial t^{2}}N^{m}\cdot\sum_{n}\overrightarrow{c^{n}}N^{n}\ dV=0.\end{multline}
+-\int_{V}\rho\sum_{m}\frac{\partial^{2}\overrightarrow{a^{m}}}{\partial t^{2}}N^{m}\cdot\sum_{n}\overrightarrow{c^{n}}N^{n}\ dV=0.
+\end{multline}
 
 \end_inset
 
@@ -1112,8 +1204,10 @@
 \end_inset
 
 , so that
-\begin_inset Formula \begin{equation}
--\int_{V}\underline{\sigma}:\nabla N^{n}\, dV+\int_{S_{T}}\vec{T}N^{n}\, dS+\int_{V}\vec{f}N^{n}\, dV-\int_{V}\rho\sum_{m}\frac{\partial^{2}\overrightarrow{a^{m}}}{\partial t^{2}}N^{m}N^{n}\, dV=\vec{0}.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\underline{\sigma}:\nabla N^{n}\, dV+\int_{S_{T}}\vec{T}N^{n}\, dS+\int_{V}\vec{f}N^{n}\, dV-\int_{V}\rho\sum_{m}\frac{\partial^{2}\overrightarrow{a^{m}}}{\partial t^{2}}N^{m}N^{n}\, dV=\vec{0}.
+\end{equation}
 
 \end_inset
 
@@ -1134,9 +1228,11 @@
 \bar no
 \noun off
 \color none
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \vec{u}=u^{o}\overrightarrow{}\text{ on }S_{u},\text{ and}\\
-\underline{R}(\overrightarrow{u^{+}}-\overrightarrow{u^{-}})=\vec{d}\text{ on }S_{f},\end{gather}
+\underline{R}(\overrightarrow{u^{+}}-\overrightarrow{u^{-}})=\vec{d}\text{ on }S_{f},
+\end{gather}
 
 \end_inset
 
@@ -1158,20 +1254,24 @@
 \end_inset
 
  reduces to
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV=\vec{0}.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV=\vec{0}.
+\end{equation}
 
 \end_inset
 
 As a result, time-dependence only enters through the constitutive relationships
  and the loading conditions.
- We considers deformation at time 
+ We consider the deformation at time 
 \begin_inset Formula $t+\Delta t$
 \end_inset
 
 ,
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}(t+\Delta t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV=\vec{0}.\label{eq:elasticity:integral:quasistatic}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}(t+\Delta t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV=\vec{0}.\label{eq:elasticity:integral:quasistatic}
+\end{equation}
 
 \end_inset
 
@@ -1190,16 +1290,20 @@
 
 ).
  The residual is simply
-\begin_inset Formula \begin{equation}
-r_{i}^{n}=-\int_{V}\sigma_{ij}(t+\Delta t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+r_{i}^{n}=-\int_{V}\sigma_{ij}(t+\Delta t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV.
+\end{equation}
 
 \end_inset
 
 We employ numerical quadrature in the finite-element discretization and
  replace the integrals with sums over the cells and quadrature points,
-\begin_inset Formula \begin{multline}
+\begin_inset Formula 
+\begin{multline}
 r_{i}^{n}=-\sum_{\text{vol cells}}\sum_{\text{quad pts}}\sigma_{ij}(x_{q},t+\Delta t)N_{,j}^{n}(x_{q})\: w_{q}|J_{cell}(x_{q})|+\sum_{\text{vol cells}}\sum_{\text{quad pt}s}f_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|\\
-+\sum_{\text{tract cells}}\sum_{\text{quad pts}}T_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|,\end{multline}
++\sum_{\text{tract cells}}\sum_{\text{quad pts}}T_{i}(x_{q},t+\Delta t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|,
+\end{multline}
 
 \end_inset
 
@@ -1239,8 +1343,10 @@
 
 \begin_layout Standard
 In order to find the Jacobian of the system, we let
-\begin_inset Formula \begin{equation}
-\sigma_{ij}(t+\Delta t)=\sigma_{ij}(t)+d\sigma_{ij}(t).\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sigma_{ij}(t+\Delta t)=\sigma_{ij}(t)+d\sigma_{ij}(t).
+\end{equation}
 
 \end_inset
 
@@ -1248,8 +1354,10 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-\int_{V}d\sigma_{ij}(t)N_{j}^{n}\ dV=-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}d\sigma_{ij}(t)N_{j}^{n}\ dV=-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t+\Delta t)N^{n}\, dS+\int_{V}f_{i}(t+\Delta t)N^{n}\, dV
+\end{equation}
 
 \end_inset
 
@@ -1260,10 +1368,12 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 d\sigma_{ij}(t)=C_{ijkl}(t)d\varepsilon_{kl}(t)\\
 d\sigma_{ij}(t)=\frac{1}{2}C_{ijkl}(t)(du_{k.l}(t)+du_{l,k}(t))\\
-d\sigma_{ij}(t)=\frac{1}{2}C_{ijkl}(t)(\sum_{m}da_{k,l}^{m}(t)N^{m}+\sum_{m}da_{l,k}^{m}(t)N^{m})\end{gather}
+d\sigma_{ij}(t)=\frac{1}{2}C_{ijkl}(t)(\sum_{m}da_{k,l}^{m}(t)N^{m}+\sum_{m}da_{l,k}^{m}(t)N^{m})
+\end{gather}
 
 \end_inset
 
@@ -1272,8 +1382,10 @@
 \end_inset
 
  is a scalar, so it is symmetric,
-\begin_inset Formula \begin{equation}
-d\sigma_{ij}\phi_{i,j}=d\sigma_{ji}\phi_{j,i},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+d\sigma_{ij}\phi_{i,j}=d\sigma_{ji}\phi_{j,i},
+\end{equation}
 
 \end_inset
 
@@ -1282,20 +1394,26 @@
 \end_inset
 
  is symmetric, so
-\begin_inset Formula \begin{equation}
-d\sigma_{ij}\phi_{i,j}=d\sigma_{ij}\phi_{j,i},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+d\sigma_{ij}\phi_{i,j}=d\sigma_{ij}\phi_{j,i},
+\end{equation}
 
 \end_inset
 
 which means
-\begin_inset Formula \begin{equation}
-\phi_{i,j}=\phi_{j,i},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\phi_{i,j}=\phi_{j,i},
+\end{equation}
 
 \end_inset
 
 which we can write as
-\begin_inset Formula \begin{equation}
-\phi_{i,j}=\frac{1}{2}(\phi_{i,j}+\phi_{j,i}).\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\phi_{i,j}=\frac{1}{2}(\phi_{i,j}+\phi_{j,i}).
+\end{equation}
 
 \end_inset
 
@@ -1303,8 +1421,10 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-\sum_{n}c_{i}^{n}N_{,j}^{n}=\frac{1}{2}(\sum_{n}c_{i}^{n}N_{,j}^{n}+\sum_{n}c_{j}^{n}N_{,i}^{n}).\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sum_{n}c_{i}^{n}N_{,j}^{n}=\frac{1}{2}(\sum_{n}c_{i}^{n}N_{,j}^{n}+\sum_{n}c_{j}^{n}N_{,i}^{n}).
+\end{equation}
 
 \end_inset
 
@@ -1314,15 +1434,19 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-A_{ij}^{nm}=\int_{V}\frac{1}{4}C_{ijkl}(N_{,l}^{m}+N_{,k}^{m})(N_{,j}^{n}+N_{,i}^{n})\ dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{ij}^{nm}=\int_{V}\frac{1}{4}C_{ijkl}(N_{,l}^{m}+N_{,k}^{m})(N_{,j}^{n}+N_{,i}^{n})\ dV.
+\end{equation}
 
 \end_inset
 
 We employ numerical quadrature in the finite-element discretization and
  replace the integral with a sum over the cells and quadrature points,
-\begin_inset Formula \begin{equation}
-A_{ij}^{nm}=\sum_{\text{vol cells}}\sum_{\text{quad pts}}\frac{1}{4}C_{ijkl}(N_{,l}^{m}(x_{q})+N_{,k}^{m}(x_{q}))(N_{,j}^{n}(x_{q})+N_{,i}^{n}(x_{q}))w_{q}|J_{cell}(x_{q}).\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{ij}^{nm}=\sum_{\text{vol cells}}\sum_{\text{quad pts}}\frac{1}{4}C_{ijkl}(N_{,l}^{m}(x_{q})+N_{,k}^{m}(x_{q}))(N_{,j}^{n}(x_{q})+N_{,i}^{n}(x_{q}))w_{q}|J_{cell}(x_{q}).
+\end{equation}
 
 \end_inset
 
@@ -1343,8 +1467,10 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{equation}
--\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t)N^{n}\, dS+\int_{V}f_{i}(t)N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ dV=\vec{0}.\label{eq:elasticity:integral:dynamic:t}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t)N^{n}\, dS+\int_{V}f_{i}(t)N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ dV=\vec{0}.\label{eq:elasticity:integral:dynamic:t}
+\end{equation}
 
 \end_inset
 
@@ -1363,16 +1489,20 @@
 
 ).
  The residual is simply
-\begin_inset Formula \begin{equation}
-r_{i}^{n}=-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t)N^{n}\, dS+\int_{V}f_{i}(t)N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+r_{i}^{n}=-\int_{V}\sigma_{ij}(t)N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}(t)N^{n}\, dS+\int_{V}f_{i}(t)N^{n}\, dV-\int_{V}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ dV.
+\end{equation}
 
 \end_inset
 
 We employ numerical quadrature in the finite-element discretization and
  replace the integrals with sums over the cells and quadrature points,
-\begin_inset Formula \begin{multline}
+\begin_inset Formula 
+\begin{multline}
 r_{i}^{n}=-\sum_{\text{vol cells}}\sum_{\text{quad pts}}\sigma_{ij}(x_{q},t)N^{n}(x_{q})\: w_{q}|J_{cell}(x_{q})|+\sum_{\text{vol cells}}\sum_{\text{quad pt}s}f_{i}(x_{q},t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|\\
-+\sum_{\text{tract cells}}\sum_{\text{quad pts}}T_{i}(x_{q},t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|-\sum_{\text{vol cells}}\sum_{\text{quad pts}}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ w_{q|J_{cell}(x_{q})},\end{multline}
++\sum_{\text{tract cells}}\sum_{\text{quad pts}}T_{i}(x_{q},t)N^{n}(x_{q})\, w_{q}|J_{cell}(x_{q})|-\sum_{\text{vol cells}}\sum_{\text{quad pts}}\rho\sum_{m}\ddot{a}_{i}^{m}(t)N^{m}N^{n}\ w_{q|J_{cell}(x_{q})},
+\end{multline}
 
 \end_inset
 
@@ -1405,9 +1535,11 @@
 .
  Using the central difference method to approximate the acceleration (and
  velocity),
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \ddot{u}_{i}(t)=\frac{1}{\Delta t^{2}}\left(u_{i}(t+\Delta t)-2u_{i}(t)+u_{i}(t-\Delta t)\right)\\
-\dot{u}_{i}(t)=\frac{1}{2\Delta t}\left(u_{i}(t+\Delta t)-u_{i}(t-\Delta t)\right)\end{gather}
+\dot{u}_{i}(t)=\frac{1}{2\Delta t}\left(u_{i}(t+\Delta t)-u_{i}(t-\Delta t)\right)
+\end{gather}
 
 \end_inset
 
@@ -1420,10 +1552,12 @@
 \end_inset
 
  (for consistency with the displacement increment quasi-static formulation),
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 u_{i}(t+\Delta t)=u_{i}(t)+du_{i}(t),\\
 \ddot{u}_{i}(t)=\frac{1}{\Delta t^{2}}\left(du_{i}(t)-u_{i}(t)+u_{i}(t-\Delta t)\right),\\
-\dot{u}_{i}(t)=\frac{1}{2\Delta t}\left(du_{i}(t)+u_{i}(t)-u_{i}(t-\Delta t)\right).\end{gather}
+\dot{u}_{i}(t)=\frac{1}{2\Delta t}\left(du_{i}(t)+u_{i}(t)-u_{i}(t-\Delta t)\right).
+\end{gather}
 
 \end_inset
 
@@ -1435,15 +1569,19 @@
 \end_inset
 
  yields
-\begin_inset Formula \begin{multline}
+\begin_inset Formula 
+\begin{multline}
 \frac{1}{\Delta t^{2}}\int_{V}\rho\sum_{m}da_{i}^{m}(t)N^{m}N^{n}\ dV=-\int_{V}\sigma_{ij}N_{,j}^{n}\: dV+\int_{S_{T}}T_{i}N^{n}\, dS+\int_{V}f_{i}N^{n}\, dV\\
--\frac{1}{\Delta t^{2}}\int_{V}\rho\sum_{m}(a_{i}^{m}(t)-a_{i}^{m}(t-\Delta t))N^{m}N^{n}\ dV.\end{multline}
+-\frac{1}{\Delta t^{2}}\int_{V}\rho\sum_{m}(a_{i}^{m}(t)-a_{i}^{m}(t-\Delta t))N^{m}N^{n}\ dV.
+\end{multline}
 
 \end_inset
 
 Thus, the Jacobian for the system is
-\begin_inset Formula \begin{equation}
-A_{ij}^{nm}=\delta_{ij}\frac{1}{\Delta t^{2}}\int_{V}\rho N^{m}N^{n}\ dV,\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{ij}^{nm}=\delta_{ij}\frac{1}{\Delta t^{2}}\int_{V}\rho N^{m}N^{n}\ dV,
+\end{equation}
 
 \end_inset
 
@@ -1452,8 +1590,10 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-A_{ij}^{nm}=\delta_{ij}\frac{1}{\Delta t^{2}}\sum_{\text{vol cells}}\sum_{\text{quad pts}}\rho(x_{q})N^{m}(x_{q})N^{n}(x_{q}),\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{ij}^{nm}=\delta_{ij}\frac{1}{\Delta t^{2}}\sum_{\text{vol cells}}\sum_{\text{quad pts}}\rho(x_{q})N^{m}(x_{q})N^{n}(x_{q}),
+\end{equation}
 
 \end_inset
 
@@ -1536,10 +1676,12 @@
 .
  The Green-Lagrange strain provides a measure of the strain relative to
  the original, undeformed configuration.
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \varepsilon_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i}+u_{k,i}u_{k,j}),\text{ or}\\
 \varepsilon_{ij}=X_{ji}X_{ij}-\delta_{ij},\text{ where}\\
-X_{ij}=x_{i,j}(t)=\frac{\partial}{\partial x_{j}}(x_{i}(0)+u_{i}(t)),\end{gather}
+X_{ij}=x_{i,j}(t)=\frac{\partial}{\partial x_{j}}(x_{i}(0)+u_{i}(t)),
+\end{gather}
 
 \end_inset
 
@@ -1557,8 +1699,10 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-S_{ij}=C_{ijkl}\varepsilon_{kl},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+S_{ij}=C_{ijkl}\varepsilon_{kl},
+\end{equation}
 
 \end_inset
 
@@ -1600,8 +1744,10 @@
 
  strain.
  Using the definition of the Green-Lagrangian strain, we have
-\begin_inset Formula \begin{equation}
-\int_{V}S_{ij}\delta\varepsilon_{ij}\: dV=\int_{V}\frac{1}{2}S_{ij}(\delta u_{i,j}+\delta u_{j,i}+u_{k,i}\delta u_{k,j}+u_{k,j}\delta u_{k,i})\: dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{V}S_{ij}\delta\varepsilon_{ij}\: dV=\int_{V}\frac{1}{2}S_{ij}(\delta u_{i,j}+\delta u_{j,i}+u_{k,i}\delta u_{k,j}+u_{k,j}\delta u_{k,i})\: dV.
+\end{equation}
 
 \end_inset
 
@@ -1615,8 +1761,10 @@
 \end_inset
 
  strain) to zero yields the elastic term in the residual,
-\begin_inset Formula \begin{equation}
-r_{i}^{n}=\int_{V}S_{ij}(N_{,i}^{n}+(\sum_{m}a_{k}^{m}N_{,j}^{m})N_{,i}^{n})\: dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+r_{i}^{n}=\int_{V}S_{ij}(N_{,i}^{n}+(\sum_{m}a_{k}^{m}N_{,j}^{m})N_{,i}^{n})\: dV.
+\end{equation}
 
 \end_inset
 
@@ -1640,14 +1788,18 @@
 \end_inset
 
  and consider the first terms of the Taylor series expansion,
-\begin_inset Formula \begin{equation}
-\int_{v}S_{ij}(t+\Delta t)\delta\varepsilon_{ij}(t+\Delta t)\: dV=\int_{V}(S_{ij}(t)\delta\varepsilon_{ij}(t)+dS_{ij}(t)\delta\varepsilon_{ij}(t)+S_{ij}(t)d\delta\varepsilon_{ij}(t))\: dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\int_{v}S_{ij}(t+\Delta t)\delta\varepsilon_{ij}(t+\Delta t)\: dV=\int_{V}(S_{ij}(t)\delta\varepsilon_{ij}(t)+dS_{ij}(t)\delta\varepsilon_{ij}(t)+S_{ij}(t)d\delta\varepsilon_{ij}(t))\: dV.
+\end{equation}
 
 \end_inset
 
 We approximate the increment in the stress tensor using the elastic constants,
-\begin_inset Formula \begin{equation}
-dS_{ij}=C_{ijkl}d\varepsilon_{kl},\end{equation}
+\begin_inset Formula 
+\begin{equation}
+dS_{ij}=C_{ijkl}d\varepsilon_{kl},
+\end{equation}
 
 \end_inset
 
@@ -1660,8 +1812,10 @@
 \end_inset
 
  strain via
-\begin_inset Formula \begin{equation}
-d\delta\varepsilon_{ij}=\frac{1}{2}(du_{k,i}\delta u_{k,j}+du_{k,j}\delta u_{k,i}).\end{equation}
+\begin_inset Formula 
+\begin{equation}
+d\delta\varepsilon_{ij}=\frac{1}{2}(du_{k,i}\delta u_{k,j}+du_{k,j}\delta u_{k,i}).
+\end{equation}
 
 \end_inset
 
@@ -1677,8 +1831,10 @@
 \end_inset
 
  strains, we have
-\begin_inset Formula \begin{equation}
-A_{ij}^{nm}=\int_{V}\frac{1}{4}C_{ijkl}(N_{,k}^{m}+(\sum_{r}a_{p}^{r}N_{,l}^{r})N_{,k}^{m})(N_{,i}^{n}+(\sum_{r}a_{p}^{r}N_{,j}^{r})N_{,i}^{n})+\frac{1}{2}S_{kl}N_{,l}^{m}N_{,l}^{n}\delta_{ij}\: dV.\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{ij}^{nm}=\int_{V}\frac{1}{4}C_{ijkl}(N_{,k}^{m}+(\sum_{r}a_{p}^{r}N_{,l}^{r})N_{,k}^{m})(N_{,i}^{n}+(\sum_{r}a_{p}^{r}N_{,j}^{r})N_{,i}^{n})+\frac{1}{2}S_{kl}N_{,l}^{m}N_{,l}^{n}\delta_{ij}\: dV.
+\end{equation}
 
 \end_inset
 

Modified: short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx	2011-06-10 01:29:23 UTC (rev 18574)
+++ short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx	2011-06-10 05:16:17 UTC (rev 18575)
@@ -1,5 +1,5 @@
-#LyX 1.6.5 created this file. For more info see http://www.lyx.org/
-\lyxformat 345
+#LyX 2.0 created this file. For more info see http://www.lyx.org/
+\lyxformat 413
 \begin_document
 \begin_header
 \textclass book
@@ -7,18 +7,26 @@
 \newcommand\Prefix[3]{\vphantom{#3}#1#2#3}
 \end_preamble
 \use_default_options false
+\maintain_unincluded_children false
 \language english
+\language_package default
 \inputencoding latin1
+\fontencoding global
 \font_roman default
 \font_sans default
 \font_typewriter default
 \font_default_family default
+\use_non_tex_fonts false
 \font_sc false
 \font_osf false
 \font_sf_scale 100
 \font_tt_scale 100
 
 \graphics default
+\default_output_format default
+\output_sync 0
+\bibtex_command default
+\index_command default
 \paperfontsize default
 \spacing single
 \use_hyperref false
@@ -26,9 +34,18 @@
 \use_geometry true
 \use_amsmath 1
 \use_esint 0
+\use_mhchem 1
+\use_mathdots 1
 \cite_engine basic
 \use_bibtopic false
+\use_indices false
 \paperorientation portrait
+\suppress_date false
+\use_refstyle 0
+\index Index
+\shortcut idx
+\color #008000
+\end_index
 \leftmargin 1in
 \topmargin 1in
 \rightmargin 1in
@@ -36,15 +53,16 @@
 \secnumdepth 3
 \tocdepth 3
 \paragraph_separation indent
-\defskip medskip
+\paragraph_indentation default
 \quotes_language english
 \papercolumns 1
 \papersides 2
 \paperpagestyle default
 \tracking_changes false
 \output_changes false
-\author "" 
-\author "" 
+\html_math_output 0
+\html_css_as_file 0
+\html_be_strict false
 \end_header
 
 \begin_body
@@ -519,7 +537,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="6" columns="4">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="1.5in">
 <column alignment="center" valignment="top" width="1.8in">
 <column alignment="center" valignment="top" width="1.5in">
@@ -857,7 +875,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="3" columns="4">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="1.25in">
 <column alignment="center" valignment="top" width="0.5in">
 <column alignment="center" valignment="top" width="1.25in">
@@ -1080,8 +1098,10 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{gather}
-\sigma_{ij}=C_{ijkl}\left(\epsilon_{kl}-\epsilon_{kl}^{I}\right)+\sigma_{ij}^{I}\,,\label{eq:1}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\sigma_{ij}=C_{ijkl}\left(\epsilon_{kl}-\epsilon_{kl}^{I}\right)+\sigma_{ij}^{I}\,,\label{eq:1}
+\end{gather}
 
 \end_inset
 
@@ -1096,8 +1116,10 @@
  of anisotropic elasticity.
  Representing the stress and strain in terms of vectors, the constitutive
  relation may be written
-\begin_inset Formula \begin{gather}
-\overrightarrow{\sigma}=\underline{C}\left(\vec{\epsilon}-\vec{\epsilon}^{I}\right)+\vec{\sigma}^{I},\label{eq:2}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\overrightarrow{\sigma}=\underline{C}\left(\vec{\epsilon}-\vec{\epsilon}^{I}\right)+\vec{\sigma}^{I},\label{eq:2}
+\end{gather}
 
 \end_inset
 
@@ -1105,14 +1127,17 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \underline{C}=\left[\begin{array}{cccccc}
 C_{1111} & C_{1122} & C_{1133} & C_{1112} & C_{1123} & C_{1113}\\
 C_{1122} & C_{2222} & C_{2233} & C_{2212} & C_{2223} & C_{2213}\\
 C_{1133} & C_{2233} & C_{3333} & C_{3312} & C_{3323} & C_{3313}\\
 C_{1112} & C_{2212} & C_{3312} & C_{1212} & C_{1223} & C_{1213}\\
 C_{1123} & C_{2223} & C_{3323} & C_{1223} & C_{2323} & C_{2313}\\
-C_{1113} & C_{2213} & C_{3313} & C_{1213} & C_{2313} & C_{1313}\end{array}\right]\:.\label{eq:3}\end{gather}
+C_{1113} & C_{2213} & C_{3313} & C_{1213} & C_{2313} & C_{1313}
+\end{array}\right]\:.\label{eq:3}
+\end{gather}
 
 \end_inset
 
@@ -1143,10 +1168,13 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
+\begin_inset Formula 
+\begin{equation}
 \begin{aligned}\mu= & \rho v_{s}^{2}\\
-\lambda= & \rho v_{p}^{2}-2\mu\end{aligned}
-\label{eq:4}\end{equation}
+\lambda= & \rho v_{p}^{2}-2\mu
+\end{aligned}
+\label{eq:4}
+\end{equation}
 
 \end_inset
 
@@ -1177,7 +1205,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="4" columns="2">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -1322,15 +1350,19 @@
 \end_inset
 
 , so we have
-\begin_inset Formula \begin{equation}
-\sigma_{11}=(\lambda+2\mu)\left(\epsilon_{11}-\epsilon_{11}^{I}\right)+\sigma_{11}^{I},\label{eq:5}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sigma_{11}=(\lambda+2\mu)\left(\epsilon_{11}-\epsilon_{11}^{I}\right)+\sigma_{11}^{I},\label{eq:5}
+\end{equation}
 
 \end_inset
 
 with
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \sigma_{22}=\sigma_{33}=\lambda\left(\epsilon_{11}-\epsilon_{11}^{I}\right)+\sigma_{22}^{I},\nonumber \\
-\sigma_{12}=\sigma_{23}=\sigma_{13}=0.\label{eq:6}\end{gather}
+\sigma_{12}=\sigma_{23}=\sigma_{13}=0.\label{eq:6}
+\end{gather}
 
 \end_inset
 
@@ -1356,15 +1388,19 @@
 \end_inset
 
 , so we have
-\begin_inset Formula \begin{equation}
-\sigma_{11}=\frac{\mu(3\lambda+2\mu)}{\lambda+\mu}\left(\epsilon_{11}-\epsilon_{11}^{I}\right)+\sigma_{11}^{I},\label{eq:7}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sigma_{11}=\frac{\mu(3\lambda+2\mu)}{\lambda+\mu}\left(\epsilon_{11}-\epsilon_{11}^{I}\right)+\sigma_{11}^{I},\label{eq:7}
+\end{equation}
 
 \end_inset
 
 with
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \epsilon_{22}=\epsilon_{33}=-\frac{\lambda}{2(\lambda+\mu)}\epsilon_{11}+\epsilon_{22}^{I},\label{eq:8}\\
-\sigma_{22}=\sigma_{33}=\sigma_{12}=\sigma_{23}=\sigma_{13}=0.\nonumber \end{gather}
+\sigma_{22}=\sigma_{33}=\sigma_{12}=\sigma_{23}=\sigma_{13}=0.\nonumber 
+\end{gather}
 
 \end_inset
 
@@ -1377,20 +1413,26 @@
 
 \begin_layout Standard
 In 2D we can write Hooke's law as
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \left[\begin{array}{c}
 \sigma_{11}\\
 \sigma_{22}\\
-\sigma_{12}\end{array}\right]=\left[\begin{array}{ccc}
+\sigma_{12}
+\end{array}\right]=\left[\begin{array}{ccc}
 C_{1111} & C_{1122} & C_{1112}\\
 C_{1122} & C_{2222} & C_{2212}\\
-C_{1112} & C_{2212} & C_{1212}\end{array}\right]\left[\begin{array}{c}
+C_{1112} & C_{2212} & C_{1212}
+\end{array}\right]\left[\begin{array}{c}
 \epsilon_{11}-\epsilon_{11}^{I}\\
 \epsilon_{22}-\epsilon_{22}^{I}\\
-\epsilon_{12}-\epsilon_{12}^{I}\end{array}\right]+\left[\begin{array}{c}
+\epsilon_{12}-\epsilon_{12}^{I}
+\end{array}\right]+\left[\begin{array}{c}
 \sigma_{11}^{I}\\
 \sigma_{22}^{I}\\
-\sigma_{12}^{I}\end{array}\right]\:.\label{eq:9}\end{gather}
+\sigma_{12}^{I}
+\end{array}\right]\:.\label{eq:9}
+\end{gather}
 
 \end_inset
 
@@ -1411,20 +1453,26 @@
 \end_inset
 
  and plane strain conditions apply, so we have 
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \left[\begin{array}{c}
 \sigma_{11}\\
 \sigma_{22}\\
-\sigma_{12}\end{array}\right]=\left[\begin{array}{ccc}
+\sigma_{12}
+\end{array}\right]=\left[\begin{array}{ccc}
 \lambda+2\mu & \lambda & 0\\
 \lambda & \lambda+2\mu & 0\\
-0 & 0 & 2\mu\end{array}\right]\left[\begin{array}{c}
+0 & 0 & 2\mu
+\end{array}\right]\left[\begin{array}{c}
 \epsilon_{11}-\epsilon_{11}^{I}\\
 \epsilon_{22}-\epsilon_{22}^{I}\\
-\epsilon_{12}-\epsilon_{12}^{I}\end{array}\right]+\left[\begin{array}{c}
+\epsilon_{12}-\epsilon_{12}^{I}
+\end{array}\right]+\left[\begin{array}{c}
 \sigma_{11}^{I}\\
 \sigma_{22}^{I}\\
-\sigma_{12}^{I}\end{array}\right]\:.\label{eq:10}\end{gather}
+\sigma_{12}^{I}
+\end{array}\right]\:.\label{eq:10}
+\end{gather}
 
 \end_inset
 
@@ -1445,28 +1493,37 @@
 \end_inset
 
  and plane stress conditions apply, so we have
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \left[\begin{array}{c}
 \sigma_{11}\\
 \sigma_{22}\\
-\sigma_{12}\end{array}\right]=\left[\begin{array}{ccc}
+\sigma_{12}
+\end{array}\right]=\left[\begin{array}{ccc}
 \frac{4\mu(\lambda+\mu)}{\lambda+2\mu} & \frac{2\mu\lambda}{\lambda+2\mu} & 0\\
 \frac{2\mu\lambda}{\lambda+2} & \frac{4\mu(\lambda+\mu)}{\lambda+2\mu} & 0\\
-0 & 0 & 2\mu\end{array}\right]\left[\begin{array}{c}
+0 & 0 & 2\mu
+\end{array}\right]\left[\begin{array}{c}
 \epsilon_{11}-\epsilon_{11}^{I}\\
 \epsilon_{22}-\epsilon_{22}^{I}\\
-\epsilon_{12}-\epsilon_{12}^{I}\end{array}\right]+\left[\begin{array}{c}
+\epsilon_{12}-\epsilon_{12}^{I}
+\end{array}\right]+\left[\begin{array}{c}
 \sigma_{11}^{I}\\
 \sigma_{22}^{I}\\
-\sigma_{12}^{I}\end{array}\right]\:,\label{eq:11}\end{gather}
+\sigma_{12}^{I}
+\end{array}\right]\:,\label{eq:11}
+\end{gather}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{equation}
+\begin_inset Formula 
+\begin{equation}
 \begin{gathered}\epsilon_{33}=-\frac{\lambda}{\lambda+2\mu}(\epsilon_{11}+\epsilon_{22})+\epsilon_{33}^{I}\\
-\epsilon_{13}=\epsilon_{23}=0\,.\end{gathered}
-\label{eq:12}\end{equation}
+\epsilon_{13}=\epsilon_{23}=0\,.
+\end{gathered}
+\label{eq:12}
+\end{equation}
 
 \end_inset
 
@@ -1490,14 +1547,17 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \underline{C}=\left[\begin{array}{cccccc}
 \lambda+2\mu & \lambda & \lambda & 0 & 0 & 0\\
 \lambda & \lambda+2\mu & \lambda & 0 & 0 & 0\\
 \lambda & \lambda & \lambda+2\mu & 0 & 0 & 0\\
 0 & 0 & 0 & 2\mu & 0 & 0\\
 0 & 0 & 0 & 0 & 2\mu & 0\\
-0 & 0 & 0 & 0 & 0 & 2\mu\end{array}\right]\:.\label{eq:13}\end{gather}
+0 & 0 & 0 & 0 & 0 & 2\mu
+\end{array}\right]\:.\label{eq:13}
+\end{gather}
 
 \end_inset
 
@@ -1505,11 +1565,17 @@
 \end_layout
 
 \begin_layout Section
+\begin_inset CommandInset label
+LatexCommand label
+name "sec:Viscoelastic-Materials"
+
+\end_inset
+
 Viscoelastic Materials
 \end_layout
 
 \begin_layout Standard
-At present, there are four viscoelastic material models available in PyLith
+At present, there are five viscoelastic material models available in PyLith
  (Table 
 \begin_inset CommandInset ref
 LatexCommand ref
@@ -1526,8 +1592,15 @@
 
 ).
  Future code versions may include alternative formulations for the various
- material models, so that users may use the most efficient formulation for
- a particular problem.
+ material models (
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "cha:Alternative-Formulations"
+
+\end_inset
+
+), so that users may use the most efficient formulation for a particular
+ problem.
  Note that both 2D and 3D viscoelastic models are described, but we present
  below only the 3D formulations.
  The 2D formulations are easily obtained from the plane strain definition.
@@ -1561,8 +1634,8 @@
 
 
 \begin_inset Tabular
-<lyxtabular version="3" rows="5" columns="2">
-<features>
+<lyxtabular version="3" rows="6" columns="2">
+<features tabularvalignment="middle">
 <column alignment="left" valignment="top" width="2.85in">
 <column alignment="center" valignment="top" width="2.47in">
 <row>
@@ -1614,6 +1687,26 @@
 \begin_inset Text
 
 \begin_layout Plain Layout
+GenMaxwellPlaneStrain
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Plane strain generalized Maxwell material (3 Maxwell models in parallel)
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
 MaxwellIsotropic3D
 \end_layout
 
@@ -1702,7 +1795,6 @@
 \begin_inset Caption
 
 \begin_layout Plain Layout
-\noindent
 \begin_inset CommandInset label
 LatexCommand label
 name "fig:Material-models"
@@ -1758,8 +1850,10 @@
  indices indicate summation over the range of the index.
  We also make frequent use of the scalar inner product.
  The scalar inner product of two second-order tensors may be written
-\begin_inset Formula \begin{gather}
-\underline{a}\cdot\underline{b}=a_{ij}b_{ij}\,.\label{eq:14}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\underline{a}\cdot\underline{b}=a_{ij}b_{ij}\,.\label{eq:14}
+\end{gather}
 
 \end_inset
 
@@ -1775,8 +1869,10 @@
 \end_inset
 
 :
-\begin_inset Formula \begin{gather}
-P=\frac{\sigma_{ii}}{3}\,,\,\,\,\,\theta=\frac{\epsilon_{ii}}{3}\,,\label{eq:15}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+P=\frac{\sigma_{ii}}{3}\,,\,\,\,\,\theta=\frac{\epsilon_{ii}}{3}\,,\label{eq:15}
+\end{gather}
 
 \end_inset
 
@@ -1790,8 +1886,10 @@
 
  represent the trace of the stress and strain tensors, respectively.
  We then define the deviatoric components of stress and strain as
-\begin_inset Formula \begin{gather}
-S_{ij}=\sigma_{ij}-P\delta_{ij}\,,\,\,\,\, e_{ij}=\epsilon_{ij}-\theta\delta_{ij}\,,\label{eq:16}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+S_{ij}=\sigma_{ij}-P\delta_{ij}\,,\,\,\,\, e_{ij}=\epsilon_{ij}-\theta\delta_{ij}\,,\label{eq:16}
+\end{gather}
 
 \end_inset
 
@@ -1817,21 +1915,25 @@
 \end_inset
 
 , as
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \overline{\sigma}=\sqrt{\frac{3}{2}\underline{S}\cdot\underline{S}}\,\,\nonumber \\
 J_{2}^{\prime}=\frac{1}{2}\underline{S}\cdot\underline{S}\,.\label{eq:17}\\
 \overline{e}=\sqrt{\frac{2}{3}\underline{e}\cdot\underline{e}}\,\,\nonumber \\
-L_{2}^{\prime}=\frac{1}{2}\underline{e}\cdot\underline{e}\,\,\nonumber \end{gather}
+L_{2}^{\prime}=\frac{1}{2}\underline{e}\cdot\underline{e}\,\,\nonumber 
+\end{gather}
 
 \end_inset
 
 Due to the symmetry of the stress and strain tensors, it is sometimes convenient
  to represent them as vectors:
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \overrightarrow{\sigma^{T}}=\left[\begin{array}{cccccc}
 \sigma_{11} & \sigma_{22} & \sigma_{33} & \sigma_{12} & \sigma_{23} & \sigma_{31}\end{array}\right]\label{eq:18}\\
 \overrightarrow{\epsilon^{T}}=\left[\begin{array}{cccccc}
-\epsilon_{11} & \epsilon_{22} & \epsilon_{33} & \epsilon_{12} & \epsilon_{23} & \epsilon_{31}\end{array}\right]\:.\nonumber \end{gather}
+\epsilon_{11} & \epsilon_{22} & \epsilon_{33} & \epsilon_{12} & \epsilon_{23} & \epsilon_{31}\end{array}\right]\:.\nonumber 
+\end{gather}
 
 \end_inset
 
@@ -1880,7 +1982,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="3" columns="3">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="0">
 <column alignment="center" valignment="top" width="0">
 <column alignment="center" valignment="top" width="0">
@@ -1996,6 +2098,121 @@
 \end_layout
 
 \begin_layout Subsection
+Linear Viscoelastic Models
+\end_layout
+
+\begin_layout Standard
+Linear viscoelastic models are obtained by various combinations of a linear
+ elastic spring and a linear viscous dashpot in series or parallel.
+ The simplest example is probably the linear Maxwell model, which consists
+ of a spring in series with a dashpot, as shown in 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "fig:Material-models"
+
+\end_inset
+
+.
+ For a one-dimensional model, the response is given by
+\begin_inset Formula 
+\begin{equation}
+\frac{d\epsilon_{Total}}{dt}=\frac{d\epsilon_{D}}{dt}+\frac{d\epsilon_{S}}{dt}=\frac{\sigma}{\eta}+\frac{1}{E}\frac{d\sigma}{dt}\:,
+\end{equation}
+
+\end_inset
+
+where 
+\begin_inset Formula $\epsilon_{Total}$
+\end_inset
+
+ is the total strain, 
+\begin_inset Formula $\epsilon_{D}$
+\end_inset
+
+ is the strain in the dashpot, 
+\begin_inset Formula $\epsilon_{S}$
+\end_inset
+
+ is the strain in the spring, 
+\begin_inset Formula $\sigma$
+\end_inset
+
+ is the stress, 
+\begin_inset Formula $\eta$
+\end_inset
+
+ is the viscosity of the dashpot, and 
+\begin_inset Formula $E$
+\end_inset
+
+ is the spring constant.
+ When a Maxwell material is subjected to constant strain, the stresses relax
+ exponentially with time.
+ When a Maxwell material is subjected to a constant stress, there is an
+ immediate elastic strain, corresponding to the response of the spring,
+ and a viscous strain that increases linearly with time.
+ Since the strain response is unbounded, the Maxwell model actually represents
+ a fluid.
+\end_layout
+
+\begin_layout Standard
+Another simple model is the Kelvin-Voigt model, which consists of a spring
+ in parallel with a dashpot.
+ In this case, the one-dimensional response is given by
+\begin_inset Formula 
+\begin{equation}
+\sigma\left(t\right)=E\epsilon\left(t\right)+\eta\frac{d\epsilon\left(t\right)}{dt}\:.
+\end{equation}
+
+\end_inset
+
+As opposed to the Maxwell model, which represents a fluid, the Kelvin-Voigt
+ model represents a solid undergoing reversible, viscoelastic strain.
+ If the material is subjected to a constant stress, it deforms at a decreasing
+ rate, gradually approaching the strain that would occur for a purely elastic
+ material.
+ When the stress is released, the material gradually relaxes back to its
+ undeformed state.
+\end_layout
+
+\begin_layout Standard
+The most general form of linear viscoelastic model is the generalized Maxwell
+ model, which consists of a spring in parallel with a number of Maxwell
+ models (last model in 
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "fig:Material-models"
+
+\end_inset
+
+).
+ Using this model, it is possible to represent a number of simpler viscoelastic
+ models.
+ For example, a simple Maxwell model is obtained by setting the elastic
+ constants of all springs to zero, with the exception of the spring contained
+ in the first Maxwell model (
+\begin_inset Formula $\mu_{1}$
+\end_inset
+
+).
+ Similarly, the Kelvin-Voigt model may be obtained by setting the elastic
+ constants 
+\begin_inset Formula $\mu_{2}=\mu_{3}=0$
+\end_inset
+
+, and setting 
+\begin_inset Formula $\mu_{1}=\infty$
+\end_inset
+
+.
+ Due to its flexibility, the generalized Maxwell model is the model used
+ in PyLith, although at present the number of Maxwell models in parallel
+ is limited to three.
+ We have also included the Maxwell model, since this is more computationally
+ efficient if the additional models are not needed.
+\end_layout
+
+\begin_layout Subsection
 Formulation for Generalized Maxwell Models
 \begin_inset CommandInset label
 LatexCommand label
@@ -2007,9 +2224,9 @@
 \end_layout
 
 \begin_layout Standard
-The generalized Maxwell viscoelastic model consists of a number of Maxwell
- linear viscoelastic models in parallel with a spring, as shown in Figure
- 
+As described above, the generalized Maxwell viscoelastic model consists
+ of a number of Maxwell linear viscoelastic models in parallel with a spring,
+ as shown in Figure 
 \begin_inset CommandInset ref
 LatexCommand ref
 reference "fig:Material-models"
@@ -2084,8 +2301,10 @@
 As for all our viscoelastic models, the volumetric strain is completely
  elastic, and the viscoelastic deformation may be expressed purely in terms
  of the deviatoric components:
-\begin_inset Formula \begin{equation}
-\underline{S}=2\mu_{tot}\left[\mu_{0}\underline{e}+\sum_{i=1}^{N}\mu_{i}\underline{q}^{i}-\underline{e}^{I}\right]+\underline{S}^{I}\,;\; P=3K\left(\theta-\theta^{I}\right)+P^{I}\,,\label{eq:19}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{S}=2\mu_{tot}\left[\mu_{0}\underline{e}+\sum_{i=1}^{N}\mu_{i}\underline{q}^{i}-\underline{e}^{I}\right]+\underline{S}^{I}\,;\; P=3K\left(\theta-\theta^{I}\right)+P^{I}\,,\label{eq:19}
+\end{equation}
 
 \end_inset
 
@@ -2102,8 +2321,10 @@
 \end_inset
 
  follows the evolution equations
-\begin_inset Formula \begin{equation}
-\underline{\dot{q}}^{i}+\frac{1}{\lambda_{i}}\underline{q}^{i}=\underline{\dot{e}}.\label{eq:20}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{\dot{q}}^{i}+\frac{1}{\lambda_{i}}\underline{q}^{i}=\underline{\dot{e}}.\label{eq:20}
+\end{equation}
 
 \end_inset
 
@@ -2132,8 +2353,10 @@
 
 , is measured.
  For a linear material we obtain:
-\begin_inset Formula \begin{equation}
-\underline{S}\left(t\right)=2\mu\left(t\right)\left(\underline{e}_{0}-\underline{e}^{I}\right)+\underline{S}^{I}\,,\label{eq:21}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{S}\left(t\right)=2\mu\left(t\right)\left(\underline{e}_{0}-\underline{e}^{I}\right)+\underline{S}^{I}\,,\label{eq:21}
+\end{equation}
 
 \end_inset
 
@@ -2144,20 +2367,26 @@
  is the shear relaxation modulus function.
  Using linearity and superposition for an arbitrary state of strain yields
  an integral equation:
-\begin_inset Formula \begin{equation}
-\underline{S}\left(t\right)=\intop_{-\infty}^{t}\mu\left(t-\tau\right)\underline{\dot{e}}\, d\tau\,.\label{eq:22}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{S}\left(t\right)=\intop_{-\infty}^{t}\mu\left(t-\tau\right)\underline{\dot{e}}\, d\tau\,.\label{eq:22}
+\end{equation}
 
 \end_inset
 
 If we assume the modulus function in Prony series form we obtain
-\begin_inset Formula \begin{equation}
-\mu\left(t\right)=\mu_{tot}\left(\mu_{0}+\sum_{i=1}^{N}\mu_{i}\exp\frac{-t}{\lambda_{i}}\right)\,,\label{eq:23}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\mu\left(t\right)=\mu_{tot}\left(\mu_{0}+\sum_{i=1}^{N}\mu_{i}\exp\frac{-t}{\lambda_{i}}\right)\,,\label{eq:23}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{equation}
-\mu_{0}+\sum_{i=1}^{N}\mu_{i}=1\,.\label{eq:24}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\mu_{0}+\sum_{i=1}^{N}\mu_{i}=1\,.\label{eq:24}
+\end{equation}
 
 \end_inset
 
@@ -2172,22 +2401,12 @@
 \end_layout
 
 \begin_layout Standard
-Considering a single Maxwell model, the relaxation modulus function is
-\begin_inset Formula \begin{equation}
-\mu\left(t\right)=\mu_{tot}\left(\mu_{0}+\mu_{1}\exp\frac{-t}{\lambda_{1}}\right)\,,\label{eq:25}\end{equation}
-
-\end_inset
-
-where
-\begin_inset Formula \begin{equation}
-\mu_{0}+\mu_{1}=1\,.\label{eq:26}\end{equation}
-
-\end_inset
-
 If we assume the material is undisturbed until a strain is suddenly applied
  at time zero, we can divide the integral into
-\begin_inset Formula \begin{equation}
-\intop_{-\infty}^{t}\left(\cdot\right)\, d\tau=\intop_{-\infty}^{0^{-}}\left(\cdot\right)\, d\tau+\intop_{0^{-}}^{0^{+}}\left(\cdot\right)\, d\tau+\intop_{0^{+}}^{t}\left(\cdot\right)\, d\tau\,.\label{eq:27}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\intop_{-\infty}^{t}\left(\cdot\right)\, d\tau=\intop_{-\infty}^{0^{-}}\left(\cdot\right)\, d\tau+\intop_{0^{-}}^{0^{+}}\left(\cdot\right)\, d\tau+\intop_{0^{+}}^{t}\left(\cdot\right)\, d\tau\,.\label{eq:27}
+\end{equation}
 
 \end_inset
 
@@ -2205,8 +2424,10 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{equation}
-\underline{S}\left(t\right)=2\mu\left(t\right)\left(\underline{e}_{0}-\underline{e}^{I}\right)+\underline{S}^{I}+2\int_{0}^{t}\mu\left(t-\tau\right)\underline{\dot{e}}\left(\tau\right)\, d\tau\,,\label{eq:28}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{S}\left(t\right)=2\mu\left(t\right)\left(\underline{e}_{0}-\underline{e}^{I}\right)+\underline{S}^{I}+2\int_{0}^{t}\mu\left(t-\tau\right)\underline{\dot{e}}\left(\tau\right)\, d\tau\,,\label{eq:28}
+\end{equation}
 
 \end_inset
 
@@ -2217,7 +2438,7 @@
 Substituting Equation 
 \begin_inset CommandInset ref
 LatexCommand ref
-reference "eq:25"
+reference "eq:23"
 
 \end_inset
 
@@ -2229,12 +2450,14 @@
 \end_inset
 
 , we obtain
-\begin_inset Formula \begin{equation}
-\underline{S}\left(t\right)=2\mu_{tot}\left[\mu_{0}\underline{e}\left(t\right)+\mu_{1}\exp\frac{-t}{\lambda_{1}}\left(\underline{e}_{0}+\intop_{0}^{t}\exp\frac{t}{\lambda_{1}}\underline{\dot{e}}\left(\tau\right)\, d\tau\right)-\underline{e}^{I}\right]+\underline{S}^{I}\,.\label{eq:29}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{S}\left(t\right)=2\mu_{tot}\left[\mu_{0}\underline{e}\left(t\right)+\sum_{i=1}^{N}\mu_{i}\exp\frac{-t}{\lambda_{i}}\left(\underline{e}_{0}+\intop_{0}^{t}\exp\frac{t}{\lambda_{i}}\underline{\dot{e}}\left(\tau\right)\, d\tau\right)-\underline{e}^{I}\right]+\underline{S}^{I}\,.\label{eq:29}
+\end{equation}
 
 \end_inset
 
-We then split the integral into two ranges: from 0 to 
+We then split each integral into two ranges: from 0 to 
 \begin_inset Formula $t_{n}$
 \end_inset
 
@@ -2246,54 +2469,70 @@
 \begin_inset Formula $t$
 \end_inset
 
-, and define the integral as
-\begin_inset Formula \begin{equation}
-\underline{i}^{1}\left(t\right)=\intop_{0}^{t}\exp\frac{\tau}{\lambda_{1}}\underline{\dot{e}}\left(\tau\right)\, d\tau\,.\label{eq:30}\end{equation}
+, and define each integral as
+\begin_inset Formula 
+\begin{equation}
+\underline{i}_{i}^{1}\left(t\right)=\intop_{0}^{t}\exp\frac{\tau}{\lambda_{i}}\underline{\dot{e}}\left(\tau\right)\, d\tau\,.\label{eq:30}
+\end{equation}
 
 \end_inset
 
 The integral then becomes
-\begin_inset Formula \begin{equation}
-\underline{i}^{1}\left(t\right)=\underline{i}^{1}\left(t_{n}\right)+\intop_{t_{n}}^{t}\exp\frac{\tau}{\lambda_{1}}\underline{\dot{e}}\left(\tau\right)\, d\tau\,.\label{eq:31}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{i}_{i}^{1}\left(t\right)=\underline{i}_{i}^{1}\left(t_{n}\right)+\intop_{t_{n}}^{t}\exp\frac{\tau}{\lambda_{i}}\underline{\dot{e}}\left(\tau\right)\, d\tau\,.\label{eq:31}
+\end{equation}
 
 \end_inset
 
 Including the negative exponential multiplier:
-\begin_inset Formula \begin{equation}
-\underline{h}^{1}\left(t\right)=\exp\frac{-t}{\lambda_{1}}\underline{i}^{1}\,.\label{eq:32}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{h}_{i}^{1}\left(t\right)=\exp\frac{-t}{\lambda_{i}}\underline{i}_{i}^{1}\,.\label{eq:32}
+\end{equation}
 
 \end_inset
 
 Then
-\begin_inset Formula \begin{equation}
-\underline{h}^{1}\left(t\right)=\exp\frac{-\Delta t}{\lambda_{1}}\underline{h}_{n}^{1}+\Delta\underline{h}\,,\label{eq:33}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{h}_{i}^{1}\left(t\right)=\exp\frac{-\Delta t}{\lambda_{i}}\underline{h}_{i}^{1}\left(t_{n}\right)+\Delta\underline{h}_{i}\,,\label{eq:33}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{equation}
-\Delta\underline{h}=\exp\frac{-t}{\lambda_{1}}\intop_{t_{n}}^{t}\exp\frac{\tau}{\lambda_{1}}\underline{\dot{e}}\left(\tau\right)\, d\tau\,.\label{eq:34}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta\underline{h}_{i}=\exp\frac{-t}{\lambda_{i}}\intop_{t_{n}}^{t}\exp\frac{\tau}{\lambda_{i}}\underline{\dot{e}}\left(\tau\right)\, d\tau\,.\label{eq:34}
+\end{equation}
 
 \end_inset
 
 Approximating the strain rate as constant over each time step, the solution
  may be found as
-\begin_inset Formula \begin{equation}
-\Delta\underline{h}=\frac{\lambda_{1}}{\Delta t}\left(1-\exp\frac{-\Delta t}{\lambda_{1}}\right)\left(\underline{e}-\underline{e}_{n}\right)=\Delta h\left(\underline{e}-\underline{e}_{n}\right)\,.\label{eq:35}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta\underline{h}_{i}=\frac{\lambda_{i}}{\Delta t}\left(1-\exp\frac{-\Delta t}{\lambda_{i}}\right)\left(\underline{e}-\underline{e}_{n}\right)=\Delta h_{i}\left(\underline{e}-\underline{e}_{n}\right)\,.\label{eq:35}
+\end{equation}
 
 \end_inset
 
 The approximation is singular for zero time steps, but a series expansion
  may be used for small time-step sizes:
-\begin_inset Formula \begin{equation}
-\Delta h\approx1-\frac{1}{2}\left(\frac{\Delta t}{\lambda_{1}}\right)+\frac{1}{3!}\left(\frac{\Delta t}{\lambda_{1}}\right)^{2}-\frac{1}{4!}\left(\frac{\Delta t}{\lambda_{1}}\right)^{3}+\cdots\,.\label{eq:36}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta h_{i}\approx1-\frac{1}{2}\left(\frac{\Delta t}{\lambda_{i}}\right)+\frac{1}{3!}\left(\frac{\Delta t}{\lambda_{i}}\right)^{2}-\frac{1}{4!}\left(\frac{\Delta t}{\lambda_{i}}\right)^{3}+\cdots\,.\label{eq:36}
+\end{equation}
 
 \end_inset
 
 This converges with only a few terms.
  With this formulation, the constitutive relation now has the simple form:
-\begin_inset Formula \begin{equation}
-\underline{S}\left(t\right)=2\mu_{tot}\left(\mu_{0}\underline{e}\left(t\right)+\mu_{1}\underline{h}^{1}\left(t\right)-\underline{e}^{I}\right)+\underline{S}^{I}\,.\label{eq:37}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\underline{S}\left(t\right)=2\mu_{tot}\left(\mu_{0}\underline{e}\left(t\right)+\sum_{i=1}^{N}\mu_{i}\underline{h}_{i}^{1}\left(t\right)-\underline{e}^{I}\right)+\underline{S}^{I}\,.\label{eq:37}
+\end{equation}
 
 \end_inset
 
@@ -2305,8 +2544,10 @@
  matrix.
  In addition to the volumetric contribution to the tangent constitutive
  matrix, we require the deviatoric part:
-\begin_inset Formula \begin{equation}
-\frac{\partial\underline{S}}{\partial\underline{\epsilon}}=\frac{\partial\underline{S}}{\partial\underline{e}}\frac{\partial\underline{e}}{\partial\underline{\epsilon}}\,,\label{eq:38}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\underline{S}}{\partial\underline{\epsilon}}=\frac{\partial\underline{S}}{\partial\underline{e}}\frac{\partial\underline{e}}{\partial\underline{\epsilon}}\,,\label{eq:38}
+\end{equation}
 
 \end_inset
 
@@ -2320,8 +2561,10 @@
 
 .
  The other derivative is given by
-\begin_inset Formula \begin{equation}
-\frac{\partial\underline{S}}{\partial\underline{e}}=2\mu_{tot}\left[\mu_{0}\underline{I}+\mu_{1}\frac{\partial\underline{h}^{1}}{\partial\underline{e}}\right]\,,\label{eq:39}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\underline{S}}{\partial\underline{e}}=2\mu_{tot}\left[\mu_{0}\underline{I}+\sum_{i=1}^{N}\mu_{i}\frac{\partial\underline{h}_{i}^{1}}{\partial\underline{e}}\right]\,,\label{eq:39}
+\end{equation}
 
 \end_inset
 
@@ -2345,14 +2588,18 @@
 \end_inset
 
 , the derivative inside the brackets is
-\begin_inset Formula \begin{equation}
-\frac{\partial\underline{h}^{1}}{\partial\underline{e}}=\Delta h\left(\Delta t\right)\underline{I}\,.\label{eq:40}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\underline{h}_{i}^{1}}{\partial\underline{e}}=\Delta h_{i}\left(\Delta t\right)\underline{I}\,.\label{eq:40}
+\end{equation}
 
 \end_inset
 
 The complete deviatoric tangent relation is then
-\begin_inset Formula \begin{equation}
-\frac{\partial\underline{S}}{\partial\underline{\epsilon}}=2\mu_{tot}\left[\mu_{0}+\mu_{1}\Delta h\left(\Delta t\right)\right]\frac{\partial\underline{e}}{\partial\underline{\epsilon}}\,.\label{eq:41}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\underline{S}}{\partial\underline{\epsilon}}=2\mu_{tot}\left[\mu_{0}+\sum_{i=1}^{N}\mu_{i}\Delta h_{i}\left(\Delta t\right)\right]\frac{\partial\underline{e}}{\partial\underline{\epsilon}}\,.\label{eq:41}
+\end{equation}
 
 \end_inset
 
@@ -2361,11 +2608,30 @@
 
 \begin_layout Standard
 We use this formulation for both our Maxwell and generalized Maxwell viscoelasti
-c models, where we include two additional Maxwell elements in the generalized
- Maxwell model.
+c models.
+ For the Maxwell model, 
+\begin_inset Formula $\mu_{0}=0$
+\end_inset
+
+ and 
+\begin_inset Formula $N=1$
+\end_inset
+
+.
+ For the generalized Maxwell model, 
+\begin_inset Formula $N=3.$
+\end_inset
+
+
 \end_layout
 
 \begin_layout Subsection
+\begin_inset CommandInset label
+LatexCommand label
+name "sub:Effective-Stress-Formulations-Viscoelastic"
+
+\end_inset
+
 Effective Stress Formulations for Viscoelastic Materials
 \end_layout
 
@@ -2381,7 +2647,14 @@
  may be employed for both a linear Maxwell model and a power-law Maxwell
  model.
  Note that this formulation is not presently employed for linear viscoelastic
- models, but it may be employed in future versions.
+ models (
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "cha:Alternative-Formulations"
+
+\end_inset
+
+), but it may be employed in future versions.
  For the viscoelastic materials considered here, the viscous volumetric
  strains are zero (incompressible flow), and it is convenient to separate
  the general stress-strain relationship at time 
@@ -2389,9 +2662,11 @@
 \end_inset
 
  into deviatoric and volumetric parts:
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \phantom{}{}^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left(^{t+\Delta t}\underline{e}-\phantom{}^{t+\Delta t}\underline{e}^{C}-\underline{e}^{I}\right)+\underline{S}^{I}=\frac{1}{a_{E}}\left(^{t+\Delta t}\underline{e}-\phantom{}^{t+\Delta t}\underline{e}^{C}-\underline{e}^{I}\right)\label{eq:42}\\
-^{t+\Delta t}P=\frac{E}{1-2\nu}\left(^{t+\Delta t}\theta-\theta^{I}\right)+P^{I}=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta-\theta^{I}\right)\:,\nonumber \end{gather}
+^{t+\Delta t}P=\frac{E}{1-2\nu}\left(^{t+\Delta t}\theta-\theta^{I}\right)+P^{I}=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta-\theta^{I}\right)\:,\nonumber 
+\end{gather}
 
 \end_inset
 
@@ -2441,20 +2716,26 @@
 \end_inset
 
  may also be written as
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{1}{a_{E}}(^{t+\Delta t}\underline{e}^{\prime}-\underline{\Delta e}^{C})+\underline{S}^{I}\,,\label{eq:43}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\underline{S}=\frac{1}{a_{E}}(^{t+\Delta t}\underline{e}^{\prime}-\underline{\Delta e}^{C})+\underline{S}^{I}\,,\label{eq:43}
+\end{gather}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{e}^{\prime}=\phantom{}^{t+\Delta t}\underline{e}-\phantom{}^{t}\underline{e}^{C}-\underline{e}^{I}\,\,,\,\,\,\underline{\Delta e}^{C}=\phantom{}^{t+\Delta t}\underline{e}^{C}-\phantom{}^{t}\underline{e}^{C}\,.\label{eq:44}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\underline{e}^{\prime}=\phantom{}^{t+\Delta t}\underline{e}-\phantom{}^{t}\underline{e}^{C}-\underline{e}^{I}\,\,,\,\,\,\underline{\Delta e}^{C}=\phantom{}^{t+\Delta t}\underline{e}^{C}-\phantom{}^{t}\underline{e}^{C}\,.\label{eq:44}
+\end{gather}
 
 \end_inset
 
 The creep strain increment is approximated using
-\begin_inset Formula \begin{gather}
-\underline{\Delta e}^{C}=\Delta t\phantom{}^{\tau}\gamma\phantom{}^{\tau}\underline{S}\,,\label{eq:45}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\underline{\Delta e}^{C}=\Delta t\phantom{}^{\tau}\gamma\phantom{}^{\tau}\underline{S}\,,\label{eq:45}
+\end{gather}
 
 \end_inset
 
@@ -2463,145 +2744,75 @@
 \end_inset
 
 -method of time integration,
-\begin_inset Formula \begin{gather}
-^{\tau}\underline{S}=(1-\alpha)_{I}^{t}\underline{S}+\alpha\phantom{}_{I}^{t+\Delta t}\underline{S}+\underline{S}^{I}=(1-\alpha)^{t}\underline{S}+\alpha\phantom{}^{t+\Delta t}\underline{S}\,\,,\label{eq:46}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{\tau}\underline{S}=(1-\alpha)_{I}^{t}\underline{S}+\alpha\phantom{}_{I}^{t+\Delta t}\underline{S}+\underline{S}^{I}=(1-\alpha)^{t}\underline{S}+\alpha\phantom{}^{t+\Delta t}\underline{S}\,\,,\label{eq:46}
+\end{gather}
 
 \end_inset
 
 and
-\begin_inset Formula \begin{gather}
-^{\tau}\gamma=\frac{3\Delta\overline{e}^{C}}{2\Delta t\phantom{}^{\tau}\overline{\sigma}}\,\,,\label{eq:47}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{\tau}\gamma=\frac{3\Delta\overline{e}^{C}}{2\Delta t\phantom{}^{\tau}\overline{\sigma}}\,\,,\label{eq:47}
+\end{gather}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{gather}
-\Delta\overline{e}^{C}=\sqrt{\frac{2}{3}\underline{\Delta e}^{C}\cdot\underline{\Delta e}^{C}}\label{eq:48}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\Delta\overline{e}^{C}=\sqrt{\frac{2}{3}\underline{\Delta e}^{C}\cdot\underline{\Delta e}^{C}}\label{eq:48}
+\end{gather}
 
 \end_inset
 
 and
-\begin_inset Formula \begin{gather}
-^{\tau}\overline{\sigma}=(1-\alpha)_{I}^{t}\overline{\sigma}+\alpha\phantom{}_{I}^{t+\Delta t}\overline{\sigma}+\overline{\sigma}^{I}=\sqrt{3\phantom{}^{\tau}J_{2}^{\prime}}\,\,.\label{eq:49}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{\tau}\overline{\sigma}=(1-\alpha)_{I}^{t}\overline{\sigma}+\alpha\phantom{}_{I}^{t+\Delta t}\overline{\sigma}+\overline{\sigma}^{I}=\sqrt{3\phantom{}^{\tau}J_{2}^{\prime}}\,\,.\label{eq:49}
+\end{gather}
 
 \end_inset
 
 
 \end_layout
 
-\begin_layout Subsubsection
-Linear Maxwell Viscoelastic Material
-\end_layout
-
 \begin_layout Standard
-A linear Maxwell viscoelastic material may be characterized by the same
- elastic parameters as an isotropic elastic material (
-\begin_inset Formula $E$
-\end_inset
-
- and 
-\begin_inset Formula $\nu$
-\end_inset
-
-), as well as the viscosity, 
-\begin_inset Formula $\eta$
-\end_inset
-
-.
- The creep strain increment is
-\begin_inset Formula \begin{gather}
-\underline{\Delta e}^{C}=\frac{\Delta t\phantom{}^{\tau}\underline{S}}{2\eta}\,\,.\label{eq:50}\end{gather}
-
-\end_inset
-
-Therefore,
-\begin_inset Formula \begin{gather}
-\Delta\overline{e}^{C}=\frac{\Delta t\sqrt{^{\tau}J_{2}^{\prime}}}{\sqrt{3\eta}}=\frac{\Delta t\phantom{}^{\tau}\overline{\sigma}}{3\eta}\,,\,\mathrm{and}\,^{\tau}\gamma=\frac{1}{2\eta}\,\,.\label{eq:51}\end{gather}
-
-\end_inset
-
-Substituting Equations 
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:46"
-
-\end_inset
-
-, 
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:50"
-
-\end_inset
-
-, and 
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:51"
-
-\end_inset
-
- into 
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:43"
-
-\end_inset
-
-, we obtain
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{1}{a_{E}}\left\{ ^{t+\Delta t}\underline{e}^{\prime}-\frac{\Delta t}{2\eta}\left[(1-\alpha)^{t}\underline{S}+\alpha\phantom{}^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,\,.\label{eq:52}\end{gather}
-
-\end_inset
-
-Solving for 
-\begin_inset Formula $^{t+\Delta t}\underline{S}$
-\end_inset
-
-,
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{1}{a_{E}+\frac{\alpha\Delta t}{2\eta}}\left[^{t+\Delta t}\underline{e}^{\prime}-\frac{\Delta t}{2\eta}(1-\alpha)^{t}\underline{S}+\frac{1+\mathrm{\nu}}{E}\underline{S}^{I}\right]\,\,.\label{eq:53}\end{gather}
-
-\end_inset
-
-In this case it is possible to solve directly for the deviatoric stresses,
- and the effective stress function approach is not needed.
- To obtain the total stress, we simply use
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\sigma_{ij}=\phantom{}^{t+\Delta t}S_{ij}+\frac{\mathit{1}}{a_{m}}\left(\,^{t+\Delta t}\theta-\theta^{I}\right)\delta_{ij}+P^{I}\delta_{ij}\,\,.\label{eq:54}\end{gather}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-It is necessary to provide a relationship for the viscoelastic tangent material
- matrix relating stress and strain.
+To form the global stiffness matrix, it is necessary to provide a relationship
+ for the viscoelastic tangent material matrix relating stress and strain.
  If we use vectors composed of the stresses and tensor strains, this relationshi
 p is
-\begin_inset Formula \begin{gather}
-\underline{C}^{VE}=\frac{\partial\phantom{}^{t+\Delta t}\overrightarrow{\sigma}}{\partial\phantom{}^{t+\Delta t}\overrightarrow{\epsilon}}\,\,.\label{eq:55}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\underline{C}^{VE}=\frac{\partial\phantom{}^{t+\Delta t}\overrightarrow{\sigma}}{\partial\phantom{}^{t+\Delta t}\overrightarrow{\epsilon}}\,\,.\label{eq:55}
+\end{gather}
 
 \end_inset
 
 In terms of the vectors, we have
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 ^{t+\Delta t}\sigma_{i}=\phantom{}^{t+\Delta t}S_{i}+\phantom{}^{t+\Delta t}P\,\,;\,\,\, i=1,2,3\label{eq:56}\\
-^{t+\Delta t}\sigma_{i}=\phantom{}^{t+\Delta t}S_{i}\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, i=4,5,6\nonumber \end{gather}
+^{t+\Delta t}\sigma_{i}=\phantom{}^{t+\Delta t}S_{i}\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, i=4,5,6\nonumber 
+\end{gather}
 
 \end_inset
 
 Therefore,
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 C_{ij}^{VE}=C_{ij}^{\prime}+\frac{1}{3a_{m}}\,;\,\,1\leq i,j\leq3\,\,.\label{eq:57}\\
-C_{ij}^{VE}=C_{ij}^{\prime}\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{otherwise}\nonumber \end{gather}
+C_{ij}^{VE}=C_{ij}^{\prime}\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{otherwise}\nonumber 
+\end{gather}
 
 \end_inset
 
 Using the chain rule,
-\begin_inset Formula \begin{gather}
-\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\frac{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}{\partial\phantom{}^{t+\Delta t}e_{l}}\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\,\,.\label{eq:58}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+C_{ij}^{\prime}=\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}\frac{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}{\partial\phantom{}^{t+\Delta t}e_{l}}\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\,\,.\label{eq:58}
+\end{gather}
 
 \end_inset
 
@@ -2613,8 +2824,10 @@
 \end_inset
 
 , we obtain
-\begin_inset Formula \begin{gather}
-\frac{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}{\partial\phantom{}^{t+\Delta t}e_{l}}=\delta_{kl}\,\,,\label{eq:59}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\frac{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}{\partial\phantom{}^{t+\Delta t}e_{l}}=\delta_{kl}\,\,,\label{eq:59}
+\end{gather}
 
 \end_inset
 
@@ -2626,85 +2839,34 @@
 \end_inset
 
 :
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{3}\left[\begin{array}{ccc}
 2 & -1 & -1\\
 -1 & 2 & -1\\
--1 & -1 & 2\end{array}\right];\,\,1\leq l,j\leq3\label{eq:60}\\
-\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\delta_{lj}\,\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{otherwise.}\nonumber \end{gather}
+-1 & -1 & 2
+\end{array}\right];\,\,1\leq l,j\leq3\label{eq:60}\\
+\frac{\partial\phantom{}^{t+\Delta t}e_{l}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\delta_{lj}\,\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{otherwise.}\nonumber 
+\end{gather}
 
 \end_inset
 
-Finally, from Equation 
+The first term of 
 \begin_inset CommandInset ref
 LatexCommand ref
-reference "eq:53"
+reference "eq:58"
 
 \end_inset
 
-, we have
-\begin_inset Formula \begin{gather}
-\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=\frac{\delta_{ik}}{a_{E}+\frac{\alpha\Delta t}{2\eta}}\,\,.\label{eq:61}\end{gather}
-
-\end_inset
-
-From Equation 
+ depends on the particular constitutive relationship, and the complete tangent
+ matrix may then be obtained from 
 \begin_inset CommandInset ref
 LatexCommand ref
 reference "eq:57"
 
 \end_inset
 
-, the final material matrix relating stress and tensor strain is
-\begin_inset Formula \begin{gather}
-C_{ij}^{VE}=\frac{1}{3a_{m}}\left[\begin{array}{cccccc}
-1 & 1 & 1 & 0 & 0 & 0\\
-1 & 1 & 1 & 0 & 0 & 0\\
-1 & 1 & 1 & 0 & 0 & 0\\
-0 & 0 & 0 & 0 & 0 & 0\\
-0 & 0 & 0 & 0 & 0 & 0\\
-0 & 0 & 0 & 0 & 0 & 0\end{array}\right]+\frac{1}{3\left(a_{E}+\frac{\alpha\Delta t}{2\eta}\right)}\left[\begin{array}{cccccc}
-2 & -1 & -1 & 0 & 0 & 0\\
--1 & 2 & -1 & 0 & 0 & 0\\
--1 & -1 & 2 & 0 & 0 & 0\\
-0 & 0 & 0 & 3 & 0 & 0\\
-0 & 0 & 0 & 0 & 3 & 0\\
-0 & 0 & 0 & 0 & 0 & 3\end{array}\right]\,.\label{eq:62}\end{gather}
-
-\end_inset
-
-Note that the coefficient of the second matrix approaches 
-\begin_inset Formula $E/3(1+\nu)=1/3a_{E}$
-\end_inset
-
- as 
-\begin_inset Formula $\eta$
-\end_inset
-
- goes to infinity.
- To check the results we make sure that the regular elastic constitutive
- matrix is obtained for selected terms in the case where 
-\begin_inset Formula $\eta$
-\end_inset
-
- goes to infinity.
-\begin_inset Formula \begin{gather}
-C_{11}^{E}=\frac{E(1-\nu)}{(1+\nu)(1-2\nu)}\,\,\nonumber \\
-C_{12}^{E}=\frac{E\nu}{(1+\nu)(1-2\nu)}\,.\label{eq:63}\\
-C_{44}^{E}=\frac{E}{1+\nu}\,\,\nonumber \end{gather}
-
-\end_inset
-
-This is consistent with the regular elasticity matrix, and Equation 
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:62"
-
-\end_inset
-
- should thus be used when forming the stiffness matrix.
- We do not presently use this formulation in PyLith 1.4, but it may be included
- in future versions.
+.
 \end_layout
 
 \begin_layout Subsubsection
@@ -2729,8 +2891,10 @@
 \end_inset
 
 ):
-\begin_inset Formula \begin{equation}
-\dot{\epsilon}_{11}^{C}=A_{E}\exp\left(\frac{-Q}{RT}\right)\left(\sigma_{1}-\sigma_{3}\right)^{n}=A_{E}\exp\left(\frac{-Q}{RT}\right)\sigma_{d}^{n}\:,\label{eq:64}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\dot{\epsilon}_{11}^{C}=A_{E}\exp\left(\frac{-Q}{RT}\right)\left(\sigma_{1}-\sigma_{3}\right)^{n}=A_{E}\exp\left(\frac{-Q}{RT}\right)\sigma_{d}^{n}\:,\label{eq:64}
+\end{equation}
 
 \end_inset
 
@@ -2783,8 +2947,10 @@
 ).
  The flow law may then be generalized in terms of the deviatoric stress
  and strain rate invariants:
-\begin_inset Formula \begin{equation}
-\sqrt{\dot{L}_{2}^{\prime C}}=A_{M}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n}\:,\label{eq:65}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sqrt{\dot{L}_{2}^{\prime C}}=A_{M}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n}\:,\label{eq:65}
+\end{equation}
 
 \end_inset
 
@@ -2796,8 +2962,10 @@
  In practice, it is necessary to compute each strain rate component using
  the flow law.
  This is accomplished using:
-\begin_inset Formula \begin{equation}
-\dot{e}_{ij}^{C}=A_{M}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n-1}S_{ij}\:.\label{eq:66}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\dot{e}_{ij}^{C}=A_{M}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n-1}S_{ij}\:.\label{eq:66}
+\end{equation}
 
 \end_inset
 
@@ -2839,10 +3007,12 @@
 \end_inset
 
 , we have
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \sigma_{2}=\sigma_{3}=P_{c}\nonumber \\
 \sigma_{1}=\sigma_{1}^{app}\label{eq:67}\\
-P=\frac{\sigma_{1}+2P_{c}}{3}\:,\nonumber \end{gather}
+P=\frac{\sigma_{1}+2P_{c}}{3}\:,\nonumber 
+\end{gather}
 
 \end_inset
 
@@ -2852,22 +3022,28 @@
 
  is the applied load.
  The deviatoric stresses are then:
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 S_{1}=\frac{2}{3}\left(\sigma_{1}-P_{c}\right)\nonumber \\
-S_{2}=S_{3}=-\frac{1}{3}\left(\sigma_{1}-P_{c}\right)\:.\label{eq:68}\end{gather}
+S_{2}=S_{3}=-\frac{1}{3}\left(\sigma_{1}-P_{c}\right)\:.\label{eq:68}
+\end{gather}
 
 \end_inset
 
 This gives
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 S_{1}=\frac{2}{3}\left(\sigma_{1}-\sigma_{3}\right)=\frac{2}{3}\sigma_{d}\nonumber \\
-S_{2}=S_{3}=-\frac{1}{3}\left(\sigma_{1}-\sigma_{3}\right)=-\frac{1}{3}\sigma_{d}\:.\label{eq:69}\end{gather}
+S_{2}=S_{3}=-\frac{1}{3}\left(\sigma_{1}-\sigma_{3}\right)=-\frac{1}{3}\sigma_{d}\:.\label{eq:69}
+\end{gather}
 
 \end_inset
 
 In terms of the second deviatoric stress invariant, we then have
-\begin_inset Formula \begin{equation}
-\sqrt{J_{2}^{\prime}}=\frac{\sigma_{d}}{\sqrt{3}}\:.\label{eq:70}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sqrt{J_{2}^{\prime}}=\frac{\sigma_{d}}{\sqrt{3}}\:.\label{eq:70}
+\end{equation}
 
 \end_inset
 
@@ -2877,15 +3053,19 @@
 \begin_layout Standard
 Under the assumption that the creep measured in the laboratory experiments
  is incompressible, we have
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \dot{e}_{11}^{C}=\dot{\epsilon}_{11}\nonumber \\
-\dot{e}_{22}^{C}=\dot{e}_{33}^{C}=-\frac{1}{2}\dot{\epsilon}_{11}\:.\label{eq:71}\end{gather}
+\dot{e}_{22}^{C}=\dot{e}_{33}^{C}=-\frac{1}{2}\dot{\epsilon}_{11}\:.\label{eq:71}
+\end{gather}
 
 \end_inset
 
 In terms of the second deviatoric strain rate invariant we then have
-\begin_inset Formula \begin{equation}
-\sqrt{\dot{L}_{2}^{\prime C}}=\frac{\sqrt{3}}{2}\dot{\epsilon}_{11}\:.\label{eq:72}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sqrt{\dot{L}_{2}^{\prime C}}=\frac{\sqrt{3}}{2}\dot{\epsilon}_{11}\:.\label{eq:72}
+\end{equation}
 
 \end_inset
 
@@ -2911,14 +3091,18 @@
 \end_inset
 
 , we obtain
-\begin_inset Formula \begin{equation}
-\sqrt{\dot{L}_{2}^{\prime C}}=A_{E}\frac{\sqrt{3}^{n+1}}{2}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n}\:,\label{eq:73}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sqrt{\dot{L}_{2}^{\prime C}}=A_{E}\frac{\sqrt{3}^{n+1}}{2}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n}\:,\label{eq:73}
+\end{equation}
 
 \end_inset
 
 and therefore,
-\begin_inset Formula \begin{equation}
-A_{M}=\frac{\sqrt{3}^{n+1}}{2}A_{E}\:.\label{eq:74}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{M}=\frac{\sqrt{3}^{n+1}}{2}A_{E}\:.\label{eq:74}
+\end{equation}
 
 \end_inset
 
@@ -2926,8 +3110,10 @@
 \end_layout
 
 \begin_layout Standard
-\begin_inset Formula \begin{equation}
-A_{T}=A_{M}\exp\left(\frac{-Q}{RT}\right)=\frac{\sqrt{3}^{n+1}}{2}A_{E}\exp\left(\frac{-Q}{RT}\right)\:.\label{eq:75}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+A_{T}=A_{M}\exp\left(\frac{-Q}{RT}\right)=\frac{\sqrt{3}^{n+1}}{2}A_{E}\exp\left(\frac{-Q}{RT}\right)\:.\label{eq:75}
+\end{equation}
 
 \end_inset
 
@@ -2961,8 +3147,10 @@
 \end_inset
 
 ): 
-\begin_inset Formula \begin{equation}
-\frac{\sqrt{\dot{L}_{2}^{\prime C}}}{\dot{e}_{0}}=\left(\frac{\sqrt{J_{2}^{\prime}}}{S_{0}}\right)^{n},\label{eq:76}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\sqrt{\dot{L}_{2}^{\prime C}}}{\dot{e}_{0}}=\left(\frac{\sqrt{J_{2}^{\prime}}}{S_{0}}\right)^{n},\label{eq:76}
+\end{equation}
 
 \end_inset
 
@@ -2976,8 +3164,10 @@
 
  are reference values for the strain rate and deviatoric stress.
  This means that
-\begin_inset Formula \begin{equation}
-\frac{\dot{e}_{0}}{S_{0}^{n}}=A_{T}\:.\label{eq:77}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\dot{e}_{0}}{S_{0}^{n}}=A_{T}\:.\label{eq:77}
+\end{equation}
 
 \end_inset
 
@@ -3042,8 +3232,10 @@
 \end_inset
 
  as
-\begin_inset Formula \begin{equation}
-S_{0}=\left(\frac{\dot{e}_{0}}{A_{T}}\right)^{\frac{1}{n}}\:.\label{eq:78}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+S_{0}=\left(\frac{\dot{e}_{0}}{A_{T}}\right)^{\frac{1}{n}}\:.\label{eq:78}
+\end{equation}
 
 \end_inset
 
@@ -3092,20 +3284,26 @@
 
 \begin_layout Standard
 The flow law in component form is 
-\begin_inset Formula \begin{equation}
-\dot{e}_{ij}^{C}=\frac{\dot{e}_{0}\sqrt{J_{2}^{\prime}}^{n-1}S_{ij}}{S_{0}^{n}}\:,\label{eq:79}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\dot{e}_{ij}^{C}=\frac{\dot{e}_{0}\sqrt{J_{2}^{\prime}}^{n-1}S_{ij}}{S_{0}^{n}}\:,\label{eq:79}
+\end{equation}
 
 \end_inset
 
 and the creep strain increment is approximated as
-\begin_inset Formula \begin{gather}
-\underline{\Delta e}^{C}\approx\frac{\Delta t\dot{e}_{0}\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}\,^{\tau}\underline{S}}{S_{0}^{n}}=\frac{\Delta t\dot{e}_{0}\phantom{}^{\tau}\overline{\sigma}^{n-1}\,^{\tau}\underline{S}}{\sqrt{3}S_{0}^{n}}\,.\label{eq:80}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\underline{\Delta e}^{C}\approx\frac{\Delta t\dot{e}_{0}\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}\,^{\tau}\underline{S}}{S_{0}^{n}}=\frac{\Delta t\dot{e}_{0}\phantom{}^{\tau}\overline{\sigma}^{n-1}\,^{\tau}\underline{S}}{\sqrt{3}S_{0}^{n}}\,.\label{eq:80}
+\end{gather}
 
 \end_inset
 
  Therefore,
-\begin_inset Formula \begin{gather}
-\Delta\bar{e}^{C}\approx\frac{2\Delta t\dot{e}_{0}\sqrt{^{\tau}J_{2}^{\prime}}^{n}}{\sqrt{3}S_{0}^{n}}=\frac{2\Delta t\dot{e}_{0}\phantom{}^{\tau}\overline{\sigma}^{n}}{\sqrt{3}^{n+1}S_{0}^{n}}\,,\,\textrm{and}\,^{\tau}\gamma=\frac{\dot{e}_{0}\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}}{S_{0}^{n}}\,.\label{eq:81}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\Delta\bar{e}^{C}\approx\frac{2\Delta t\dot{e}_{0}\sqrt{^{\tau}J_{2}^{\prime}}^{n}}{\sqrt{3}S_{0}^{n}}=\frac{2\Delta t\dot{e}_{0}\phantom{}^{\tau}\overline{\sigma}^{n}}{\sqrt{3}^{n+1}S_{0}^{n}}\,,\,\textrm{and}\,^{\tau}\gamma=\frac{\dot{e}_{0}\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}}{S_{0}^{n}}\,.\label{eq:81}
+\end{gather}
 
 \end_inset
 
@@ -3138,29 +3336,37 @@
 \end_inset
 
 , we obtain:
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{1}{a_{E}}\left\{ ^{t+\Delta t}\underline{e}^{\prime}-\Delta t\phantom{}^{\tau}\gamma\left[\left(1-\alpha\right)^{t}\underline{S}+\alpha{}^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,,\label{eq:82}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\underline{S}=\frac{1}{a_{E}}\left\{ ^{t+\Delta t}\underline{e}^{\prime}-\Delta t\phantom{}^{\tau}\gamma\left[\left(1-\alpha\right)^{t}\underline{S}+\alpha{}^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,,\label{eq:82}
+\end{gather}
 
 \end_inset
 
 which may be rewritten:
-\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}\left(a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma\right)={}^{t+\Delta t}\underline{e}^{\prime}-\Delta t\phantom{}^{\tau}\gamma\left(1-\alpha\right)^{t}\underline{S}+a_{E}\underline{S}^{I}\,.\label{eq:83}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{t+\Delta t}\underline{S}\left(a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma\right)={}^{t+\Delta t}\underline{e}^{\prime}-\Delta t\phantom{}^{\tau}\gamma\left(1-\alpha\right)^{t}\underline{S}+a_{E}\underline{S}^{I}\,.\label{eq:83}
+\end{gather}
 
 \end_inset
 
 Taking the scalar inner product of both sides we obtain:
-\begin_inset Formula \begin{gather}
-a^{2}\,\,{}^{t+\Delta t}J_{2}^{\prime}-b+c\phantom{}^{\tau}\gamma-d^{2}\,^{\tau}\gamma^{2}=F=0\,,\label{eq:84}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+a^{2}\,\,{}^{t+\Delta t}J_{2}^{\prime}-b+c\phantom{}^{\tau}\gamma-d^{2}\,^{\tau}\gamma^{2}=F=0\,,\label{eq:84}
+\end{gather}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 a=a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma\,\,\nonumber \\
 b=\frac{1}{2}{}^{t+\Delta t}\underline{e}^{\prime}\cdot{}^{t+\Delta t}\underline{e}^{\prime}+a_{E}{}^{t+\Delta t}\underline{e}^{\prime}\cdot\underline{S}^{I}+a_{E}^{2}\,^{I}J_{2}^{\prime}\,.\label{eq:85}\\
 c=\Delta t\left(1-\alpha\right){}^{t+\Delta t}\underline{e}^{\prime}\cdot^{t}\underline{S}+\Delta t\left(1-\alpha\right)a_{E}\,^{t}\underline{S}\cdot\underline{S}^{I}\,\,\nonumber \\
-d=\Delta t\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\,\,\nonumber \end{gather}
+d=\Delta t\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\,\,\nonumber 
+\end{gather}
 
 \end_inset
 
@@ -3199,10 +3405,11 @@
 
 \end_inset
 
-, and the total stresses may be found from Equation 
+, and the total stresses may be found by combining the deviatoric and volumetric
+ components from 
 \begin_inset CommandInset ref
 LatexCommand ref
-reference "eq:54"
+reference "eq:42"
 
 \end_inset
 
@@ -3210,15 +3417,15 @@
 \end_layout
 
 \begin_layout Standard
-To compute the tangent stress-strain relation, we proceed as for the linear
- case, but the evaluation of the first derivative in Equation 
+To compute the tangent stress-strain relation, we need to compute the first
+ term in 
 \begin_inset CommandInset ref
 LatexCommand ref
 reference "eq:58"
 
 \end_inset
 
- is more complicated.
+.
  We begin by rewriting Equation 
 \begin_inset CommandInset ref
 LatexCommand ref
@@ -3227,8 +3434,10 @@
 \end_inset
 
  as
-\begin_inset Formula \begin{gather}
-F=^{t+\Delta t}S_{i}\left(a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma\right)-\phantom{}^{t+\Delta t}e_{i}^{\prime}+\Delta t\phantom{}^{\tau}\gamma\left(1-\alpha\right)^{t}S_{i}-a_{E}S_{i}^{I}=0\:.\label{eq:86}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+F=^{t+\Delta t}S_{i}\left(a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma\right)-\phantom{}^{t+\Delta t}e_{i}^{\prime}+\Delta t\phantom{}^{\tau}\gamma\left(1-\alpha\right)^{t}S_{i}-a_{E}S_{i}^{I}=0\:.\label{eq:86}
+\end{gather}
 
 \end_inset
 
@@ -3237,8 +3446,10 @@
 \end_inset
 
  is
-\begin_inset Formula \begin{gather}
-\frac{\partial F}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=-\delta_{ik}\:,\label{eq:87}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\frac{\partial F}{\partial\phantom{}^{t+\Delta t}e_{k}^{\prime}}=-\delta_{ik}\:,\label{eq:87}
+\end{gather}
 
 \end_inset
 
@@ -3247,8 +3458,10 @@
 \end_inset
 
  is
-\begin_inset Formula \begin{gather}
-\frac{\partial F}{\partial\phantom{}^{t+\Delta t}S_{i}}=a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma+\frac{\partial\phantom{}^{\tau}\gamma}{\partial\phantom{}^{t+\Delta t}S_{i}}\Delta t\left[\alpha\phantom{}^{t+\Delta t}S_{i}+\left(1-\alpha\right)^{t}S_{i}\right]\:.\label{eq:88}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\frac{\partial F}{\partial\phantom{}^{t+\Delta t}S_{i}}=a_{E}+\alpha\Delta t\phantom{}^{\tau}\gamma+\frac{\partial\phantom{}^{\tau}\gamma}{\partial\phantom{}^{t+\Delta t}S_{i}}\Delta t\left[\alpha\phantom{}^{t+\Delta t}S_{i}+\left(1-\alpha\right)^{t}S_{i}\right]\:.\label{eq:88}
+\end{gather}
 
 \end_inset
 
@@ -3267,22 +3480,28 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{gather}
-^{\tau}\gamma=\frac{\dot{e}_{0}}{S_{0}^{n}}\left[\alpha\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\right]^{n-1}\:.\label{eq:89}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+^{\tau}\gamma=\frac{\dot{e}_{0}}{S_{0}^{n}}\left[\alpha\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\right]^{n-1}\:.\label{eq:89}
+\end{gather}
 
 \end_inset
 
 Then
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 \frac{\partial\phantom{}^{\tau}\gamma}{\partial{}^{t+\Delta t}S_{i}}=\frac{\partial\phantom{}^{\tau}\gamma}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\frac{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}{\partial\phantom{}^{t+\Delta t}S_{l}}\label{eq:90}\\
-=\frac{\dot{e}_{0}\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}{}^{t+\Delta t}T_{i}}{2S_{0}^{n}}\,,\nonumber \end{gather}
+=\frac{\dot{e}_{0}\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}{}^{t+\Delta t}T_{i}}{2S_{0}^{n}}\,,\nonumber 
+\end{gather}
 
 \end_inset
 
 Where
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 ^{t+\Delta t}T_{i}=\phantom{}^{t+\Delta t}S_{i}\:;\:\:1\leq i\leq3\label{eq:91}\\
-^{t+\Delta t}T_{i}=2\phantom{}^{t+\Delta t}S_{i}\:;\:\:\textrm{otherwise.}\nonumber \end{gather}
+^{t+\Delta t}T_{i}=2\phantom{}^{t+\Delta t}S_{i}\:;\:\:\textrm{otherwise.}\nonumber 
+\end{gather}
 
 \end_inset
 
@@ -3308,8 +3527,10 @@
 \end_inset
 
 , and the quotient rule for derivatives of an implicit function,
-\begin_inset Formula \begin{gather}
-\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime}}=\frac{\delta_{ik}}{a_{E}+\alpha\Delta t\left[^{\tau}\gamma+\frac{\dot{e}_{0}{}^{\tau}S_{i}\left(n-1\right){}^{t+\Delta t}T_{i}\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{2\sqrt{^{t+\Delta t}J_{2}^{\prime}}S_{0}^{n}}\right]}\,.\label{eq:92}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\frac{\partial\phantom{}^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime}}=\frac{\delta_{ik}}{a_{E}+\alpha\Delta t\left[^{\tau}\gamma+\frac{\dot{e}_{0}{}^{\tau}S_{i}\left(n-1\right){}^{t+\Delta t}T_{i}\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{2\sqrt{^{t+\Delta t}J_{2}^{\prime}}S_{0}^{n}}\right]}\,.\label{eq:92}
+\end{gather}
 
 \end_inset
 
@@ -3317,10 +3538,10 @@
 \begin_inset Formula $\left(n=1\right)$
 \end_inset
 
-, this equation is identical to Equation 
+, this equation is identical to the linear formulation in 
 \begin_inset CommandInset ref
 LatexCommand ref
-reference "eq:61"
+reference "sub:Effective-Stress-Formulation-Maxwell"
 
 \end_inset
 
@@ -3344,20 +3565,24 @@
 \end_inset
 
 ,
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 C_{ij}^{VE}=\frac{1}{3a_{m}}\left[\begin{array}{cccccc}
 1 & 1 & 1 & 0 & 0 & 0\\
 1 & 1 & 1 & 0 & 0 & 0\\
 1 & 1 & 1 & 0 & 0 & 0\\
 0 & 0 & 0 & 0 & 0 & 0\\
 0 & 0 & 0 & 0 & 0 & 0\\
-0 & 0 & 0 & 0 & 0 & 0\end{array}\right]+\frac{1}{3}\frac{\partial{}^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime}}\left[\begin{array}{cccccc}
+0 & 0 & 0 & 0 & 0 & 0
+\end{array}\right]+\frac{1}{3}\frac{\partial{}^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime}}\left[\begin{array}{cccccc}
 2 & -1 & -1 & 0 & 0 & 0\\
 -1 & 2 & -1 & 0 & 0 & 0\\
 -1 & -1 & 2 & 0 & 0 & 0\\
 0 & 0 & 0 & 3 & 0 & 0\\
 0 & 0 & 0 & 0 & 3 & 0\\
-0 & 0 & 0 & 0 & 0 & 3\end{array}\right]\,.\label{eq:93}\end{gather}
+0 & 0 & 0 & 0 & 0 & 3
+\end{array}\right]\,.\label{eq:93}
+\end{gather}
 
 \end_inset
 
@@ -3393,8 +3618,10 @@
 \end_inset
 
 , which may be written:
-\begin_inset Formula \begin{gather}
-\frac{\partial F}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}=2a^{2}\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\frac{\dot{e}_{0}\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{S_{0}^{n}}\left(2a\alpha\Delta t{}^{t+\Delta t}J_{2}^{\prime}+c-2d^{2}\,^{\tau}\gamma\right)\,.\label{eq:94}\end{gather}
+\begin_inset Formula 
+\begin{gather}
+\frac{\partial F}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}=2a^{2}\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\frac{\dot{e}_{0}\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{S_{0}^{n}}\left(2a\alpha\Delta t{}^{t+\Delta t}J_{2}^{\prime}+c-2d^{2}\,^{\tau}\gamma\right)\,.\label{eq:94}
+\end{gather}
 
 \end_inset
 
@@ -3402,9 +3629,8 @@
 \end_layout
 
 \begin_layout Standard
-The effective stress formulation is the method used in PyLith 1.4 for power-law
- viscoelasticity.
- We may add additional formulations in the future.
+The effective stress formulation is the method used in PyLith 1.4 and above
+ for power-law viscoelasticity.
 \end_layout
 
 \begin_layout Standard
@@ -3437,7 +3663,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="5" columns="2">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -3607,7 +3833,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="10" columns="2">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -3902,7 +4128,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="7" columns="2">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>
@@ -4111,14 +4337,18 @@
 \begin_layout Standard
 The elastoplasticity formulation in PyLith is based on an additive decomposition
  of the total strain into elastic and plastic parts:
-\begin_inset Formula \begin{equation}
-d\epsilon_{ij}=d\epsilon_{ij}^{E}+d\epsilon_{ij}^{P}\:.\label{eq:95}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+d\epsilon_{ij}=d\epsilon_{ij}^{E}+d\epsilon_{ij}^{P}\:.\label{eq:95}
+\end{equation}
 
 \end_inset
 
 The stress increment is then given by
-\begin_inset Formula \begin{equation}
-d\sigma_{ij}=C_{ijrs}^{E}\left(d\epsilon_{rs}-d\epsilon_{rs}^{P}\right)\:,\label{eq:96}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+d\sigma_{ij}=C_{ijrs}^{E}\left(d\epsilon_{rs}-d\epsilon_{rs}^{P}\right)\:,\label{eq:96}
+\end{equation}
 
 \end_inset
 
@@ -4131,8 +4361,10 @@
  We first require a yield condition, which specifies the state of stress
  at which plastic flow initiates.
  This is generally given in the form:
-\begin_inset Formula \begin{equation}
-f\left(\underline{\sigma},k\right)=0\:,\label{eq:97}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+f\left(\underline{\sigma},k\right)=0\:,\label{eq:97}
+\end{equation}
 
 \end_inset
 
@@ -4144,14 +4376,18 @@
  It is then necessary to specify a flow rule, which describes the relationship
  between plastic strain and stress.
  The flow rule is given in the form:
-\begin_inset Formula \begin{equation}
-g\left(\underline{\sigma},k\right)=0\:.\label{eq:98}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+g\left(\underline{\sigma},k\right)=0\:.\label{eq:98}
+\end{equation}
 
 \end_inset
 
 The plastic strain increment is then given as
-\begin_inset Formula \begin{equation}
-d\epsilon_{ij}^{P}=d\lambda\frac{\partial g}{\partial\sigma_{ij}}\:,\label{eq:99}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+d\epsilon_{ij}^{P}=d\lambda\frac{\partial g}{\partial\sigma_{ij}}\:,\label{eq:99}
+\end{equation}
 
 \end_inset
 
@@ -4212,20 +4448,26 @@
  which allows control over the unreasonable amounts of dilatation that are
  sometimes predicted by the associated model.
  The model is described by the following yield condition:
-\begin_inset Formula \begin{equation}
-f\left(\underline{\sigma},k\right)=\alpha_{f}I_{1}+\sqrt{J_{2}^{\prime}}-\beta\:,\label{eq:100}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+f\left(\underline{\sigma},k\right)=\alpha_{f}I_{1}+\sqrt{J_{2}^{\prime}}-\beta\:,\label{eq:100}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{equation}
-\alpha_{f}=\frac{2\sin\phi\left(k\right)}{\sqrt{3}\left(3-\sin\phi\left(k\right)\right)}\:,\label{eq:101}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\alpha_{f}=\frac{2\sin\phi\left(k\right)}{\sqrt{3}\left(3-\sin\phi\left(k\right)\right)}\:,\label{eq:101}
+\end{equation}
 
 \end_inset
 
 and
-\begin_inset Formula \begin{equation}
-\beta=\frac{6\bar{c}\left(k\right)\cos\phi_{0}}{\sqrt{3}\left(3-\sin\phi_{0}\right)}\:.\label{eq:102}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\beta=\frac{6\bar{c}\left(k\right)\cos\phi_{0}}{\sqrt{3}\left(3-\sin\phi_{0}\right)}\:.\label{eq:102}
+\end{equation}
 
 \end_inset
 
@@ -4268,14 +4510,18 @@
  represents a circular cone in principal stress space that is coincident
  with the outer apices of the corresponding Mohr-Coulomb yield surface.
  The flow rule is given by:
-\begin_inset Formula \begin{equation}
-g\left(\underline{\sigma},k\right)=\sqrt{J_{2}^{\prime}}+\alpha_{g}I_{1}\:,\label{eq:103}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+g\left(\underline{\sigma},k\right)=\sqrt{J_{2}^{\prime}}+\alpha_{g}I_{1}\:,\label{eq:103}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{equation}
-\alpha_{g}=\frac{2\sin\psi(k)}{\sqrt{3}\left(3-\sin\psi\left(k\right)\right)}\:.\label{eq:104}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\alpha_{g}=\frac{2\sin\psi(k)}{\sqrt{3}\left(3-\sin\psi\left(k\right)\right)}\:.\label{eq:104}
+\end{equation}
 
 \end_inset
 
@@ -4289,18 +4535,22 @@
 \begin_layout Standard
 As for the viscoelastic models, it is convenient to separate the deformation
  into deviatoric and volumetric parts:
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 ^{t+\Delta t}S_{ij}=\frac{1}{a_{E}}\left(^{t+\Delta t}e_{ij}-\phantom{}^{t+\Delta t}e_{ij}^{P}-e_{ij}^{I}\right)+S_{ij}^{I}=\frac{1}{a_{E}}\left(^{t+\Delta t}e_{ij}^{\prime}-\Delta e_{ij}^{P}\right)+S_{ij}^{I}\label{eq:105}\\
-^{t+\Delta t}P=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta-\phantom{}^{t+\Delta t}\theta^{P}-\theta^{I}\right)+P^{I}=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\Delta\theta^{P}\right)+P^{I}\:,\nonumber \end{gather}
+^{t+\Delta t}P=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta-\phantom{}^{t+\Delta t}\theta^{P}-\theta^{I}\right)+P^{I}=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\Delta\theta^{P}\right)+P^{I}\:,\nonumber 
+\end{gather}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 ^{t+\Delta t}e_{ij}^{\prime}=\phantom{}^{t+\Delta t}e_{ij}-\phantom{}^{t}e_{ij}^{P}-e_{ij}^{I}\nonumber \\
 \Delta e_{ij}^{P}=\phantom{}^{t+\Delta t}e_{ij}^{P}-\phantom{}^{t}e_{ij}^{P}\nonumber \\
 ^{t+\Delta t}\theta^{\prime}=\phantom{}^{t+\Delta t}\theta-\phantom{}^{t}\theta^{P}-\theta^{I}\nonumber \\
-\Delta\theta^{P}=\phantom{}^{t+\Delta t}\theta^{P}-\phantom{}^{t}\theta^{P}\:.\label{eq:106}\end{gather}
+\Delta\theta^{P}=\phantom{}^{t+\Delta t}\theta^{P}-\phantom{}^{t}\theta^{P}\:.\label{eq:106}
+\end{gather}
 
 \end_inset
 
@@ -4314,20 +4564,26 @@
 \end_inset
 
 , the plastic strain increment is
-\begin_inset Formula \begin{equation}
-\Delta\epsilon_{ij}^{P}=\lambda\frac{\partial\phantom{}^{t+\Delta t}g}{\partial\phantom{}^{t+\Delta t}\sigma_{ij}}=\lambda\alpha_{g}\delta_{ij}+\lambda\frac{^{t+\Delta t}S_{ij}}{2\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\:.\label{eq:107}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta\epsilon_{ij}^{P}=\lambda\frac{\partial\phantom{}^{t+\Delta t}g}{\partial\phantom{}^{t+\Delta t}\sigma_{ij}}=\lambda\alpha_{g}\delta_{ij}+\lambda\frac{^{t+\Delta t}S_{ij}}{2\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\:.\label{eq:107}
+\end{equation}
 
 \end_inset
 
 The volumetric part is
-\begin_inset Formula \begin{equation}
-\Delta\theta^{P}=\frac{1}{3}\Delta\epsilon_{ii}^{P}=\lambda\alpha_{g}\:,\label{eq:108}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta\theta^{P}=\frac{1}{3}\Delta\epsilon_{ii}^{P}=\lambda\alpha_{g}\:,\label{eq:108}
+\end{equation}
 
 \end_inset
 
 and the deviatoric part is
-\begin_inset Formula \begin{equation}
-\Delta e_{ij}^{P}=\Delta\epsilon_{ij}^{P}-\Delta\epsilon_{m}^{P}\delta_{ij}=\lambda\frac{^{t+\Delta t}S_{ij}}{2\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\:.\label{eq:109}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta e_{ij}^{P}=\Delta\epsilon_{ij}^{P}-\Delta\epsilon_{m}^{P}\delta_{ij}=\lambda\frac{^{t+\Delta t}S_{ij}}{2\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\:.\label{eq:109}
+\end{equation}
 
 \end_inset
 
@@ -4443,20 +4699,26 @@
 \end_inset
 
  and taking the scalar product of both sides:
-\begin_inset Formula \begin{equation}
-\lambda=\sqrt{2}\,\phantom{}^{t+\Delta t}d-2a_{E}\sqrt{^{t+\Delta t}J_{2}^{\prime}}\:,\label{eq:110}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\lambda=\sqrt{2}\,\phantom{}^{t+\Delta t}d-2a_{E}\sqrt{^{t+\Delta t}J_{2}^{\prime}}\:,\label{eq:110}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{equation}
-^{t+\Delta t}d^{2}=2a_{E}^{2}J_{2}^{\prime I}+2a_{E}S_{ij}^{I}\,\phantom{}^{t+\Delta t}e_{ij}^{\prime}+\phantom{}^{t+\Delta t}e_{ij}^{\prime}\,\phantom{}^{t+\Delta t}e_{ij}^{\prime}\:.\label{eq:111}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+^{t+\Delta t}d^{2}=2a_{E}^{2}J_{2}^{\prime I}+2a_{E}S_{ij}^{I}\,\phantom{}^{t+\Delta t}e_{ij}^{\prime}+\phantom{}^{t+\Delta t}e_{ij}^{\prime}\,\phantom{}^{t+\Delta t}e_{ij}^{\prime}\:.\label{eq:111}
+\end{equation}
 
 \end_inset
 
 The second deviatoric stress invariant is therefore
-\begin_inset Formula \begin{equation}
-\sqrt{^{t+\Delta t}J_{2}^{\prime}}=\frac{\sqrt{2}\,\phantom{}^{t+\Delta t}d-\lambda}{2a_{E}}\:,\label{eq:112}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\sqrt{^{t+\Delta t}J_{2}^{\prime}}=\frac{\sqrt{2}\,\phantom{}^{t+\Delta t}d-\lambda}{2a_{E}}\:,\label{eq:112}
+\end{equation}
 
 \end_inset
 
@@ -4475,8 +4737,10 @@
 \end_inset
 
  as:
-\begin_inset Formula \begin{equation}
-^{t+\Delta t}P=\frac{^{t+\Delta t}I_{1}}{3}=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\lambda\alpha_{g}\right)\:.\label{eq:113}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+^{t+\Delta t}P=\frac{^{t+\Delta t}I_{1}}{3}=\frac{1}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\lambda\alpha_{g}\right)\:.\label{eq:113}
+\end{equation}
 
 \end_inset
 
@@ -4489,8 +4753,10 @@
 \end_inset
 
  to obtain:
-\begin_inset Formula \begin{equation}
-\lambda=\frac{2a_{E}a_{m}\left(\frac{3\alpha_{f}}{a_{m}}\phantom{}^{t+\Delta t}\theta^{\prime}+\frac{^{t+\Delta t}d}{\sqrt{2}a_{E}}-\beta\bar{c}\right)}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\:.\label{eq:114}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\lambda=\frac{2a_{E}a_{m}\left(\frac{3\alpha_{f}}{a_{m}}\phantom{}^{t+\Delta t}\theta^{\prime}+\frac{^{t+\Delta t}d}{\sqrt{2}a_{E}}-\beta\bar{c}\right)}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\:.\label{eq:114}
+\end{equation}
 
 \end_inset
 
@@ -4513,8 +4779,10 @@
 \end_inset
 
  to obtain
-\begin_inset Formula \begin{equation}
-^{t+\Delta t}S_{ij}=\frac{\Delta e_{ij}^{P}\left(\sqrt{2}\,\phantom{\,}^{t+\Delta t}d-\lambda\right)}{\lambda a_{E}}\:.\label{eq:115}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+^{t+\Delta t}S_{ij}=\frac{\Delta e_{ij}^{P}\left(\sqrt{2}\,\phantom{\,}^{t+\Delta t}d-\lambda\right)}{\lambda a_{E}}\:.\label{eq:115}
+\end{equation}
 
 \end_inset
 
@@ -4526,8 +4794,10 @@
 \end_inset
 
 , we obtain the deviatoric plastic strain increment:
-\begin_inset Formula \begin{equation}
-\Delta e_{ij}^{P}=\frac{\lambda}{\sqrt{2}\,\phantom{}^{t+\Delta t}d}\left(^{t+\Delta t}e_{ij}^{\prime}+a_{E}S_{ij}^{I}\right)\:.\label{eq:116}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\Delta e_{ij}^{P}=\frac{\lambda}{\sqrt{2}\,\phantom{}^{t+\Delta t}d}\left(^{t+\Delta t}e_{ij}^{\prime}+a_{E}S_{ij}^{I}\right)\:.\label{eq:116}
+\end{equation}
 
 \end_inset
 
@@ -4571,21 +4841,27 @@
 \end_inset
 
  as a single equation in terms of stress and strain vectors:
-\begin_inset Formula \begin{equation}
-^{t+\Delta t}\sigma_{i}=\frac{1}{a_{E}}\left(^{t+\Delta t}e_{i}^{\prime}-\Delta e_{i}^{P}\right)+S_{i}^{I}+\frac{R_{i}}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\Delta\theta^{P}\right)+P^{I}\label{eq:117}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+^{t+\Delta t}\sigma_{i}=\frac{1}{a_{E}}\left(^{t+\Delta t}e_{i}^{\prime}-\Delta e_{i}^{P}\right)+S_{i}^{I}+\frac{R_{i}}{a_{m}}\left(^{t+\Delta t}\theta^{\prime}-\Delta\theta^{P}\right)+P^{I}\label{eq:117}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{gather}
+\begin_inset Formula 
+\begin{gather}
 R_{i}=1\:;\; i=1,2,3\label{eq:118}\\
-R_{i}=0\:;\; i=4,5,6\:.\nonumber \end{gather}
+R_{i}=0\:;\; i=4,5,6\:.\nonumber 
+\end{gather}
 
 \end_inset
 
 The elastoplastic tangent matrix is then given by
-\begin_inset Formula \begin{equation}
-C_{ij}^{EP}=\frac{\partial\phantom{}^{t+\Delta t}\sigma_{i}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{a_{E}}\left(\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)+\frac{R_{i}}{a_{m}}\left(\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\partial\Delta\theta^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\:.\label{eq:119}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+C_{ij}^{EP}=\frac{\partial\phantom{}^{t+\Delta t}\sigma_{i}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{a_{E}}\left(\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)+\frac{R_{i}}{a_{m}}\left(\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\partial\Delta\theta^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\:.\label{eq:119}
+\end{equation}
 
 \end_inset
 
@@ -4608,14 +4884,17 @@
 \end_inset
 
 , we have
-\begin_inset Formula \begin{equation}
+\begin_inset Formula 
+\begin{equation}
 \frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{3}\left[\begin{array}{cccccc}
 2 & -1 & -1 & 0 & 0 & 0\\
 -1 & 2 & -1 & 0 & 0 & 0\\
 -1 & -1 & 2 & 0 & 0 & 0\\
 0 & 0 & 0 & 3 & 0 & 0\\
 0 & 0 & 0 & 0 & 3 & 0\\
-0 & 0 & 0 & 0 & 0 & 3\end{array}\right]\:,\label{eq:120}\end{equation}
+0 & 0 & 0 & 0 & 0 & 3
+\end{array}\right]\:,\label{eq:120}
+\end{equation}
 
 \end_inset
 
@@ -4634,8 +4913,10 @@
 \end_inset
 
  we have
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{R_{j}}{3}\:.\label{eq:121}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{R_{j}}{3}\:.\label{eq:121}
+\end{equation}
 
 \end_inset
 
@@ -4647,8 +4928,10 @@
 \end_inset
 
  we have
-\begin_inset Formula \begin{equation}
-\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{\sqrt{2}\,\phantom{}^{t+\Delta t}d}\left[\left(^{t+\Delta t}e_{i}^{\prime}+a_{E}S_{i}^{I}\right)\left(\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\lambda}{\phantom{}^{t+\Delta t}d}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)+\lambda\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right]\:.\label{eq:122}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\Delta e_{i}^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{1}{\sqrt{2}\,\phantom{}^{t+\Delta t}d}\left[\left(^{t+\Delta t}e_{i}^{\prime}+a_{E}S_{i}^{I}\right)\left(\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}-\frac{\lambda}{\phantom{}^{t+\Delta t}d}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)+\lambda\frac{\partial\phantom{}^{t+\Delta t}e_{i}^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right]\:.\label{eq:122}
+\end{equation}
 
 \end_inset
 
@@ -4657,15 +4940,19 @@
 \end_inset
 
  is
-\begin_inset Formula \begin{equation}
-\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\phantom{}^{t+\Delta t}d}\:,\label{eq:123}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\phantom{}^{t+\Delta t}d}\:,\label{eq:123}
+\end{equation}
 
 \end_inset
 
 where
-\begin_inset Formula \begin{align}
+\begin_inset Formula 
+\begin{align}
 T_{j}^{I} & =S_{j}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{j}=\phantom{}^{t+\Delta t}e_{j}^{\prime}\:;\; j=1,2,3\nonumber \\
-T_{j}^{I} & =2S_{j}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{j}=2\phantom{}^{t+\Delta t}e_{j}^{\prime}\:;\; j=4,5,6\:.\label{eq:124}\end{align}
+T_{j}^{I} & =2S_{j}^{I}\;\mathrm{and}\;\phantom{}^{t+\Delta t}E_{j}=2\phantom{}^{t+\Delta t}e_{j}^{\prime}\:;\; j=4,5,6\:.\label{eq:124}
+\end{align}
 
 \end_inset
 
@@ -4674,9 +4961,11 @@
 \end_inset
 
  is a function of derivatives already computed:
-\begin_inset Formula \begin{align}
+\begin_inset Formula 
+\begin{align}
 \frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}} & =\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{3\alpha_{f}}{a_{m}}\frac{\partial\phantom{}^{t+\Delta t}\theta^{\prime}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}+\frac{1}{\sqrt{2}a_{E}}\frac{\partial\phantom{}^{t+\Delta t}d}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\right)\nonumber \\
- & =\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{\alpha_{f}R_{j}}{a_{m}}+\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\sqrt{2}a_{E}\phantom{}^{t+\Delta t}d}\right)\:.\label{eq:125}\end{align}
+ & =\frac{2a_{E}a_{m}}{6\alpha_{f}\alpha_{g}a_{E}+a_{m}}\left(\frac{\alpha_{f}R_{j}}{a_{m}}+\frac{a_{E}T_{j}^{I}+\phantom{}^{t+\Delta t}E_{j}}{\sqrt{2}a_{E}\phantom{}^{t+\Delta t}d}\right)\:.\label{eq:125}
+\end{align}
 
 \end_inset
 
@@ -4688,8 +4977,10 @@
 \end_inset
 
 , the derivative of the volumetric plastic strain increment is:
-\begin_inset Formula \begin{equation}
-\frac{\partial\Delta\theta^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\alpha_{g}\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\:.\label{eq:126}\end{equation}
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\Delta\theta^{P}}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}=\alpha_{g}\frac{\partial\lambda}{\partial\phantom{}^{t+\Delta t}\epsilon_{j}}\:.\label{eq:126}
+\end{equation}
 
 \end_inset
 
@@ -4719,7 +5010,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="7" columns="2">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="0">
 <column alignment="center" valignment="top" width="0">
 <row>
@@ -4995,7 +5286,7 @@
 
 \begin_inset Tabular
 <lyxtabular version="3" rows="3" columns="2">
-<features>
+<features tabularvalignment="middle">
 <column alignment="center" valignment="middle" width="0.85in">
 <column alignment="center" valignment="middle" width="2.47in">
 <row>

Modified: short/3D/PyLith/trunk/doc/userguide/userguide.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/userguide.lyx	2011-06-10 01:29:23 UTC (rev 18574)
+++ short/3D/PyLith/trunk/doc/userguide/userguide.lyx	2011-06-10 05:16:17 UTC (rev 18575)
@@ -339,6 +339,13 @@
 
 
 \begin_inset CommandInset include
+LatexCommand include
+filename "alternativeformul/alternativeformul.lyx"
+
+\end_inset
+
+
+\begin_inset CommandInset include
 LatexCommand input
 filename "analyticalsolns/analyticalsolns.lyx"
 



More information about the CIG-COMMITS mailing list