[cig-commits] r12498 - doc/cigma/manual

[luis at geodynamics.org]

Author: luis
Date: 2008-07-31 10:58:23 -0700 (Thu, 31 Jul 2008)
New Revision: 12498

Added:
   doc/cigma/manual/cigma.lyx~
Modified:
   doc/cigma/manual/cigma.lyx
Log:
Manual update

Modified: doc/cigma/manual/cigma.lyx
===================================================================
--- doc/cigma/manual/cigma.lyx	2008-07-30 14:49:37 UTC (rev 12497)
+++ doc/cigma/manual/cigma.lyx	2008-07-31 17:58:23 UTC (rev 12498)
@@ -1,2505 +1,2509 @@
-#LyX 1.5.1 created this file. For more info see http://www.lyx.org/
-\lyxformat 276
-\begin_document
-\begin_header
-\textclass book
-\begin_preamble
-\usepackage{hyperref}
-
-\let\myUrl\url
-\renewcommand{\url}[1]{(\myUrl{#1})}
-\end_preamble
-\language english
-\inputencoding auto
-\font_roman default
-\font_sans default
-\font_typewriter default
-\font_default_family default
-\font_sc false
-\font_osf false
-\font_sf_scale 100
-\font_tt_scale 100
-\graphics default
-\paperfontsize default
-\spacing single
-\papersize default
-\use_geometry true
-\use_amsmath 1
-\use_esint 0
-\cite_engine basic
-\use_bibtopic false
-\paperorientation portrait
-\leftmargin 1in
-\topmargin 1in
-\rightmargin 1in
-\bottommargin 1in
-\secnumdepth 3
-\tocdepth 3
-\paragraph_separation indent
-\defskip medskip
-\quotes_language swedish
-\papercolumns 1
-\papersides 2
-\paperpagestyle headings
-\tracking_changes false
-\output_changes false
-\author "" 
-\author "" 
-\end_header
-
-\begin_body
-
-\begin_layout Standard
-\noindent
-\align center
-
-\color none
-\begin_inset ERT
-status collapsed
-
-\begin_layout Standard
-
-
-\backslash
-thispagestyle{empty}
-\end_layout
-
-\end_inset
-
-
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status collapsed
-
-\begin_layout Standard
-\begin_inset Graphics
-	filename figures/cigma_cover.pdf
-	display color
-	width 75page%
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\noindent
-\align center
-
-\color none
-\begin_inset ERT
-status collapsed
-
-\begin_layout Standard
-
-
-\backslash
-thispagestyle{empty}
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Title
-Cigma
-\end_layout
-
-\begin_layout Author
-© California Institute of Technology
-\newline
-Version 0.9.1
-\end_layout
-
-\begin_layout Date
-\begin_inset ERT
-status collapsed
-
-\begin_layout Standard
-
-
-\backslash
-today
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset LatexCommand tableofcontents
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset FloatList figure
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset ERT
-status collapsed
-
-\begin_layout Standard
-
-
-\backslash
-raggedbottom
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-
-\newpage
-
-\end_layout
-
-\begin_layout Chapter
-Introduction
-\end_layout
-
-\begin_layout Section
-About Cigma
-\end_layout
-
-\begin_layout Standard
-The CIG Model Analyzer (Cigma) consists of a suite of tools for comparing
- numerical models.
- CIG has developed Cigma in response to demand from the short-term tectonics
- community for a simple tool that can perform rigorous error analysis on
- their finite element codes.
- The long-term goal for Cigma, however, is for it to be used for nearly
- all geodynamics modeling codes.
-\end_layout
-
-\begin_layout Standard
-The current version of Cigma is intended for the calculation of 
-\begin_inset Formula $L_{2}$
-\end_inset
-
- residuals for finite element models.
- This software was designed for performing the following three main tasks:
- (1) error analysis, (2) benchmarking, and (3) code verification.
- In error analysis, Cigma can calculate a global measure of the error between
- two solutions by performing an integration over a discretized version of
- the common domain.
- This comparison can take place even when the underlying discretizations
- do not overlap.
- In benchmarking, Cigma can help the geodynamics community agree on a standard
- solution to specific problems by facilitating the process of comparing
- different numerical codes against each other.
- Lastly, as an automated tool, Cigma can help application developers create
- regression tests to ensure that software changes do not affect the consistency
- of the results.
-\end_layout
-
-\begin_layout Standard
-At its core, Cigma draws from a variety of libraries.
- The 
-\begin_inset LatexCommand htmlurl
-name "FInite element Automatic Tabulator (FIAT)"
-target "www.fenics.org/wiki/FIAT"
-
-\end_inset
-
- Python Library supports generation of arbitrary order instances of the
- Lagrange elements on lines, triangles, and tetrahedra.
- It can also generate higher order instances of Jacobi-type quadrature rules
- on the same element shapes.
- The 
-\begin_inset LatexCommand htmlurl
-name "Approximate Nearest Neighbor (ANN)"
-target "www.cs.umd.edu/~mount/ANN/"
-
-\end_inset
-
- Library, written in C++, supports data structures and algorithms for both
- exact and nearest-neighbor searching in arbitrarily high dimensions.
-\end_layout
-
-\begin_layout Standard
-Both of these libraries extend and generalize Cigma's functionality so it
- can handle other types of elements, and provide the ability to compare
- vector fields.
-\end_layout
-
-\begin_layout Section
-Citation
-\end_layout
-
-\begin_layout Standard
-Computational Infrastructure for Geodynamics (CIG) is making this source
- code available to you in the hope that the software will enhance your research
- in geophysics.
- This is a brand-new code and at present no papers are published or press.
- Please cite this manual as follows:
-\end_layout
-
-\begin_layout Itemize
-Armendariz, L., and S.
- Kientz.
- 
-\emph on
-Cigma User Manual.
-
-\emph default
- Pasadena, CA: Computational Infrastructure of Geodynamics, 2008.
- URL: geodynamics.org/cig/software/cs/cigma/cigma.pdf
-\end_layout
-
-\begin_layout Standard
-CIG requests that in your oral presentations and in your papers that you
- indicate your use of this code and acknowledge the author of the code and
- 
-\begin_inset LatexCommand htmlurl
-name "CIG"
-target "geodynamics.org"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Section
-Support
-\end_layout
-
-\begin_layout Standard
-Cigma development is based upon work supported by the National Science Foundatio
-n under Grant No.
- EAR-0406751.
- Any opinions, findings, and conclusions or recommendations expressed in
- this material are those of the authors and do not necessarily reflect the
- views of the National Science Foundation.
- The code is being released under the GNU General Public License.
-\end_layout
-
-\begin_layout Chapter
-Installation and Getting Help
-\end_layout
-
-\begin_layout Section
-Getting Help
-\end_layout
-
-\begin_layout Standard
-For help, send an e-mail to the 
-\begin_inset LatexCommand url
-name "CIG Computational Science Mailing List"
-target "cig-cs at geodynamics.org"
-
-\end_inset
-
-.
- You can subscribe to the 
-\family typewriter
-cig-cs
-\family default
- mailing list and view archived discussions at the 
-\begin_inset LatexCommand htmlurl
-name "CIG Mail Lists web page"
-target "geodynamics.org/cig/lists"
-
-\end_inset
-
-.
- If you encounter any bugs or have problems installing Cigma, please submit
- a report to the 
-\begin_inset LatexCommand htmlurl
-name "CIG Bug Tracker"
-target "geodynamics.org/bugs"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Section
-Installating from Source
-\end_layout
-
-\begin_layout Standard
-To install Cigma, download the source package from the 
-\begin_inset LatexCommand htmlurl
-name "CIG Cigma web page"
-target "geodynamics.org/cig/software/packages/cs/cigma"
-
-\end_inset
-
-.
- You will need the GNU C++ compiler for this step.
- The HDF5 and VTK libraries, described later in this section, are also required.
-\end_layout
-
-\begin_layout Standard
-After installing the necessary dependencies, simply unpack the source tar
- file and issue the following commands
-\end_layout
-
-\begin_layout LyX-Code
-$ tar xvfz cigma-0.9.1.tar.gz
-\end_layout
-
-\begin_layout LyX-Code
-$ cd cigma-0.9.1
-\end_layout
-
-\begin_layout LyX-Code
-$ ./configure --with-hdf5=$HDF5_PREFIX --with-vtk=$VTK_PREFIX
-\end_layout
-
-\begin_layout LyX-Code
-$ make
-\end_layout
-
-\begin_layout LyX-Code
-$ make install
-\end_layout
-
-\begin_layout Standard
-Both HDF5 and VTK prefixes are required if you have installed these libraries
- in a custom location.
-\end_layout
-
-\begin_layout Subsection
-HDF5
-\end_layout
-
-\begin_layout Standard
-The Hierarchical Data Format (HDF5) library is available for download from
- 
-\begin_inset LatexCommand htmlurl
-name "The HDF Group"
-target "hdfgroup.org/HDF5"
-
-\end_inset
-
-.
- Binaries can be obtained at 
-\begin_inset LatexCommand htmlurl
-name "hdfgroup.org/HDF5/release/obtain5.html"
-target "hdfgroup.org/HDF5/release/obtain5.html"
-
-\end_inset
-
-.
- To install from source, download the latest stable 1.6 version of this library
- (currently 1.6.7) and issue the following commands
-\end_layout
-
-\begin_layout LyX-Code
-$ tar xvfz hdf5-1.6.7
-\end_layout
-
-\begin_layout LyX-Code
-$ cd hdf5-1.6.7
-\end_layout
-
-\begin_layout LyX-Code
-$ ./configure --prefix=$HOME/opt/hdf5/1.6.7
-\end_layout
-
-\begin_layout LyX-Code
-$ make
-\end_layout
-
-\begin_layout LyX-Code
-$ make install
-\end_layout
-
-\begin_layout Standard
-Without the custom prefix, this installation procedure will target your
- 
-\family typewriter
-/usr
-\family default
- directory.
-\end_layout
-
-\begin_layout Subsection
-VTK
-\end_layout
-
-\begin_layout Standard
-The Visualization Toolkit (VTK) library is a popular open source graphics
- library for scientific visualizations.
- We recommend that you obtain the binaries from a package manager, making
- sure to install the associated development headers along with the library.
- The source for the VTK library is available from 
-\begin_inset LatexCommand htmlurl
-name "Kitware, Inc."
-target "www.vtk.org/get-software.php"
-
-\end_inset
-
-.
- If you decide to compile VTK from source, you will also need the CMake
- build environment, available from 
-\begin_inset LatexCommand htmlurl
-name "CMake"
-target "cmake.org"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Standard
-After downloading the source package, you can issue the following commands
-\end_layout
-
-\begin_layout LyX-Code
-$ tar xvfz vtk-5.0.4.tar.gz
-\end_layout
-
-\begin_layout LyX-Code
-$ cd VTK
-\end_layout
-
-\begin_layout LyX-Code
-$ mkdir build
-\end_layout
-
-\begin_layout LyX-Code
-$ cd build
-\end_layout
-
-\begin_layout LyX-Code
-$ ccmake ..
-\end_layout
-
-\begin_layout LyX-Code
-$ make
-\end_layout
-
-\begin_layout LyX-Code
-$ make install
-\end_layout
-
-\begin_layout Subsection
-NumPy
-\end_layout
-
-\begin_layout Standard
-NumPy is a Python module which adds support for large multi-dimensional
- arrays and matrices, and is available for download from the 
-\begin_inset LatexCommand htmlurl
-name "NumPy Home Page"
-target "numpy.scipy.org "
-
-\end_inset
-
-.
- To install this extension module from source, download it and issue the
- following command in the NumPy source directory
-\end_layout
-
-\begin_layout LyX-Code
-$ python setup.py install
-\end_layout
-
-\begin_layout Standard
-To verify the installation, the command `
-\family typewriter
-import numpy
-\family default
-' should run successfully from within your standard Python interpreter shell.
-\end_layout
-
-\begin_layout Subsection
-PyTables
-\end_layout
-
-\begin_layout Standard
-PyTables is a Python extension module that builds on top of the HDF5 library.
- It provides a convenient scripting interface to manipulate HDF5 files,
- which can be used to manipulate the input/output files created by Cigma.
- PyTables is available for download from 
-\begin_inset LatexCommand htmlurl
-name "PyTables"
-target "www.pytables.org"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Standard
-To install this extension from source, download the latest stable version
- (currently 2.0.2) and issue the following commands
-\end_layout
-
-\begin_layout LyX-Code
-$ tar xvfz pytables-2.0.2
-\end_layout
-
-\begin_layout LyX-Code
-$ cd pytables-2.0.2
-\end_layout
-
-\begin_layout LyX-Code
-$ python setup.py install
-\end_layout
-
-\begin_layout Standard
-To verify the installation, the command `
-\family typewriter
-import tables
-\family default
-' should run successfully from within your standard Python interpreter shell.
-\end_layout
-
-\begin_layout Subsection
-FIAT
-\end_layout
-
-\begin_layout Standard
-The Finite Element Automatic Tabulator (FIAT) Python module supports generation
- of arbitrary order instances of the Lagrange elements, as well as arbitrary
- order instances of the Jacobi-type quadrature rules on the same element
- shapes.
- It is available from 
-\begin_inset LatexCommand htmlurl
-name "The FEniCS Project"
-target "www.fenics.org/wiki/FIAT"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Standard
-To install this module, download it and issue the following command inside
- the FIAT source directory:
-\end_layout
-
-\begin_layout LyX-Code
-$ sudo python setup.py install
-\end_layout
-
-\begin_layout Standard
-You should now be able to run `
-\family typewriter
-import FIAT
-\family default
-' from within your standard Python interpreter shell.
-\end_layout
-
-\begin_layout Subsection
-HDFView (optional)
-\end_layout
-
-\begin_layout Standard
-NCSA HDFView is a graphical user interface tool for accessing data in your
- HDF5 files.
- You can use it for viewing the internal group hierarchy (see Figure 
-\begin_inset LatexCommand ref
-reference "fig:With-HFDFView,-you"
-
-\end_inset
-
-), adding new datasets, and modifying or deleting existing datasets, or
- altering metadata on groups and datasets.
- You can download it from the 
-\begin_inset LatexCommand htmlurl
-name "HDFView home page"
-target "hdf.ncsa.uiuc.edu/hdf-java-html/hdfview"
-
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Standard
-\noindent
-\align center
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status collapsed
-
-\begin_layout Standard
-\noindent
-\align center
-\begin_inset Graphics
-	filename figures/hdfview-bmrsnog.png
-	lyxscale 70
-	width 60page%
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset Caption
-
-\begin_layout Standard
-\begin_inset LatexCommand label
-name "fig:With-HFDFView,-you"
-
-\end_inset
-
-With HFDFView, you can view the internal file hierarchy in a tree structure,
- add new datasets, and modify or delete existing datasets
-\end_layout
-
-\end_inset
-
- 
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Chapter
-\begin_inset LatexCommand label
-name "cha:Error-Analysis"
-
-\end_inset
-
-Error Analysis
-\end_layout
-
-\begin_layout Section
-Introduction
-\end_layout
-
-\begin_layout Standard
-When studying differential equations that represent physical systems we
- often obtain solutions by using a variety of techniques.
- Solving these equations using classical analytical methods is almost impossible
-, so we often resort to numerical algorithms such as the finite element
- method.
-\end_layout
-
-\begin_layout Standard
-For assessing the quality of the solution to our equations, it is important
- to quantify the error in our numerical solutions.
- To calculate this error, we will have to compare these solutions against
- each other through the use of a distance measure between two functions.
-\end_layout
-
-\begin_layout Section
-Distance Measures
-\end_layout
-
-\begin_layout Standard
-The simplest possible quantitative measure of the difference between two
- distinct functions consists of taking the pointwise difference at a common
- set of points.
- While no finite sample of points can perfectly represent a continuum of
- values, valuable information can be inferred from the statistics of the
- resulting set of differences.
-\end_layout
-
-\begin_layout Standard
-XXX: motivate this measure..validation
-\end_layout
-
-\begin_layout Standard
-A very useful distance measure we can use is the 
-\begin_inset Formula $L_{2}$
-\end_inset
-
- norm, defined by the following integral
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\varepsilon=||u-v||_{L_{2}}=\sqrt{\int_{\Omega}||u(\vec{x})-v(\vec{x})||^{2}d\vec{x}}\]
-
-\end_inset
-
-where 
-\begin_inset Formula $u$
-\end_inset
-
- and 
-\begin_inset Formula $v$
-\end_inset
-
- are the two functions on a global coordinate system.
- This gives us a single global estimate 
-\begin_inset Formula $\varepsilon$
-\end_inset
-
- representing the distance between the two fields 
-\begin_inset Formula $u(\vec{x})$
-\end_inset
-
- and 
-\begin_inset Formula $v(\vec{x})$
-\end_inset
-
-.
- Alternatively, you may think of this as the size, or norm, of the residual
- field 
-\begin_inset Formula $\rho(\vec{x})=u(\vec{x})-v(\vec{x})$
-\end_inset
-
-.
- 
-\end_layout
-
-\begin_layout Section
-Global Error Approximation
-\end_layout
-
-\begin_layout Standard
-If we discretize the domain 
-\begin_inset Formula $\Omega$
-\end_inset
-
- into an appropriate set of cells 
-\begin_inset Formula $\Omega_{1},\Omega_{2},\ldots,\Omega_{n_{el}}$
-\end_inset
-
-, we can compute the global error 
-\begin_inset Formula $\varepsilon^{2}$
-\end_inset
-
- as a sum over localized cell contributions, 
-\begin_inset Formula $\varepsilon^{2}=\sum_{e}\varepsilon_{e}^{2}$
-\end_inset
-
-, where each cell residual 
-\begin_inset Formula $\varepsilon_{e}^{2}$
-\end_inset
-
- is given by
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \begin{eqnarray*}
-\varepsilon_{e}^{2} & = & \int_{\Omega_{e}}||u(\vec{x})-v(\vec{x})||^{2}d\vec{x}\end{eqnarray*}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-In general, we won't be able to integrate each cell residual 
-\begin_inset Formula $\varepsilon_{e}^{2}$
-\end_inset
-
- exactly since either of the functions 
-\begin_inset Formula $u$
-\end_inset
-
- and 
-\begin_inset Formula $v$
-\end_inset
-
- may have an incompatible representation relative to the finite element
- space on each domain 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
-.
- However, we can still calculate an approximation of each cell residual
- by applying an appropriate quadrature rule with a tolerable truncation
- error 
-\begin_inset LatexCommand cite
-key "Encyclopaedia of Cubature Formulas 2005"
-
-\end_inset
-
-.
- 
-\end_layout
-
-\begin_layout Standard
-To obtain an approximation to the integral of a function 
-\begin_inset Formula $F(\vec{x})$
-\end_inset
-
- over a cell 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
-, we simply apply the quadrature rule with weights 
-\begin_inset Formula $w_{e,1},w_{e,2},\ldots,w_{e,n_{Q}}$
-\end_inset
-
- and integration points 
-\begin_inset Formula $\vec{x}_{e,1},\vec{x}_{e,2},\ldots,\vec{x}_{e,n_{Q}}$
-\end_inset
-
- appropriate for the physical element 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
-:
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\int_{\Omega_{e}}\ F(\vec{x})\ d\vec{x}=\sum_{q=1}^{n_{Q}}w_{e,q}F(\vec{x}_{e,q})\]
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-Applying this quadrature rule directly over the entire physical domain 
-\begin_inset Formula $\Omega=\cup\Omega_{e}$
-\end_inset
-
- gives us 
-\begin_inset Formula \[
-\int_{\Omega}\ F(\vec{x})\ d\vec{x}=\sum_{e=1}^{n_{el}}\sum_{q=1}^{n_{Q}}w_{e,q}F(\vec{x}_{e,q})\]
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-For efficiency reasons, it is undesirable in finite element applications
- to perform calculations in a global coordinate system.
- To avoid duplication of work, shape function evaluations may be performed
- once on a reference cell 
-\begin_inset Formula $\hat{\Omega}$
-\end_inset
-
- and then transformed back into the corresponding physical cell 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
- as needed.
-\end_layout
-
-\begin_layout Standard
-To compute integrals of 
-\begin_inset Formula $F$
-\end_inset
-
- in a reference coordinate system, we need to apply a change of variables:
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\int_{\Omega_{e}}\ F(\vec{x})\ d\vec{x}=\int_{\hat{\Omega}}\ F(\vec{x}_{e}(\vec{\xi}))J_{e}(\vec{\xi})\ d\vec{\xi}\]
-
-\end_inset
-
-where 
-\color none
-the additional factor 
-\begin_inset Formula $J_{e}(\vec{\xi})=\det\left[\frac{d\vec{x}_{e}}{d\vec{\xi}}\right]$
-\end_inset
-
- is the Jacobian determinant of the reference map 
-\begin_inset Formula $\vec{x}_{e}(\vec{\xi}):\hat{\Omega}\to\Omega_{e}$
-\end_inset
-
-.
- This map describes how the physical points 
-\begin_inset Formula $\vec{x}\in\Omega_{e}$
-\end_inset
-
- are transformed from the reference points 
-\begin_inset Formula $\vec{\xi}\in\hat{\Omega}$
-\end_inset
-
-.
- Put another way, the inverse reference map 
-\color inherit
-
-\begin_inset Formula $\vec{\xi}=\vec{x}_{e}^{-1}(\vec{x})$
-\end_inset
-
- tells us how the physical domain 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
- maps into the reference cell 
-\begin_inset Formula $\hat{\Omega}$
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Standard
-At this point, we can assume without loss of generality that every physical
- cell 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
- can be derived from a single reference cell 
-\begin_inset Formula $\hat{\Omega}$
-\end_inset
-
-, so that our quadrature rule becomes simply a set of weights 
-\begin_inset Formula $w_{1},\ldots,w_{n_{Q}}$
-\end_inset
-
- and points 
-\begin_inset Formula $\vec{\xi}_{1},\ldots,\vec{\xi}_{n_{Q}}$
-\end_inset
-
- over 
-\begin_inset Formula $\hat{\Omega}$
-\end_inset
-
-.
- After changing variables, we end up with the final form of the quadrature
- rule that we can use to integrate the global residual field:
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\int_{\hat{\Omega}}\ F(\vec{x}_{e}(\vec{\xi}))J_{e}(\vec{\xi})\ d\vec{\xi}=\sum_{q=1}^{n_{Q}}w_{q}F(\vec{x}_{e}(\vec{\xi}_{q}))J_{e}(\vec{\xi}_{q})\]
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-If we let 
-\begin_inset Formula $F(\vec{x})=||u(\vec{x})-v(\vec{x})||^{2}$
-\end_inset
-
- in the above expression, we can find the cell residual contribution over
- 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
- from
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \begin{eqnarray*}
-\varepsilon_{e}^{2} & = & \int_{\Omega_{e}}\ ||u(\vec{x})-v(\vec{x})||^{2}\ d\vec{x}\\
- & = & \int_{\hat{\Omega}}\ ||u(\vec{x}_{e}(\vec{\xi}))-v(\vec{x}_{e}(\vec{\xi}))||^{2}J_{e}(\vec{\xi})\ d\vec{\xi}\\
- & = & \sum_{q=1}^{\mathrm{n_{Q}}}w_{q}||u(\vec{x}_{e}(\vec{\xi}_{q}))-v(\vec{x}_{e}(\vec{\xi}_{q}))||^{2}J_{e}(\vec{\xi}_{q})\end{eqnarray*}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-The global error 
-\begin_inset Formula $\varepsilon=\sqrt{\sum_{e}\varepsilon_{e}^{2}}$
-\end_inset
-
- is then approximated by the expression
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \begin{eqnarray*}
-\varepsilon & = & \sqrt{\sum_{e=1}^{\mathrm{n_{el}}}\sum_{q=1}^{\mathrm{n_{Q}}}w_{q}||u(\vec{x}_{e}(\vec{\xi}_{q}))-v(\vec{x}_{e}(\vec{\xi}_{q}))||^{2}J_{e}(\vec{\xi}_{q})}\end{eqnarray*}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-We can make a number of observations based on our final expression for the
- global error, (XXX: add observations here).
-\end_layout
-
-\begin_layout Section
-Comparing Residuals
-\end_layout
-
-\begin_layout Standard
-For comparing errors of different solutions, you can normalize the global
- error by the norm of the exact solution.
- For a family of solutions 
-\begin_inset Formula $u^{h}(\vec{x})$
-\end_inset
-
-, parametrized by the the maximum element size 
-\begin_inset Formula $h$
-\end_inset
-
- of the underlying discretization, and an exact solution 
-\begin_inset Formula $u(\vec{x})$
-\end_inset
-
-, this normalized error is given by
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \begin{eqnarray*}
-\varepsilon_{rel} & = & \frac{||u-u^{h}||_{L_{2}}}{||u||_{L_{2}}}\\
- & = & \frac{\int_{\Omega}||u(\vec{x})-u^{h}(\vec{x})||^{2}d\vec{x}}{\int_{\Omega}||u(\vec{x})||^{2}d\vec{x}}\end{eqnarray*}
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-This normalized error can be interpreted as the average error in the physical
- quantity being evaluated, so that a value of 0.01 corresponds to a 1% averaged
- error.
-\end_layout
-
-\begin_layout Standard
-Even if the exact solution is not currently known, this normalized error
- may be used to test the accuracy between two or more numerical solutions,
- defined on successively refined meshes.
-\end_layout
-
-\begin_layout Chapter
-\begin_inset LatexCommand label
-name "cha:Running-Cigma"
-
-\end_inset
-
-Running Cigma
-\end_layout
-
-\begin_layout Standard
-Cigma is primarily designed for calculating error estimates between arbitrary
- fields, so its primary usage is centered on the following operation
-\end_layout
-
-\begin_layout LyX-Code
-$ 
-\series bold
-cigma compare -a 
-\emph on
-\bar under
-field1
-\emph default
-\bar default
- -b 
-\emph on
-\bar under
-field2
-\emph default
-\bar default
- -o residuals.vtk
-\end_layout
-
-\begin_layout Standard
-You will need to provide two datasets describing each of the two fields,
- along with an integration rule and a discretization over which to integrate,
- although these last two will have reasonable defaults if they are not specified.
- This chapter aims to summarize the information you'll need to run Cigma.
-\end_layout
-
-\begin_layout Section
-Command Line Interface
-\end_layout
-
-\begin_layout Standard
-Cigma is designed to be scriptable.
- Thus, all operations can be specified as command-line arguments given to
- a single executable called 
-\family typewriter
-cigma
-\family default
-.
- A set of available commands can be obtained by typing
-\end_layout
-
-\begin_layout LyX-Code
-
-\family typewriter
-$ 
-\family default
-\series bold
-cigma help
-\end_layout
-
-\begin_layout LyX-Code
-
-\end_layout
-
-\begin_layout LyX-Code
-Usage: cigma <subcommand> [options]
-\end_layout
-
-\begin_layout LyX-Code
-Type 'cigma help <subcommand>' for help on a specific subcommand.
-\end_layout
-
-\begin_layout LyX-Code
-Available subcommands
-\end_layout
-
-\begin_layout LyX-Code
-   help
-\end_layout
-
-\begin_layout LyX-Code
-   list
-\end_layout
-
-\begin_layout LyX-Code
-   extract
-\end_layout
-
-\begin_layout LyX-Code
-   eval
-\end_layout
-
-\begin_layout LyX-Code
-   compare 
-\end_layout
-
-\begin_layout Standard
-XXX Quick description of each command
-\end_layout
-
-\begin_layout Standard
-XXX: Also add the same description to the help page of the respective command
-\end_layout
-
-\begin_layout Section
-Input and Output Formats
-\end_layout
-
-\begin_layout Standard
-XXX Reformat from newfile1.lyx
-\end_layout
-
-\begin_layout Standard
-The default data storage format for Cigma is the Hierarchical Data Format
- (HDF5), a portable file format developed at the 
-\begin_inset LatexCommand htmlurl
-name "National Center for Supercomputing Applications (NCSA)"
-target "hdf.ncsa.uiuc.edu/HDF5"
-
-\end_inset
-
-.
- The HDF5 is designed for storing multi-dimensional arrays together with
- meta-data in a portable self-describing format.
-\end_layout
-
-\begin_layout Standard
-The underlying data storage format for Cigma is the HDF5 format, due to
- its flexibility for storing and organizing large amounts of data.
- The Hierarchical Data Format (HDF) is designed for storing, retrieving,
- analyzing, visualizing, and converting scientific data.
- It uses a hierarchical structure that provides application programmers
- with a host of options for organizing how data is stored in HDF5 files.
- Using HDF5 datasets in Cigma allows us to avoid having to convert between
- too many distinct formats.
- Moreover, due to the amount of disk I/O, large finite element meshes can
- be handled more efficiently in binary format.
- 
-\end_layout
-
-\begin_layout Standard
-Another VTK format for providing mesh and field inputs.
-
-\color none
- 
-\color inherit
-[<-- XXX something wrong with that fragment.
- Did you mean 
-\begin_inset Quotes sld
-\end_inset
-
-Another format for providing mesh and field inputs is the VTK format
-\begin_inset Quotes srd
-\end_inset
-
-?] You can easily examine the structure of an input file by using the 
-\family typewriter
-cigma list
-\family default
- command, which will simply reveal the names and dimensions of all datasets
- inside the specified file.
-\end_layout
-
-\begin_layout Standard
-Inputs can be either HDF5 or VTK.
- Specifying the complete path to a dataset consists of the special form
- 
-\family typewriter
-\series bold
-filepath:dataset
-\family default
-\series default
-, a colon-delimited pair of file path and dataset path.
-\end_layout
-
-\begin_layout Standard
-Mesh inputs are currently only unstructured grids.
-\end_layout
-
-\begin_layout Standard
-Only exception are the residuals, which are written in legacy VTK format
- as scalar cell data over an unstructured grid.
-\end_layout
-
-\begin_layout Standard
-XXX New section?
-\end_layout
-
-\begin_layout Standard
-Because Cigma relies on the ability of the user to specify dataset paths,
- we have provided a command called 
-\family typewriter
-list
-\family default
- for viewing the structure of an input file.
- Its usage is very simple.
- 
-\end_layout
-
-\begin_layout Standard
-To view the structure of an HDF5 file:
-\end_layout
-
-\begin_layout LyX-Code
-$ 
-\series bold
-cigma list file.h5
-\end_layout
-
-\begin_deeper
-\begin_layout LyX-Code
-/mesh/coordinates           Dataset {119827, 3}
-\end_layout
-
-\begin_layout LyX-Code
-/mesh/connectivity          Dataset {661929, 4}
-\end_layout
-
-\begin_layout LyX-Code
-/vars/displacement/step0    Dataset {119827, 3}
-\end_layout
-
-\end_deeper
-\begin_layout Standard
-You can also view the structure of a VTK file with this command:
-\end_layout
-
-\begin_layout LyX-Code
-$ 
-\series bold
-cigma list file.vtk
-\end_layout
-
-\begin_deeper
-\begin_layout LyX-Code
-Reading file.vtk
-\end_layout
-
-\begin_layout LyX-Code
-Points = 119827
-\end_layout
-
-\begin_layout LyX-Code
-Cells = 661929
-\end_layout
-
-\begin_layout LyX-Code
-PointDataArray[0] = displacements_t0 (119827 x 3)
-\end_layout
-
-\end_deeper
-\begin_layout LyX-Code
-
-\end_layout
-
-\begin_layout Section
-Cigma Datasets
-\end_layout
-
-\begin_layout Standard
-Here we describe in a bit more detail the kind of datasets that Cigma expects.
- 
-\end_layout
-
-\begin_layout Standard
-The first kind of dataset required are meshes.
- In this version of Cigma, meshes are expected to be specifed as an unstructured
- grid.
-\end_layout
-
-\begin_layout Standard
-A field is a function which assigns a physical quantity to every point in
- a domain 
-\begin_inset Formula $\Omega$
-\end_inset
-
-.
- This quantity may correspond to a scalar, a vector, or even a tensor.
- As in Chapter 
-\begin_inset LatexCommand ref
-reference "cha:Error-Analysis"
-
-\end_inset
-
-, we shall use the same discretization 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
-.
- 
-\end_layout
-
-\begin_layout Standard
-A finite element approximation to an solution field 
-\begin_inset Formula $\phi(\vec{x})$
-\end_inset
-
- on an element 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
- is essentially specified by the weighed sum over a fixed set of local shape
- functions 
-\begin_inset Formula $\phi_{e,1}(\vec{\xi}),\phi_{e,2}(\vec{\xi}),\ldots,\phi_{e,N}(\vec{\xi})$
-\end_inset
-
- defined on the appropriate reference cell,
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\phi(\vec{x})=\sum_{n=1}^{N}d_{e,n}\phi_{e,n}(\vec{\xi})\]
-
-\end_inset
-
-where the global point 
-\begin_inset Formula $\vec{x}\in\Omega_{e}$
-\end_inset
-
- corresponds to the local reference point 
-\begin_inset Formula $\vec{\xi}=\vec{x}_{e}^{-1}(\vec{x})\in\hat{\Omega}$
-\end_inset
-
-, and the weights 
-\begin_inset Formula $d_{e,n}$
-\end_inset
-
-, also known as degrees of freedom, are given on.
- [<-- XXX does this sentence end correctly?] These shape functions define
- a basis for the function space on the finite element 
-\begin_inset Formula $\Omega_{e}$
-\end_inset
-
-, [<-- XXX more to this?]
-\end_layout
-
-\begin_layout Standard
-For a mesh consisting of a single element type, we need to specify the 
-\begin_inset Formula $(x_{n},y_{n},z_{n})$
-\end_inset
-
- of its degrees of freedom, and the connectivity relations 
-\begin_inset Formula $\Omega_{e}=\{n_{1},n_{2},\ldots\}$
-\end_inset
-
- among them which define each individual element in the corresponding discretiza
-tion.
-\end_layout
-
-\begin_layout Standard
-The most common element types used in the finite element method are given
- by the following shape functions, which define the interior of a tetrahedron
- with vertices at 
-\begin_inset Formula $a=(0,0,0)$
-\end_inset
-
-, 
-\begin_inset Formula $b=(1,0,0)$
-\end_inset
-
-, 
-\begin_inset Formula $c=(0,1,0)$
-\end_inset
-
-, 
-\begin_inset Formula $d=(0,0,1)$
-\end_inset
-
-, 
-\end_layout
-
-\begin_layout Standard
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status open
-
-\begin_layout Standard
-\begin_inset Box Frameless
-position "c"
-hor_pos "c"
-has_inner_box 1
-inner_pos "b"
-use_parbox 0
-width "50col%"
-special "none"
-height "1in"
-height_special "totalheight"
-status open
-
-\begin_layout Standard
-\begin_inset Formula \begin{eqnarray*}
-N_{a} & = & \frac{1}{2}(-1-x-y-z)\\
-N_{b} & = & \frac{1}{2}(1+x)\\
-N_{c} & = & \frac{1}{2}(1+y)\\
-N_{d} & = & \frac{1}{2}(1+z)\end{eqnarray*}
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\begin_inset Box Frameless
-position "c"
-hor_pos "c"
-has_inner_box 1
-inner_pos "c"
-use_parbox 0
-width "50col%"
-special "none"
-height "1in"
-height_special "totalheight"
-status collapsed
-
-\begin_layout Standard
-\begin_inset Graphics
-	filename figures/reference-tet4.png
-	width 5cm
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset Caption
-
-\begin_layout Standard
-Reference Tetrahedron
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-Likewise, a hexahedral element with vertices located at 
-\begin_inset Formula $a=(-1,-1,-1)$
-\end_inset
-
-, 
-\begin_inset Formula $b=(1,-1,-1)$
-\end_inset
-
-, 
-\begin_inset Formula $c=(1,1,-1)$
-\end_inset
-
-, 
-\begin_inset Formula $d=(-1,1,-1)$
-\end_inset
-
-, 
-\begin_inset Formula $e=(-1,-1,1)$
-\end_inset
-
-, 
-\begin_inset Formula $f=(1,-1,1)$
-\end_inset
-
-, 
-\begin_inset Formula $g=(1,1,1)$
-\end_inset
-
-, 
-\begin_inset Formula $h=(-1,1,1)$
-\end_inset
-
-.
-\end_layout
-
-\begin_layout Standard
-\begin_inset Float figure
-placement H
-wide false
-sideways false
-status open
-
-\begin_layout Standard
-\begin_inset Box Frameless
-position "c"
-hor_pos "c"
-has_inner_box 1
-inner_pos "c"
-use_parbox 0
-width "50col%"
-special "none"
-height "1in"
-height_special "totalheight"
-status collapsed
-
-\begin_layout Standard
-\begin_inset Graphics
-	filename figures/reference-hex8.png
-	width 5cm
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\begin_inset Box Frameless
-position "c"
-hor_pos "c"
-has_inner_box 1
-inner_pos "c"
-use_parbox 0
-width "50col%"
-special "none"
-height "1in"
-height_special "totalheight"
-status open
-
-\begin_layout Standard
-\begin_inset Formula \begin{eqnarray*}
-N_{a} & = & \frac{1}{8}\left(1-x\right)\left(1-y\right)\left(1-z\right)\\
-N_{b} & = & \frac{1}{8}\left(1+x\right)\left(1-y\right)\left(1-z\right)\\
-N_{c} & = & \frac{1}{8}\left(1-x\right)\left(1+y\right)\left(1-z\right)\\
-N_{d} & = & \frac{1}{8}\left(1+x\right)\left(1+y\right)\left(1-z\right)\\
-N_{e} & = & \frac{1}{8}\left(1-x\right)\left(1-y\right)\left(1+z\right)\\
-N_{f} & = & \frac{1}{8}\left(1+x\right)\left(1-y\right)\left(1+z\right)\\
-N_{g} & = & \frac{1}{8}\left(1-x\right)\left(1+y\right)\left(1+z\right)\\
-N_{h} & = & \frac{1}{8}\left(1+x\right)\left(1+y\right)\left(1+z\right)\end{eqnarray*}
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset Caption
-
-\begin_layout Standard
-Reference Hexahedron
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-XXX: Take description of higher order elements from Karniadakis book
-\end_layout
-
-\begin_layout Subsection
-Comparing Two Finite Element Fields
-\end_layout
-
-\begin_layout Standard
-Comparing two arbitrary finite element fields can be accomplished with the
- 
-\family typewriter
-cigma compare
-\family default
- command-line utility.
- By default, the comparison will take place over the elements in the mesh
- of the first field.
- Finally, the square of each of the local residual values are written to
- the specified VTK output file as cell-based scalars.
-\end_layout
-
-\begin_layout Standard
-A basic comparison can be as simple as specifying the following arguments:
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare --output=residuals.vtk
-\end_layout
-
-\begin_layout LyX-Code
-
-\series bold
- 
-\series default
-             --first=field1.h5:/field1/stepN
-\end_layout
-
-\begin_layout LyX-Code
-
-\series bold
- 
-\series default
-             --second=field2.h5:/field2/stepN
-\end_layout
-
-\begin_layout LyX-Code
-
-\end_layout
-
-\begin_layout Subsection
-Comparing against Known Values
-\end_layout
-
-\begin_layout Standard
-A finite element description might not always be available for one of the
- fields.
- However, you can break the comparison into several steps if you have a
- means to compute that field on any of the required points.
-\end_layout
-
-\begin_layout Standard
-First, extract the global coordinates of the integration points.
- This will result in an explicit list of points over which to evaluate your
- field.
-\end_layout
-
-\begin_layout LyX-Code
-cigma extract --mesh=field1.h5:/model/mesh/
-\end_layout
-
-\begin_layout LyX-Code
-              --output=points.h5:/projected_points
-\end_layout
-
-\begin_layout Standard
-Now you can evaluate your function at the designated points.
- You can provide the path to the explicit set of values with
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare --output=residuals.vtk
-\end_layout
-
-\begin_layout LyX-Code
-              --first=field1.h5:/field1/stepN
-\end_layout
-
-\begin_layout LyX-Code
-              --second-points=points.h5:/projected_points
-\end_layout
-
-\begin_layout LyX-Code
-              --second-values=values.h5:/projected_values
-\end_layout
-
-\begin_layout Subsection
-\begin_inset LatexCommand label
-name "sub:Comparing-against-a"
-
-\end_inset
-
-Comparing against a Known Function
-\end_layout
-
-\begin_layout Standard
-If one of your fields is easily described by an analytic expression, then
- you also have the option to compile your analytic function into Cigma,
- which will enable the 
-\family typewriter
-compare
-\family default
- command to reference your function by name.
- For example,
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare --output=residuals.vtk
-\end_layout
-
-\begin_layout LyX-Code
-
-\series bold
- 
-\series default
-             --first=field1.h5:/vars/displacement/step0
-\end_layout
-
-\begin_layout LyX-Code
-
-\series bold
- 
-\series default
-             --second=
-\bar under
-disloc3d
-\series bold
-\bar default
- 
-\series default
- 
-\end_layout
-
-\begin_layout Standard
-You may also interact with your analytic function by using the 
-\family typewriter
-cigma
-\family default
- 
-\family typewriter
-eval
-\family default
- command, and obtain a set of values which may then be passed back to the
- 
-\family typewriter
-compare
-\family default
- command.
- 
-\end_layout
-
-\begin_layout LyX-Code
-cigma eval --function=
-\bar under
-disloc3d
-\end_layout
-
-\begin_layout LyX-Code
-           --points=points.h5:/projected_points
-\end_layout
-
-\begin_layout LyX-Code
-           --values=values.h5:/disloc3d_values
-\end_layout
-
-\begin_layout Standard
-Once these values have been obtained, we can use it to resume the calculation
- of 
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare --output=residuals.vtk
-\end_layout
-
-\begin_layout LyX-Code
-              --first=field1.h5:/vars/displacement/step0
-\end_layout
-
-\begin_layout LyX-Code
-              --second-points=points.h5:/projected_points
-\end_layout
-
-\begin_layout LyX-Code
-              --second-values=values.h5:/disloc3d_values
-\end_layout
-
-\begin_layout Standard
-Another built-in function you might find useful is the 
-\family typewriter
-\bar under
-zero
-\family default
-\bar default
- function, which is defined to return 0 everywhere.
- You can use this function for estimating the norm of your field.
-\end_layout
-
-\begin_layout Subsection
-Specifying an Integration Mesh
-\end_layout
-
-\begin_layout Standard
-To override the mesh used in the integration, you can specify an extra argument
- providing the location of the mesh,
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare [...]
-\end_layout
-
-\begin_layout LyX-Code
-  --quadrature-mesh=mesh.h5:/model/mesh/
-\end_layout
-
-\begin_layout Standard
-Alternatively, you can also specify the coordinates and connectivity arrays
- separately, in case they reside in separate files.
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare [...]
-\end_layout
-
-\begin_layout LyX-Code
-  --quadrature-mesh-coordinates=file1.h5:/model/mesh/coordinates
-\end_layout
-
-\begin_layout LyX-Code
-  --quadrature-mesh-connectivity=file2.h5:/model/mesh/connectivity
-\end_layout
-
-\begin_layout Subsection
-Specifying a Quadrature Rule
-\end_layout
-
-\begin_layout Standard
-To specify a quadrature rule, you will have to provide the quadrature weights
- and points in the appropriate reference element.
- This can be done with the following additional argument:
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare [...]
-\end_layout
-
-\begin_layout LyX-Code
-  --quadrature-rule=quadrature-rules.h5:/path/to/rule
-\end_layout
-
-\begin_layout Standard
-You may also specify the location of the points and weights separately:
-\end_layout
-
-\begin_layout LyX-Code
-cigma compare [...]
-\end_layout
-
-\begin_layout LyX-Code
-  --quadrature-rule-points=file.h5:/path/to/rule/points
-\end_layout
-
-\begin_layout LyX-Code
-  --quadrature-rule-weights=file.h5:/path/to/rule/weights
-\end_layout
-
-\begin_layout LyX-Code
-
-\end_layout
-
-\begin_layout Standard
-In the 
-\family typewriter
-src/
-\family default
- directory of the Cigma distribution, you can find a Python script called
- 
-\family typewriter
-rules.py
-\family default
- that uses FIAT for calculating an appropriate set of quadrature rules for
- use in Cigma.
-\end_layout
-
-\begin_layout Chapter
-Examples
-\end_layout
-
-\begin_layout Standard
-In this chapter we will use sample datasets from the strike-slip benchmark
- case defined by the CIG Short-Term Tectonics working group.
- In this benchmark problem, we solve for the viscoelastic relaxation of
- stresses from a single finite earthquake while ignoring gravity.
- The problem is defined on a cube domain with sides of 24 km consisting
- of two layers of different material types.
- The top layer of the cube is nearly elastic, while the bottom layer is
- viscoelastic.
- (XXX: Anything else from Pylith manual or source?)
-\end_layout
-
-\begin_layout Standard
-\begin_inset Float figure
-wide false
-sideways false
-status open
-
-\begin_layout Standard
-\begin_inset Graphics
-	filename figures/strike-slip-fault.png
-	width 5cm
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-\begin_inset Caption
-
-\begin_layout Standard
-Geometry of Strike-Slip Benchmark
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Section
-Convergence
-\end_layout
-
-\begin_layout Standard
-When the exact solution of your equations is known, the norm of the error
- for a finite element solution is easily computed through the procedure
- given in Chapter 
-\begin_inset LatexCommand ref
-reference "cha:Error-Analysis"
-
-\end_inset
-
-.
- 
-\end_layout
-
-\begin_layout Standard
-Another way of checking the accuracy of a solution is by rerunning the same
- problem with a finer mesh.
- Recall that in the finite element method, the weight functions and trial
- solutions are chosen so that the finite element method is convergent.
- In other words, as the element size 
-\begin_inset Formula $h$
-\end_inset
-
- decreases, the solution tends to the correct solution.
- Thus, if the difference between the coarse and fine mesh solutions is small,
- we can conclude that the coarse mesh solution is quite accurate.
- However, if a solution changes noticeably with refinement of the mesh,
- the coarse mesh solution is inaccurate, and even the finer mesh may be
- inadequate as well.
-\end_layout
-
-\begin_layout Standard
-As can be seen from these results, the log of the error varies linearly
- with element size and the slope depends on the order of the element.
- Denoting the slope by 
-\begin_inset Formula $\alpha$
-\end_inset
-
-, then the error in the function can be expressed as
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\log\varepsilon=\log C+\alpha\log h\]
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Standard
-where 
-\begin_inset Formula $\log C$
-\end_inset
-
- is an arbitrary constant, the y-intercept of the curve.
- The slope 
-\begin_inset Formula $\alpha$
-\end_inset
-
- is the rate of convergence of the element.
- We can rewrite this as
-\end_layout
-
-\begin_layout Standard
-\begin_inset Formula \[
-\varepsilon=Ch^{\alpha}\]
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Section
-Visualization
-\end_layout
-
-\begin_layout Standard
-Using MayaVi2, we may obtain a visualization of the Cigma datasets 
-\end_layout
-
-\begin_layout Chapter
-\start_of_appendix
-File Formats
-\end_layout
-
-\begin_layout Standard
-To differentiate between these formats, you will need to provide an appropriate
- extension for your input and output files.
- The supported extensions are currently 
-\family typewriter
-.h5
-\family default
- for HDF5, 
-\family typewriter
-.vtk
-\family default
- for VTK files, and 
-\family typewriter
-.txt
-\family default
- for simple text files.
-\end_layout
-
-\begin_layout Section
-Input File Format
-\end_layout
-
-\begin_layout Standard
-There are five different objects you will want to use as input to Cigma,
- all of which can be represented as two-dimensional datasets.
- These are the 
-\end_layout
-
-\begin_layout Itemize
-mesh coordinates
-\end_layout
-
-\begin_layout Itemize
-mesh connectivity
-\end_layout
-
-\begin_layout Itemize
-field coefficients (degrees of freedom)
-\end_layout
-
-\begin_layout Itemize
-quadrature rule points
-\end_layout
-
-\begin_layout Itemize
-quadrature rule weights
-\end_layout
-
-\begin_layout Subsection
-HDF5
-\end_layout
-
-\begin_layout Standard
-HDF5 files are organized in a hierarchical structure, similar to a UNIX
- file system.
- Two types of primary objects, groups and datasets, are stored in this structure.
- A group contains instances of zero or more groups or datasets, while a
- dataset stores a multi-dimensional array of data elements.
- Both kinds of objects are accompanied by supporting metadata known as attribute
-s.
-\end_layout
-
-\begin_layout Standard
-A dataset is physically stored in two parts: a header and a data array.
- The header contains miscellaneous metadata describing the dataset as well
- as information that is needed to interpret the array portion of the dataset.
- Essentially, it includes the name, datatype, dataspace, and storage layout
- of the dataset.
- The name is a text string identifying the dataset.
- The datatype describes the type of the data array elements.
- The dataspace defines the dimensionality of the dataset, i.e., the size and
- shape of the multi-dimensional array.
-\end_layout
-
-\begin_layout Standard
-In Cigma you always provide an explicit path to every dataset, so you have
- a fair amount of flexibility for how you organize your datasets inside
- HDF5 files.
- For example, a typical Cigma HDF5 file could have the following structure.
-\end_layout
-
-\begin_layout LyX-Code
-model.h5
-\end_layout
-
-\begin_layout LyX-Code
-
-\backslash
-__ model
-\end_layout
-
-\begin_layout LyX-Code
-    |__ mesh
-\end_layout
-
-\begin_layout LyX-Code
-    |   |__ coordinates    [nno x nsd]
-\end_layout
-
-\begin_layout LyX-Code
-    |   
-\backslash
-__ connectivity   [nel x ndof]
-\end_layout
-
-\begin_layout LyX-Code
-    
-\backslash
-__ variables
-\end_layout
-
-\begin_layout LyX-Code
-        |__ temperature
-\end_layout
-
-\begin_layout LyX-Code
-        |   |__ step00000  [nno x 1]
-\end_layout
-
-\begin_layout LyX-Code
-        |   |__ step00010  [nno x 1]
-\end_layout
-
-\begin_layout LyX-Code
-        |   
-\backslash
-...
-\end_layout
-
-\begin_layout LyX-Code
-        |__ displacement
-\end_layout
-
-\begin_layout LyX-Code
-        |   |__ step00000  [nno x 3]
-\end_layout
-
-\begin_layout LyX-Code
-        |   |__ step00010  [nno x 3]
-\end_layout
-
-\begin_layout LyX-Code
-        |   
-\backslash
-...
-\end_layout
-
-\begin_layout LyX-Code
-        
-\backslash
-__ velocity
-\end_layout
-
-\begin_layout LyX-Code
-            |__ step00000  [nno x 3]
-\end_layout
-
-\begin_layout LyX-Code
-            |__ step00010  [nno x 3]
-\end_layout
-
-\begin_layout LyX-Code
-            
-\backslash
-...
-\end_layout
-
-\begin_layout Standard
-Generally, Cigma will only require you to specify the path to a specific
- field dataset.
- If a small amount of metadata is present in your field dataset, the rest
- of the required information, such as which mesh and finite elements to
- use, will be deduced from that metadata.
-\end_layout
-
-\begin_layout Description
-MeshID an identifier assigned for use in linking child datasets to their
- parent mesh.
-\end_layout
-
-\begin_layout Description
-MeshLocation points to the HDF5 group which contains the appropriate coordinates
- and connectivity datasets.
-\end_layout
-
-\begin_layout Description
-FunctionSpace string identifier to determine which shape functions to use
- for interpolating values inside the element (e.g., tet4, hex8, quad4, tri3,
- ...).
-\end_layout
-
-\begin_layout Description
-DatasetType string identifier for classifying the type of data contained
- in the dataset (e.g., points, connectivity, degrees of freedom, quadrature
- rules, global quadrature points, global field values).
-\end_layout
-
-\begin_layout Subsection
-VTK
-\end_layout
-
-\begin_layout Standard
-The current version of Cigma only supports VTK unstructured grid datasets
- of a single element type.
- You can still compare various unstructured grids with different element
- types to each other.
-\end_layout
-
-\begin_layout Standard
-Note that while you typically provide a path (or name) for every dataset,
- this is not necessary when specifying a VTK mesh, since this data is taken
- from the special Points and Cells arrays, which you cannot rename.
- However, you will still need to provide a name when referring to the field
- coefficients, which are assumed to be stored as Point Data in the input
- VTK file.
-\end_layout
-
-\begin_layout Standard
-For more information, you may want to refer to Visualization ToolKit's 
-\begin_inset LatexCommand htmlurl
-name "File Formats"
-target "www.vtk.org/pdf/file-formats.pdf"
-
-\end_inset
-
- document.
- 
-\end_layout
-
-\begin_layout Subsection
-Text
-\end_layout
-
-\begin_layout Standard
-The text format is always in table form, with the dimensions of the table
- specified in the first line.
- For example, mesh coordinates can be specified in the following format
-\end_layout
-
-\begin_layout LyX-Code
-nno nsd
-\end_layout
-
-\begin_layout LyX-Code
-x1 y1 z1
-\end_layout
-
-\begin_layout LyX-Code
-x2 y2 z2
-\end_layout
-
-\begin_layout LyX-Code
-x3 y3 z3
-\end_layout
-
-\begin_layout LyX-Code
-...
-\end_layout
-
-\begin_layout Standard
-Mesh connectivity with 
-\end_layout
-
-\begin_layout LyX-Code
-nel ndofs
-\end_layout
-
-\begin_layout LyX-Code
-node_1 node_2 node_3 node_4 ...
-\end_layout
-
-\begin_layout LyX-Code
-node_1 node_2 node_3 node_4 ...
-\end_layout
-
-\begin_layout LyX-Code
-node_1 node_2 node_3 node_4 ...
-\end_layout
-
-\begin_layout LyX-Code
-...
-\end_layout
-
-\begin_layout Standard
-A generic field with 
-\family typewriter
-ndim
-\family default
- components (i.e., scalar, vector, or tensor) is specified by
-\end_layout
-
-\begin_layout LyX-Code
-num ndim
-\end_layout
-
-\begin_layout LyX-Code
-f1 f2 f3 ..
-\end_layout
-
-\begin_layout LyX-Code
-f1 f2 f3 ..
-\end_layout
-
-\begin_layout LyX-Code
-...
-\end_layout
-
-\begin_layout Standard
-In this last case, the number of rows could refer to either 
-\family typewriter
-nno
-\family default
- or 
-\family typewriter
-nel
-\family default
-, depending on whether the field is node-based or cell-based.
-\end_layout
-
-\begin_layout Section
-Output File Format
-\end_layout
-
-\begin_layout Standard
-The output format for residuals consists simply of a list of scalars for
- each element in the discretization.
- If you specify a 
-\family typewriter
-.vtk
-\family default
- extension for your output file, this will result in an scalar dataset named
- epsilon stored using the legacy VTK file format in the Cell Data section
- of the output file.
- Note that this output consists of the squared values of the local residuals,
- so further post-processing will be necessary.
-\end_layout
-
-\begin_layout Standard
-\begin_inset Include \include{gpl-license.lyx}
-preview false
-
-\end_inset
-
-
-\end_layout
-
-\begin_layout Bibliography
-\begin_inset LatexCommand bibitem
-label "1"
-key "Akin 2005"
-
-\end_inset
-
-Akin, J.E.
- (2005), 
-\emph on
-Finite Element Analysis with Error Estimators,
-\emph default
- Butterworth-Heinemann, Oxford, 447 pp.
- 
-\end_layout
-
-\begin_layout Bibliography
-\begin_inset LatexCommand bibitem
-label "2"
-key "Encyclopaedia of Cubature Formulas 2005"
-
-\end_inset
-
-Encyclopaedia of Cubature Formulas (1998-2005), Katholieke Universiteit
- Leuven, Dept.
- Computerwetenschappen, www.cs.kuleuven.be/~nines/research/ecf/
-\end_layout
-
-\begin_layout Bibliography
-\begin_inset LatexCommand bibitem
-label "3"
-key "Hughes 2000"
-
-\end_inset
-
-Hughes, Thomas J.R.
- (2000), 
-\emph on
-The Finite Element Method: Linear Static and Dynamic Finite Element Analysis
-\emph default
-, Dover, 672 pp.
-\end_layout
-
-\begin_layout Bibliography
-\begin_inset LatexCommand bibitem
-label "4"
-key "KarniadakisSherwin2005"
-
-\end_inset
-
-Karniadakis, George E.
- and Spencer J.
- Sherwin (2005), 
-\emph on
-Spectral/
-\emph default
-hp 
-\emph on
-Element Methods for Computational Fluid Dynamics, Second Edition
-\emph default
-, Numerical Mathematics and Scientific Computation, Oxford, 657 pp.
-\end_layout
-
-\begin_layout Bibliography
-\begin_inset LatexCommand bibitem
-label "5"
-key "Uesu-etal2005"
-
-\end_inset
-
-Uesu, D., L.
- Bavoil, S.
- Fleishman, J.
- Shepherd, and C.
- Silva (2005), Simplification of Unstructured Tetrahedral Meshes by Point-Sampli
-ng, 
-\emph on
-Proceedings of the 2005 International Workshop on Volume Graphics,
-\emph default
- 157-165, doi: 10.2312/VG/VG05/157-165.
-\end_layout
-
-\end_body
-\end_document
+#LyX 1.5.4 created this file. For more info see http://www.lyx.org/
+\lyxformat 276
+\begin_document
+\begin_header
+\textclass book
+\begin_preamble
+\usepackage{hyperref}
+
+\let\myUrl\url
+\renewcommand{\url}[1]{(\myUrl{#1})}
+\end_preamble
+\language english
+\inputencoding auto
+\font_roman default
+\font_sans default
+\font_typewriter default
+\font_default_family default
+\font_sc false
+\font_osf false
+\font_sf_scale 100
+\font_tt_scale 100
+\graphics default
+\paperfontsize default
+\spacing single
+\papersize default
+\use_geometry true
+\use_amsmath 1
+\use_esint 0
+\cite_engine basic
+\use_bibtopic false
+\paperorientation portrait
+\leftmargin 1in
+\topmargin 1in
+\rightmargin 1in
+\bottommargin 1in
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\defskip medskip
+\quotes_language swedish
+\papercolumns 1
+\papersides 2
+\paperpagestyle headings
+\tracking_changes false
+\output_changes false
+\author "" 
+\author "" 
+\end_header
+
+\begin_body
+
+\begin_layout Standard
+\noindent
+\align center
+
+\color none
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+thispagestyle{empty}
+\end_layout
+
+\end_inset
+
+
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/cigma_cover.pdf
+	display color
+	width 75page%
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\noindent
+\align center
+
+\color none
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+thispagestyle{empty}
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Title
+Cigma
+\end_layout
+
+\begin_layout Author
+© California Institute of Technology
+\newline
+Version 1.0.0
+\end_layout
+
+\begin_layout Date
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+today
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset LatexCommand tableofcontents
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset FloatList figure
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+raggedbottom
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+
+\newpage
+
+\end_layout
+
+\begin_layout Chapter
+Introduction
+\end_layout
+
+\begin_layout Section
+About Cigma
+\end_layout
+
+\begin_layout Standard
+The CIG Model Analyzer (Cigma) consists of a suite of tools for comparing
+ numerical models.
+ CIG has developed Cigma in response to demand from the short-term tectonics
+ community for a simple tool that can perform rigorous error analysis on
+ their finite element codes.
+ The long-term goal for Cigma, however, is for it to be used for nearly
+ all geodynamics modeling codes.
+\end_layout
+
+\begin_layout Standard
+The current version of Cigma is intended for the calculation of 
+\begin_inset Formula $L_{2}$
+\end_inset
+
+ residuals for finite element models.
+ This software was designed for performing the following three main tasks:
+ (1) error analysis, (2) benchmarking, and (3) code verification.
+ In error analysis, Cigma can calculate a global measure of the error between
+ two solutions by performing an integration over a discretized version of
+ the common domain.
+ This comparison can take place even when the underlying discretizations
+ do not overlap.
+ In benchmarking, Cigma can help the geodynamics community agree on a standard
+ solution to specific problems by facilitating the process of comparing
+ different numerical codes against each other.
+ Lastly, as an automated tool, Cigma can help application developers create
+ regression tests to ensure that software changes do not affect the consistency
+ of the results.
+\end_layout
+
+\begin_layout Standard
+At its core, Cigma draws from a variety of libraries.
+ The 
+\begin_inset LatexCommand htmlurl
+name "FInite element Automatic Tabulator (FIAT)"
+target "www.fenics.org/wiki/FIAT"
+
+\end_inset
+
+ Python Library supports generation of arbitrary order instances of the
+ Lagrange elements on lines, triangles, and tetrahedra.
+ It can also generate higher order instances of Jacobi-type quadrature rules
+ on the same element shapes.
+ The 
+\begin_inset LatexCommand htmlurl
+name "Approximate Nearest Neighbor (ANN)"
+target "www.cs.umd.edu/~mount/ANN/"
+
+\end_inset
+
+ Library, written in C++, supports data structures and algorithms for both
+ exact and nearest-neighbor searching in arbitrarily high dimensions.
+\end_layout
+
+\begin_layout Standard
+Both of these libraries extend and generalize Cigma's functionality so it
+ can handle other types of elements, and provide the ability to compare
+ vector fields.
+\end_layout
+
+\begin_layout Section
+Citation
+\end_layout
+
+\begin_layout Standard
+Computational Infrastructure for Geodynamics (CIG) is making this source
+ code available to you in the hope that the software will enhance your research
+ in geophysics.
+ This is a brand-new code and at present no papers are published or press.
+ Please cite this manual as follows:
+\end_layout
+
+\begin_layout Itemize
+Armendariz, L., and S.
+ Kientz.
+ 
+\emph on
+Cigma User Manual.
+
+\emph default
+ Pasadena, CA: Computational Infrastructure of Geodynamics, 2008.
+ URL: geodynamics.org/cig/software/cs/cigma/cigma.pdf
+\end_layout
+
+\begin_layout Standard
+CIG requests that in your oral presentations and in your papers that you
+ indicate your use of this code and acknowledge the author of the code and
+ 
+\begin_inset LatexCommand htmlurl
+name "CIG"
+target "geodynamics.org"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Section
+Support
+\end_layout
+
+\begin_layout Standard
+Cigma development is based upon work supported by the National Science Foundatio
+n under Grant No.
+ EAR-0406751.
+ Any opinions, findings, and conclusions or recommendations expressed in
+ this material are those of the authors and do not necessarily reflect the
+ views of the National Science Foundation.
+ The code is being released under the GNU General Public License.
+\end_layout
+
+\begin_layout Chapter
+Installation and Getting Help
+\end_layout
+
+\begin_layout Section
+Getting Help
+\end_layout
+
+\begin_layout Standard
+For help, send an e-mail to the 
+\begin_inset LatexCommand url
+name "CIG Computational Science Mailing List"
+target "cig-cs at geodynamics.org"
+
+\end_inset
+
+.
+ You can subscribe to the 
+\family typewriter
+cig-cs
+\family default
+ mailing list and view archived discussions at the 
+\begin_inset LatexCommand htmlurl
+name "CIG Mail Lists web page"
+target "geodynamics.org/cig/lists"
+
+\end_inset
+
+.
+ If you encounter any bugs or have problems installing Cigma, please submit
+ a report to the 
+\begin_inset LatexCommand htmlurl
+name "CIG Bug Tracker"
+target "geodynamics.org/bugs"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Section
+Installating from Source
+\end_layout
+
+\begin_layout Standard
+To install Cigma, download the source package from the 
+\begin_inset LatexCommand htmlurl
+name "CIG Cigma web page"
+target "geodynamics.org/cig/software/packages/cs/cigma"
+
+\end_inset
+
+.
+ You will need the GNU C++ compiler for this step.
+ The Boost, HDF5, and VTK libraries, described later in this section, are
+ also required.
+\end_layout
+
+\begin_layout Standard
+After installing the necessary dependencies, simply unpack the source tar
+ file and issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz cigma-1.0.0.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd cigma-1.0.0
+\end_layout
+
+\begin_layout LyX-Code
+$ ./configure --with-boost=$BOOST_PREFIX
+\end_layout
+
+\begin_layout LyX-Code
+              --with-hdf5=$HDF5_PREFIX 
+\end_layout
+
+\begin_layout LyX-Code
+              --with-vtk=$VTK_PREFIX 
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+Both HDF5 and VTK prefixes are required if you have installed these libraries
+ in a custom location.
+\end_layout
+
+\begin_layout Subsection
+Boost
+\end_layout
+
+\begin_layout Standard
+The Boost C++ libraries are used in Cigma for parsing command-line options,
+ and ensuring cross-platform access to the filesystem.
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz boost_1_35_0.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd boost_1_35_0
+\end_layout
+
+\begin_layout LyX-Code
+$ ./configure
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+Without a custom installation prefix, Boost will target your 
+\family typewriter
+/usr/local
+\family default
+ directory.
+\end_layout
+
+\begin_layout Subsection
+HDF5
+\end_layout
+
+\begin_layout Standard
+The Hierarchical Data Format (HDF5) library is available for download from
+ 
+\begin_inset LatexCommand htmlurl
+name "The HDF Group"
+target "hdfgroup.org/HDF5"
+
+\end_inset
+
+.
+ Binaries can be obtained at 
+\begin_inset LatexCommand htmlurl
+name "hdfgroup.org/HDF5/release/obtain5.html"
+target "hdfgroup.org/HDF5/release/obtain5.html"
+
+\end_inset
+
+.
+ To install from source, download the latest stable 1.6 version of this library
+ (1.6.7, for example) and issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz hdf5-1.6.7.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd hdf5-1.6.7
+\end_layout
+
+\begin_layout LyX-Code
+$ ./configure --prefix=$HOME/opt/hdf5/1.6.7
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+Without the custom prefix, this installation procedure will target your
+ 
+\family typewriter
+/usr
+\family default
+ directory.
+\end_layout
+
+\begin_layout Subsection
+VTK
+\end_layout
+
+\begin_layout Standard
+The Visualization Toolkit (VTK) library is a popular open source graphics
+ library for scientific visualizations.
+ We recommend that you obtain the binaries from a package manager, making
+ sure to install the associated development headers along with the library.
+ The source for the VTK library is available from 
+\begin_inset LatexCommand htmlurl
+name "Kitware, Inc."
+target "www.vtk.org/get-software.php"
+
+\end_inset
+
+.
+ If you decide to compile VTK from source, you will also need the CMake
+ build environment, available from 
+\begin_inset LatexCommand htmlurl
+name "CMake"
+target "cmake.org"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+After downloading the source package, you can issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz vtk-5.0.4.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd VTK
+\end_layout
+
+\begin_layout LyX-Code
+$ mkdir build
+\end_layout
+
+\begin_layout LyX-Code
+$ cd build
+\end_layout
+
+\begin_layout LyX-Code
+$ ccmake ..
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+You can set the installation prefix during the 
+\family typewriter
+ccmake
+\family default
+ step above, which defaults to your 
+\family typewriter
+/usr/local
+\family default
+ directory.
+\end_layout
+
+\begin_layout Section
+Additional Software
+\end_layout
+
+\begin_layout Subsection
+NumPy
+\end_layout
+
+\begin_layout Standard
+NumPy is a Python module which adds support for large multi-dimensional
+ arrays and matrices, and is available for download from the 
+\begin_inset LatexCommand htmlurl
+name "NumPy Home Page"
+target "numpy.scipy.org "
+
+\end_inset
+
+.
+ To install this extension module from source, download it and issue the
+ following command in the NumPy source directory
+\end_layout
+
+\begin_layout LyX-Code
+$ python setup.py install
+\end_layout
+
+\begin_layout Standard
+To verify the installation, the command `
+\family typewriter
+import numpy
+\family default
+' should run successfully from within your standard Python interpreter shell.
+\end_layout
+
+\begin_layout Subsection
+PyTables
+\end_layout
+
+\begin_layout Standard
+PyTables is a Python extension module that builds on top of the HDF5 library.
+ It provides a convenient scripting interface to manipulate HDF5 files,
+ which can be used to manipulate the input/output files created by Cigma.
+ PyTables is available for download from 
+\begin_inset LatexCommand htmlurl
+name "PyTables"
+target "www.pytables.org"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+To install this extension from source, download the latest stable version
+ (currently 2.0.2) and issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz pytables-2.0.2
+\end_layout
+
+\begin_layout LyX-Code
+$ cd pytables-2.0.2
+\end_layout
+
+\begin_layout LyX-Code
+$ python setup.py install
+\end_layout
+
+\begin_layout Standard
+To verify the installation, the command `
+\family typewriter
+import tables
+\family default
+' should run successfully from within your standard Python interpreter shell.
+\end_layout
+
+\begin_layout Subsection
+FIAT
+\end_layout
+
+\begin_layout Standard
+The Finite Element Automatic Tabulator (FIAT) Python module supports generation
+ of arbitrary order instances of the Lagrange elements, as well as arbitrary
+ order instances of the Jacobi-type quadrature rules on the same element
+ shapes.
+ It is available from 
+\begin_inset LatexCommand htmlurl
+name "The FEniCS Project"
+target "www.fenics.org/wiki/FIAT"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+To install this module, download it and issue the following command inside
+ the FIAT source directory:
+\end_layout
+
+\begin_layout LyX-Code
+$ sudo python setup.py install
+\end_layout
+
+\begin_layout Standard
+You should now be able to run `
+\family typewriter
+import FIAT
+\family default
+' from within your standard Python interpreter shell.
+\end_layout
+
+\begin_layout Subsection
+HDFView
+\end_layout
+
+\begin_layout Standard
+NCSA HDFView is a graphical user interface tool for accessing data in your
+ HDF5 files.
+ You can use it for viewing the internal group hierarchy, adding new datasets,
+ and modifying or deleting existing datasets, or altering metadata on groups
+ and datasets.
+ You can download it from the 
+\begin_inset LatexCommand htmlurl
+name "HDFView home page"
+target "hdf.ncsa.uiuc.edu/hdf-java-html/hdfview"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Chapter
+\begin_inset LatexCommand label
+name "cha:Error-Analysis"
+
+\end_inset
+
+Error Analysis
+\end_layout
+
+\begin_layout Section
+Introduction
+\end_layout
+
+\begin_layout Standard
+When studying differential equations that represent physical systems we
+ often obtain solutions by using a variety of techniques.
+ Solving these equations using classical analytical methods is almost impossible
+, so we often resort to numerical algorithms such as the finite element
+ method.
+\end_layout
+
+\begin_layout Standard
+For assessing the quality of the solution to our equations, it is important
+ to quantify the error in our numerical solutions.
+ To calculate this error, we will have to compare these solutions against
+ each other through the use of a distance measure between two functions.
+\end_layout
+
+\begin_layout Standard
+The simplest possible quantitative measure of the difference between two
+ distinct functions consists of taking the pointwise difference at a common
+ set of points.
+ While no finite sample of points can perfectly represent a continuum of
+ values, valuable information can be inferred from the statistics of the
+ resulting set of differences.
+ However, the fields we want to compare may not be defined at a common set
+ of points.
+ This simple measure becomes unapplicable.
+\end_layout
+
+\begin_layout Section
+Interpolation Functions
+\end_layout
+
+\begin_layout Standard
+We can apply a number of interpolation
+\end_layout
+
+\begin_layout Section
+Distance Measures
+\end_layout
+
+\begin_layout Standard
+A very useful distance measure we can use is the 
+\begin_inset Formula $L_{2}$
+\end_inset
+
+ norm, defined by the following integral
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\varepsilon=||u-v||_{L_{2}}=\sqrt{\int_{\Omega}||u(\vec{x})-v(\vec{x})||^{2}d\vec{x}}\]
+
+\end_inset
+
+where 
+\begin_inset Formula $u$
+\end_inset
+
+ and 
+\begin_inset Formula $v$
+\end_inset
+
+ are the two functions on a global coordinate system.
+ This gives us a single global estimate 
+\begin_inset Formula $\varepsilon$
+\end_inset
+
+ representing the distance between the two fields 
+\begin_inset Formula $u(\vec{x})$
+\end_inset
+
+ and 
+\begin_inset Formula $v(\vec{x})$
+\end_inset
+
+.
+ Alternatively, you may think of this as the size, or norm, of the residual
+ field 
+\begin_inset Formula $\rho(\vec{x})=u(\vec{x})-v(\vec{x})$
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Section
+Global Error Approximation
+\end_layout
+
+\begin_layout Standard
+If we discretize the domain 
+\begin_inset Formula $\Omega$
+\end_inset
+
+ into an appropriate set of cells 
+\begin_inset Formula $\Omega_{1},\Omega_{2},\ldots,\Omega_{n_{el}}$
+\end_inset
+
+, we can compute the global error 
+\begin_inset Formula $\varepsilon^{2}$
+\end_inset
+
+ as a sum over localized cell contributions, 
+\begin_inset Formula $\varepsilon^{2}=\sum_{e}\varepsilon_{e}^{2}$
+\end_inset
+
+, where each cell residual 
+\begin_inset Formula $\varepsilon_{e}^{2}$
+\end_inset
+
+ is given by
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon_{e}^{2} & = & \int_{\Omega_{e}}||u(\vec{x})-v(\vec{x})||^{2}d\vec{x}\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+In general, we won't be able to integrate each cell residual 
+\begin_inset Formula $\varepsilon_{e}^{2}$
+\end_inset
+
+ exactly since either of the functions 
+\begin_inset Formula $u$
+\end_inset
+
+ and 
+\begin_inset Formula $v$
+\end_inset
+
+ may have an incompatible representation relative to the finite element
+ space on each domain 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+.
+ However, we can still calculate an approximation of each cell residual
+ by applying an appropriate quadrature rule with a tolerable truncation
+ error 
+\begin_inset LatexCommand cite
+key "Encyclopaedia of Cubature Formulas 2005"
+
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+To obtain an approximation to the integral of a function 
+\begin_inset Formula $F(\vec{x})$
+\end_inset
+
+ over a cell 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+, we simply apply the quadrature rule with weights 
+\begin_inset Formula $w_{e,1},w_{e,2},\ldots,w_{e,n_{Q}}$
+\end_inset
+
+ and integration points 
+\begin_inset Formula $\vec{x}_{e,1},\vec{x}_{e,2},\ldots,\vec{x}_{e,n_{Q}}$
+\end_inset
+
+ appropriate for the physical element 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\int_{\Omega_{e}}\ F(\vec{x})\ d\vec{x}=\sum_{q=1}^{n_{Q}}w_{e,q}F(\vec{x}_{e,q})\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Applying this quadrature rule directly over the entire physical domain 
+\begin_inset Formula $\Omega=\cup\Omega_{e}$
+\end_inset
+
+ gives us 
+\begin_inset Formula \[
+\int_{\Omega}\ F(\vec{x})\ d\vec{x}=\sum_{e=1}^{n_{el}}\sum_{q=1}^{n_{Q}}w_{e,q}F(\vec{x}_{e,q})\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+For efficiency reasons, it is undesirable in finite element applications
+ to perform calculations in a global coordinate system.
+ To avoid duplication of work, shape function evaluations may be performed
+ once on a reference cell 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+ and then transformed back into the corresponding physical cell 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ as needed.
+\end_layout
+
+\begin_layout Standard
+To compute integrals of 
+\begin_inset Formula $F$
+\end_inset
+
+ in a reference coordinate system, we need to apply a change of variables:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\int_{\Omega_{e}}\ F(\vec{x})\ d\vec{x}=\int_{\hat{\Omega}}\ F(\vec{x}_{e}(\vec{\xi}))J_{e}(\vec{\xi})\ d\vec{\xi}\]
+
+\end_inset
+
+where 
+\color none
+the additional factor 
+\begin_inset Formula $J_{e}(\vec{\xi})=\det\left[\frac{d\vec{x}_{e}}{d\vec{\xi}}\right]$
+\end_inset
+
+ is the Jacobian determinant of the reference map 
+\begin_inset Formula $\vec{x}_{e}(\vec{\xi}):\hat{\Omega}\to\Omega_{e}$
+\end_inset
+
+.
+ This map describes how the physical points 
+\begin_inset Formula $\vec{x}\in\Omega_{e}$
+\end_inset
+
+ are transformed from the reference points 
+\begin_inset Formula $\vec{\xi}\in\hat{\Omega}$
+\end_inset
+
+.
+ Put another way, the inverse reference map 
+\color inherit
+
+\begin_inset Formula $\vec{\xi}=\vec{x}_{e}^{-1}(\vec{x})$
+\end_inset
+
+ tells us how the physical domain 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ maps into the reference cell 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+At this point, we can assume without loss of generality that every physical
+ cell 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ can be derived from a single reference cell 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+, so that our quadrature rule becomes simply a set of weights 
+\begin_inset Formula $w_{1},\ldots,w_{n_{Q}}$
+\end_inset
+
+ and points 
+\begin_inset Formula $\vec{\xi}_{1},\ldots,\vec{\xi}_{n_{Q}}$
+\end_inset
+
+ over 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+.
+ After changing variables, we end up with the final form of the quadrature
+ rule that we can use to integrate the global residual field:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\int_{\hat{\Omega}}\ F(\vec{x}_{e}(\vec{\xi}))J_{e}(\vec{\xi})\ d\vec{\xi}=\sum_{q=1}^{n_{Q}}w_{q}F(\vec{x}_{e}(\vec{\xi}_{q}))J_{e}(\vec{\xi}_{q})\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+If we let 
+\begin_inset Formula $F(\vec{x})=||u(\vec{x})-v(\vec{x})||^{2}$
+\end_inset
+
+ in the above expression, we can find the cell residual contribution over
+ 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ from
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon_{e}^{2} & = & \int_{\Omega_{e}}\ ||u(\vec{x})-v(\vec{x})||^{2}\ d\vec{x}\\
+ & = & \int_{\hat{\Omega}}\ ||u(\vec{x}_{e}(\vec{\xi}))-v(\vec{x}_{e}(\vec{\xi}))||^{2}J_{e}(\vec{\xi})\ d\vec{\xi}\\
+ & = & \sum_{q=1}^{\mathrm{n_{Q}}}w_{q}||u(\vec{x}_{e}(\vec{\xi}_{q}))-v(\vec{x}_{e}(\vec{\xi}_{q}))||^{2}J_{e}(\vec{\xi}_{q})\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+The global error 
+\begin_inset Formula $\varepsilon=\sqrt{\sum_{e}\varepsilon_{e}^{2}}$
+\end_inset
+
+ is then approximated by the expression
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon & = & \sqrt{\sum_{e=1}^{\mathrm{n_{el}}}\sum_{q=1}^{\mathrm{n_{Q}}}w_{q}||u(\vec{x}_{e}(\vec{\xi}_{q}))-v(\vec{x}_{e}(\vec{\xi}_{q}))||^{2}J_{e}(\vec{\xi}_{q})}\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+We can make a number of observations based on our final expression for the
+ global error,
+\end_layout
+
+\begin_layout Section
+Comparing Residuals
+\end_layout
+
+\begin_layout Standard
+For comparing errors of different solutions, you can normalize the global
+ error by the norm of the exact solution.
+ For a family of solutions 
+\begin_inset Formula $u^{h}(\vec{x})$
+\end_inset
+
+, parametrized by the the maximum element size 
+\begin_inset Formula $h$
+\end_inset
+
+ of the underlying discretization, and an exact solution 
+\begin_inset Formula $u(\vec{x})$
+\end_inset
+
+, this normalized error is given by
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon_{rel} & = & \frac{||u-u^{h}||_{L_{2}}}{||u||_{L_{2}}}\\
+ & = & \frac{\int_{\Omega}||u(\vec{x})-u^{h}(\vec{x})||^{2}d\vec{x}}{\int_{\Omega}||u(\vec{x})||^{2}d\vec{x}}\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+This normalized error can be interpreted as the average error in the physical
+ quantity being evaluated, so that a value of 0.01 corresponds to a 1% averaged
+ error.
+\end_layout
+
+\begin_layout Standard
+Even if the exact solution is not currently known, this normalized error
+ may be used to test the accuracy between two or more numerical solutions,
+ defined on successively refined meshes.
+\end_layout
+
+\begin_layout Chapter
+\begin_inset LatexCommand label
+name "cha:Running-Cigma"
+
+\end_inset
+
+Running Cigma
+\end_layout
+
+\begin_layout Standard
+Cigma is primarily designed for calculating error estimates between arbitrary
+ fields, so its primary usage is centered on the following operation
+\end_layout
+
+\begin_layout LyX-Code
+$ 
+\series bold
+cigma compare 
+\emph on
+\bar under
+field1
+\emph default
+\bar default
+ 
+\emph on
+\bar under
+field2
+\emph default
+\bar default
+ -o residuals.vtk
+\end_layout
+
+\begin_layout Standard
+You will need to provide two datasets describing each of the two fields,
+ specify an integration rule and a domain discretization over which to integrate
+, although these last two will have reasonable defaults if they are not
+ specified.
+\end_layout
+
+\begin_layout Standard
+The default data storage format for Cigma is the Hierarchical Data Format
+ (HDF5), a portable file format developed at the 
+\begin_inset LatexCommand htmlurl
+name "National Center for Supercomputing Applications (NCSA)"
+target "hdf.ncsa.uiuc.edu/HDF5"
+
+\end_inset
+
+.
+ The HDF5 is designed for storing multi-dimensional arrays together with
+ meta-data in a portable self-describing format.
+\end_layout
+
+\begin_layout Section
+Command Line Interface
+\end_layout
+
+\begin_layout Standard
+Cigma is designed to be scriptable.
+ Thus, all operations can be specified as command-line arguments given to
+ a single executable called 
+\family typewriter
+cigma
+\family default
+.
+ A set of available commands can be obtained by typing
+\end_layout
+
+\begin_layout LyX-Code
+
+\family typewriter
+$ 
+\family default
+\series bold
+cigma help
+\end_layout
+
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout LyX-Code
+Usage: cigma <subcommand> [options]
+\end_layout
+
+\begin_layout LyX-Code
+Type 'cigma help <subcommand>' for help on a specific subcommand.
+\end_layout
+
+\begin_layout LyX-Code
+Available subcommands
+\end_layout
+
+\begin_layout LyX-Code
+   help
+\end_layout
+
+\begin_layout LyX-Code
+   list
+\end_layout
+
+\begin_layout LyX-Code
+   extract
+\end_layout
+
+\begin_layout LyX-Code
+   eval
+\end_layout
+
+\begin_layout LyX-Code
+   compare 
+\end_layout
+
+\begin_layout Standard
+Examples for each of these commands are given in the sections below.
+\end_layout
+
+\begin_layout Section
+Input and Output Formats
+\end_layout
+
+\begin_layout Standard
+The underlying data storage format for Cigma is the HDF5 format, due to
+ its flexibility for storing and organizing large amounts of data.
+ The Hierarchical Data Format (HDF) is designed for storing, retrieving,
+ analyzing, visualizing, and converting scientific data.
+ It uses a hierarchical structure that provides users a host of options
+ for organizing how their data is stored in HDF5 files.
+ Using HDF5 datasets in Cigma allows us to avoid having to convert between
+ too many distinct formats.
+ Moreover, due to the amount of disk I/O, large finite element meshes can
+ be handled more efficiently in binary format.
+ 
+\end_layout
+
+\begin_layout Standard
+Another popular format for providing mesh and field inputs.You can easily
+ examine the structure of an input file by using the 
+\family typewriter
+cigma list
+\family default
+ command, which will simply reveal the names and dimensions of all datasets
+ inside the specified file.
+\end_layout
+
+\begin_layout Standard
+Inputs can be either HDF5 or VTK.
+ Specifying the complete path to a dataset consists of the special form
+ 
+\family typewriter
+\series bold
+filepath:dataset
+\family default
+\series default
+, a colon-delimited pair of file path and dataset path.
+\end_layout
+
+\begin_layout Standard
+Mesh inputs are currently only unstructured grids.
+\end_layout
+
+\begin_layout Standard
+Only exception are the residuals, which are written in legacy VTK format
+ as scalar cell data over an unstructured grid.
+\end_layout
+
+\begin_layout Standard
+XXX New section?
+\end_layout
+
+\begin_layout Standard
+Because Cigma relies on the ability of the user to specify dataset paths,
+ we have provided a command called 
+\family typewriter
+list
+\family default
+ for viewing the structure of an input file.
+ Its usage is very simple.
+ 
+\end_layout
+
+\begin_layout Standard
+To view the structure of an HDF5 file:
+\end_layout
+
+\begin_layout LyX-Code
+$ 
+\series bold
+cigma list file.h5
+\end_layout
+
+\begin_deeper
+\begin_layout LyX-Code
+/mesh/coordinates           Dataset {119827, 3}
+\end_layout
+
+\begin_layout LyX-Code
+/mesh/connectivity          Dataset {661929, 4}
+\end_layout
+
+\begin_layout LyX-Code
+/vars/displacement/step0    Dataset {119827, 3}
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+You can also view the structure of a VTK file with this command:
+\end_layout
+
+\begin_layout LyX-Code
+$ 
+\series bold
+cigma list file.vtk
+\end_layout
+
+\begin_deeper
+\begin_layout LyX-Code
+Reading file.vtk
+\end_layout
+
+\begin_layout LyX-Code
+Points = 119827
+\end_layout
+
+\begin_layout LyX-Code
+Cells = 661929
+\end_layout
+
+\begin_layout LyX-Code
+PointDataArray[0] = displacements_t0 (119827 x 3)
+\end_layout
+
+\end_deeper
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout Section
+Cigma Datasets
+\end_layout
+
+\begin_layout Standard
+A field is a function which assigns a physical quantity to every point in
+ a domain 
+\begin_inset Formula $\Omega$
+\end_inset
+
+.
+ This quantity may correspond to a scalar, a vector, or even a tensor.
+ As in Chapter 
+\begin_inset LatexCommand ref
+reference "cha:Error-Analysis"
+
+\end_inset
+
+, we shall use the same discretization 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+A finite element approximation to an solution field 
+\begin_inset Formula $\phi(\vec{x})$
+\end_inset
+
+ on an element 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ is essentially specified by the weighed sum over a fixed set of local shape
+ functions 
+\begin_inset Formula $\phi_{e,1}(\vec{\xi}),\phi_{e,2}(\vec{\xi}),\ldots,\phi_{e,N}(\vec{\xi})$
+\end_inset
+
+ defined on the appropriate reference cell,
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\phi(\vec{x})=\sum_{n=1}^{N}d_{e,n}\phi_{e,n}(\vec{\xi})\]
+
+\end_inset
+
+where the global point 
+\begin_inset Formula $\vec{x}\in\Omega_{e}$
+\end_inset
+
+ corresponds to the local reference point 
+\begin_inset Formula $\vec{\xi}=\vec{x}_{e}^{-1}(\vec{x})\in\hat{\Omega}$
+\end_inset
+
+, and the weights 
+\begin_inset Formula $d_{e,n}$
+\end_inset
+
+, also known as degrees of freedom, are given on each global node.
+ These shape functions define a basis for the function space on the finite
+ element 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+For a mesh consisting of a single element type, we need to specify the 
+\begin_inset Formula $(x_{n},y_{n},z_{n})$
+\end_inset
+
+ of its degrees of freedom, and the connectivity relations 
+\begin_inset Formula $\Omega_{e}=\{n_{1},n_{2},\ldots\}$
+\end_inset
+
+ among them which define each individual element in the corresponding discretiza
+tion.
+\end_layout
+
+\begin_layout Standard
+The most common element types used in the finite element method are given
+ by the following shape functions, which define the interior of a tetrahedron
+ with vertices at 
+\begin_inset Formula $a=(0,0,0)$
+\end_inset
+
+, 
+\begin_inset Formula $b=(1,0,0)$
+\end_inset
+
+, 
+\begin_inset Formula $c=(0,1,0)$
+\end_inset
+
+, 
+\begin_inset Formula $d=(0,0,1)$
+\end_inset
+
+, 
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Standard
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "b"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status open
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+N_{a} & = & \frac{1}{2}(-1-x-y-z)\\
+N_{b} & = & \frac{1}{2}(1+x)\\
+N_{c} & = & \frac{1}{2}(1+y)\\
+N_{d} & = & \frac{1}{2}(1+z)\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "c"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/reference-tet4.png
+	width 5cm
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Caption
+
+\begin_layout Standard
+Reference Tetrahedron
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Likewise, a hexahedral element with vertices located at 
+\begin_inset Formula $a=(-1,-1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $b=(1,-1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $c=(1,1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $d=(-1,1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $e=(-1,-1,1)$
+\end_inset
+
+, 
+\begin_inset Formula $f=(1,-1,1)$
+\end_inset
+
+, 
+\begin_inset Formula $g=(1,1,1)$
+\end_inset
+
+, 
+\begin_inset Formula $h=(-1,1,1)$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Standard
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "c"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/reference-hex8.png
+	width 5cm
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "c"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status open
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+N_{a} & = & \frac{1}{8}\left(1-x\right)\left(1-y\right)\left(1-z\right)\\
+N_{b} & = & \frac{1}{8}\left(1+x\right)\left(1-y\right)\left(1-z\right)\\
+N_{c} & = & \frac{1}{8}\left(1-x\right)\left(1+y\right)\left(1-z\right)\\
+N_{d} & = & \frac{1}{8}\left(1+x\right)\left(1+y\right)\left(1-z\right)\\
+N_{e} & = & \frac{1}{8}\left(1-x\right)\left(1-y\right)\left(1+z\right)\\
+N_{f} & = & \frac{1}{8}\left(1+x\right)\left(1-y\right)\left(1+z\right)\\
+N_{g} & = & \frac{1}{8}\left(1-x\right)\left(1+y\right)\left(1+z\right)\\
+N_{h} & = & \frac{1}{8}\left(1+x\right)\left(1+y\right)\left(1+z\right)\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Caption
+
+\begin_layout Standard
+Reference Hexahedron
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+XXX: Take description of higher order elements from Karniadakis book
+\end_layout
+
+\begin_layout Subsection
+Comparing Two Finite Element Fields
+\end_layout
+
+\begin_layout Standard
+Comparing two arbitrary finite element fields can be accomplished with the
+ 
+\family typewriter
+cigma compare
+\family default
+ command-line utility.
+ By default, the comparison will involve a numerical integration over each
+ of the elements in the mesh associated with the first field.
+ If this discretization is inadequate, you may also specify an alternative
+ discretization over which to perform the integration.
+ The comparison operation will output the local residual values are into
+ the specified output file.
+\end_layout
+
+\begin_layout Standard
+A basic comparison can be as simple as specifying the following arguments:
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --first=field1.h5:/field1/stepN
+\end_layout
+
+\begin_layout LyX-Code
+
+\series bold
+ 
+\series default
+             --second=field2.h5:/field2/stepN
+\end_layout
+
+\begin_layout LyX-Code
+              --output=residuals.vtk
+\end_layout
+
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout Subsection
+Comparing against Known Values
+\end_layout
+
+\begin_layout Standard
+A finite element description might not always be available for one of the
+ fields.
+ However, you can break the comparison into several steps if you have a
+ means to compute that field on any of the required points.
+\end_layout
+
+\begin_layout Standard
+First, extract the global coordinates of the integration points.
+ This will result in an explicit list of points over which to evaluate your
+ field.
+\end_layout
+
+\begin_layout LyX-Code
+cigma extract field1.h5:/model/mesh/ -o points.h5:/projected_points
+\end_layout
+
+\begin_layout Standard
+Now you can evaluate your function at the designated points.
+ You can provide the path to the explicit set of values with
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --first=field1.h5:/field1/stepN
+\end_layout
+
+\begin_layout LyX-Code
+              --second-points=points.h5:/projected_points
+\end_layout
+
+\begin_layout LyX-Code
+              --second-values=values.h5:/projected_values
+\end_layout
+
+\begin_layout LyX-Code
+              --output=residuals.vtk
+\end_layout
+
+\begin_layout Subsection
+\begin_inset LatexCommand label
+name "sub:Comparing-against-a"
+
+\end_inset
+
+Comparing against a Known Function
+\end_layout
+
+\begin_layout Standard
+If one of your fields is easily described by an analytic expression, then
+ you also have the option to compile your analytic function as a builtin
+ Cigma function.
+ This will enable you to reference your function by name when using the
+ 
+\family typewriter
+compare
+\family default
+ command.
+ For example,
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --first=field1.h5:/vars/displacement/step0
+\end_layout
+
+\begin_layout LyX-Code
+
+\series bold
+ 
+\series default
+             --second=
+\bar under
+disloc3d
+\end_layout
+
+\begin_layout LyX-Code
+              --output=residuals.vtk
+\end_layout
+
+\begin_layout Standard
+You may also interact with your analytic function by using the 
+\family typewriter
+cigma
+\family default
+ 
+\family typewriter
+eval
+\family default
+ command, and obtain a set of values which may then be passed back to the
+ 
+\family typewriter
+compare
+\family default
+ command.
+ 
+\end_layout
+
+\begin_layout LyX-Code
+cigma eval --function=
+\bar under
+disloc3d
+\end_layout
+
+\begin_layout LyX-Code
+           --points=points.h5:/projected_points
+\end_layout
+
+\begin_layout LyX-Code
+           --values=values.h5:/disloc3d_values
+\end_layout
+
+\begin_layout Standard
+Once these values have been obtained, we can use it to resume the calculation
+ of 
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --output=residuals.vtk
+\end_layout
+
+\begin_layout LyX-Code
+              --first=field1.h5:/vars/displacement/step0
+\end_layout
+
+\begin_layout LyX-Code
+              --second-points=points.h5:/projected_points
+\end_layout
+
+\begin_layout LyX-Code
+              --second-values=values.h5:/disloc3d_values
+\end_layout
+
+\begin_layout Standard
+Another built-in function you might find useful is the 
+\family typewriter
+\bar under
+zero
+\family default
+\bar default
+ function, which is defined to return 0 everywhere.
+ You can use this function for estimating the norm of your field.
+\end_layout
+
+\begin_layout Subsection
+Specifying a Integration Mesh
+\end_layout
+
+\begin_layout Standard
+To override the mesh used in the integration, you can specify an extra argument
+ providing the location of the mesh,
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-mesh=mesh.h5:/model/mesh/
+\end_layout
+
+\begin_layout Standard
+Alternatively, you can also specify the coordinates and connectivity arrays
+ separately, in case they reside in separate files.
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-mesh-coordinates=file1.h5:/model/mesh/coordinates
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-mesh-connectivity=file2.h5:/model/mesh/connectivity
+\end_layout
+
+\begin_layout Subsection
+Specifying a Quadrature Rule
+\end_layout
+
+\begin_layout Standard
+To specify a quadrature rule, you will have to provide the quadrature weights
+ and points in the appropriate reference element.
+ This can be done with the following additional argument:
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-rule=quadrature-rules.h5:/path/to/rule
+\end_layout
+
+\begin_layout Standard
+You may also specify the location of the points and weights separately:
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-rule-points=file.h5:/path/to/rule/points
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-rule-weights=file.h5:/path/to/rule/weights
+\end_layout
+
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout Standard
+In the 
+\family typewriter
+src/
+\family default
+ directory of the Cigma distribution, you can find a Python script called
+ 
+\family typewriter
+rules.py
+\family default
+ that uses FIAT for calculating an appropriate set of quadrature rules for
+ use in Cigma.
+\end_layout
+
+\begin_layout Chapter
+Examples
+\end_layout
+
+\begin_layout Standard
+In this chapter we will use sample datasets from the strike-slip benchmark
+ case defined by the CIG Short-Term Tectonics working group.
+ In this benchmark problem, we solve for the viscoelastic relaxation of
+ stresses from a single finite earthquake while ignoring gravity.
+ The problem is defined on a cube domain with sides of 24 km consisting
+ of two layers of different material types.
+ The top layer of the cube is nearly elastic, while the bottom layer is
+ viscoelastic.
+ (XXX: Anything else from Pylith manual or source?)
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+wide false
+sideways false
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/strike-slip-fault.png
+	width 5cm
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Caption
+
+\begin_layout Standard
+Geometry of Strike-Slip Benchmark
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Section
+Convergence
+\end_layout
+
+\begin_layout Standard
+When the exact solution of your equations is known, the norm of the error
+ for a finite element solution is easily computed through the procedure
+ given in Chapter 
+\begin_inset LatexCommand ref
+reference "cha:Error-Analysis"
+
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+Another way of checking the accuracy of a solution is by rerunning the same
+ problem with a finer mesh.
+ Recall that in the finite element method, the weight functions and trial
+ solutions are chosen so that the finite element method is convergent.
+ In other words, as the element size 
+\begin_inset Formula $h$
+\end_inset
+
+ decreases, the solution tends to the correct solution.
+ Thus, if the difference between the coarse and fine mesh solutions is small,
+ we can conclude that the coarse mesh solution is quite accurate.
+ However, if a solution changes noticeably with refinement of the mesh,
+ the coarse mesh solution is inaccurate, and even the finer mesh may be
+ inadequate as well.
+\end_layout
+
+\begin_layout Standard
+As can be seen from these results, the log of the error varies linearly
+ with element size and the slope depends on the order of the element.
+ Denoting the slope by 
+\begin_inset Formula $\alpha$
+\end_inset
+
+, then the error in the function can be expressed as
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\log\varepsilon=\log C+\alpha\log h\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+where 
+\begin_inset Formula $\log C$
+\end_inset
+
+ is an arbitrary constant, the y-intercept of the curve.
+ The slope 
+\begin_inset Formula $\alpha$
+\end_inset
+
+ is the rate of convergence of the element.
+ We can rewrite this as
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\varepsilon=Ch^{\alpha}\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+XXX: list of examples to discuss here
+\end_layout
+
+\begin_layout Standard
+Viscoelastic strikeslip example (tet vs.
+ hex elts) -- pylith-tet / pylith-hex
+\end_layout
+
+\begin_layout Standard
+Cylinder extension for maxwell material (newtonian) -- pylith / geofest
+ / analytic
+\end_layout
+
+\begin_layout Standard
+Cylinder extension for maxwell material (non-newtonian) -- pylith / geofest
+ / analytic
+\end_layout
+
+\begin_layout Standard
+Cylinder gravitation -- pylith / geofest / analytic
+\end_layout
+
+\begin_layout Standard
+Laplace example (linear, quadratic elements) -- deal.II / libmesh / analytic
+\end_layout
+
+\begin_layout Standard
+Citcom example (structured datasets) -- citcom (3d) / analytic
+\end_layout
+
+\begin_layout Standard
+Gale example (structured datasets) -- gale (2d) / analytic
+\end_layout
+
+\begin_layout Standard
+MADDs example -- 
+\end_layout
+
+\begin_layout Section
+Visualization
+\end_layout
+
+\begin_layout Standard
+The residual field obtained in the 
+\end_layout
+
+\begin_layout Chapter
+\start_of_appendix
+File Formats
+\end_layout
+
+\begin_layout Standard
+To differentiate between these input and output file formats, you will need
+ to provide an appropriate file extension.
+ Currently recognized file extensions include 
+\family typewriter
+*.h5
+\family default
+ for HDF5 files, 
+\family typewriter
+*.vtk
+\family default
+ for legacy VTK files, and 
+\family typewriter
+*.dat
+\family default
+ for simple text files.
+\end_layout
+
+\begin_layout Section
+Input File Format
+\end_layout
+
+\begin_layout Standard
+There are five different objects you will want to use as input to Cigma,
+ all of which can be represented as two-dimensional datasets.
+ These are the 
+\end_layout
+
+\begin_layout Itemize
+mesh coordinates
+\end_layout
+
+\begin_layout Itemize
+mesh connectivity
+\end_layout
+
+\begin_layout Itemize
+field coefficients (degrees of freedom)
+\end_layout
+
+\begin_layout Itemize
+quadrature rule points
+\end_layout
+
+\begin_layout Itemize
+quadrature rule weights
+\end_layout
+
+\begin_layout Subsection
+HDF5
+\end_layout
+
+\begin_layout Standard
+HDF5 files are organized in a hierarchical structure, similar to a UNIX
+ file system.
+ Two types of primary objects, groups and datasets, are stored in this structure.
+ A group contains instances of zero or more groups or datasets, while a
+ dataset stores a multi-dimensional array of data elements.
+ Both kinds of objects are accompanied by supporting metadata known as attribute
+s.
+\end_layout
+
+\begin_layout Standard
+A dataset is physically stored in two parts: a header and a data array.
+ The header contains miscellaneous metadata describing the dataset as well
+ as information that is needed to interpret the array portion of the dataset.
+ Essentially, it includes the name, datatype, dataspace, and storage layout
+ of the dataset.
+ The name is a text string identifying the dataset.
+ The datatype describes the type of the data array elements.
+ The dataspace defines the dimensionality of the dataset, i.e., the size and
+ shape of the multi-dimensional array.
+\end_layout
+
+\begin_layout Standard
+In Cigma you always provide an explicit path to every dataset, so you have
+ a fair amount of flexibility for how you organize your datasets inside
+ HDF5 files.
+ For example, a typical Cigma HDF5 file could have the following structure.
+\end_layout
+
+\begin_layout LyX-Code
+model.h5
+\end_layout
+
+\begin_layout LyX-Code
+
+\backslash
+__ model
+\end_layout
+
+\begin_layout LyX-Code
+    |__ mesh
+\end_layout
+
+\begin_layout LyX-Code
+    |   |__ coordinates    [nno x nsd]
+\end_layout
+
+\begin_layout LyX-Code
+    |   
+\backslash
+__ connectivity   [nel x ndof]
+\end_layout
+
+\begin_layout LyX-Code
+    
+\backslash
+__ variables
+\end_layout
+
+\begin_layout LyX-Code
+        |__ temperature
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00000  [nno x 1]
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00010  [nno x 1]
+\end_layout
+
+\begin_layout LyX-Code
+        |   
+\backslash
+...
+\end_layout
+
+\begin_layout LyX-Code
+        |__ displacement
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00000  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00010  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+        |   
+\backslash
+...
+\end_layout
+
+\begin_layout LyX-Code
+        
+\backslash
+__ velocity
+\end_layout
+
+\begin_layout LyX-Code
+            |__ step00000  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+            |__ step00010  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+            
+\backslash
+...
+\end_layout
+
+\begin_layout Standard
+Generally, Cigma will only require you to specify the path to a specific
+ field dataset.
+ If a small amount of metadata is present in your field dataset, the rest
+ of the required information, such as which mesh and finite elements to
+ use, will be deduced from that metadata.
+\end_layout
+
+\begin_layout Description
+MeshID an identifier assigned for use in linking child datasets to their
+ parent mesh.
+\end_layout
+
+\begin_layout Description
+MeshLocation points to the HDF5 group which contains the appropriate coordinates
+ and connectivity datasets.
+\end_layout
+
+\begin_layout Description
+FunctionSpace string identifier to determine which shape functions to use
+ for interpolating values inside the element (e.g., tet4, hex8, quad4, tri3,
+ ...).
+\end_layout
+
+\begin_layout Description
+DatasetType string identifier for classifying the type of data contained
+ in the dataset (e.g., points, connectivity, degrees of freedom, quadrature
+ rules, global quadrature points, global field values).
+\end_layout
+
+\begin_layout Subsection
+VTK
+\end_layout
+
+\begin_layout Standard
+The current version of Cigma only supports VTK unstructured grid datasets
+ of a single element type.
+ You can still compare various unstructured grids with different element
+ types to each other.
+\end_layout
+
+\begin_layout Standard
+Note that while you typically provide a path (or name) for every dataset,
+ this is not necessary when specifying a VTK mesh, since this data is taken
+ from the special Points and Cells arrays, which you cannot rename.
+ However, you will still need to provide a name when referring to the field
+ coefficients, which are assumed to be stored as Point Data in the input
+ VTK file.
+\end_layout
+
+\begin_layout Standard
+For more information, you may want to refer to Visualization ToolKit's 
+\begin_inset LatexCommand htmlurl
+name "File Formats"
+target "www.vtk.org/pdf/file-formats.pdf"
+
+\end_inset
+
+ document.
+ 
+\end_layout
+
+\begin_layout Subsection
+Text
+\end_layout
+
+\begin_layout Standard
+The text format is always in table form, with the dimensions of the table
+ specified in the first line.
+ For example, mesh coordinates can be specified in the following format
+\end_layout
+
+\begin_layout LyX-Code
+<nno> <nsd>
+\end_layout
+
+\begin_layout LyX-Code
+x1 y1 z1
+\end_layout
+
+\begin_layout LyX-Code
+x2 y2 z2
+\end_layout
+
+\begin_layout LyX-Code
+x3 y3 z3
+\end_layout
+
+\begin_layout LyX-Code
+...
+\end_layout
+
+\begin_layout Standard
+Mesh connectivity with 
+\end_layout
+
+\begin_layout LyX-Code
+nel ndofs
+\end_layout
+
+\begin_layout LyX-Code
+node_1 node_2 node_3 node_4 ...
+\end_layout
+
+\begin_layout LyX-Code
+node_1 node_2 node_3 node_4 ...
+\end_layout
+
+\begin_layout LyX-Code
+node_1 node_2 node_3 node_4 ...
+\end_layout
+
+\begin_layout LyX-Code
+...
+\end_layout
+
+\begin_layout Standard
+A generic field with 
+\family typewriter
+ndim
+\family default
+ components (i.e., scalar, vector, or tensor) is specified by
+\end_layout
+
+\begin_layout LyX-Code
+num ndim
+\end_layout
+
+\begin_layout LyX-Code
+f1 f2 f3 ..
+\end_layout
+
+\begin_layout LyX-Code
+f1 f2 f3 ..
+\end_layout
+
+\begin_layout LyX-Code
+...
+\end_layout
+
+\begin_layout Standard
+In this last case, the number of rows could refer to either 
+\family typewriter
+nno
+\family default
+ or 
+\family typewriter
+nel
+\family default
+, depending on whether the field is node-based or cell-based.
+\end_layout
+
+\begin_layout Section
+Output File Format
+\end_layout
+
+\begin_layout Standard
+The output format for residuals consists simply of a list of scalars for
+ each element in the discretization.
+ If you specify a 
+\family typewriter
+.vtk
+\family default
+ extension for your output file, this will result in an scalar dataset named
+ epsilon stored using the legacy VTK file format in the Cell Data section
+ of the output file.
+ Note that this output consists of the squared values of the local residuals,
+ so further post-processing will be necessary.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "1"
+key "Akin 2005"
+
+\end_inset
+
+Akin, J.E.
+ (2005), 
+\emph on
+Finite Element Analysis with Error Estimators,
+\emph default
+ Butterworth-Heinemann, Oxford, 447 pp.
+ 
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "2"
+key "Encyclopaedia of Cubature Formulas 2005"
+
+\end_inset
+
+Encyclopaedia of Cubature Formulas (1998-2005), Katholieke Universiteit
+ Leuven, Dept.
+ Computerwetenschappen, www.cs.kuleuven.be/~nines/research/ecf/
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "3"
+key "Hughes 2000"
+
+\end_inset
+
+Hughes, Thomas J.R.
+ (2000), 
+\emph on
+The Finite Element Method: Linear Static and Dynamic Finite Element Analysis
+\emph default
+, Dover, 672 pp.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "4"
+key "KarniadakisSherwin2005"
+
+\end_inset
+
+Karniadakis, George E.
+ and Spencer J.
+ Sherwin (2005), 
+\emph on
+Spectral/
+\emph default
+hp 
+\emph on
+Element Methods for Computational Fluid Dynamics, Second Edition
+\emph default
+, Numerical Mathematics and Scientific Computation, Oxford, 657 pp.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "5"
+key "Uesu-etal2005"
+
+\end_inset
+
+Uesu, D., L.
+ Bavoil, S.
+ Fleishman, J.
+ Shepherd, and C.
+ Silva (2005), Simplification of Unstructured Tetrahedral Meshes by Point-Sampli
+ng, 
+\emph on
+Proceedings of the 2005 International Workshop on Volume Graphics,
+\emph default
+ 157-165, doi: 10.2312/VG/VG05/157-165.
+\end_layout
+
+\end_body
+\end_document

Copied: doc/cigma/manual/cigma.lyx~ (from rev 11788, doc/cigma/manual/cigma.lyx)
===================================================================
--- doc/cigma/manual/cigma.lyx~	                        (rev 0)
+++ doc/cigma/manual/cigma.lyx~	2008-07-31 17:58:23 UTC (rev 12498)
@@ -0,0 +1,2505 @@
+#LyX 1.5.4 created this file. For more info see http://www.lyx.org/
+\lyxformat 276
+\begin_document
+\begin_header
+\textclass book
+\begin_preamble
+\usepackage{hyperref}
+
+\let\myUrl\url
+\renewcommand{\url}[1]{(\myUrl{#1})}
+\end_preamble
+\language english
+\inputencoding auto
+\font_roman default
+\font_sans default
+\font_typewriter default
+\font_default_family default
+\font_sc false
+\font_osf false
+\font_sf_scale 100
+\font_tt_scale 100
+\graphics default
+\paperfontsize default
+\spacing single
+\papersize default
+\use_geometry true
+\use_amsmath 1
+\use_esint 0
+\cite_engine basic
+\use_bibtopic false
+\paperorientation portrait
+\leftmargin 1in
+\topmargin 1in
+\rightmargin 1in
+\bottommargin 1in
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\defskip medskip
+\quotes_language swedish
+\papercolumns 1
+\papersides 2
+\paperpagestyle headings
+\tracking_changes false
+\output_changes false
+\author "" 
+\author "" 
+\end_header
+
+\begin_body
+
+\begin_layout Standard
+\noindent
+\align center
+
+\color none
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+thispagestyle{empty}
+\end_layout
+
+\end_inset
+
+
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/cigma_cover.pdf
+	display color
+	width 75page%
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\noindent
+\align center
+
+\color none
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+thispagestyle{empty}
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Title
+Cigma
+\end_layout
+
+\begin_layout Author
+© California Institute of Technology
+\newline
+Version 1.0.0
+\end_layout
+
+\begin_layout Date
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+today
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset LatexCommand tableofcontents
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset FloatList figure
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+raggedbottom
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+
+\newpage
+
+\end_layout
+
+\begin_layout Chapter
+Introduction
+\end_layout
+
+\begin_layout Section
+About Cigma
+\end_layout
+
+\begin_layout Standard
+The CIG Model Analyzer (Cigma) consists of a suite of tools for comparing
+ numerical models.
+ CIG has developed Cigma in response to demand from the short-term tectonics
+ community for a simple tool that can perform rigorous error analysis on
+ their finite element codes.
+ The long-term goal for Cigma, however, is for it to be used for nearly
+ all geodynamics modeling codes.
+\end_layout
+
+\begin_layout Standard
+The current version of Cigma is intended for the calculation of 
+\begin_inset Formula $L_{2}$
+\end_inset
+
+ residuals for finite element models.
+ This software was designed for performing the following three main tasks:
+ (1) error analysis, (2) benchmarking, and (3) code verification.
+ In error analysis, Cigma can calculate a global measure of the error between
+ two solutions by performing an integration over a discretized version of
+ the common domain.
+ This comparison can take place even when the underlying discretizations
+ do not overlap.
+ In benchmarking, Cigma can help the geodynamics community agree on a standard
+ solution to specific problems by facilitating the process of comparing
+ different numerical codes against each other.
+ Lastly, as an automated tool, Cigma can help application developers create
+ regression tests to ensure that software changes do not affect the consistency
+ of the results.
+\end_layout
+
+\begin_layout Standard
+At its core, Cigma draws from a variety of libraries.
+ The 
+\begin_inset LatexCommand htmlurl
+name "FInite element Automatic Tabulator (FIAT)"
+target "www.fenics.org/wiki/FIAT"
+
+\end_inset
+
+ Python Library supports generation of arbitrary order instances of the
+ Lagrange elements on lines, triangles, and tetrahedra.
+ It can also generate higher order instances of Jacobi-type quadrature rules
+ on the same element shapes.
+ The 
+\begin_inset LatexCommand htmlurl
+name "Approximate Nearest Neighbor (ANN)"
+target "www.cs.umd.edu/~mount/ANN/"
+
+\end_inset
+
+ Library, written in C++, supports data structures and algorithms for both
+ exact and nearest-neighbor searching in arbitrarily high dimensions.
+\end_layout
+
+\begin_layout Standard
+Both of these libraries extend and generalize Cigma's functionality so it
+ can handle other types of elements, and provide the ability to compare
+ vector fields.
+\end_layout
+
+\begin_layout Section
+Citation
+\end_layout
+
+\begin_layout Standard
+Computational Infrastructure for Geodynamics (CIG) is making this source
+ code available to you in the hope that the software will enhance your research
+ in geophysics.
+ This is a brand-new code and at present no papers are published or press.
+ Please cite this manual as follows:
+\end_layout
+
+\begin_layout Itemize
+Armendariz, L., and S.
+ Kientz.
+ 
+\emph on
+Cigma User Manual.
+
+\emph default
+ Pasadena, CA: Computational Infrastructure of Geodynamics, 2008.
+ URL: geodynamics.org/cig/software/cs/cigma/cigma.pdf
+\end_layout
+
+\begin_layout Standard
+CIG requests that in your oral presentations and in your papers that you
+ indicate your use of this code and acknowledge the author of the code and
+ 
+\begin_inset LatexCommand htmlurl
+name "CIG"
+target "geodynamics.org"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Section
+Support
+\end_layout
+
+\begin_layout Standard
+Cigma development is based upon work supported by the National Science Foundatio
+n under Grant No.
+ EAR-0406751.
+ Any opinions, findings, and conclusions or recommendations expressed in
+ this material are those of the authors and do not necessarily reflect the
+ views of the National Science Foundation.
+ The code is being released under the GNU General Public License.
+\end_layout
+
+\begin_layout Chapter
+Installation and Getting Help
+\end_layout
+
+\begin_layout Section
+Getting Help
+\end_layout
+
+\begin_layout Standard
+For help, send an e-mail to the 
+\begin_inset LatexCommand url
+name "CIG Computational Science Mailing List"
+target "cig-cs at geodynamics.org"
+
+\end_inset
+
+.
+ You can subscribe to the 
+\family typewriter
+cig-cs
+\family default
+ mailing list and view archived discussions at the 
+\begin_inset LatexCommand htmlurl
+name "CIG Mail Lists web page"
+target "geodynamics.org/cig/lists"
+
+\end_inset
+
+.
+ If you encounter any bugs or have problems installing Cigma, please submit
+ a report to the 
+\begin_inset LatexCommand htmlurl
+name "CIG Bug Tracker"
+target "geodynamics.org/bugs"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Section
+Installating from Source
+\end_layout
+
+\begin_layout Standard
+To install Cigma, download the source package from the 
+\begin_inset LatexCommand htmlurl
+name "CIG Cigma web page"
+target "geodynamics.org/cig/software/packages/cs/cigma"
+
+\end_inset
+
+.
+ You will need the GNU C++ compiler for this step.
+ The Boost, HDF5, and VTK libraries, described later in this section, are
+ also required.
+\end_layout
+
+\begin_layout Standard
+After installing the necessary dependencies, simply unpack the source tar
+ file and issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz cigma-1.0.0.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd cigma-1.0.0
+\end_layout
+
+\begin_layout LyX-Code
+$ ./configure --with-boost=$BOOST_PREFIX
+\end_layout
+
+\begin_layout LyX-Code
+              --with-hdf5=$HDF5_PREFIX 
+\end_layout
+
+\begin_layout LyX-Code
+              --with-vtk=$VTK_PREFIX 
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+Both HDF5 and VTK prefixes are required if you have installed these libraries
+ in a custom location.
+\end_layout
+
+\begin_layout Subsection
+Boost
+\end_layout
+
+\begin_layout Standard
+The Boost C++ libraries are used in Cigma for parsing command-line options,
+ and ensuring cross-platform access to the filesystem.
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz boost_1_35_0.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd boost_1_35_0
+\end_layout
+
+\begin_layout LyX-Code
+$ ./configure
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+Without a custom installation prefix, Boost will target your 
+\family typewriter
+/usr/local
+\family default
+ directory.
+\end_layout
+
+\begin_layout Subsection
+HDF5
+\end_layout
+
+\begin_layout Standard
+The Hierarchical Data Format (HDF5) library is available for download from
+ 
+\begin_inset LatexCommand htmlurl
+name "The HDF Group"
+target "hdfgroup.org/HDF5"
+
+\end_inset
+
+.
+ Binaries can be obtained at 
+\begin_inset LatexCommand htmlurl
+name "hdfgroup.org/HDF5/release/obtain5.html"
+target "hdfgroup.org/HDF5/release/obtain5.html"
+
+\end_inset
+
+.
+ To install from source, download the latest stable 1.6 version of this library
+ (1.6.7, for example) and issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz hdf5-1.6.7.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd hdf5-1.6.7
+\end_layout
+
+\begin_layout LyX-Code
+$ ./configure --prefix=$HOME/opt/hdf5/1.6.7
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+Without the custom prefix, this installation procedure will target your
+ 
+\family typewriter
+/usr
+\family default
+ directory.
+\end_layout
+
+\begin_layout Subsection
+VTK
+\end_layout
+
+\begin_layout Standard
+The Visualization Toolkit (VTK) library is a popular open source graphics
+ library for scientific visualizations.
+ We recommend that you obtain the binaries from a package manager, making
+ sure to install the associated development headers along with the library.
+ The source for the VTK library is available from 
+\begin_inset LatexCommand htmlurl
+name "Kitware, Inc."
+target "www.vtk.org/get-software.php"
+
+\end_inset
+
+.
+ If you decide to compile VTK from source, you will also need the CMake
+ build environment, available from 
+\begin_inset LatexCommand htmlurl
+name "CMake"
+target "cmake.org"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+After downloading the source package, you can issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz vtk-5.0.4.tar.gz
+\end_layout
+
+\begin_layout LyX-Code
+$ cd VTK
+\end_layout
+
+\begin_layout LyX-Code
+$ mkdir build
+\end_layout
+
+\begin_layout LyX-Code
+$ cd build
+\end_layout
+
+\begin_layout LyX-Code
+$ ccmake ..
+\end_layout
+
+\begin_layout LyX-Code
+$ make
+\end_layout
+
+\begin_layout LyX-Code
+$ make install
+\end_layout
+
+\begin_layout Standard
+You can set the installation prefix during the 
+\family typewriter
+ccmake
+\family default
+ step above, which defaults to your 
+\family typewriter
+/usr/local
+\family default
+ directory.
+\end_layout
+
+\begin_layout Section
+Additional Software
+\end_layout
+
+\begin_layout Subsection
+NumPy
+\end_layout
+
+\begin_layout Standard
+NumPy is a Python module which adds support for large multi-dimensional
+ arrays and matrices, and is available for download from the 
+\begin_inset LatexCommand htmlurl
+name "NumPy Home Page"
+target "numpy.scipy.org "
+
+\end_inset
+
+.
+ To install this extension module from source, download it and issue the
+ following command in the NumPy source directory
+\end_layout
+
+\begin_layout LyX-Code
+$ python setup.py install
+\end_layout
+
+\begin_layout Standard
+To verify the installation, the command `
+\family typewriter
+import numpy
+\family default
+' should run successfully from within your standard Python interpreter shell.
+\end_layout
+
+\begin_layout Subsection
+PyTables
+\end_layout
+
+\begin_layout Standard
+PyTables is a Python extension module that builds on top of the HDF5 library.
+ It provides a convenient scripting interface to manipulate HDF5 files,
+ which can be used to manipulate the input/output files created by Cigma.
+ PyTables is available for download from 
+\begin_inset LatexCommand htmlurl
+name "PyTables"
+target "www.pytables.org"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+To install this extension from source, download the latest stable version
+ (currently 2.0.2) and issue the following commands
+\end_layout
+
+\begin_layout LyX-Code
+$ tar xvfz pytables-2.0.2
+\end_layout
+
+\begin_layout LyX-Code
+$ cd pytables-2.0.2
+\end_layout
+
+\begin_layout LyX-Code
+$ python setup.py install
+\end_layout
+
+\begin_layout Standard
+To verify the installation, the command `
+\family typewriter
+import tables
+\family default
+' should run successfully from within your standard Python interpreter shell.
+\end_layout
+
+\begin_layout Subsection
+FIAT
+\end_layout
+
+\begin_layout Standard
+The Finite Element Automatic Tabulator (FIAT) Python module supports generation
+ of arbitrary order instances of the Lagrange elements, as well as arbitrary
+ order instances of the Jacobi-type quadrature rules on the same element
+ shapes.
+ It is available from 
+\begin_inset LatexCommand htmlurl
+name "The FEniCS Project"
+target "www.fenics.org/wiki/FIAT"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+To install this module, download it and issue the following command inside
+ the FIAT source directory:
+\end_layout
+
+\begin_layout LyX-Code
+$ sudo python setup.py install
+\end_layout
+
+\begin_layout Standard
+You should now be able to run `
+\family typewriter
+import FIAT
+\family default
+' from within your standard Python interpreter shell.
+\end_layout
+
+\begin_layout Subsection
+HDFView
+\end_layout
+
+\begin_layout Standard
+NCSA HDFView is a graphical user interface tool for accessing data in your
+ HDF5 files.
+ You can use it for viewing the internal group hierarchy, adding new datasets,
+ and modifying or deleting existing datasets, or altering metadata on groups
+ and datasets.
+ You can download it from the 
+\begin_inset LatexCommand htmlurl
+name "HDFView home page"
+target "hdf.ncsa.uiuc.edu/hdf-java-html/hdfview"
+
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Chapter
+\begin_inset LatexCommand label
+name "cha:Error-Analysis"
+
+\end_inset
+
+Error Analysis
+\end_layout
+
+\begin_layout Section
+Introduction
+\end_layout
+
+\begin_layout Standard
+When studying differential equations that represent physical systems we
+ often obtain solutions by using a variety of techniques.
+ Solving these equations using classical analytical methods is almost impossible
+, so we often resort to numerical algorithms such as the finite element
+ method.
+\end_layout
+
+\begin_layout Standard
+For assessing the quality of the solution to our equations, it is important
+ to quantify the error in our numerical solutions.
+ To calculate this error, we will have to compare these solutions against
+ each other through the use of a distance measure between two functions.
+\end_layout
+
+\begin_layout Standard
+The simplest possible quantitative measure of the difference between two
+ distinct functions consists of taking the pointwise difference at a common
+ set of points.
+ While no finite sample of points can perfectly represent a continuum of
+ values, valuable information can be inferred from the statistics of the
+ resulting set of differences.
+ However, the fields we want to compare may not be defined at a common set
+ of points.
+ This simple measure becomes unapplicable.
+\end_layout
+
+\begin_layout Section
+Interpolation Functions
+\end_layout
+
+\begin_layout Section
+Distance Measures
+\end_layout
+
+\begin_layout Standard
+A very useful distance measure we can use is the 
+\begin_inset Formula $L_{2}$
+\end_inset
+
+ norm, defined by the following integral
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\varepsilon=||u-v||_{L_{2}}=\sqrt{\int_{\Omega}||u(\vec{x})-v(\vec{x})||^{2}d\vec{x}}\]
+
+\end_inset
+
+where 
+\begin_inset Formula $u$
+\end_inset
+
+ and 
+\begin_inset Formula $v$
+\end_inset
+
+ are the two functions on a global coordinate system.
+ This gives us a single global estimate 
+\begin_inset Formula $\varepsilon$
+\end_inset
+
+ representing the distance between the two fields 
+\begin_inset Formula $u(\vec{x})$
+\end_inset
+
+ and 
+\begin_inset Formula $v(\vec{x})$
+\end_inset
+
+.
+ Alternatively, you may think of this as the size, or norm, of the residual
+ field 
+\begin_inset Formula $\rho(\vec{x})=u(\vec{x})-v(\vec{x})$
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Section
+Global Error Approximation
+\end_layout
+
+\begin_layout Standard
+If we discretize the domain 
+\begin_inset Formula $\Omega$
+\end_inset
+
+ into an appropriate set of cells 
+\begin_inset Formula $\Omega_{1},\Omega_{2},\ldots,\Omega_{n_{el}}$
+\end_inset
+
+, we can compute the global error 
+\begin_inset Formula $\varepsilon^{2}$
+\end_inset
+
+ as a sum over localized cell contributions, 
+\begin_inset Formula $\varepsilon^{2}=\sum_{e}\varepsilon_{e}^{2}$
+\end_inset
+
+, where each cell residual 
+\begin_inset Formula $\varepsilon_{e}^{2}$
+\end_inset
+
+ is given by
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon_{e}^{2} & = & \int_{\Omega_{e}}||u(\vec{x})-v(\vec{x})||^{2}d\vec{x}\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+In general, we won't be able to integrate each cell residual 
+\begin_inset Formula $\varepsilon_{e}^{2}$
+\end_inset
+
+ exactly since either of the functions 
+\begin_inset Formula $u$
+\end_inset
+
+ and 
+\begin_inset Formula $v$
+\end_inset
+
+ may have an incompatible representation relative to the finite element
+ space on each domain 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+.
+ However, we can still calculate an approximation of each cell residual
+ by applying an appropriate quadrature rule with a tolerable truncation
+ error 
+\begin_inset LatexCommand cite
+key "Encyclopaedia of Cubature Formulas 2005"
+
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+To obtain an approximation to the integral of a function 
+\begin_inset Formula $F(\vec{x})$
+\end_inset
+
+ over a cell 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+, we simply apply the quadrature rule with weights 
+\begin_inset Formula $w_{e,1},w_{e,2},\ldots,w_{e,n_{Q}}$
+\end_inset
+
+ and integration points 
+\begin_inset Formula $\vec{x}_{e,1},\vec{x}_{e,2},\ldots,\vec{x}_{e,n_{Q}}$
+\end_inset
+
+ appropriate for the physical element 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\int_{\Omega_{e}}\ F(\vec{x})\ d\vec{x}=\sum_{q=1}^{n_{Q}}w_{e,q}F(\vec{x}_{e,q})\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Applying this quadrature rule directly over the entire physical domain 
+\begin_inset Formula $\Omega=\cup\Omega_{e}$
+\end_inset
+
+ gives us 
+\begin_inset Formula \[
+\int_{\Omega}\ F(\vec{x})\ d\vec{x}=\sum_{e=1}^{n_{el}}\sum_{q=1}^{n_{Q}}w_{e,q}F(\vec{x}_{e,q})\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+For efficiency reasons, it is undesirable in finite element applications
+ to perform calculations in a global coordinate system.
+ To avoid duplication of work, shape function evaluations may be performed
+ once on a reference cell 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+ and then transformed back into the corresponding physical cell 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ as needed.
+\end_layout
+
+\begin_layout Standard
+To compute integrals of 
+\begin_inset Formula $F$
+\end_inset
+
+ in a reference coordinate system, we need to apply a change of variables:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\int_{\Omega_{e}}\ F(\vec{x})\ d\vec{x}=\int_{\hat{\Omega}}\ F(\vec{x}_{e}(\vec{\xi}))J_{e}(\vec{\xi})\ d\vec{\xi}\]
+
+\end_inset
+
+where 
+\color none
+the additional factor 
+\begin_inset Formula $J_{e}(\vec{\xi})=\det\left[\frac{d\vec{x}_{e}}{d\vec{\xi}}\right]$
+\end_inset
+
+ is the Jacobian determinant of the reference map 
+\begin_inset Formula $\vec{x}_{e}(\vec{\xi}):\hat{\Omega}\to\Omega_{e}$
+\end_inset
+
+.
+ This map describes how the physical points 
+\begin_inset Formula $\vec{x}\in\Omega_{e}$
+\end_inset
+
+ are transformed from the reference points 
+\begin_inset Formula $\vec{\xi}\in\hat{\Omega}$
+\end_inset
+
+.
+ Put another way, the inverse reference map 
+\color inherit
+
+\begin_inset Formula $\vec{\xi}=\vec{x}_{e}^{-1}(\vec{x})$
+\end_inset
+
+ tells us how the physical domain 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ maps into the reference cell 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+At this point, we can assume without loss of generality that every physical
+ cell 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ can be derived from a single reference cell 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+, so that our quadrature rule becomes simply a set of weights 
+\begin_inset Formula $w_{1},\ldots,w_{n_{Q}}$
+\end_inset
+
+ and points 
+\begin_inset Formula $\vec{\xi}_{1},\ldots,\vec{\xi}_{n_{Q}}$
+\end_inset
+
+ over 
+\begin_inset Formula $\hat{\Omega}$
+\end_inset
+
+.
+ After changing variables, we end up with the final form of the quadrature
+ rule that we can use to integrate the global residual field:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\int_{\hat{\Omega}}\ F(\vec{x}_{e}(\vec{\xi}))J_{e}(\vec{\xi})\ d\vec{\xi}=\sum_{q=1}^{n_{Q}}w_{q}F(\vec{x}_{e}(\vec{\xi}_{q}))J_{e}(\vec{\xi}_{q})\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+If we let 
+\begin_inset Formula $F(\vec{x})=||u(\vec{x})-v(\vec{x})||^{2}$
+\end_inset
+
+ in the above expression, we can find the cell residual contribution over
+ 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ from
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon_{e}^{2} & = & \int_{\Omega_{e}}\ ||u(\vec{x})-v(\vec{x})||^{2}\ d\vec{x}\\
+ & = & \int_{\hat{\Omega}}\ ||u(\vec{x}_{e}(\vec{\xi}))-v(\vec{x}_{e}(\vec{\xi}))||^{2}J_{e}(\vec{\xi})\ d\vec{\xi}\\
+ & = & \sum_{q=1}^{\mathrm{n_{Q}}}w_{q}||u(\vec{x}_{e}(\vec{\xi}_{q}))-v(\vec{x}_{e}(\vec{\xi}_{q}))||^{2}J_{e}(\vec{\xi}_{q})\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+The global error 
+\begin_inset Formula $\varepsilon=\sqrt{\sum_{e}\varepsilon_{e}^{2}}$
+\end_inset
+
+ is then approximated by the expression
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon & = & \sqrt{\sum_{e=1}^{\mathrm{n_{el}}}\sum_{q=1}^{\mathrm{n_{Q}}}w_{q}||u(\vec{x}_{e}(\vec{\xi}_{q}))-v(\vec{x}_{e}(\vec{\xi}_{q}))||^{2}J_{e}(\vec{\xi}_{q})}\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+We can make a number of observations based on our final expression for the
+ global error, (XXX: add observations here).
+\end_layout
+
+\begin_layout Section
+Comparing Residuals
+\end_layout
+
+\begin_layout Standard
+For comparing errors of different solutions, you can normalize the global
+ error by the norm of the exact solution.
+ For a family of solutions 
+\begin_inset Formula $u^{h}(\vec{x})$
+\end_inset
+
+, parametrized by the the maximum element size 
+\begin_inset Formula $h$
+\end_inset
+
+ of the underlying discretization, and an exact solution 
+\begin_inset Formula $u(\vec{x})$
+\end_inset
+
+, this normalized error is given by
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+\varepsilon_{rel} & = & \frac{||u-u^{h}||_{L_{2}}}{||u||_{L_{2}}}\\
+ & = & \frac{\int_{\Omega}||u(\vec{x})-u^{h}(\vec{x})||^{2}d\vec{x}}{\int_{\Omega}||u(\vec{x})||^{2}d\vec{x}}\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+This normalized error can be interpreted as the average error in the physical
+ quantity being evaluated, so that a value of 0.01 corresponds to a 1% averaged
+ error.
+\end_layout
+
+\begin_layout Standard
+Even if the exact solution is not currently known, this normalized error
+ may be used to test the accuracy between two or more numerical solutions,
+ defined on successively refined meshes.
+\end_layout
+
+\begin_layout Chapter
+\begin_inset LatexCommand label
+name "cha:Running-Cigma"
+
+\end_inset
+
+Running Cigma
+\end_layout
+
+\begin_layout Standard
+Cigma is primarily designed for calculating error estimates between arbitrary
+ fields, so its primary usage is centered on the following operation
+\end_layout
+
+\begin_layout LyX-Code
+$ 
+\series bold
+cigma compare 
+\emph on
+\bar under
+field1
+\emph default
+\bar default
+ 
+\emph on
+\bar under
+field2
+\emph default
+\bar default
+ -o residuals.vtk
+\end_layout
+
+\begin_layout Standard
+You will need to provide two datasets describing each of the two fields,
+ specify an integration rule and a domain discretization over which to integrate
+, although these last two will have reasonable defaults if they are not
+ specified.
+\end_layout
+
+\begin_layout Standard
+The default data storage format for Cigma is the Hierarchical Data Format
+ (HDF5), a portable file format developed at the 
+\begin_inset LatexCommand htmlurl
+name "National Center for Supercomputing Applications (NCSA)"
+target "hdf.ncsa.uiuc.edu/HDF5"
+
+\end_inset
+
+.
+ The HDF5 is designed for storing multi-dimensional arrays together with
+ meta-data in a portable self-describing format.
+\end_layout
+
+\begin_layout Section
+Command Line Interface
+\end_layout
+
+\begin_layout Standard
+Cigma is designed to be scriptable.
+ Thus, all operations can be specified as command-line arguments given to
+ a single executable called 
+\family typewriter
+cigma
+\family default
+.
+ A set of available commands can be obtained by typing
+\end_layout
+
+\begin_layout LyX-Code
+
+\family typewriter
+$ 
+\family default
+\series bold
+cigma help
+\end_layout
+
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout LyX-Code
+Usage: cigma <subcommand> [options]
+\end_layout
+
+\begin_layout LyX-Code
+Type 'cigma help <subcommand>' for help on a specific subcommand.
+\end_layout
+
+\begin_layout LyX-Code
+Available subcommands
+\end_layout
+
+\begin_layout LyX-Code
+   help
+\end_layout
+
+\begin_layout LyX-Code
+   list
+\end_layout
+
+\begin_layout LyX-Code
+   extract
+\end_layout
+
+\begin_layout LyX-Code
+   eval
+\end_layout
+
+\begin_layout LyX-Code
+   compare 
+\end_layout
+
+\begin_layout Standard
+Examples for each of these commands are given in the sections below.
+\end_layout
+
+\begin_layout Section
+Input and Output Formats
+\end_layout
+
+\begin_layout Standard
+The underlying data storage format for Cigma is the HDF5 format, due to
+ its flexibility for storing and organizing large amounts of data.
+ The Hierarchical Data Format (HDF) is designed for storing, retrieving,
+ analyzing, visualizing, and converting scientific data.
+ It uses a hierarchical structure that provides users a host of options
+ for organizing how their data is stored in HDF5 files.
+ Using HDF5 datasets in Cigma allows us to avoid having to convert between
+ too many distinct formats.
+ Moreover, due to the amount of disk I/O, large finite element meshes can
+ be handled more efficiently in binary format.
+ 
+\end_layout
+
+\begin_layout Standard
+Another popular format for providing mesh and field inputs.You can easily
+ examine the structure of an input file by using the 
+\family typewriter
+cigma list
+\family default
+ command, which will simply reveal the names and dimensions of all datasets
+ inside the specified file.
+\end_layout
+
+\begin_layout Standard
+Inputs can be either HDF5 or VTK.
+ Specifying the complete path to a dataset consists of the special form
+ 
+\family typewriter
+\series bold
+filepath:dataset
+\family default
+\series default
+, a colon-delimited pair of file path and dataset path.
+\end_layout
+
+\begin_layout Standard
+Mesh inputs are currently only unstructured grids.
+\end_layout
+
+\begin_layout Standard
+Only exception are the residuals, which are written in legacy VTK format
+ as scalar cell data over an unstructured grid.
+\end_layout
+
+\begin_layout Standard
+XXX New section?
+\end_layout
+
+\begin_layout Standard
+Because Cigma relies on the ability of the user to specify dataset paths,
+ we have provided a command called 
+\family typewriter
+list
+\family default
+ for viewing the structure of an input file.
+ Its usage is very simple.
+ 
+\end_layout
+
+\begin_layout Standard
+To view the structure of an HDF5 file:
+\end_layout
+
+\begin_layout LyX-Code
+$ 
+\series bold
+cigma list file.h5
+\end_layout
+
+\begin_deeper
+\begin_layout LyX-Code
+/mesh/coordinates           Dataset {119827, 3}
+\end_layout
+
+\begin_layout LyX-Code
+/mesh/connectivity          Dataset {661929, 4}
+\end_layout
+
+\begin_layout LyX-Code
+/vars/displacement/step0    Dataset {119827, 3}
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+You can also view the structure of a VTK file with this command:
+\end_layout
+
+\begin_layout LyX-Code
+$ 
+\series bold
+cigma list file.vtk
+\end_layout
+
+\begin_deeper
+\begin_layout LyX-Code
+Reading file.vtk
+\end_layout
+
+\begin_layout LyX-Code
+Points = 119827
+\end_layout
+
+\begin_layout LyX-Code
+Cells = 661929
+\end_layout
+
+\begin_layout LyX-Code
+PointDataArray[0] = displacements_t0 (119827 x 3)
+\end_layout
+
+\end_deeper
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout Section
+Cigma Datasets
+\end_layout
+
+\begin_layout Standard
+A field is a function which assigns a physical quantity to every point in
+ a domain 
+\begin_inset Formula $\Omega$
+\end_inset
+
+.
+ This quantity may correspond to a scalar, a vector, or even a tensor.
+ As in Chapter 
+\begin_inset LatexCommand ref
+reference "cha:Error-Analysis"
+
+\end_inset
+
+, we shall use the same discretization 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+A finite element approximation to an solution field 
+\begin_inset Formula $\phi(\vec{x})$
+\end_inset
+
+ on an element 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+ is essentially specified by the weighed sum over a fixed set of local shape
+ functions 
+\begin_inset Formula $\phi_{e,1}(\vec{\xi}),\phi_{e,2}(\vec{\xi}),\ldots,\phi_{e,N}(\vec{\xi})$
+\end_inset
+
+ defined on the appropriate reference cell,
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\phi(\vec{x})=\sum_{n=1}^{N}d_{e,n}\phi_{e,n}(\vec{\xi})\]
+
+\end_inset
+
+where the global point 
+\begin_inset Formula $\vec{x}\in\Omega_{e}$
+\end_inset
+
+ corresponds to the local reference point 
+\begin_inset Formula $\vec{\xi}=\vec{x}_{e}^{-1}(\vec{x})\in\hat{\Omega}$
+\end_inset
+
+, and the weights 
+\begin_inset Formula $d_{e,n}$
+\end_inset
+
+, also known as degrees of freedom, are given on each global node.
+ These shape functions define a basis for the function space on the finite
+ element 
+\begin_inset Formula $\Omega_{e}$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+For a mesh consisting of a single element type, we need to specify the 
+\begin_inset Formula $(x_{n},y_{n},z_{n})$
+\end_inset
+
+ of its degrees of freedom, and the connectivity relations 
+\begin_inset Formula $\Omega_{e}=\{n_{1},n_{2},\ldots\}$
+\end_inset
+
+ among them which define each individual element in the corresponding discretiza
+tion.
+\end_layout
+
+\begin_layout Standard
+The most common element types used in the finite element method are given
+ by the following shape functions, which define the interior of a tetrahedron
+ with vertices at 
+\begin_inset Formula $a=(0,0,0)$
+\end_inset
+
+, 
+\begin_inset Formula $b=(1,0,0)$
+\end_inset
+
+, 
+\begin_inset Formula $c=(0,1,0)$
+\end_inset
+
+, 
+\begin_inset Formula $d=(0,0,1)$
+\end_inset
+
+, 
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Standard
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "b"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status open
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+N_{a} & = & \frac{1}{2}(-1-x-y-z)\\
+N_{b} & = & \frac{1}{2}(1+x)\\
+N_{c} & = & \frac{1}{2}(1+y)\\
+N_{d} & = & \frac{1}{2}(1+z)\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "c"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/reference-tet4.png
+	width 5cm
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Caption
+
+\begin_layout Standard
+Reference Tetrahedron
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Likewise, a hexahedral element with vertices located at 
+\begin_inset Formula $a=(-1,-1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $b=(1,-1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $c=(1,1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $d=(-1,1,-1)$
+\end_inset
+
+, 
+\begin_inset Formula $e=(-1,-1,1)$
+\end_inset
+
+, 
+\begin_inset Formula $f=(1,-1,1)$
+\end_inset
+
+, 
+\begin_inset Formula $g=(1,1,1)$
+\end_inset
+
+, 
+\begin_inset Formula $h=(-1,1,1)$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+placement H
+wide false
+sideways false
+status open
+
+\begin_layout Standard
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "c"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/reference-hex8.png
+	width 5cm
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\begin_inset Box Frameless
+position "c"
+hor_pos "c"
+has_inner_box 1
+inner_pos "c"
+use_parbox 0
+width "50col%"
+special "none"
+height "1in"
+height_special "totalheight"
+status open
+
+\begin_layout Standard
+\begin_inset Formula \begin{eqnarray*}
+N_{a} & = & \frac{1}{8}\left(1-x\right)\left(1-y\right)\left(1-z\right)\\
+N_{b} & = & \frac{1}{8}\left(1+x\right)\left(1-y\right)\left(1-z\right)\\
+N_{c} & = & \frac{1}{8}\left(1-x\right)\left(1+y\right)\left(1-z\right)\\
+N_{d} & = & \frac{1}{8}\left(1+x\right)\left(1+y\right)\left(1-z\right)\\
+N_{e} & = & \frac{1}{8}\left(1-x\right)\left(1-y\right)\left(1+z\right)\\
+N_{f} & = & \frac{1}{8}\left(1+x\right)\left(1-y\right)\left(1+z\right)\\
+N_{g} & = & \frac{1}{8}\left(1-x\right)\left(1+y\right)\left(1+z\right)\\
+N_{h} & = & \frac{1}{8}\left(1+x\right)\left(1+y\right)\left(1+z\right)\end{eqnarray*}
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Caption
+
+\begin_layout Standard
+Reference Hexahedron
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+XXX: Take description of higher order elements from Karniadakis book
+\end_layout
+
+\begin_layout Subsection
+Comparing Two Finite Element Fields
+\end_layout
+
+\begin_layout Standard
+Comparing two arbitrary finite element fields can be accomplished with the
+ 
+\family typewriter
+cigma compare
+\family default
+ command-line utility.
+ By default, the comparison will involve a numerical integration over each
+ of the elements in the mesh associated with the first field.
+ If this discretization is inadequate, you may also specify an alternative
+ discretization over which to perform the integration.
+ The comparison operation will output the local residual values are into
+ the specified output file.
+\end_layout
+
+\begin_layout Standard
+A basic comparison can be as simple as specifying the following arguments:
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --first=field1.h5:/field1/stepN
+\end_layout
+
+\begin_layout LyX-Code
+
+\series bold
+ 
+\series default
+             --second=field2.h5:/field2/stepN
+\end_layout
+
+\begin_layout LyX-Code
+              --output=residuals.vtk
+\end_layout
+
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout Subsection
+Comparing against Known Values
+\end_layout
+
+\begin_layout Standard
+A finite element description might not always be available for one of the
+ fields.
+ However, you can break the comparison into several steps if you have a
+ means to compute that field on any of the required points.
+\end_layout
+
+\begin_layout Standard
+First, extract the global coordinates of the integration points.
+ This will result in an explicit list of points over which to evaluate your
+ field.
+\end_layout
+
+\begin_layout LyX-Code
+cigma extract field1.h5:/model/mesh/ -o points.h5:/projected_points
+\end_layout
+
+\begin_layout Standard
+Now you can evaluate your function at the designated points.
+ You can provide the path to the explicit set of values with
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --first=field1.h5:/field1/stepN
+\end_layout
+
+\begin_layout LyX-Code
+              --second-points=points.h5:/projected_points
+\end_layout
+
+\begin_layout LyX-Code
+              --second-values=values.h5:/projected_values
+\end_layout
+
+\begin_layout LyX-Code
+              --output=residuals.vtk
+\end_layout
+
+\begin_layout Subsection
+\begin_inset LatexCommand label
+name "sub:Comparing-against-a"
+
+\end_inset
+
+Comparing against a Known Function
+\end_layout
+
+\begin_layout Standard
+If one of your fields is easily described by an analytic expression, then
+ you also have the option to compile your analytic function as a builtin
+ Cigma function.
+ This will enable you to reference your function by name when using the
+ 
+\family typewriter
+compare
+\family default
+ command.
+ For example,
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --first=field1.h5:/vars/displacement/step0
+\end_layout
+
+\begin_layout LyX-Code
+
+\series bold
+ 
+\series default
+             --second=
+\bar under
+disloc3d
+\end_layout
+
+\begin_layout LyX-Code
+              --output=residuals.vtk
+\end_layout
+
+\begin_layout Standard
+You may also interact with your analytic function by using the 
+\family typewriter
+cigma
+\family default
+ 
+\family typewriter
+eval
+\family default
+ command, and obtain a set of values which may then be passed back to the
+ 
+\family typewriter
+compare
+\family default
+ command.
+ 
+\end_layout
+
+\begin_layout LyX-Code
+cigma eval --function=
+\bar under
+disloc3d
+\end_layout
+
+\begin_layout LyX-Code
+           --points=points.h5:/projected_points
+\end_layout
+
+\begin_layout LyX-Code
+           --values=values.h5:/disloc3d_values
+\end_layout
+
+\begin_layout Standard
+Once these values have been obtained, we can use it to resume the calculation
+ of 
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare --output=residuals.vtk
+\end_layout
+
+\begin_layout LyX-Code
+              --first=field1.h5:/vars/displacement/step0
+\end_layout
+
+\begin_layout LyX-Code
+              --second-points=points.h5:/projected_points
+\end_layout
+
+\begin_layout LyX-Code
+              --second-values=values.h5:/disloc3d_values
+\end_layout
+
+\begin_layout Standard
+Another built-in function you might find useful is the 
+\family typewriter
+\bar under
+zero
+\family default
+\bar default
+ function, which is defined to return 0 everywhere.
+ You can use this function for estimating the norm of your field.
+\end_layout
+
+\begin_layout Subsection
+Specifying a Integration Mesh
+\end_layout
+
+\begin_layout Standard
+To override the mesh used in the integration, you can specify an extra argument
+ providing the location of the mesh,
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-mesh=mesh.h5:/model/mesh/
+\end_layout
+
+\begin_layout Standard
+Alternatively, you can also specify the coordinates and connectivity arrays
+ separately, in case they reside in separate files.
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-mesh-coordinates=file1.h5:/model/mesh/coordinates
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-mesh-connectivity=file2.h5:/model/mesh/connectivity
+\end_layout
+
+\begin_layout Subsection
+Specifying a Quadrature Rule
+\end_layout
+
+\begin_layout Standard
+To specify a quadrature rule, you will have to provide the quadrature weights
+ and points in the appropriate reference element.
+ This can be done with the following additional argument:
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-rule=quadrature-rules.h5:/path/to/rule
+\end_layout
+
+\begin_layout Standard
+You may also specify the location of the points and weights separately:
+\end_layout
+
+\begin_layout LyX-Code
+cigma compare [...]
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-rule-points=file.h5:/path/to/rule/points
+\end_layout
+
+\begin_layout LyX-Code
+  --quadrature-rule-weights=file.h5:/path/to/rule/weights
+\end_layout
+
+\begin_layout LyX-Code
+
+\end_layout
+
+\begin_layout Standard
+In the 
+\family typewriter
+src/
+\family default
+ directory of the Cigma distribution, you can find a Python script called
+ 
+\family typewriter
+rules.py
+\family default
+ that uses FIAT for calculating an appropriate set of quadrature rules for
+ use in Cigma.
+\end_layout
+
+\begin_layout Chapter
+Examples
+\end_layout
+
+\begin_layout Standard
+In this chapter we will use sample datasets from the strike-slip benchmark
+ case defined by the CIG Short-Term Tectonics working group.
+ In this benchmark problem, we solve for the viscoelastic relaxation of
+ stresses from a single finite earthquake while ignoring gravity.
+ The problem is defined on a cube domain with sides of 24 km consisting
+ of two layers of different material types.
+ The top layer of the cube is nearly elastic, while the bottom layer is
+ viscoelastic.
+ (XXX: Anything else from Pylith manual or source?)
+\end_layout
+
+\begin_layout Standard
+\begin_inset Float figure
+wide false
+sideways false
+status collapsed
+
+\begin_layout Standard
+\begin_inset Graphics
+	filename figures/strike-slip-fault.png
+	width 5cm
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Caption
+
+\begin_layout Standard
+Geometry of Strike-Slip Benchmark
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Section
+Convergence
+\end_layout
+
+\begin_layout Standard
+When the exact solution of your equations is known, the norm of the error
+ for a finite element solution is easily computed through the procedure
+ given in Chapter 
+\begin_inset LatexCommand ref
+reference "cha:Error-Analysis"
+
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+Another way of checking the accuracy of a solution is by rerunning the same
+ problem with a finer mesh.
+ Recall that in the finite element method, the weight functions and trial
+ solutions are chosen so that the finite element method is convergent.
+ In other words, as the element size 
+\begin_inset Formula $h$
+\end_inset
+
+ decreases, the solution tends to the correct solution.
+ Thus, if the difference between the coarse and fine mesh solutions is small,
+ we can conclude that the coarse mesh solution is quite accurate.
+ However, if a solution changes noticeably with refinement of the mesh,
+ the coarse mesh solution is inaccurate, and even the finer mesh may be
+ inadequate as well.
+\end_layout
+
+\begin_layout Standard
+As can be seen from these results, the log of the error varies linearly
+ with element size and the slope depends on the order of the element.
+ Denoting the slope by 
+\begin_inset Formula $\alpha$
+\end_inset
+
+, then the error in the function can be expressed as
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\log\varepsilon=\log C+\alpha\log h\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+where 
+\begin_inset Formula $\log C$
+\end_inset
+
+ is an arbitrary constant, the y-intercept of the curve.
+ The slope 
+\begin_inset Formula $\alpha$
+\end_inset
+
+ is the rate of convergence of the element.
+ We can rewrite this as
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula \[
+\varepsilon=Ch^{\alpha}\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+XXX: list of examples to discuss here
+\end_layout
+
+\begin_layout Standard
+Viscoelastic strikeslip example (tet vs.
+ hex elts) -- pylith-tet / pylith-hex
+\end_layout
+
+\begin_layout Standard
+Cylinder extension for maxwell material (newtonian) -- pylith / geofest
+ / analytic
+\end_layout
+
+\begin_layout Standard
+Cylinder extension for maxwell material (non-newtonian) -- pylith / geofest
+ / analytic
+\end_layout
+
+\begin_layout Standard
+Cylinder gravitation -- pylith / geofest / analytic
+\end_layout
+
+\begin_layout Standard
+Laplace example (linear, quadratic elements) -- deal.II / libmesh / analytic
+\end_layout
+
+\begin_layout Standard
+Citcom example (structured datasets) -- citcom (3d) / analytic
+\end_layout
+
+\begin_layout Standard
+Gale example (structured datasets) -- gale (2d) / analytic
+\end_layout
+
+\begin_layout Standard
+MADDs example -- 
+\end_layout
+
+\begin_layout Section
+Visualization
+\end_layout
+
+\begin_layout Standard
+The residual field obtained in the 
+\end_layout
+
+\begin_layout Chapter
+\start_of_appendix
+File Formats
+\end_layout
+
+\begin_layout Standard
+To differentiate between these input and output file formats, you will need
+ to provide an appropriate file extension.
+ Currently recognized file extensions include 
+\family typewriter
+*.h5
+\family default
+ for HDF5 files, 
+\family typewriter
+*.vtk
+\family default
+ for legacy VTK files, and 
+\family typewriter
+*.dat
+\family default
+ for simple text files.
+\end_layout
+
+\begin_layout Section
+Input File Format
+\end_layout
+
+\begin_layout Standard
+There are five different objects you will want to use as input to Cigma,
+ all of which can be represented as two-dimensional datasets.
+ These are the 
+\end_layout
+
+\begin_layout Itemize
+mesh coordinates
+\end_layout
+
+\begin_layout Itemize
+mesh connectivity
+\end_layout
+
+\begin_layout Itemize
+field coefficients (degrees of freedom)
+\end_layout
+
+\begin_layout Itemize
+quadrature rule points
+\end_layout
+
+\begin_layout Itemize
+quadrature rule weights
+\end_layout
+
+\begin_layout Subsection
+HDF5
+\end_layout
+
+\begin_layout Standard
+HDF5 files are organized in a hierarchical structure, similar to a UNIX
+ file system.
+ Two types of primary objects, groups and datasets, are stored in this structure.
+ A group contains instances of zero or more groups or datasets, while a
+ dataset stores a multi-dimensional array of data elements.
+ Both kinds of objects are accompanied by supporting metadata known as attribute
+s.
+\end_layout
+
+\begin_layout Standard
+A dataset is physically stored in two parts: a header and a data array.
+ The header contains miscellaneous metadata describing the dataset as well
+ as information that is needed to interpret the array portion of the dataset.
+ Essentially, it includes the name, datatype, dataspace, and storage layout
+ of the dataset.
+ The name is a text string identifying the dataset.
+ The datatype describes the type of the data array elements.
+ The dataspace defines the dimensionality of the dataset, i.e., the size and
+ shape of the multi-dimensional array.
+\end_layout
+
+\begin_layout Standard
+In Cigma you always provide an explicit path to every dataset, so you have
+ a fair amount of flexibility for how you organize your datasets inside
+ HDF5 files.
+ For example, a typical Cigma HDF5 file could have the following structure.
+\end_layout
+
+\begin_layout LyX-Code
+model.h5
+\end_layout
+
+\begin_layout LyX-Code
+
+\backslash
+__ model
+\end_layout
+
+\begin_layout LyX-Code
+    |__ mesh
+\end_layout
+
+\begin_layout LyX-Code
+    |   |__ coordinates    [nno x nsd]
+\end_layout
+
+\begin_layout LyX-Code
+    |   
+\backslash
+__ connectivity   [nel x ndof]
+\end_layout
+
+\begin_layout LyX-Code
+    
+\backslash
+__ variables
+\end_layout
+
+\begin_layout LyX-Code
+        |__ temperature
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00000  [nno x 1]
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00010  [nno x 1]
+\end_layout
+
+\begin_layout LyX-Code
+        |   
+\backslash
+...
+\end_layout
+
+\begin_layout LyX-Code
+        |__ displacement
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00000  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+        |   |__ step00010  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+        |   
+\backslash
+...
+\end_layout
+
+\begin_layout LyX-Code
+        
+\backslash
+__ velocity
+\end_layout
+
+\begin_layout LyX-Code
+            |__ step00000  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+            |__ step00010  [nno x 3]
+\end_layout
+
+\begin_layout LyX-Code
+            
+\backslash
+...
+\end_layout
+
+\begin_layout Standard
+Generally, Cigma will only require you to specify the path to a specific
+ field dataset.
+ If a small amount of metadata is present in your field dataset, the rest
+ of the required information, such as which mesh and finite elements to
+ use, will be deduced from that metadata.
+\end_layout
+
+\begin_layout Description
+MeshID an identifier assigned for use in linking child datasets to their
+ parent mesh.
+\end_layout
+
+\begin_layout Description
+MeshLocation points to the HDF5 group which contains the appropriate coordinates
+ and connectivity datasets.
+\end_layout
+
+\begin_layout Description
+FunctionSpace string identifier to determine which shape functions to use
+ for interpolating values inside the element (e.g., tet4, hex8, quad4, tri3,
+ ...).
+\end_layout
+
+\begin_layout Description
+DatasetType string identifier for classifying the type of data contained
+ in the dataset (e.g., points, connectivity, degrees of freedom, quadrature
+ rules, global quadrature points, global field values).
+\end_layout
+
+\begin_layout Subsection
+VTK
+\end_layout
+
+\begin_layout Standard
+The current version of Cigma only supports VTK unstructured grid datasets
+ of a single element type.
+ You can still compare various unstructured grids with different element
+ types to each other.
+\end_layout
+
+\begin_layout Standard
+Note that while you typically provide a path (or name) for every dataset,
+ this is not necessary when specifying a VTK mesh, since this data is taken
+ from the special Points and Cells arrays, which you cannot rename.
+ However, you will still need to provide a name when referring to the field
+ coefficients, which are assumed to be stored as Point Data in the input
+ VTK file.
+\end_layout
+
+\begin_layout Standard
+For more information, you may want to refer to Visualization ToolKit's 
+\begin_inset LatexCommand htmlurl
+name "File Formats"
+target "www.vtk.org/pdf/file-formats.pdf"
+
+\end_inset
+
+ document.
+ 
+\end_layout
+
+\begin_layout Subsection
+Text
+\end_layout
+
+\begin_layout Standard
+The text format is always in table form, with the dimensions of the table
+ specified in the first line.
+ For example, mesh coordinates can be specified in the following format
+\end_layout
+
+\begin_layout LyX-Code
+<nno> <nsd>
+\end_layout
+
+\begin_layout LyX-Code
+x1 y1 z1
+\end_layout
+
+\begin_layout LyX-Code
+x2 y2 z2
+\end_layout
+
+\begin_layout LyX-Code
+x3 y3 z3
+\end_layout
+
+\begin_layout LyX-Code
+...
+\end_layout
+
+\begin_layout Standard
+Mesh connectivity with 
+\end_layout
+
+\begin_layout LyX-Code
+nel ndofs
+\end_layout
+
+\begin_layout LyX-Code
+node_1 node_2 node_3 node_4 ...
+\end_layout
+
+\begin_layout LyX-Code
+node_1 node_2 node_3 node_4 ...
+\end_layout
+
+\begin_layout LyX-Code
+node_1 node_2 node_3 node_4 ...
+\end_layout
+
+\begin_layout LyX-Code
+...
+\end_layout
+
+\begin_layout Standard
+A generic field with 
+\family typewriter
+ndim
+\family default
+ components (i.e., scalar, vector, or tensor) is specified by
+\end_layout
+
+\begin_layout LyX-Code
+num ndim
+\end_layout
+
+\begin_layout LyX-Code
+f1 f2 f3 ..
+\end_layout
+
+\begin_layout LyX-Code
+f1 f2 f3 ..
+\end_layout
+
+\begin_layout LyX-Code
+...
+\end_layout
+
+\begin_layout Standard
+In this last case, the number of rows could refer to either 
+\family typewriter
+nno
+\family default
+ or 
+\family typewriter
+nel
+\family default
+, depending on whether the field is node-based or cell-based.
+\end_layout
+
+\begin_layout Section
+Output File Format
+\end_layout
+
+\begin_layout Standard
+The output format for residuals consists simply of a list of scalars for
+ each element in the discretization.
+ If you specify a 
+\family typewriter
+.vtk
+\family default
+ extension for your output file, this will result in an scalar dataset named
+ epsilon stored using the legacy VTK file format in the Cell Data section
+ of the output file.
+ Note that this output consists of the squared values of the local residuals,
+ so further post-processing will be necessary.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "1"
+key "Akin 2005"
+
+\end_inset
+
+Akin, J.E.
+ (2005), 
+\emph on
+Finite Element Analysis with Error Estimators,
+\emph default
+ Butterworth-Heinemann, Oxford, 447 pp.
+ 
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "2"
+key "Encyclopaedia of Cubature Formulas 2005"
+
+\end_inset
+
+Encyclopaedia of Cubature Formulas (1998-2005), Katholieke Universiteit
+ Leuven, Dept.
+ Computerwetenschappen, www.cs.kuleuven.be/~nines/research/ecf/
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "3"
+key "Hughes 2000"
+
+\end_inset
+
+Hughes, Thomas J.R.
+ (2000), 
+\emph on
+The Finite Element Method: Linear Static and Dynamic Finite Element Analysis
+\emph default
+, Dover, 672 pp.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "4"
+key "KarniadakisSherwin2005"
+
+\end_inset
+
+Karniadakis, George E.
+ and Spencer J.
+ Sherwin (2005), 
+\emph on
+Spectral/
+\emph default
+hp 
+\emph on
+Element Methods for Computational Fluid Dynamics, Second Edition
+\emph default
+, Numerical Mathematics and Scientific Computation, Oxford, 657 pp.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset LatexCommand bibitem
+label "5"
+key "Uesu-etal2005"
+
+\end_inset
+
+Uesu, D., L.
+ Bavoil, S.
+ Fleishman, J.
+ Shepherd, and C.
+ Silva (2005), Simplification of Unstructured Tetrahedral Meshes by Point-Sampli
+ng, 
+\emph on
+Proceedings of the 2005 International Workshop on Volume Graphics,
+\emph default
+ 157-165, doi: 10.2312/VG/VG05/157-165.
+\end_layout
+
+\end_body
+\end_document