Reorganizing theory and user manuals.

Refs #3 and idaholab#156.

doc/content/bib/mastodon.bib

@@ -19,16 +19,6 @@ @article{housner1952hsi
   year={1952}
 }
 
-@article{lsymer1969finite,
-  author={J.~Lsymer and R.~L.~Kuhlemeyer},
-  title={Finite dynamic model for infinite media},
-  journal={Journal of the Engineering Mechanics Division},
-  volume={95},
-  number={4},
-  pages={859-878},
-  year={1969}
-}
-
 @article{chiang1994anew,
   author={D.~Y.~Chiang and J.~L.~Beck},
   title={A new class of distributed-element models for cyclic plasticity -I. Theory and application.},

doc/content/manuals/include/bcs/bcs-user.md

@@ -0,0 +1,62 @@
+## Boundary Conditions
+
+Boundary conditions are required by MASTODON finite element analysis framework to be able to run
+simple quasi-static and dynamic analysis. Herein, the basic boundary conditions are given and other
+constraints that are used to solve for geotechnical earthquake engineering problems are presented
+separately (e.g. periodic boundary conditions.) Following input block creates a fully fixed boundary condition at the bottom of a single element.
+
+!listing test/tests/materials/isoil/HYS_darendeli.i
+         start=BCs
+         end=Periodic
+
+In the above input, `type = PresetBC` sets the variable = disp_”related degree of freedom” with the
+value = 0 which provides the fixity by defining a zero displacement at the node. boundary = 0 command
+selects the nodes at the bottom surface of the element (labeled as surface 0) and assigns the
+boundary conditions.
+
+#### Prescribed Displacement
+
+The preset displacement boundary condition can be used to apply a displacement time history to a
+boundary (at the nodes). The displacement boundary condition first converts the user defined
+displacement time history to an acceleration time history using Backward Euler finite difference
+scheme. This acceleration is then integrated to get displacement using Newmark-beta method.  The
+resulting displacement is then applied as a kinematic displacement boundary condition. The following
+command can be used to apply the preset displacement boundary condition:
+
+!listing test/tests/materials/isoil/HYS_darendeli.i
+         start=top_x
+         end=Functions
+
+The above command should be embedded inside the BCs command block. “boundary = 5” assigns the preset
+displacement to boundary 5 which, in this case, is a predefined boundary of a single element as
+described in single element problem above. Alternatively, the boundary number can be identified using [Cubit](https://cubit.sandia.gov/) or [Trelis](https://www.csimsoft.com/trelis.jsp).  “variable = disp_x” imposes the boundary condition on the x
+direction. “beta” is the Newmark-beta integration parameter. The “function = top_disp” specifies the
+function that defines the loading time history. It is defined in the “Functions” block as follows:
+
+!listing test/tests/materials/isoil/HYS_darendeli.i
+         start=Functions
+         end=Materials
+
+Displacement2.csv is the file, located in the same directory of the input file, containing the
+displacement time history. The first column of this file should contain the time vector starting at
+0.0. The second column should contain the displacement values. “type = PiecewiseLinear” defines the
+type of the function which is in this case piecewise-linear. “format” specifies the format of the
+data file, i.e. whether the data is in columns or rows.
+
+#### Prescribed Acceleration
+
+The preset acceleration boundary condition can be used to apply an acceleration time history to a
+boundary. The preset acceleration boundary condition integrates the given acceleration time history
+to get the displacement using Newmark-beta method. This displacement is then applied as a kinematic
+displacement boundary condition. Syntax is the same as prescribing a displacement boundary condition
+but with type = PresetAcceleration and the function describing time vs acceleration data instead of
+time vs displacement.
+
+#### Periodic Boundary Conditions
+
+Periodic boundary conditions are used to constrain the nodes to move together in the specified
+directions. The following input is an example applied on the single element problem above and should be embedded into the BCs block segment as:
+
+!listing test/tests/materials/isoil/HYS_darendeli.i
+         start=Periodic
+         end=top_x

doc/content/manuals/include/bcs/drm-theory.md

@@ -0,0 +1,34 @@
+#### Domain reduction method (DRM)
+
+Earthquake 'source-to-site' simulations require simulating a huge soil domain (order of many
+kilometers) with a earthquake fault. The nuclear power plant structure, which is usually less than
+100 m wide, is located very far from the earthquake fault, and the presence of the structure only
+affects the response of the soil in the vicinity of the structure. In most of these situations, where
+a localized feature such as a structure is present in a huge soil domain, the problem can be divided
+into two parts: (i) a free-field 'source-to-site' simulation is run on the huge soil domain (
+[fig:DRM](a)) that does not contain the localized feature, and (ii) the forces from the free-field
+simulation at one element layer, which is the element layer separating the bigger and smaller soil
+domain, can be transferred to a much smaller domain containing the localized feature (
+[fig:DRM](b)). This method of reducing the domain is called the domain reduction method (DRM)
+[citep:bielak2003domain]. [fig:DRM] is reproduced from [citet:bielak2003domain].
+
+!media media/theory/DRM.png
+       style=width:100%;float:center;
+       id=fig:DRM
+       caption=Domain reduction method summary: (a) Big soil domain containing the earthquake fault
+               but not the localized feature. The displacements are obtained at the boundaries
+               $\Gamma$ and $\Gamma_e$ and are converted to equivalent forces. (b) Smaller soil
+               domain containing the localized feature but not the earthquake fault. The equivalent
+               forced calculated in (a) are applied at the boundaries $\Gamma$ and $\Gamma_e$.
+
+To convert the displacements at $\Gamma$ and $\Gamma_e$ from part (i) to equivalent forces, a finite
+element model of the one element layer between $\Gamma$ and $\Gamma_e$ is simulated in two
+steps. First, the boundary $\Gamma_e$ is fixed and the boundary $\Gamma$ is moved with the
+displacements recorded at $\Gamma$. This step gives the equivalent forces at $\Gamma_e$. Second, the
+boundary $\Gamma$ is fixed and the boundary $\Gamma_e$ is moved with the displacements recorded at
+$\Gamma_e$. This steps gives the equivalent forces at $\Gamma$.
+
+Note: The meshes for the bigger soil domain and smaller soil domain need not align between $\Gamma$
+and $\Gamma_e$. The equivalent forces can be applied as point forces at the same coordinate location
+at which nodes are present in the bigger model, irrespective of whether these locations correspond to
+nodal locations in the smaller model.

doc/content/manuals/include/bcs/intro_bcs-theory.md

@@ -0,0 +1 @@
+## Boundary Conditions

doc/content/manuals/include/bcs/non_reflecting-theory.md

@@ -0,0 +1,27 @@
+#### Non-reflecting boundary
+
+This boundary condition applies a Lysmer damper [citep:lysmer1969finite] on a given boundary to
+absorb the waves hitting the boundary. To understand Lysmer dampers, let us consider an uniform
+linear elastic soil column and say a 1D vertically propagating P wave is traveling through this soil
+column. Then the normal stress at any location in the soil column is given by:
+
+\begin{equation}
+\label{eqn:normal_stress}
+\sigma = E \epsilon = E \frac{du}{dx} = \frac{E}{V_p} \frac{du}{dt}= \rho V_p \frac{du}{dt}
+\end{equation}
+
+where, $E$ is the Young's modulus, $\sigma$ is the normal stress, $\epsilon$ is the normal strain,
+$\rho$ is the density, $V_p = \sqrt{\frac{E}{\rho}}$ is the P-wave speed and $\frac{du}{dt}$ is the
+particle velocity. Note that for a 3D problem, the P-wave speed is $V_p = \sqrt{\frac{E(1-\
+nu)}{(1+\nu)(1-2\nu)}}$.
+
+The stress in the above equation is directly proportional to the particle velocity which makes this
+boundary condition analogous to a viscous damper with damping coefficient of $\rho V_p$. So
+truncating the soil domain and placing this damper at the end of the domain is equivalent to
+simulating wave propagation in an infinite soil column provided the soil is made of linear elastic
+material and the wave is vertically incident on the boundary.
+
+If the soil is not linear elastic or if the wave is incident at an angle on the boundary, the waves
+are not completely absorbed by the Lysmer damper. However, if the non-reflecting boundary is placed
+sufficiently far from the region of interest, any reflected waves will get damped out by Rayleigh
+damping or hysteretic material behavior before it reaches the region of interest.

doc/content/manuals/include/bcs/preset_acceleration-theory.md

@@ -0,0 +1,10 @@
+#### Preset acceleration
+
+If the ground excitation was measured at a depth within the soil by placing an accelerometer at that
+location, then it is termed as a within-soil input. This time history contains the wave that was
+generated by the earthquake (incoming wave) and the wave that is reflected off the free surface. This
+ground excitation time history is usually available in the form of a acceleration time history. Since
+MASTODON is a displacement controlled algorithm, i.e., displacements are the primary unknown
+variables, the acceleration time history is first converted to the corresponding displacement time
+history using Newmark time integration equation. This displacement time
+history is then prescribed to the boundary.

doc/content/manuals/include/bcs/seismic_force-theory.md

@@ -0,0 +1,27 @@
+#### Seismic force
+
+In some cases, the ground excitation is measured at a rock outcrop (where rock is found at surface
+level and there is no soil above it). To apply this to a location where rock is say $10$m deep and
+there is soil above it, a sideset is created at $10$m depth and the ground excitation (converted into
+a stress) is applied at this depth. To apply ground excitation as a stress, the input function should
+be given as ground velocity.
+
+To convert a velocity applied normal to the boundary into a normal stress, the normal stress equation above can be used. Using a similar argument as discussed in the section above, to
+convert a velocity applied tangential to the boundary into a shear stress, Equation
+[eqn:shear_stress] can be used.
+
+\begin{equation}
+\label{eqn:shear_stress}
+\tau = \rho  V_s  \frac{du}{dt}
+\end{equation}
+where, $V_s$ is the shear wave speed and $\tau$ is the shear stress.
+
+In some situations, the ground motion measured at a depth within the soil is available. This ground
+motion is the summation of the wave that enters and exits the soil domain. MASTODON has the
+capability to extract the incoming wave from the within soil ground motion. To calculate the incoming
+wave velocity, an iterative procedure is used. The initial guess for the incoming wave velocity
+($v_i$) at time t is taken to be the same as the within soil velocity measured at that location. The
+velocity at this boundary obtained from MASTODON ($v_{mastodon}$) is now going to be different from
+the measured within soil velocity ($v_{measured}$) at time t. Half the difference between
+$v_{mastodon}$ and $v_{measured}$ is added to $v_o$ and the iterations are continued until $v_i$
+converges (within a numerical tolerance).

doc/content/manuals/include/contact/intro_contact-theory.md

@@ -0,0 +1,5 @@
+## Foundation-soil interface models
+
+The foundation-soil interface is an important aspect of NLSSI modeling. The foundation-soil interface
+simulates geometric nonlinearities in the soil-structure system: gapping (opening and closing of gaps
+between the soil and the foundation), sliding, and uplift.

doc/content/manuals/include/contact/thin_layer-theory.md

@@ -0,0 +1,57 @@
+#### Thin-layer method
+
+An efficient approach to modeling the foundation-soil interface is to create a thin layer of the
+I-Soil material at the interface, as illustrated in [fig:thin_layer] below.
+
+!media media/theory/thin_layer.png
+       style=width:60%;margin-left:100px;float:center;
+       id=fig:thin_layer
+       caption=Modeling the foundation-soil interface as a thin layer for a sample surface foundation.
+
+The red layer between the foundation (green) and soil (yellow) is the thin layer of I-Soil. The
+properties of this thin layer are then adjusted to simulate Coulomb friction between the
+surfaces. The Coulomb-friction-type behavior can be achieved by modeling the material of the thin
+soil layer as follows:
+
+1. Define an I-Soil material with a user-defined stress-strain curve.
+
+2. Calculate the shear strength of the thin layer as $\tau_{max}=\mu \sigma_N$ , where $\tau_{max}$
+   is the shear strength, $\mu$ is the friction coefficient, and $\sigma_N$ is the normal stress on
+   the contact surface. The shear strength of the thin layer is the stress at which sliding starts at
+   the interface. Therefore, this shear strength should be proportional to the normal stress to
+   simulate Coulomb friction. This can be achieved by setting the initial shear strength equal to the
+   reference pressure, $p_{ref}$. The reference pressure can then be set to an arbitrary positive
+   value, such as the average normal stress at the interface from gravity loads. The shear strength
+   will then follow the equation
+
+   \begin{equation}
+   \tau_{max} = \mu p_{ref}
+   \end{equation}
+
+3. Define the stress-strain curve to be almost elastic-perfectly-plastic, and such that the shear
+   modulus of the thin layer is equal to the shear modulus of the surrounding soil, in case of an
+   embedded foundation. If the foundation is resting on the surface such as in [fig:thin_layer]
+   above, the shear modulus of the thin layer soil should be as high as possible, such that the
+   linear horizontal foundation stiffness is not reduced due to the presence of the thin layer. A
+   sample stress-strain curve is shown in [fig:thin_layer_stress_strain] below. The sample curve in
+   the figure shows an almost bilinear shear behavior with gradual yielding and strain hardening,
+   both of which, are provided to reduce possible high-frequency response. High-frequency response is
+   likely to occur if a pure Coulomb friction model (elastic-perfectly-plastic shear behavior at the
+   interface) is employed, due to the sudden change in the interface shear stiffness to zero.
+
+!media media/theory/thin_layer_stress_strain.png
+       style=width:60%;margin-left:150px;float:center;
+       id=fig:thin_layer_stress_strain
+       caption=Sample shear-stress shear-strain curve for modeling the thin-layer interface using I-Soil.
+
+4. Turn on pressure dependency of the soil stress-strain curve and set $a_0$, $a_1$ and $a_2$ to 0, 0
+   and 1, respectively. This ensures that the stress-strain curve scales linearly with the normal
+   pressure on the interface, thereby also increasing the shear strength in the above equation
+   linearly with the normal pressure, similar to Coulomb friction.
+
+5. Use a large value for the Poisson’s ratio, in order to avoid sudden changes in the volume of the
+   thin-layer elements after the yield point is reached. A suitable value for the Poisson’s ratio can
+   be calculated by trial and error.
+
+Following the above steps should enable the user to reasonably simulate geometric
+nonlinearities. These steps will be automated in MASTODON in the near future.

doc/content/manuals/include/damping/frequency_independent-theory.md

@@ -0,0 +1,3 @@
+#### Frequency-independent damping
+
+As seen in the previous sub-section, the damping ratio using Rayleigh damping varies with frequency. Although the parameters $\eta$ and $\zeta$ can be tuned to arrive at a constant damping ratio for a short frequency range, as the frequency range increases, the damping ratio no longer remains constant. For scenarios like these, where a constant damping ratio is required over a large frequency range, frequency independent damping formulations work better. This formulations is under consideration for adding to MASTODON.

doc/content/manuals/include/damping/intro_damping-theory.md

@@ -0,0 +1,10 @@
+## Damping
+
+When the soil-structure system (including both soil and concrete)
+responds to an earthquake excitation, energy is dissipated in two primary
+ways: (1)small-strain and hysteretic material damping, and (2) damping due to gapping,
+sliding and uplift at the soil-foundation interface. Dissipation of
+energy due to item (1) is modeled (approximately) using following methods: (i) viscous damping for small strain damping experienced at very small strain
+levels ($\gamma$ $\leq 0.001 \%$) where the material behavior is largely linear viscoelastic; (ii)
+hysteretic damping due to nonlinear hysteretic behavior of the material.
+Dissipation of energy due to (2) is discussed in [foundation-soil interface models](#Foundation-soil interface models). This section discusses the damping that is present at small strain levels.