[cig-commits] r20414 - in short/3D/PyLith/branches/v1.7-stable/doc/userguide: governingeqns intro runpylith
sue at geodynamics.org
sue at geodynamics.org
Tue Jun 26 09:08:27 PDT 2012
Author: sue
Date: 2012-06-26 09:08:27 -0700 (Tue, 26 Jun 2012)
New Revision: 20414
Modified:
short/3D/PyLith/branches/v1.7-stable/doc/userguide/governingeqns/governingeqns.lyx
short/3D/PyLith/branches/v1.7-stable/doc/userguide/intro/intro.lyx
short/3D/PyLith/branches/v1.7-stable/doc/userguide/runpylith/runpylith.lyx
Log:
fixed typos and format problems
Modified: short/3D/PyLith/branches/v1.7-stable/doc/userguide/governingeqns/governingeqns.lyx
===================================================================
--- short/3D/PyLith/branches/v1.7-stable/doc/userguide/governingeqns/governingeqns.lyx 2012-06-26 16:08:18 UTC (rev 20413)
+++ short/3D/PyLith/branches/v1.7-stable/doc/userguide/governingeqns/governingeqns.lyx 2012-06-26 16:08:27 UTC (rev 20414)
@@ -1,1859 +1,1859 @@
-#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
-\papersize default
-\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
-\bottommargin 1in
-\secnumdepth 3
-\tocdepth 3
-\paragraph_separation indent
-\paragraph_indentation default
-\quotes_language english
-\papercolumns 1
-\papersides 1
-\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:Governing-Equations"
-
-\end_inset
-
-Governing Equations
-\end_layout
-
-\begin_layout Standard
-We present here a brief derivation of the equations for both quasi-static
- and dynamic computations.
- Since the general equations are the same (except for the absence of inertial
- terms in the quasi-static case), we first derive these equations.
- We then present solution methods for each specific case.
- In all of our derivations, we use the notation described in Table
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "tab:Mathematical-notation"
-
-\end_inset
-
- for both index and vector notation.
- When using index notation, we use the common convention that repeated indices
- indicate summation over the range of the index.
-\end_layout
-
-\begin_layout Standard
-\noindent
-\align center
-\begin_inset Float table
-placement H
-wide false
-sideways false
-status open
-
-\begin_layout Plain Layout
-\begin_inset Caption
-
-\begin_layout Plain Layout
-\begin_inset CommandInset label
-LatexCommand label
-name "tab:Mathematical-notation"
-
-\end_inset
-
-Mathematical notation
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Plain Layout
-\noindent
-\align center
-\begin_inset VSpace defskip
-\end_inset
-
-
-\begin_inset Tabular
-<lyxtabular version="3" rows="11" columns="3">
-<features tabularvalignment="middle">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
-<column alignment="center" valignment="top" width="0">
-<row>
-<cell multicolumn="1" alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Symbol
-\end_layout
-
-\end_inset
-</cell>
-<cell multicolumn="2" alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-
-\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
-Description
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Index notation
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Vector Notation
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" rightline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $a_{i}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\overrightarrow{a}$
-\end_inset
-
-
-\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
-Vector field a
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $a_{ij}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\underline{a}$
-\end_inset
-
-
-\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
-Second order tensor field a
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $u_{i}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\overrightarrow{u}$
-\end_inset
-
-
-\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
-Displacement vector field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $d_{i}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\vec{{d}}$
-\end_inset
-
-
-\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
-Fault slip vector field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $f_{i}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\overrightarrow{f}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Body force vector field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $T_{i}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\overrightarrow{T}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" leftline="true" rightline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Traction vector field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\sigma_{ij}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\underline{\sigma}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Stress tensor field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $n_{i}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\overrightarrow{n}$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" rightline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Normal vector field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-<row>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\rho$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-\begin_inset Formula $\rho$
-\end_inset
-
-
-\end_layout
-
-\end_inset
-</cell>
-<cell alignment="center" valignment="top" bottomline="true" leftline="true" rightline="true" usebox="none">
-\begin_inset Text
-
-\begin_layout Plain Layout
-Mass density scalar field
-\end_layout
-
-\end_inset
-</cell>
-</row>
-</lyxtabular>
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Section
-Derivation of Elasticity Equation
-\end_layout
-
-\begin_layout Subsection
-Index Notation
-\end_layout
-
-\begin_layout Standard
-Consider volume
-\begin_inset Formula $V$
-\end_inset
-
- bounded by surface
-\begin_inset Formula $S$
-\end_inset
-
-.
- Applying a Lagrangian description of the conservation of momentum gives
-\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}
-
-\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}
-
-\end_inset
-
-where
-\begin_inset Formula $n_{j}$
-\end_inset
-
- is the vector normal to
-\begin_inset Formula $S$
-\end_inset
-
-.
- Substituting into equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eqn:momentum:index"
-
-\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}
-
-\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}
-
-\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}
-
-\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}
-
-\end_inset
-
-Because the volume
-\begin_inset Formula $V$
-\end_inset
-
- is arbitrary, the integrand must be zero at every location in the volume,
- so that we end up with
-\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}
-
-\end_inset
-
-We specify tractions,
-\begin_inset Formula $T_{i}$
-\end_inset
-
-, on surface
-\begin_inset Formula $S_{f}$
-\end_inset
-
-, displacements,
-\begin_inset Formula $u_{i}^{o}$
-\end_inset
-
-, on surface
-\begin_inset Formula $S_{u}$
-\end_inset
-
-, and slip,
-\begin_inset Formula $d_{k}$
-\end_inset
-
-, on fault surface
-\begin_inset Formula $S_{f}$
-\end_inset
-
- (we will consider the case of fault constitutive models in Section
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "sec:fault"
-
-\end_inset
-
-).
- The rotation matrix
-\begin_inset Formula $R_{ki}$
-\end_inset
-
- transforms vectors from the global coordinate system to the fault coordinate
- system.
- Note that since both
-\begin_inset Formula $T_{i}$
-\end_inset
-
- and
-\begin_inset Formula $u_{i}$
-\end_inset
-
- are vector quantities, there can be some spatial overlap of the surfaces
-
-\begin_inset Formula $S_{T}$
-\end_inset
-
- and
-\begin_inset Formula $S_{u}$
-\end_inset
-
-; however, the same degree of freedom cannot simultaneously have both types
- of boundary conditions.
-\end_layout
-
-\begin_layout Subsection
-Vector Notation
-\end_layout
-
-\begin_layout Standard
-Consider volume
-\begin_inset Formula $V$
-\end_inset
-
- bounded by surface
-\begin_inset Formula $S$
-\end_inset
-
-.
- Applying a Lagrangian description of the conservation of momentum gives
-\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}
-
-\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}
-
-\end_inset
-
-where
-\begin_inset Formula $\overrightarrow{n}$
-\end_inset
-
- is the vector normal to
-\begin_inset Formula $S$
-\end_inset
-
-.
- Substituting into equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eqn:momentum:vec"
-
-\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}
-
-\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}
-
-\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}
-
-\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}
-
-\end_inset
-
-Because the volume
-\begin_inset Formula $V$
-\end_inset
-
- 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}
-\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}
-
-\end_inset
-
-We specify tractions,
-\begin_inset Formula $\vec{T}$
-\end_inset
-
-, on surface
-\begin_inset Formula $S_{f}$
-\end_inset
-
-, displacements,
-\begin_inset Formula $\overrightarrow{u^{o}}$
-\end_inset
-
-, on surface
-\begin_inset Formula $S_{u}$
-\end_inset
-
-, and slip,
-\begin_inset Formula $\vec{d}$
-\end_inset
-
-, on fault surface
-\begin_inset Formula $S_{f}$
-\end_inset
-
- (we will consider the case of fault constitutive models in Section
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "sec:fault"
-
-\end_inset
-
-).
- The rotation matrix
-\begin_inset Formula $\underline{R}$
-\end_inset
-
- transforms vectors from the global coordinate system to the fault coordinate
- system.
- Note that since both
-\begin_inset Formula $\overrightarrow{T}$
-\end_inset
-
- and
-\begin_inset Formula $\overrightarrow{u}$
-\end_inset
-
- are vector quantities, there can be some spatial overlap of the surfaces
-
-\begin_inset Formula $S_{T}$
-\end_inset
-
- and
-\begin_inset Formula $S_{u}$
-\end_inset
-
-; however, the same degree of freedom cannot simultaneously have both types
- of boundary conditions.
-\end_layout
-
-\begin_layout Section
-Finite-Element Formulation of Elasticity Equation
-\end_layout
-
-\begin_layout Standard
-We formulate a set of algebraic equations using Galerkin's method.
- We consider (1) a trial solution,
-\begin_inset Formula $\vec{u}$
-\end_inset
-
-, that is a piecewise differentiable vector field and satisfies the Dirichlet
- boundary conditions on
-\begin_inset Formula $S_{u}$
-\end_inset
-
-, and (2) a weighting function,
-\begin_inset Formula $\vec{\phi}$
-\end_inset
-
-, that is a piecewise differentiable vector field and is zero on
-\begin_inset Formula $S_{u}$
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Subsection
-Index Notation
-\end_layout
-
-\begin_layout Standard
-We start with the wave equation (strong form),
-\end_layout
-
-\begin_layout Standard
-\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}
-
-\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}
-\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}
-
-\end_inset
-
- Consider the divergence theorem applied to the dot product of the stress
- tensor and the weighting function,
-\begin_inset Formula $\sigma_{ij}\phi_{i}$
-\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}
-
-\end_inset
-
-Expanding the left-hand side yields
-\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}
-
-\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}
-
-\end_inset
-
-Turning our attention to the second term, we separate the integration over
-
-\begin_inset Formula $S$
-\end_inset
-
- into integration over
-\begin_inset Formula $S_{T}$
-\end_inset
-
- and
-\begin_inset Formula $S_{u}$
-\end_inset
-
- (we will consider tractions over the fault surface,
-\begin_inset Formula $S_{f}$
-\end_inset
-
-, associated with the fault constitutive model in Section
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "sec:fault"
-
-\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}
-
-\end_inset
-
-and recognize that
-\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}
-
-\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}
-
-\end_inset
-
-We express the trial solution and weighting function as linear combinations
- of basis functions,
-\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}
-
-\end_inset
-
-Note that because the trial solution satisfies the Dirichlet boundary condition,
- the number of basis functions for
-\begin_inset Formula $u$
-\end_inset
-
- is generally greater than the number of basis functions for
-\begin_inset Formula $\phi$
-\end_inset
-
-, i.e.,
-\begin_inset Formula $m>n$
-\end_inset
-
-.
- Substituting in the expressions for the trial solution and weighting function
- yields
-\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}
-
-\end_inset
-
- Because the weighting function is arbitrary, this equation must hold for
- all
-\begin_inset Formula $c_{i}^{n}$
-\end_inset
-
-, so that the quantity in parenthesis is zero for each
-\begin_inset Formula $c_{i}^{n}$
-\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}
-
-\end_inset
-
-We want to solve this equation for the unknown coefficients
-\begin_inset Formula $a_{i}^{m}$
-\end_inset
-
- subject to
-\end_layout
-
-\begin_layout Standard
-
-\family roman
-\series medium
-\shape up
-\size normal
-\emph off
-\bar no
-\noun off
-\color none
-\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}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Subsection
-Vector Notation
-\end_layout
-
-\begin_layout Standard
-We start with the wave equation (strong form),
-\end_layout
-
-\begin_layout Standard
-\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}
-
-\end_inset
-
-We construct the weak form by multiplying the wave equation by a weighting
- function and setting the integral over the domain to zero.
- The weighting function is a piecewise differential vector field,
-\begin_inset Formula $\overrightarrow{\phi}$
-\end_inset
-
-, where
-\begin_inset Formula $\overrightarrow{\phi}=0$
-\end_inset
-
- on
-\begin_inset Formula $S_{u}.$
-\end_inset
-
- Hence our weak form is
-\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}
-
-\end_inset
-
- Consider the divergence theorem applied to the dot product of the stress
- tensor and the trial function,
-\begin_inset Formula $\underline{\sigma}\cdot\overrightarrow{\phi}$
-\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}
-
-\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}
-
-\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}
-
-\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}
-
-\end_inset
-
-We separate the integration over
-\begin_inset Formula $S$
-\end_inset
-
- into integration over
-\begin_inset Formula $S_{T}$
-\end_inset
-
- and
-\begin_inset Formula $S_{u}$
-\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}
-
-\end_inset
-
-and recognize that
-\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}
-
-\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}
-
-\end_inset
-
-We express the trial solution and weighting function as linear combinations
- of basis functions,
-\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}
-
-\end_inset
-
-Note that because the weighting function is zero on
-\begin_inset Formula $S_{u}$
-\end_inset
-
-, the number of basis functions for
-\begin_inset Formula $\vec{u}$
-\end_inset
-
- is generally greater than the number of basis functions for
-\begin_inset Formula $\vec{\phi}$
-\end_inset
-
-, i.e.,
-\begin_inset Formula $m>n$
-\end_inset
-
-.
- Substituting in the expressions for the trial solution and weighting function
- yields
-\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}
-
-\end_inset
-
- Because the weighting function is arbitrary, this equation must hold for
- all
-\begin_inset Formula $\overrightarrow{c^{n}}$
-\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}
-
-\end_inset
-
-We want to solve this equation for the unknown coefficients
-\begin_inset Formula $\overrightarrow{a^{m}}$
-\end_inset
-
- subject to
-\end_layout
-
-\begin_layout Standard
-
-\family roman
-\series medium
-\shape up
-\size normal
-\emph off
-\bar no
-\noun off
-\color none
-\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}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Section
-Solution Method for Quasi-Static Problems
-\end_layout
-
-\begin_layout Standard
-For brevity we outline the solution method for quasi-static problems using
- only index notation.
- In quasi-static problems we neglect the inertial terms, so equation
-\begin_inset CommandInset ref
-LatexCommand eqref
-reference "eq:elasticity:integral:discretized"
-
-\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}
-
-\end_inset
-
-As a result, time-dependence only enters through the constitutive relationships
- and the loading conditions.
- 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}
-
-\end_inset
-
-We solve this equation through formulation of a linear algebraic system
- of equations (
-\begin_inset Formula $Au=b$
-\end_inset
-
-), involving the residual (
-\begin_inset Formula $r=b-Au$
-\end_inset
-
-) and Jacobian (
-\begin_inset Formula $A$
-\end_inset
-
-).
- 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}
-
-\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}
-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}
-
-\end_inset
-
-where
-\begin_inset Formula $r_{i}^{n}$
-\end_inset
-
- is an
-\begin_inset Formula $nd$
-\end_inset
-
- vector (
-\begin_inset Formula $d$
-\end_inset
-
- is the dimension of the vector space) and
-\begin_inset Formula $i$
-\end_inset
-
- is a vector space component,
-\begin_inset Formula $x_{q}$
-\end_inset
-
- are the coordinates of the quadrature points,
-\begin_inset Formula $w_{q}$
-\end_inset
-
- are the weights of the quadrature points, and
-\begin_inset Formula $|J_{cell}(x_{q})|$
-\end_inset
-
- is the determinant of the Jacobian matrix evaluated at the quadrature points
- associated with mapping the reference cell to the actual cell.
- The quadrature scheme for the integral over the tractions is one dimension
- lower than the one used in integrating the terms for the volume cells.
-\end_layout
-
-\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}
-
-\end_inset
-
-Isolating the term associated with the increment in stresses yields
-\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}
-
-\end_inset
-
-We associate the term on the left-hand-side with the action of the system
- Jacobian on the increment of the displacement field.
- We approximate the increment in stresses using linear elasticity and infinitesi
-mal strains,
-\end_layout
-
-\begin_layout Standard
-\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}
-
-\end_inset
-
-Now,
-\begin_inset Formula $d\sigma_{ij}\phi_{i,j}$
-\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}
-
-\end_inset
-
-and we know that
-\begin_inset Formula $d\sigma_{ij}$
-\end_inset
-
- is symmetric, so
-\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}
-
-\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}
-
-\end_inset
-
-In terms of the basis functions, we have
-\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}
-
-\end_inset
-
-Combining these expressions for the increment in stresses and making use
- of the symmetry of the weighting functions, we find the system Jacobian
- is
-\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}
-
-\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}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Section
-Solution Method for Dynamic Problems
-\end_layout
-
-\begin_layout Standard
-For brevity we outline the solution method for dynamic problems using only
- index notation.
- Time-dependence enters through the constitutive relationships, loading
- conditions, and the inertial terms.
- We consider the deformation at time
-\begin_inset Formula $t$
-\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}
-
-\end_inset
-
-We solve this equation through formulation of a linear algebraic system
- of equations (
-\begin_inset Formula $Au=b$
-\end_inset
-
-), involving the residual (
-\begin_inset Formula $r=b-Au$
-\end_inset
-
-) and Jacobian (
-\begin_inset Formula $A$
-\end_inset
-
-).
- 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}
-
-\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}
-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}
-
-\end_inset
-
-where
-\begin_inset Formula $x_{q}$
-\end_inset
-
- are the coordinates of the quadrature points,
-\begin_inset Formula $w_{q}$
-\end_inset
-
- are the weights of the quadrature points, and
-\begin_inset Formula $|J_{cell}(x_{q})|$
-\end_inset
-
- is the determinant of the Jacobian matrix evaluated at the quadrature points
- associated with mapping the reference cell to the actual cell.
- The quadrature scheme for the integral over the tractions is one dimension
- lower than the one used in integrating the terms for the volume cells.
-
-\end_layout
-
-\begin_layout Standard
-We find the system Jacobian matrix by making use of the temporal discretization
- and isolating the term for the increment in the displacement field at time
-
-\begin_inset Formula $t$
-\end_inset
-
-.
- Using the central difference method to approximate the acceleration (and
- velocity),
-\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}
-
-\end_inset
-
-and writing the displacement at time
-\begin_inset Formula $t+\Delta t$
-\end_inset
-
- in terms of the displacement at
-\begin_inset Formula $t$
-\end_inset
-
- (for consistency with the displacement increment quasi-static formulation),
-\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}
-
-\end_inset
-
-Substituting into equation
-\begin_inset CommandInset ref
-LatexCommand eqref
-reference "eq:elasticity:integral:dynamic:t"
-
-\end_inset
-
- yields
-\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}
-
-\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}
-
-\end_inset
-
-and using numerical quadrature in the finite-element discretization to replace
- the integrals with sums over the cells and quadrature points,
-\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}
-
-\end_inset
-
-where
-\begin_inset Formula $A_{ij}^{mn}$
-\end_inset
-
- is a
-\begin_inset Formula $nd$
-\end_inset
-
- by
-\begin_inset Formula $md$
-\end_inset
-
- matrix (
-\begin_inset Formula $d$
-\end_inset
-
- is the dimension of the vector space),
-\begin_inset Formula $m$
-\end_inset
-
- and
-\begin_inset Formula $n$
-\end_inset
-
- refer to the basis functions and
-\begin_inset Formula $i$
-\end_inset
-
- and
-\begin_inset Formula $j$
-\end_inset
-
- are vector space components.
- We consider the contributions associated with the fault in section
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "sec:fault"
-
-\end_inset
-
- and with absorbing boundaries is section
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "sec:absorbing:boundaries"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Section
-Small Strain Formulation
-\begin_inset CommandInset label
-LatexCommand label
-name "sec:Small-Strain-Formulation"
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-In some crustal deformation problems sufficient deformation may occur that
- the assumptions associated with infinitesimal strains no longer hold.
- This is often the case for problems when one wants to include the effects
- of gravitational body forces and deformation on the overburden pressure.
- In such cases we want to account for both rigid body motion and small strains.
- The elasticity formulation in PyLith for small strains uses the Green-Lagrange
- strain tensor and the Second Piola-Kirchhoff stress tensor as is based
- on the one presented by Bathe
-\begin_inset CommandInset citation
-LatexCommand cite
-key "Bathe:1995"
-
-\end_inset
-
-.
- The Green-Lagrange strain provides a measure of the strain relative to
- the original, undeformed configuration.
-\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}
-
-\end_inset
-
-and
-\begin_inset Formula $X_{ij}$
-\end_inset
-
- is the deformation tensor.
- The Second Piola-Kirchhoff stress tensor,
-\begin_inset Formula $S_{ij}$
-\end_inset
-
-, is related to the Green-Lagrange strain tensor through the elasticity
- constants,
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula
-\begin{equation}
-S_{ij}=C_{ijkl}\varepsilon_{kl},
-\end{equation}
-
-\end_inset
-
-in the same manner as in the infinitesimal strain formulation.
-\end_layout
-
-\begin_layout Standard
-The elasticity integral in the finite-element formulation includes additional
- terms when we account for small strains.
- Recognizing the similarity between the weighting function and an increment
- in strain in the infinitesimal formulation (many finite-element texts derive
- the finite-element formulation for elasticity using the Principle of Virtual
- Work), we replace
-\begin_inset Formula $\int_{V}\sigma_{ij}\phi_{i,j}\: dV$
-\end_inset
-
- with
-\begin_inset Formula $\int_{V}S_{ij}\delta\varepsilon_{ij}\: dV$
-\end_inset
-
- in equation
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "eq:elasticity:integral"
-
-\end_inset
-
-, where
-\begin_inset Formula $\delta\varepsilon_{ij}$
-\end_inset
-
- is the
-\begin_inset Quotes eld
-\end_inset
-
-virtual
-\begin_inset Quotes erd
-\end_inset
-
- 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}
-
-\end_inset
-
-Writing the displacements in terms of the basis functions and forcing the
- terms associated with the arbitrary weighting function (
-\begin_inset Quotes eld
-\end_inset
-
-virtual
-\begin_inset Quotes erd
-\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}
-
-\end_inset
-
-Thus, we have one additional term (the second term) compared with the residual
- for infinitesimal strains.
- Just as in the infinitesimal formulation, we evaluate the integral over
- the volume using numerical quadrature with sums over the quadrature points
- of each cell.
-\end_layout
-
-\begin_layout Subsection
-Quasi-static Problems
-\end_layout
-
-\begin_layout Standard
-The system Jacobian for quasi-static problems includes terms associated
- with elasticity.
- For the small strain formulation, we write the elasticity term at time
-
-\begin_inset Formula $t+\Delta t$
-\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}
-
-\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}
-
-\end_inset
-
-and the increment in the
-\begin_inset Quotes eld
-\end_inset
-
-virtual
-\begin_inset Quotes erd
-\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}
-
-\end_inset
-
-We associate the system Jacobian with the terms involving the increment
- in displacements.
- After substituting in the expressions for the increment in the stresses
- and the increment in the
-\begin_inset Quotes eld
-\end_inset
-
-virtual
-\begin_inset Quotes erd
-\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}
-
-\end_inset
-
-The small strain formulation produces additional terms associated with the
- elastic constants and new a new term associated with the stress tensor.
-\end_layout
-
-\begin_layout Subsection
-Dynamic Problems
-\end_layout
-
-\begin_layout Standard
-The system Jacobian matrix in dynamic problems does not include any terms
- associated with elasticity, so the system Jacobian matrix in the small
- strain formulation matches the one used in the infinitesimal strain formulation.
-\end_layout
-
-\end_body
-\end_document
+#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
+\papersize default
+\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
+\bottommargin 1in
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\paragraph_indentation default
+\quotes_language english
+\papercolumns 1
+\papersides 1
+\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:Governing-Equations"
+
+\end_inset
+
+Governing Equations
+\end_layout
+
+\begin_layout Standard
+We present here a brief derivation of the equations for both quasi-static
+ and dynamic computations.
+ Since the general equations are the same (except for the absence of inertial
+ terms in the quasi-static case), we first derive these equations.
+ We then present solution methods for each specific case.
+ In all of our derivations, we use the notation described in Table
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "tab:Mathematical-notation"
+
+\end_inset
+
+ for both index and vector notation.
+ When using index notation, we use the common convention that repeated indices
+ indicate summation over the range of the index.
+\end_layout
+
+\begin_layout Standard
+\noindent
+\align center
+\begin_inset Float table
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Plain Layout
+\begin_inset Caption
+
+\begin_layout Plain Layout
+\begin_inset CommandInset label
+LatexCommand label
+name "tab:Mathematical-notation"
+
+\end_inset
+
+Mathematical notation
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\noindent
+\align center
+\begin_inset VSpace defskip
+\end_inset
+
+
+\begin_inset Tabular
+<lyxtabular version="3" rows="11" columns="3">
+<features tabularvalignment="middle">
+<column alignment="center" valignment="top" width="0">
+<column alignment="center" valignment="top" width="0">
+<column alignment="center" valignment="top" width="0">
+<row>
+<cell multicolumn="1" alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Symbol
+\end_layout
+
+\end_inset
+</cell>
+<cell multicolumn="2" alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+
+\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
+Description
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Index notation
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Vector Notation
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $a_{i}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\overrightarrow{a}$
+\end_inset
+
+
+\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
+Vector field a
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $a_{ij}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\underline{a}$
+\end_inset
+
+
+\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
+Second order tensor field a
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $u_{i}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\overrightarrow{u}$
+\end_inset
+
+
+\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
+Displacement vector field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $d_{i}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\vec{{d}}$
+\end_inset
+
+
+\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
+Fault slip vector field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $f_{i}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\overrightarrow{f}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Body force vector field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $T_{i}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\overrightarrow{T}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Traction vector field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\sigma_{ij}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\underline{\sigma}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Stress tensor field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $n_{i}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\overrightarrow{n}$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Normal vector field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\rho$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+\begin_inset Formula $\rho$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" bottomline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Mass density scalar field
+\end_layout
+
+\end_inset
+</cell>
+</row>
+</lyxtabular>
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Section
+Derivation of Elasticity Equation
+\end_layout
+
+\begin_layout Subsection
+Index Notation
+\end_layout
+
+\begin_layout Standard
+Consider volume
+\begin_inset Formula $V$
+\end_inset
+
+ bounded by surface
+\begin_inset Formula $S$
+\end_inset
+
+.
+ Applying a Lagrangian description of the conservation of momentum gives
+\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}
+
+\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}
+
+\end_inset
+
+where
+\begin_inset Formula $n_{j}$
+\end_inset
+
+ is the vector normal to
+\begin_inset Formula $S$
+\end_inset
+
+.
+ Substituting into Equation
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eqn:momentum:index"
+
+\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}
+
+\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}
+
+\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}
+
+\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}
+
+\end_inset
+
+Because the volume
+\begin_inset Formula $V$
+\end_inset
+
+ is arbitrary, the integrand must be zero at every location in the volume,
+ so that we end up with
+\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}
+
+\end_inset
+
+We specify tractions,
+\begin_inset Formula $T_{i}$
+\end_inset
+
+, on surface
+\begin_inset Formula $S_{f}$
+\end_inset
+
+, displacements,
+\begin_inset Formula $u_{i}^{o}$
+\end_inset
+
+, on surface
+\begin_inset Formula $S_{u}$
+\end_inset
+
+, and slip,
+\begin_inset Formula $d_{k}$
+\end_inset
+
+, on fault surface
+\begin_inset Formula $S_{f}$
+\end_inset
+
+ (we will consider the case of fault constitutive models in Section
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sec:fault"
+
+\end_inset
+
+).
+ The rotation matrix
+\begin_inset Formula $R_{ki}$
+\end_inset
+
+ transforms vectors from the global coordinate system to the fault coordinate
+ system.
+ Note that since both
+\begin_inset Formula $T_{i}$
+\end_inset
+
+ and
+\begin_inset Formula $u_{i}$
+\end_inset
+
+ are vector quantities, there can be some spatial overlap of the surfaces
+
+\begin_inset Formula $S_{T}$
+\end_inset
+
+ and
+\begin_inset Formula $S_{u}$
+\end_inset
+
+; however, the same degree of freedom cannot simultaneously have both types
+ of boundary conditions.
+\end_layout
+
+\begin_layout Subsection
+Vector Notation
+\end_layout
+
+\begin_layout Standard
+Consider volume
+\begin_inset Formula $V$
+\end_inset
+
+ bounded by surface
+\begin_inset Formula $S$
+\end_inset
+
+.
+ Applying a Lagrangian description of the conservation of momentum gives
+\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}
+
+\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}
+
+\end_inset
+
+where
+\begin_inset Formula $\overrightarrow{n}$
+\end_inset
+
+ is the vector normal to
+\begin_inset Formula $S$
+\end_inset
+
+.
+ Substituting into Equation
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eqn:momentum:vec"
+
+\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}
+
+\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}
+
+\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}
+
+\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}
+
+\end_inset
+
+Because the volume
+\begin_inset Formula $V$
+\end_inset
+
+ 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}
+\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}
+
+\end_inset
+
+We specify tractions,
+\begin_inset Formula $\vec{T}$
+\end_inset
+
+, on surface
+\begin_inset Formula $S_{f}$
+\end_inset
+
+, displacements,
+\begin_inset Formula $\overrightarrow{u^{o}}$
+\end_inset
+
+, on surface
+\begin_inset Formula $S_{u}$
+\end_inset
+
+, and slip,
+\begin_inset Formula $\vec{d}$
+\end_inset
+
+, on fault surface
+\begin_inset Formula $S_{f}$
+\end_inset
+
+ (we will consider the case of fault constitutive models in Section
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sec:fault"
+
+\end_inset
+
+).
+ The rotation matrix
+\begin_inset Formula $\underline{R}$
+\end_inset
+
+ transforms vectors from the global coordinate system to the fault coordinate
+ system.
+ Note that since both
+\begin_inset Formula $\overrightarrow{T}$
+\end_inset
+
+ and
+\begin_inset Formula $\overrightarrow{u}$
+\end_inset
+
+ are vector quantities, there can be some spatial overlap of the surfaces
+
+\begin_inset Formula $S_{T}$
+\end_inset
+
+ and
+\begin_inset Formula $S_{u}$
+\end_inset
+
+; however, the same degree of freedom cannot simultaneously have both types
+ of boundary conditions.
+\end_layout
+
+\begin_layout Section
+Finite-Element Formulation of Elasticity Equation
+\end_layout
+
+\begin_layout Standard
+We formulate a set of algebraic equations using Galerkin's method.
+ We consider (1) a trial solution,
+\begin_inset Formula $\vec{u}$
+\end_inset
+
+, that is a piecewise differentiable vector field and satisfies the Dirichlet
+ boundary conditions on
+\begin_inset Formula $S_{u}$
+\end_inset
+
+, and (2) a weighting function,
+\begin_inset Formula $\vec{\phi}$
+\end_inset
+
+, that is a piecewise differentiable vector field and is zero on
+\begin_inset Formula $S_{u}$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Subsection
+Index Notation
+\end_layout
+
+\begin_layout Standard
+We start with the wave equation (strong form),
+\end_layout
+
+\begin_layout Standard
+\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}
+
+\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}
+\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}
+
+\end_inset
+
+ Consider the divergence theorem applied to the dot product of the stress
+ tensor and the weighting function,
+\begin_inset Formula $\sigma_{ij}\phi_{i}$
+\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}
+
+\end_inset
+
+Expanding the left-hand side yields
+\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}
+
+\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}
+
+\end_inset
+
+Turning our attention to the second term, we separate the integration over
+
+\begin_inset Formula $S$
+\end_inset
+
+ into integration over
+\begin_inset Formula $S_{T}$
+\end_inset
+
+ and
+\begin_inset Formula $S_{u}$
+\end_inset
+
+ (we will consider tractions over the fault surface,
+\begin_inset Formula $S_{f}$
+\end_inset
+
+, associated with the fault constitutive model in Section
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sec:fault"
+
+\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}
+
+\end_inset
+
+and recognize that
+\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}
+
+\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}
+
+\end_inset
+
+We express the trial solution and weighting function as linear combinations
+ of basis functions,
+\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}
+
+\end_inset
+
+Note that because the trial solution satisfies the Dirichlet boundary condition,
+ the number of basis functions for
+\begin_inset Formula $u$
+\end_inset
+
+ is generally greater than the number of basis functions for
+\begin_inset Formula $\phi$
+\end_inset
+
+, i.e.,
+\begin_inset Formula $m>n$
+\end_inset
+
+.
+ Substituting in the expressions for the trial solution and weighting function
+ yields
+\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}
+
+\end_inset
+
+ Because the weighting function is arbitrary, this equation must hold for
+ all
+\begin_inset Formula $c_{i}^{n}$
+\end_inset
+
+, so that the quantity in parenthesis is zero for each
+\begin_inset Formula $c_{i}^{n}$
+\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}
+
+\end_inset
+
+We want to solve this equation for the unknown coefficients
+\begin_inset Formula $a_{i}^{m}$
+\end_inset
+
+ subject to
+\end_layout
+
+\begin_layout Standard
+
+\family roman
+\series medium
+\shape up
+\size normal
+\emph off
+\bar no
+\noun off
+\color none
+\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}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Subsection
+Vector Notation
+\end_layout
+
+\begin_layout Standard
+We start with the wave equation (strong form),
+\end_layout
+
+\begin_layout Standard
+\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}
+
+\end_inset
+
+We construct the weak form by multiplying the wave equation by a weighting
+ function and setting the integral over the domain to zero.
+ The weighting function is a piecewise differential vector field,
+\begin_inset Formula $\overrightarrow{\phi}$
+\end_inset
+
+, where
+\begin_inset Formula $\overrightarrow{\phi}=0$
+\end_inset
+
+ on
+\begin_inset Formula $S_{u}.$
+\end_inset
+
+ Hence our weak form is
+\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}
+
+\end_inset
+
+ Consider the divergence theorem applied to the dot product of the stress
+ tensor and the trial function,
+\begin_inset Formula $\underline{\sigma}\cdot\overrightarrow{\phi}$
+\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}
+
+\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}
+
+\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}
+
+\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}
+
+\end_inset
+
+We separate the integration over
+\begin_inset Formula $S$
+\end_inset
+
+ into integration over
+\begin_inset Formula $S_{T}$
+\end_inset
+
+ and
+\begin_inset Formula $S_{u}$
+\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}
+
+\end_inset
+
+and recognize that
+\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}
+
+\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}
+
+\end_inset
+
+We express the trial solution and weighting function as linear combinations
+ of basis functions,
+\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}
+
+\end_inset
+
+Note that because the weighting function is zero on
+\begin_inset Formula $S_{u}$
+\end_inset
+
+, the number of basis functions for
+\begin_inset Formula $\vec{u}$
+\end_inset
+
+ is generally greater than the number of basis functions for
+\begin_inset Formula $\vec{\phi}$
+\end_inset
+
+, i.e.,
+\begin_inset Formula $m>n$
+\end_inset
+
+.
+ Substituting in the expressions for the trial solution and weighting function
+ yields
+\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}
+
+\end_inset
+
+ Because the weighting function is arbitrary, this equation must hold for
+ all
+\begin_inset Formula $\overrightarrow{c^{n}}$
+\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}
+
+\end_inset
+
+We want to solve this equation for the unknown coefficients
+\begin_inset Formula $\overrightarrow{a^{m}}$
+\end_inset
+
+ subject to
+\end_layout
+
+\begin_layout Standard
+
+\family roman
+\series medium
+\shape up
+\size normal
+\emph off
+\bar no
+\noun off
+\color none
+\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}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Section
+Solution Method for Quasi-Static Problems
+\end_layout
+
+\begin_layout Standard
+For brevity we outline the solution method for quasi-static problems using
+ only index notation.
+ In quasi-static problems we neglect the inertial terms, so Equation
+\begin_inset CommandInset ref
+LatexCommand eqref
+reference "eq:elasticity:integral:discretized"
+
+\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}
+
+\end_inset
+
+As a result, time-dependence only enters through the constitutive relationships
+ and the loading conditions.
+ 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}
+
+\end_inset
+
+We solve this equation through formulation of a linear algebraic system
+ of equations (
+\begin_inset Formula $Au=b$
+\end_inset
+
+), involving the residual (
+\begin_inset Formula $r=b-Au$
+\end_inset
+
+) and Jacobian (
+\begin_inset Formula $A$
+\end_inset
+
+).
+ 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}
+
+\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}
+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}
+
+\end_inset
+
+where
+\begin_inset Formula $r_{i}^{n}$
+\end_inset
+
+ is an
+\begin_inset Formula $nd$
+\end_inset
+
+ vector (
+\begin_inset Formula $d$
+\end_inset
+
+ is the dimension of the vector space) and
+\begin_inset Formula $i$
+\end_inset
+
+ is a vector space component,
+\begin_inset Formula $x_{q}$
+\end_inset
+
+ are the coordinates of the quadrature points,
+\begin_inset Formula $w_{q}$
+\end_inset
+
+ are the weights of the quadrature points, and
+\begin_inset Formula $|J_{cell}(x_{q})|$
+\end_inset
+
+ is the determinant of the Jacobian matrix evaluated at the quadrature points
+ associated with mapping the reference cell to the actual cell.
+ The quadrature scheme for the integral over the tractions is one dimension
+ lower than the one used in integrating the terms for the volume cells.
+\end_layout
+
+\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}
+
+\end_inset
+
+Isolating the term associated with the increment in stresses yields
+\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}
+
+\end_inset
+
+We associate the term on the left-hand-side with the action of the system
+ Jacobian on the increment of the displacement field.
+ We approximate the increment in stresses using linear elasticity and infinitesi
+mal strains,
+\end_layout
+
+\begin_layout Standard
+\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}
+
+\end_inset
+
+Now,
+\begin_inset Formula $d\sigma_{ij}\phi_{i,j}$
+\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}
+
+\end_inset
+
+and we know that
+\begin_inset Formula $d\sigma_{ij}$
+\end_inset
+
+ is symmetric, so
+\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}
+
+\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}
+
+\end_inset
+
+In terms of the basis functions, we have
+\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}
+
+\end_inset
+
+Combining these expressions for the increment in stresses and making use
+ of the symmetry of the weighting functions, we find the system Jacobian
+ is
+\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}
+
+\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}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Section
+Solution Method for Dynamic Problems
+\end_layout
+
+\begin_layout Standard
+For brevity we outline the solution method for dynamic problems using only
+ index notation.
+ Time-dependence enters through the constitutive relationships, loading
+ conditions, and the inertial terms.
+ We consider the deformation at time
+\begin_inset Formula $t$
+\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}
+
+\end_inset
+
+We solve this equation through formulation of a linear algebraic system
+ of equations (
+\begin_inset Formula $Au=b$
+\end_inset
+
+), involving the residual (
+\begin_inset Formula $r=b-Au$
+\end_inset
+
+) and Jacobian (
+\begin_inset Formula $A$
+\end_inset
+
+).
+ 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}
+
+\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}
+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}
+
+\end_inset
+
+where
+\begin_inset Formula $x_{q}$
+\end_inset
+
+ are the coordinates of the quadrature points,
+\begin_inset Formula $w_{q}$
+\end_inset
+
+ are the weights of the quadrature points, and
+\begin_inset Formula $|J_{cell}(x_{q})|$
+\end_inset
+
+ is the determinant of the Jacobian matrix evaluated at the quadrature points
+ associated with mapping the reference cell to the actual cell.
+ The quadrature scheme for the integral over the tractions is one dimension
+ lower than the one used in integrating the terms for the volume cells.
+
+\end_layout
+
+\begin_layout Standard
+We find the system Jacobian matrix by making use of the temporal discretization
+ and isolating the term for the increment in the displacement field at time
+
+\begin_inset Formula $t$
+\end_inset
+
+.
+ Using the central difference method to approximate the acceleration (and
+ velocity),
+\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}
+
+\end_inset
+
+and writing the displacement at time
+\begin_inset Formula $t+\Delta t$
+\end_inset
+
+ in terms of the displacement at
+\begin_inset Formula $t$
+\end_inset
+
+ (for consistency with the displacement increment quasi-static formulation),
+\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}
+
+\end_inset
+
+Substituting into Equation
+\begin_inset CommandInset ref
+LatexCommand eqref
+reference "eq:elasticity:integral:dynamic:t"
+
+\end_inset
+
+ yields
+\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}
+
+\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}
+
+\end_inset
+
+and using numerical quadrature in the finite-element discretization to replace
+ the integrals with sums over the cells and quadrature points,
+\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}
+
+\end_inset
+
+where
+\begin_inset Formula $A_{ij}^{mn}$
+\end_inset
+
+ is an
+\begin_inset Formula $nd$
+\end_inset
+
+ by
+\begin_inset Formula $md$
+\end_inset
+
+ matrix (
+\begin_inset Formula $d$
+\end_inset
+
+ is the dimension of the vector space),
+\begin_inset Formula $m$
+\end_inset
+
+ and
+\begin_inset Formula $n$
+\end_inset
+
+ refer to the basis functions and
+\begin_inset Formula $i$
+\end_inset
+
+ and
+\begin_inset Formula $j$
+\end_inset
+
+ are vector space components.
+ We consider the contributions associated with the fault in Section
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sec:fault"
+
+\end_inset
+
+ and with absorbing boundaries is Section
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "sec:absorbing:boundaries"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Section
+Small Strain Formulation
+\begin_inset CommandInset label
+LatexCommand label
+name "sec:Small-Strain-Formulation"
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+In some crustal deformation problems sufficient deformation may occur that
+ the assumptions associated with infinitesimal strains no longer hold.
+ This is often the case for problems when one wants to include the effects
+ of gravitational body forces and deformation on the overburden pressure.
+ In such cases we want to account for both rigid body motion and small strains.
+ The elasticity formulation in PyLith for small strains uses the Green-Lagrange
+ strain tensor and the Second Piola-Kirchhoff stress tensor as is based
+ on the one presented by Bathe
+\begin_inset CommandInset citation
+LatexCommand cite
+key "Bathe:1995"
+
+\end_inset
+
+.
+ The Green-Lagrange strain provides a measure of the strain relative to
+ the original, undeformed configuration.
+\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}
+
+\end_inset
+
+and
+\begin_inset Formula $X_{ij}$
+\end_inset
+
+ is the deformation tensor.
+ The Second Piola-Kirchhoff stress tensor,
+\begin_inset Formula $S_{ij}$
+\end_inset
+
+, is related to the Green-Lagrange strain tensor through the elasticity
+ constants,
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula
+\begin{equation}
+S_{ij}=C_{ijkl}\varepsilon_{kl},
+\end{equation}
+
+\end_inset
+
+in the same manner as in the infinitesimal strain formulation.
+\end_layout
+
+\begin_layout Standard
+The elasticity integral in the finite-element formulation includes additional
+ terms when we account for small strains.
+ Recognizing the similarity between the weighting function and an increment
+ in strain in the infinitesimal formulation (many finite-element texts derive
+ the finite-element formulation for elasticity using the Principle of Virtual
+ Work), we replace
+\begin_inset Formula $\int_{V}\sigma_{ij}\phi_{i,j}\: dV$
+\end_inset
+
+ with
+\begin_inset Formula $\int_{V}S_{ij}\delta\varepsilon_{ij}\: dV$
+\end_inset
+
+ in Equation
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:elasticity:integral"
+
+\end_inset
+
+, where
+\begin_inset Formula $\delta\varepsilon_{ij}$
+\end_inset
+
+ is the
+\begin_inset Quotes eld
+\end_inset
+
+virtual
+\begin_inset Quotes erd
+\end_inset
+
+ 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}
+
+\end_inset
+
+Writing the displacements in terms of the basis functions and forcing the
+ terms associated with the arbitrary weighting function (
+\begin_inset Quotes eld
+\end_inset
+
+virtual
+\begin_inset Quotes erd
+\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}
+
+\end_inset
+
+Thus, we have one additional term (the second term) compared with the residual
+ for infinitesimal strains.
+ Just as in the infinitesimal formulation, we evaluate the integral over
+ the volume using numerical quadrature with sums over the quadrature points
+ of each cell.
+\end_layout
+
+\begin_layout Subsection
+Quasi-static Problems
+\end_layout
+
+\begin_layout Standard
+The system Jacobian for quasi-static problems includes terms associated
+ with elasticity.
+ For the small strain formulation, we write the elasticity term at time
+
+\begin_inset Formula $t+\Delta t$
+\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}
+
+\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}
+
+\end_inset
+
+and the increment in the
+\begin_inset Quotes eld
+\end_inset
+
+virtual
+\begin_inset Quotes erd
+\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}
+
+\end_inset
+
+We associate the system Jacobian with the terms involving the increment
+ in displacements.
+ After substituting in the expressions for the increment in the stresses
+ and the increment in the
+\begin_inset Quotes eld
+\end_inset
+
+virtual
+\begin_inset Quotes erd
+\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}
+
+\end_inset
+
+The small strain formulation produces additional terms associated with the
+ elastic constants and new a new term associated with the stress tensor.
+\end_layout
+
+\begin_layout Subsection
+Dynamic Problems
+\end_layout
+
+\begin_layout Standard
+The system Jacobian matrix in dynamic problems does not include any terms
+ associated with elasticity, so the system Jacobian matrix in the small
+ strain formulation matches the one used in the infinitesimal strain formulation.
+\end_layout
+
+\end_body
+\end_document
Modified: short/3D/PyLith/branches/v1.7-stable/doc/userguide/intro/intro.lyx
===================================================================
--- short/3D/PyLith/branches/v1.7-stable/doc/userguide/intro/intro.lyx 2012-06-26 16:08:18 UTC (rev 20413)
+++ short/3D/PyLith/branches/v1.7-stable/doc/userguide/intro/intro.lyx 2012-06-26 16:08:27 UTC (rev 20414)
@@ -1,542 +1,542 @@
-#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 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
-\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
-Introduction
-\end_layout
-
-\begin_layout Section
-Overview
-\end_layout
-
-\begin_layout Standard
-PyLith is a multi-scale simulation software package for earthquake physics.
- It is portable, scalable software for simulation of crustal deformation
- across spatial scales ranging from meters to hundreds of kilometers and
- temporal scales ranging from milliseconds to thousands of years.
-\end_layout
-
-\begin_layout Section
-History
-\end_layout
-
-\begin_layout Standard
-PyLith 1.0 was the first version to allow the solution of both implicit (quasi-st
-atic) and explicit (dynamic) problems and was a complete rewrite of the
- original PyLith (version 0.8).
- PyLith 1.0 combines the functionality of EqSim
-\begin_inset CommandInset citation
-LatexCommand cite
-key "Aagaard:etal:2001a,Aagaard:etal:2001b"
-
-\end_inset
-
- and PyLith 0.8.
- PyLith 0.8 was a direct descendant of LithoMop and was the first version
- that ran in parallel, as well as providing several other improvements over
- LithoMop.
- LithoMop was the product of major reengineering of Tecton, a finite-element
- code for simulating static and quasi-static crustal deformation.
- The major new features present in LithoMop included dynamic memory allocation
- and the use of the Pyre simulation framework and PETSc solvers.
- EqSim was written by Brad Aagaard to solve problems in earthquake dynamics,
- including rupture propagation and seismic wave propagation.
-\end_layout
-
-\begin_layout Standard
-The release of PyLith 1.0 has been followed by additional releases that expand
- the number of features as well as improve performance.
- The PyLith 1.x series of releases allows the solution of both quasi-static
- and dynamic problems in one, two, or three dimensions.
- The code runs in either serial or parallel, and the design allows for relativel
-y easy scripting using the Python programming language.
- Material properties and values for boundary and fault conditions are specified
- using spatial databases, which permit easy prescription of complex spatial
- variations of properties and parameters.
- Simulation parameters are generally specified through the use of simple
- ASCII files or the command line.
- At present, mesh information may be provided using a simple ASCII file
- (PyLith mesh ASCII format) or imported from CUBIT or LaGriT, two widely-used
- meshing packages.
- The elements currently available include a linear bar in 1-D, linear triangles
- and quadrilaterals in 2-D, and linear tetrahedra and hexahedra in 3-D.
- Materials presently available include isotropic elastic, linear Maxwell
- viscoelastic, generalized Maxwell viscoelastic, power-law viscoelastic,
- and Drucker-Prager elastoplastic.
- Boundary conditions include Dirichlet (prescribed displacements and velocities)
-, Neumann (traction), point forces, and absorbing boundaries.
- Cohesive elements are used to implement slip across interior surfaces (faults)
- with both kinematically-specified fault slip and slip governed by fault
- constitutive models.
- PyLith also includes an interface for computing static Green's functions
- for fault slip.
-\end_layout
-
-\begin_layout Standard
-PyLith is under active development and we expect a number of additions and
- improvements in the near future.
- Likely enhancements will include additional bulk and fault constitutive
- models, coupled quasi-static and dynamic simulations for earthquake cycle
- modeling, and coupling between elasticity, heat flow, and/or fluid flow.
-\end_layout
-
-\begin_layout Section
-PyLith Workflow
-\end_layout
-
-\begin_layout Standard
-PyLith is one component in the process of investigating problems in tectonics
- (Figure
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "fig:Workflow-summary"
-
-\end_inset
-
-).
- Given a geological problem of interest, a scientist must first provide
- a geometrical representation of the desired structure.
- Once the structure has been defined, a computational mesh must be created.
- PyLith presently provides three mesh importing options: CUBIT Exodus format,
- LaGriT GMV and Pset files, and PyLith mesh ASCII format.
- The modeling of the physical processes of interest is performed by a code
- such as PyLith.
- Present output consists of VTK or HDF5/Xdmf files which can be used by
- a number of visualization codes (e.g., ParaView, Visit, MayaVi, and Matlab).
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status open
-
-\begin_layout Plain Layout
-\align center
-\begin_inset Graphics
- filename figs/workflow.eps
- scale 67
- keepAspectRatio
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Plain Layout
-\begin_inset Caption
-
-\begin_layout Plain Layout
-\begin_inset CommandInset label
-LatexCommand label
-name "fig:Workflow-summary"
-
-\end_inset
-
-Workflow involved in going from geologic structure to problem analysis.
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Section
-PyLith Design
-\end_layout
-
-\begin_layout Standard
-PyLith is separated into modules to encapsulate behavior and facilitate
- use across multiple applications.
- This allows expert users to replace functionality of a wide variety of
- components without recompiling or polluting the main code.
- PyLith employs external packages (see Figure
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "fig:pylith-dependencies"
-
-\end_inset
-
-) to reduce development time and enhance computational efficiency; for example,
- PyLith 0.8 ran two times faster when the PETSc linear solver was used.
-\end_layout
-
-\begin_layout Standard
-\noindent
-\align center
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status open
-
-\begin_layout Plain Layout
-\align center
-\begin_inset Graphics
- filename figs/packages.eps
- scale 40
-
-\end_inset
-
-
-\begin_inset Caption
-
-\begin_layout Plain Layout
-\begin_inset CommandInset label
-LatexCommand label
-name "fig:pylith-dependencies"
-
-\end_inset
-
-PyLith dependencies.
- PyLith makes direct use of several other packages, some of which have their
- own dependencies.
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-PyLith is written in two programming languages.
- High-level code is written in Python; this rich, expressive interpreted
- language with dynamic typing reduces development time and permits flexible
- addition of user-contributed modules.
- This high-level code makes use of Pyre, a science-neutral simulation framework
- developed at Caltech, to link the modules together at runtime and gather
- user-input.
- Low-level code is written in C++, providing fast execution while still
- allowing an object-oriented implementation.
- This low-level code relies on PETSc to perform operations on matrices and
- vectors in parallel.
- We also make extensive use of two Python packages.
- SWIG is a package that simplifies the task of adding C++ extensions to
- Python code, and FIAT provides tabulated basis functions and numerical
- quadrature points.
-
-\end_layout
-
-\begin_layout Standard
-In writing PyLith 1.0, the code was designed to be object-oriented and modular.
- Each type of module is accessed through a specified interface (set of functions
-).
- This permits adding, replacing, and rewriting modules without affecting
- other parts of the code.
- This code structure simplifies code maintenance and development.
- Extending the set of code features is also easier, since developers can
- create new modules derived from the existing ones.
-\end_layout
-
-\begin_layout Standard
-The new code design leverages Pyre, Sieve, and PETSc much more extensively
- than the previous version.
- Pyre is used to glue together the various modules used to construct a
- simulation and specify the parameters.
- Sieve is used for all finite-element storage and manipulation and handles
- the creation of the PETSc matrices and vectors.
- As a result, most of the PyLith source code pertains to implementing the
- geodynamics, such as bulk rheology, boundary conditions, and slip on faults.
-
-\end_layout
-
-\begin_layout Standard
-PyLith also uses FIAT to tabulate the finite-element basis functions at
- the numerical integration (quadrature) points.
- Nemesis allows PyLith to run Python using the Message Passing Interface
- (MPI) for parallel processing.
- Additional, indirect dependencies (see Figure
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "fig:pylith-dependencies"
-
-\end_inset
-
-) include numpy (efficient operations on numerical arrays in Python), Proj.4
- (geographic projections), and SWIG (calling C++ functions from Python).
-\end_layout
-
-\begin_layout Standard
-During development, tests were constructed for nearly every module function.
- These unit tests are distributed with the source code.
- These tests are run throughout the development cycle to expose bugs and
- isolate their origin.
- As additional changes are made to the code, the tests are rerun to help
- prevent introduction of new bugs.
- A number of simple, full-scale tests, such as axial compression and extension,
- simple shear, and slip on through-going faults, have been used to test
- the code.
- Additionally, we have run the Southern California Earthquake Center crustal
- deformation and several of the spontaneous rupture benchmarks for strike-slip
- and reverse-slip to determine the relative local and global error (see
- Chapter
-\begin_inset CommandInset ref
-LatexCommand ref
-reference "cha:Benchmarks"
-
-\end_inset
-
-).
-\end_layout
-
-\begin_layout Subsection
-Pyre
-\end_layout
-
-\begin_layout Standard
-Pyre is an object-oriented environment capable of specifying and launching
- numerical simulations on multiple platforms, including Beowulf-class parallel
- computers and grid computing systems.
- Pyre allows the binding of multiple components such as solid and fluid
- models used in Earth science simulations, and different meshers.
- The Pyre framework enables the elegant setup, modification and launching
- of massively parallel solver applications.
-\end_layout
-
-\begin_layout Standard
-\noindent
-\align center
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status open
-
-\begin_layout Plain Layout
-\align center
-\begin_inset Graphics
- filename figs/pyre_overview.png
- width 4in
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Plain Layout
-\begin_inset Caption
-
-\begin_layout Plain Layout
-\begin_inset CommandInset label
-LatexCommand label
-name "fig:Pyre:Architecture"
-
-\end_inset
-
-Pyre Architecture.
- The integration framework is a set of cooperating abstract services.
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-Pyre is a framework, a combination of software and design philosophy that
- promotes the reuse of code.
- In their canonical software design book,
-\emph on
-Design Patterns
-\emph default
-, Erich Gamma
-\shape italic
-et al
-\shape default
-.
- condense the concept of a framework concept down to, ``When you use a framework
-, you reuse the main body and write the code it calls.'' In the context of
- frameworks and object-oriented programming, Pyre can be thought of as a
- collection of classes and the way their instances interact.
- Programming applications based on Pyre will look similar to those written
- in any other object-oriented language.
- The Pyre framework contains a subset of parts that make up the overall
- framework.
- Each of those parts is designed to solve a specific problem.
-\end_layout
-
-\begin_layout Standard
-The framework approach to computation offers many advantages.
- It permits the exchange of codes and promotes the reuse of standardized
- software while preserving efficiency.
- Frameworks are also an efficient way to handle changes in computer architecture.
- They present programmers and scientists with a unified and well-defined
- task and allow for shared costs of the housekeeping aspects of software
- development.
- They provide greater institutional continuity to model development than
- piecemeal approaches.
-\end_layout
-
-\begin_layout Standard
-The Pyre framework incorporates features aimed at enabling the scientific
- non-expert to perform tasks easily without hindering the expert.
- Target features for end users allow complete and intuitive simulation specifica
-tion, reasonable defaults, consistency checks of input, good diagnostics,
- easy access to remote facilities, and status monitoring.
- Target features for developers include easy access to user input, a shorter
- development cycle, and good debugging support.
-\end_layout
-
-\begin_layout Subsection
-Sieve and PETSc
-\end_layout
-
-\begin_layout Standard
-PyLith 1.x makes use of a set of data structures and routines in PETSc called
-
-\family typewriter
-Sieve
-\family default
-, which is still under active development.
-
-\family typewriter
-Sieve
-\family default
- provides data structures and routines for for representing and manipulating
- computational meshes, and it greatly simplifies finite-element computations.
-
-\family typewriter
- Sieve
-\family default
- represents the topology of the domain.
- Zero volume elements are inserted along all fault surfaces to implement
- kinematic (prescribed) or dynamic (constitutive model) implementations
- of fault slip.
- Material properties and other parameters are represented as sections (scalar
- and vector fields) over the mesh, and values for a vertex or cell can be
- retrieved by restricting the section to the vertex or cell.
- For each problem, functions are provided to calculate the residual and
- its Jacobian.
- All numerical integration is done in these functions, and parallel assembly
- is accomplished using the restrict/update paradigm of the
-\family typewriter
-Sieve
-\family default
- framework.
- We assemble into PETSc linear algebra objects and then call PETSc solvers.
- The solution is mapped back into a section, which can be output in VTK
- format.
-\end_layout
-
-\begin_layout Standard
-PETSc
-\begin_inset Flex URL
-status collapsed
-
-\begin_layout Plain Layout
-
-www-unix.mcs.anl.gov/petsc/petsc-as
-\end_layout
-
-\end_inset
-
-, the Portable, Extensible Toolkit for Scientific computation, provides
- a suite of routines for parallel, numerical solution of partial differential
- equations for linear and nonlinear systems with large, sparse systems of
- equations.
- PETSc includes solvers that implement a variety of Newton and Krylov subspace
- methods.
- It can also interface with many external packages, including ESSL, MUMPS,
- Matlab, ParMETIS, PVODE, and Hypre, thereby providing additional solvers
- and interaction with other software packages.
-\end_layout
-
-\begin_layout Standard
-PETSc includes interfaces for FORTRAN 77/90, C, C++, and Python for nearly
- all of the routines, and PETSc can be installed on most Unix systems.
- PETSc can be built with user-supplied, highly optimized linear algebra
- routines (e.g., ATLAS and commercial versions of BLAS/LAPACK), thereby improving
- application performance.
- Users can use PETSc parallel matrices, vectors, and other data structures
- for most parallel operations, eliminating the need for explicit calls to
- Message Passing Interface (MPI) routines.
- Many settings and options can be controlled with PETSc-specific command-line
- arguments, including selection of preconditions, solvers, and generation
- of performance logs.
-\end_layout
-
-\end_body
-\end_document
+#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 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
+\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
+Introduction
+\end_layout
+
+\begin_layout Section
+Overview
+\end_layout
+
+\begin_layout Standard
+PyLith is a multi-scale simulation software package for earthquake physics.
+ It is portable, scalable software for simulation of crustal deformation
+ across spatial scales ranging from meters to hundreds of kilometers and
+ temporal scales ranging from milliseconds to thousands of years.
+\end_layout
+
+\begin_layout Section
+History
+\end_layout
+
+\begin_layout Standard
+PyLith 1.0 was the first version to allow the solution of both implicit (quasi-st
+atic) and explicit (dynamic) problems and was a complete rewrite of the
+ original PyLith (version 0.8).
+ PyLith 1.0 combines the functionality of EqSim
+\begin_inset CommandInset citation
+LatexCommand cite
+key "Aagaard:etal:2001a,Aagaard:etal:2001b"
+
+\end_inset
+
+ and PyLith 0.8.
+ PyLith 0.8 was a direct descendant of LithoMop and was the first version
+ that ran in parallel, as well as providing several other improvements over
+ LithoMop.
+ LithoMop was the product of major reengineering of Tecton, a finite-element
+ code for simulating static and quasi-static crustal deformation.
+ The major new features present in LithoMop included dynamic memory allocation
+ and the use of the Pyre simulation framework and PETSc solvers.
+ EqSim was written by Brad Aagaard to solve problems in earthquake dynamics,
+ including rupture propagation and seismic wave propagation.
+\end_layout
+
+\begin_layout Standard
+The release of PyLith 1.0 has been followed by additional releases that expand
+ the number of features as well as improve performance.
+ The PyLith 1.x series of releases allows the solution of both quasi-static
+ and dynamic problems in one, two, or three dimensions.
+ The code runs in either serial or parallel, and the design allows for relativel
+y easy scripting using the Python programming language.
+ Material properties and values for boundary and fault conditions are specified
+ using spatial databases, which permit easy prescription of complex spatial
+ variations of properties and parameters.
+ Simulation parameters are generally specified through the use of simple
+ ASCII files or the command line.
+ At present, mesh information may be provided using a simple ASCII file
+ (PyLith mesh ASCII format) or imported from CUBIT or LaGriT, two widely-used
+ meshing packages.
+ The elements currently available include a linear bar in 1D, linear triangles
+ and quadrilaterals in 2D, and linear tetrahedra and hexahedra in 3D.
+ Materials presently available include isotropic elastic, linear Maxwell
+ viscoelastic, generalized Maxwell viscoelastic, power-law viscoelastic,
+ and Drucker-Prager elastoplastic.
+ Boundary conditions include Dirichlet (prescribed displacements and velocities)
+, Neumann (traction), point forces, and absorbing boundaries.
+ Cohesive elements are used to implement slip across interior surfaces (faults)
+ with both kinematically-specified fault slip and slip governed by fault
+ constitutive models.
+ PyLith also includes an interface for computing static Green's functions
+ for fault slip.
+\end_layout
+
+\begin_layout Standard
+PyLith is under active development and we expect a number of additions and
+ improvements in the near future.
+ Likely enhancements will include additional bulk and fault constitutive
+ models, coupled quasi-static and dynamic simulations for earthquake cycle
+ modeling, and coupling between elasticity, heat flow, and/or fluid flow.
+\end_layout
+
+\begin_layout Section
+PyLith Workflow
+\end_layout
+
+\begin_layout Standard
+PyLith is one component in the process of investigating problems in tectonics
+ (Figure
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "fig:Workflow-summary"
+
+\end_inset
+
+).
+ Given a geological problem of interest, a scientist must first provide
+ a geometrical representation of the desired structure.
+ Once the structure has been defined, a computational mesh must be created.
+ PyLith presently provides three mesh importing options: CUBIT Exodus format,
+ LaGriT GMV and Pset files, and PyLith mesh ASCII format.
+ The modeling of the physical processes of interest is performed by a code
+ such as PyLith.
+ Present output consists of VTK or HDF5/Xdmf files which can be used by
+ a number of visualization codes (e.g., ParaView, Visit, MayaVi, and Matlab).
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Plain Layout
+\align center
+\begin_inset Graphics
+ filename figs/workflow.eps
+ scale 67
+ keepAspectRatio
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Caption
+
+\begin_layout Plain Layout
+\begin_inset CommandInset label
+LatexCommand label
+name "fig:Workflow-summary"
+
+\end_inset
+
+Workflow involved in going from geologic structure to problem analysis.
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Section
+PyLith Design
+\end_layout
+
+\begin_layout Standard
+PyLith is separated into modules to encapsulate behavior and facilitate
+ use across multiple applications.
+ This allows expert users to replace functionality of a wide variety of
+ components without recompiling or polluting the main code.
+ PyLith employs external packages (see Figure
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "fig:pylith-dependencies"
+
+\end_inset
+
+) to reduce development time and enhance computational efficiency; for example,
+ PyLith 0.8 ran two times faster when the PETSc linear solver was used.
+\end_layout
+
+\begin_layout Standard
+\noindent
+\align center
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Plain Layout
+\align center
+\begin_inset Graphics
+ filename figs/packages.eps
+ scale 40
+
+\end_inset
+
+
+\begin_inset Caption
+
+\begin_layout Plain Layout
+\begin_inset CommandInset label
+LatexCommand label
+name "fig:pylith-dependencies"
+
+\end_inset
+
+PyLith dependencies.
+ PyLith makes direct use of several other packages, some of which have their
+ own dependencies.
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+PyLith is written in two programming languages.
+ High-level code is written in Python; this rich, expressive interpreted
+ language with dynamic typing reduces development time and permits flexible
+ addition of user-contributed modules.
+ This high-level code makes use of Pyre, a science-neutral simulation framework
+ developed at Caltech, to link the modules together at runtime and gather
+ user-input.
+ Low-level code is written in C++, providing fast execution while still
+ allowing an object-oriented implementation.
+ This low-level code relies on PETSc to perform operations on matrices and
+ vectors in parallel.
+ We also make extensive use of two Python packages.
+ SWIG is a package that simplifies the task of adding C++ extensions to
+ Python code, and FIAT provides tabulated basis functions and numerical
+ quadrature points.
+
+\end_layout
+
+\begin_layout Standard
+In writing PyLith 1.0, the code was designed to be object-oriented and modular.
+ Each type of module is accessed through a specified interface (set of functions
+).
+ This permits adding, replacing, and rewriting modules without affecting
+ other parts of the code.
+ This code structure simplifies code maintenance and development.
+ Extending the set of code features is also easier, since developers can
+ create new modules derived from the existing ones.
+\end_layout
+
+\begin_layout Standard
+The new code design leverages Pyre, Sieve, and PETSc much more extensively
+ than the previous version.
+ Pyre is used to glue together the various modules used to construct a
+ simulation and specify the parameters.
+ Sieve is used for all finite-element storage and manipulation and handles
+ the creation of the PETSc matrices and vectors.
+ As a result, most of the PyLith source code pertains to implementing the
+ geodynamics, such as bulk rheology, boundary conditions, and slip on faults.
+
+\end_layout
+
+\begin_layout Standard
+PyLith also uses FIAT to tabulate the finite-element basis functions at
+ the numerical integration (quadrature) points.
+ Nemesis allows PyLith to run Python using the Message Passing Interface
+ (MPI) for parallel processing.
+ Additional, indirect dependencies (see Figure
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "fig:pylith-dependencies"
+
+\end_inset
+
+) include numpy (efficient operations on numerical arrays in Python), Proj.4
+ (geographic projections), and SWIG (calling C++ functions from Python).
+\end_layout
+
+\begin_layout Standard
+During development, tests were constructed for nearly every module function.
+ These unit tests are distributed with the source code.
+ These tests are run throughout the development cycle to expose bugs and
+ isolate their origin.
+ As additional changes are made to the code, the tests are rerun to help
+ prevent introduction of new bugs.
+ A number of simple, full-scale tests, such as axial compression and extension,
+ simple shear, and slip on through-going faults, have been used to test
+ the code.
+ Additionally, we have run the Southern California Earthquake Center crustal
+ deformation and several of the spontaneous rupture benchmarks for strike-slip
+ and reverse-slip to determine the relative local and global error (see
+ Chapter
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "cha:Benchmarks"
+
+\end_inset
+
+).
+\end_layout
+
+\begin_layout Subsection
+Pyre
+\end_layout
+
+\begin_layout Standard
+Pyre is an object-oriented environment capable of specifying and launching
+ numerical simulations on multiple platforms, including Beowulf-class parallel
+ computers and grid computing systems.
+ Pyre allows the binding of multiple components such as solid and fluid
+ models used in Earth science simulations, and different meshers.
+ The Pyre framework enables the elegant setup, modification and launching
+ of massively parallel solver applications.
+\end_layout
+
+\begin_layout Standard
+\noindent
+\align center
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Plain Layout
+\align center
+\begin_inset Graphics
+ filename figs/pyre_overview.png
+ width 4in
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Caption
+
+\begin_layout Plain Layout
+\begin_inset CommandInset label
+LatexCommand label
+name "fig:Pyre:Architecture"
+
+\end_inset
+
+Pyre Architecture.
+ The integration framework is a set of cooperating abstract services.
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Pyre is a framework, a combination of software and design philosophy that
+ promotes the reuse of code.
+ In their canonical software design book,
+\emph on
+Design Patterns
+\emph default
+, Erich Gamma
+\shape italic
+et al
+\shape default
+.
+ condense the concept of a framework concept down to, ``When you use a framework
+, you reuse the main body and write the code it calls.'' In the context of
+ frameworks and object-oriented programming, Pyre can be thought of as a
+ collection of classes and the way their instances interact.
+ Programming applications based on Pyre will look similar to those written
+ in any other object-oriented language.
+ The Pyre framework contains a subset of parts that make up the overall
+ framework.
+ Each of those parts is designed to solve a specific problem.
+\end_layout
+
+\begin_layout Standard
+The framework approach to computation offers many advantages.
+ It permits the exchange of codes and promotes the reuse of standardized
+ software while preserving efficiency.
+ Frameworks are also an efficient way to handle changes in computer architecture.
+ They present programmers and scientists with a unified and well-defined
+ task and allow for shared costs of the housekeeping aspects of software
+ development.
+ They provide greater institutional continuity to model development than
+ piecemeal approaches.
+\end_layout
+
+\begin_layout Standard
+The Pyre framework incorporates features aimed at enabling the scientific
+ non-expert to perform tasks easily without hindering the expert.
+ Target features for end users allow complete and intuitive simulation specifica
+tion, reasonable defaults, consistency checks of input, good diagnostics,
+ easy access to remote facilities, and status monitoring.
+ Target features for developers include easy access to user input, a shorter
+ development cycle, and good debugging support.
+\end_layout
+
+\begin_layout Subsection
+Sieve and PETSc
+\end_layout
+
+\begin_layout Standard
+PyLith 1.x makes use of a set of data structures and routines in PETSc called
+
+\family typewriter
+Sieve
+\family default
+, which is still under active development.
+
+\family typewriter
+Sieve
+\family default
+ provides data structures and routines for for representing and manipulating
+ computational meshes, and it greatly simplifies finite-element computations.
+
+\family typewriter
+ Sieve
+\family default
+ represents the topology of the domain.
+ Zero volume elements are inserted along all fault surfaces to implement
+ kinematic (prescribed) or dynamic (constitutive model) implementations
+ of fault slip.
+ Material properties and other parameters are represented as sections (scalar
+ and vector fields) over the mesh, and values for a vertex or cell can be
+ retrieved by restricting the section to the vertex or cell.
+ For each problem, functions are provided to calculate the residual and
+ its Jacobian.
+ All numerical integration is done in these functions, and parallel assembly
+ is accomplished using the restrict/update paradigm of the
+\family typewriter
+Sieve
+\family default
+ framework.
+ We assemble into PETSc linear algebra objects and then call PETSc solvers.
+ The solution is mapped back into a section, which can be output in VTK
+ format.
+\end_layout
+
+\begin_layout Standard
+PETSc
+\begin_inset Flex URL
+status collapsed
+
+\begin_layout Plain Layout
+
+www-unix.mcs.anl.gov/petsc/petsc-as
+\end_layout
+
+\end_inset
+
+, the Portable, Extensible Toolkit for Scientific computation, provides
+ a suite of routines for parallel, numerical solution of partial differential
+ equations for linear and nonlinear systems with large, sparse systems of
+ equations.
+ PETSc includes solvers that implement a variety of Newton and Krylov subspace
+ methods.
+ It can also interface with many external packages, including ESSL, MUMPS,
+ Matlab, ParMETIS, PVODE, and Hypre, thereby providing additional solvers
+ and interaction with other software packages.
+\end_layout
+
+\begin_layout Standard
+PETSc includes interfaces for FORTRAN 77/90, C, C++, and Python for nearly
+ all of the routines, and PETSc can be installed on most Unix systems.
+ PETSc can be built with user-supplied, highly optimized linear algebra
+ routines (e.g., ATLAS and commercial versions of BLAS/LAPACK), thereby improving
+ application performance.
+ Users can use PETSc parallel matrices, vectors, and other data structures
+ for most parallel operations, eliminating the need for explicit calls to
+ Message Passing Interface (MPI) routines.
+ Many settings and options can be controlled with PETSc-specific command-line
+ arguments, including selection of preconditions, solvers, and generation
+ of performance logs.
+\end_layout
+
+\end_body
+\end_document
Modified: short/3D/PyLith/branches/v1.7-stable/doc/userguide/runpylith/runpylith.lyx
===================================================================
--- short/3D/PyLith/branches/v1.7-stable/doc/userguide/runpylith/runpylith.lyx 2012-06-26 16:08:18 UTC (rev 20413)
+++ short/3D/PyLith/branches/v1.7-stable/doc/userguide/runpylith/runpylith.lyx 2012-06-26 16:08:27 UTC (rev 20414)
@@ -79,7 +79,7 @@
A set of parameters describing the problem.
These parameters describe the type of problem to be run, solver information,
time-stepping information, boundary conditions, materials, etc.
- This information can be provided from the command-line or by using a
+ This information can be provided from the command line or by using a
\family typewriter
.cfg
\family default
@@ -704,7 +704,7 @@
\begin_layout Standard
All of the example problems are set up using configuration files in the
example directory, and specific problems are defined by including the appropria
-te configuration file on the command-line.
+te configuration file on the command line.
Referring to the directory
\family typewriter
examples/twocells/twohex8
@@ -729,8 +729,12 @@
\end_layout
\begin_layout Standard
-The settings in pylithapp.cfg will be read automatically, and additional
- settings are included by specifying one of the other files on the command-line:
+The settings in
+\family typewriter
+pylithapp.cfg
+\family default
+ will be read automatically, and additional settings are included by specifying
+ one of the other files on the command line:
\end_layout
\begin_layout LyX-Code
@@ -1184,7 +1188,7 @@
\begin_layout Standard
The MeshIOCubit object reads the NetCDF Exodus II files output from CUBIT.
Beginning with CUBIT 11.0, the names of the nodesets are included in the
- Exodus II files and PyLith can use these nodeset names or revert to using
+ Exodus II files, and PyLith can use these nodeset names or revert to using
the nodeset ids.
The properties and facilities associated with the MeshIOCubit object are:
\end_layout
@@ -1222,7 +1226,7 @@
\end_layout
\begin_layout Description
-filename_gmv Name of GMV file.
+filename_gmv Name of the GMV file.
\end_layout
\begin_layout Description
@@ -1241,7 +1245,7 @@
\begin_layout Description
record_header_32bt Flag indicating FORTRAN record header is 32-bit (default
- is True)
+ is True).
\end_layout
\begin_layout Description
@@ -1253,7 +1257,7 @@
\end_layout
\begin_layout Standard
-The distributor users a partitioner to compute which cells should be placed
+The distributor uses a partitioner to compute which cells should be placed
on each processor, computes the overlap among the processors, and then
distributes the mesh among the processors.
The properties and facilities of the Distributor include:
@@ -1276,7 +1280,7 @@
\begin_inset Quotes erd
\end_inset
-, default is
+; default is
\begin_inset Quotes eld
\end_inset
@@ -1380,10 +1384,10 @@
This allows one to run much larger simulations by (1) permitting the mesh
generator to construct a mesh with a node spacing twice as large as that
needed in the simulation and (2) operations performed in serial during
- the simulation setup phase, such as, adjusting the topology to insert cohesive
- cells and distribution of the mesh among processors uses this much smaller
+ the simulation setup phase, such as adjusting the topology to insert cohesive
+ cells and distribution of the mesh among processors, uses this much smaller
coarse mesh.
- For 2D problems the global mesh refinement increases the maximum problem
+ For 2D problems, the global mesh refinement increases the maximum problem
size by a factor of four, and for 3D problems it increases the maximum
problem size by a factor of eight.
\end_layout
@@ -1404,7 +1408,7 @@
\family typewriter
TimeDependent
\family default
- for use in static, quasi-static, and dynamic simulations and
+ for use in static, quasi-static, and dynamic simulations, and
\family typewriter
GreensFns
\family default
@@ -1442,7 +1446,7 @@
\end_layout
\begin_layout Description
-dimension Spatial dimension of the problem (default is 3)
+dimension Spatial dimension of the problem (default is 3).
\end_layout
\begin_layout Standard
@@ -1611,8 +1615,9 @@
specified using the FIATLagrange object, whereas the parameters for Simplex
cells (lines, triangles, tetrahedra) are specified using the FIATSimplex
object.
- Both objects use the same set of parameters and PyLith will setup the basis
- functions and quadrature scheme appropriately for the two families of cells.
+ Both objects use the same set of parameters, and PyLith will set up the
+ basis functions and quadrature scheme appropriately for the two families
+ of cells.
The quadrature scheme and basis functions must be set for each material
and boundary condition involving finite-element integrations (Dirichlet
boundary conditions are constraints and do not involve integrations).
@@ -1730,7 +1735,7 @@
).
A more advanced set of solver settings that may provide better performance
- in many elasticity problems are given in Table
+ in many elasticity problems is given in Table
\begin_inset CommandInset ref
LatexCommand ref
reference "tab:petsc:options:advanced"
@@ -1768,7 +1773,7 @@
\series default
\color none
-These settings are only available if you build PETSc with Fortran enabled
+These settings are only available if you build PETSc with FORTRAN enabled
and the ML package.
These features are included in the PyLith binary packages.
\end_layout
@@ -2309,6 +2314,8 @@
\begin_inset Text
\begin_layout Plain Layout
+
+\family typewriter
ksp_type
\end_layout
@@ -2342,6 +2349,7 @@
\begin_layout Plain Layout
\family typewriter
+\size small
sub_pc_factor_shift_type
\end_layout
@@ -2486,7 +2494,6 @@
\begin_layout Plain Layout
\shape italic
-\emph on
1.0e-12
\end_layout
@@ -2789,6 +2796,7 @@
\begin_layout Plain Layout
\family typewriter
+\size small
fs_pc_fieldsplit_real_diagonal
\end_layout
@@ -3280,27 +3288,39 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp.timedependent.formulation]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
time_step = pylith.problems.TimeStepUniform
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
solver = pylith.problems.SolverLinear ; Nonlinear solver is pylith.problems.SolverNo
nlinear
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
output = [domain,ground_surface]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
matrix_type = sbaij ; To use a non-symmetric sparse matrix, set it to aij
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
view_jacobian = false
\end_layout
@@ -3312,7 +3332,7 @@
In explicit time-stepping formulations for elasticity, boundary conditions
and fault slip can excite short waveform elastic waves that are not accurately
resolved by the discretization.
- We use numerical damping via an artificial viscosity
+ We use numerical damping via an artificial viscosity
\begin_inset CommandInset citation
LatexCommand cite
key "Knopoff:Ni:2001,Day:Ely:2002"
@@ -3470,7 +3490,7 @@
\end_layout
\begin_layout Description
-start_time Start time for simulation (default is 0.0 s)
+start_time Start time for simulation (default is 0.0 s).
\end_layout
\begin_layout Description
@@ -3514,7 +3534,7 @@
\end_layout
\begin_layout Standard
-The nonuniform, user-specified, time-step implementation allows the user
+The nonuniform, user-specified time-step implementation allows the user
to specify the time steps in an ASCII file (see Section
\begin_inset CommandInset ref
LatexCommand ref
@@ -3538,7 +3558,7 @@
\end_layout
\begin_layout Description
-loop_steps If true, cycle through time steps, otherwise keep using last
+loop_steps If true, cycle through time steps; otherwise keep using last
time-step size for any time remaining.
\end_layout
@@ -3552,10 +3572,14 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp.problem.formulation]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
time_step = pylith.problems.TimeStepUser ; Change the time step algorithm
\end_layout
@@ -3564,18 +3588,26 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp.problem.formulation.time_step]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
total_time = 1000.0*year
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
filename = timesteps.txt
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
loop_steps = false ; Default value
\end_layout
@@ -3622,10 +3654,14 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp.problem.formulation]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
time_step = pylith.problems.TimeStepAdapt ; Change the time step algorithm
\end_layout
@@ -3634,22 +3670,32 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp.problem.formulation.time_step]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
total_time = 1000.0*year
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
max_dt = 10.0*year
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
adapt_skip = 10 ; Default value
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
stability_factor = 2.0 ; Default value
\end_layout
@@ -3715,10 +3761,14 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
problem = pylith.problems.GreensFns ; Change problem type from the default
\end_layout
@@ -3727,10 +3777,14 @@
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
[pylithapp.greensfns]
\end_layout
\begin_layout LyX-Code
+
+\size footnotesize
fault_id = 100 ; Default value
\end_layout
@@ -4069,26 +4123,26 @@
\begin_layout LyX-Code
-\size small
+\size footnotesize
[pylithapp.timedependent.materials.material]
\end_layout
\begin_layout LyX-Code
-\size small
+\size footnotesize
db_properties = spatialdata.spatialdb.UniformDB ; Set the db to a UniformDB
\end_layout
\begin_layout LyX-Code
-\size small
+\size footnotesize
db_properties.values = [vp,vs,density] ; Set the names of the values in the
database
\end_layout
\begin_layout LyX-Code
-\size small
+\size footnotesize
db_properties.data = [5773.5*m/s, 3333.3*m/s, 2700.0*kg/m**3] ; Set the values
in the database
\end_layout
@@ -4410,7 +4464,7 @@
\begin_layout Standard
As in the other Pyre objects, spatial database objects contain parameters
- that can be set from the command line or using
+ that can be set from the command line or by using
\family typewriter
.cfg or .pml
\family default
@@ -4706,7 +4760,7 @@
OutputSolnPoints
\family default
, to interpolate the solution to arbitrary points.
- By default, the output manager will include the displaceent time histories
+ By default, the output manager will include the displacement time histories
in the output.
The locations are specified in a text file.
In addition to the
More information about the CIG-COMMITS
mailing list