| 1FlagsDefine(iComputeDeriv) | If set to 1, you must provide the computation of the amplitude derivative. The default is 0, which means that Abaqus computes the derivative automatically. |
| 1FlagsDefine(iComputeSecDeriv) | If set to 1, you must provide the computation of the amplitude second derivative. The default is 0, which means that Abaqus computes the second derivative automatically. |
| 1FlagsDefine(iComputeInteg) | If set to 1, you must provide the computation of the amplitude incremental integral. The default is 0, which means that Abaqus computes the incremental integral automatically. |
| 1FlagsDefine(iComputeDoubleInteg) | If set to 1, you must provide the computation of the amplitude incremental double integral. The default is 0, which means that Abaqus computes the incremental integral automatically. |
| 1FlagsDefine(iStopAnalysis) | If set to 1, the analysis will be stopped and an error message will be issued. The default is 0, which means that Abaqus will not stop the analysis. |
| 1FlagsDefine(iConcludeStep) | If set to 1, Abaqus will conclude the step execution and advance to the next step (if a next step exists). The default is 0. |
+
+# svars
+
+An array containing the values of the solution-dependent state variables associated with this amplitude definition. The number of such variables is nsvars (see above). You define the meaning of these variables.
+
+This array is passed into UAMP containing the values of these variables at the start of the current increment. In most cases they should be updated to be the values at the end of the increment.
+
+# AmpDerivative
+
+Current value of the amplitude derivative.
+
+# AmpSecDerivative
+
+Current value of the amplitude second derivative.
+
+
+
+# AmpIncIntegral
+
+Current value of the amplitude incremental integral.
+
+# AmpDoubleIntegral
+
+Current value of the amplitude incremental double integral.
+
+# Variables passed in for information
+
+# ampName
+
+User-specified amplitude name, left justified.
+
+# time(iStepTime)
+
+Current value of step time or frequency.
+
+# time(iTotalTime)
+
+Current value of total time.
+
+# ampValueOld
+
+Old value of the amplitude from the previous increment.
+
+# dt
+
+Time increment.
+
+# props
+
+User-specified array of material constants associated with this amplitude definition.
+
+# nProps
+
+User-defined number of material constants associated with this amplitude definition.
+
+# nSvars
+
+User-defined number of solution-dependent state variables associated with this amplitude definition.
+
+# lFlagsInfo
+
+Integer flag array with information regrading the current call to UAMP.
+
+# lFlagsInfo(iInitialization)
+
+This flag is equal to 1 if UAMP is called from the initialization phase of the first analysis step and is set to 0 otherwise.
+
+# lFlagsInfo(iRegularInc)
+
+This flag is equal to 1 if UAMP is called from a regular increment and is set to 0 if called from the initialization phase of the first analysis step.
+
+# lFlagsInfo(iCuts)
+
+Number of cutbacks in this increment.
+
+# lFlagsInfo(ikStep)
+
+Step number.
+
+
+
+# nSensor
+
+Total number of sensors in the model.
+
+# sensorValues
+
+Array with sensor values at the end of the previous increment. Each sensor value corresponds to a history output variable associated with the output database request defining the sensor.
+
+# sensorNames
+
+Array with user-defined sensor names in the entire model, left justified. Each sensor name corresponds to a sensor value provided with the output database request. All names will be converted to uppercase characters if lowercase or mixed-case characters were used in their definition.
+
+# jSensorLookUpTable
+
+Variable that must be passed into the utility functions IGETSENSORID and GETSENSORVALUE.
+
+Example: Amplitude definition using sensor and state variables
+```python
+c user amplitude subroutine
+Subroutine UAMP(
+C passed in for information and state variables
+* ampName, time, ampValueOld, dt, nProps, props, nSvars,
+* svars, lFlagsInfo,
+* nSensor, sensorValues, sensorNames,
+* jSensorLookUpTable,
+C to be defined
+* ampValueNew,
+* lFlagsDefine,
+* AmpDerivative, AmpSecDerivative, AmpIncIntegral,
+* AmpDoubleIntegral)
+include 'aba_param.inc'
+C svars - additional state variables, similar to (V)UEL
+ dimension sensorValues(nSensor), svars(nSvars),
+* props(nProps)
+ character*80 sensorNames(nSensor)
+ character*80 ampName
+C time indices
+ parameter( iStepTime = 1,
+* iTotalTime = 2,
+* nTime = 2)
+C flags passed in for information
+```
+
+
+
+```txt
+parameter( iInitialization = 1,
+* iRegularInc = 2,
+* iCuts = 3,
+* ikStep = 4,
+* nFlagsInfo = 4)
+
+C optional flags to be defined
+parameter( iComputeDeriv = 1,
+* iComputeSecDeriv = 2,
+* iComputeInteg = 3,
+* iComputeDoubleInteg = 4,
+* iStopAnalysis = 5,
+* iConcludeStep = 6,
+* nFlagsDefine = 6)
+
+parameter( tStep=0.18d0, tAccelerateMotor = .00375d0,
+* omegaFinal=23.26d0,
+* zero=0.0d0, one=1.0d0, two=2.0d0, four=4.0d0)
+
+dimension time(nTime), lFlagsInfo(nFlagsInfo),
+* lFlagsDefine(nFlagsDefine)
+dimension jSensorLookUpTable(*)
+
+lFlagsDefine(iComputeDeriv) = 1
+lFlagsDefine(iComputeSecDeriv) = 1
+lFlagsDefine(iComputeInteg) = 1
+lFlagsDefine(iComputeDoubleInteg) = 1
+
+c get sensor value
+vTrans_CU1 = GetSensorValue('HORIZ_TRANSL_MOTION',
+* jSensorLookUpTable,
+* sensorValues)
+
+if (ampName(1:22).eq. 'MOTOR_WITH_STOP_SENSOR') then
+ if (lFlagsInfo(iInitialization).eq.1) then
+ AmpSecDerivative = zero
+ AmpDerivative = omegaFinal/tAccelerateMotor
+ ampValueNew = zero
+ AmpIncIntegral = zero
+ AmpDoubleIntegral = zero
+
+ svars(1) = zero
+ svars(2) = zero
+```
+
+
+
+```vba
+else
+ tim = time(iStepTime)
+
+c ramp up the angular rot velocity of the
+c electric motor
+c after which hold constant
+ if (tim .le. tAccelerateMotor) then
+ AmpSecDerivative = zero
+ AmpDerivative = omegaFinal/tAccelerateMotor
+ ampValueNew = omegaFinal*tim/tAccelerateMotor
+ AmpIncIntegral = dt*(ampValueOld+ampValueNew)/two
+ AmpDoubleIntegral = dt**2*(ampValueOld+ampValueNew)/four
+ else
+ AmpSecDerivative = zero
+ AmpDerivative = zero
+ ampValueNew = omegaFinal
+ AmpIncIntegral = dt*(ampValueOld+ampValueNew)/two
+ AmpDoubleIntegral = dt**2*(ampValueOld+ampValueNew)/four
+ end if
+
+c retrieve old sensor value
+ vTrans_CU1_old = svars(1)
+
+c detect a zero crossing and count the number of
+c crossings
+ if (vTrans_CU1_old*vTrans_CU1 .le. zero .and.
+ * tim .gt. tAccelerateMotor ) then
+ svars(2) = svars(2) + one
+ end if
+ nrCrossings = int(svars(2))
+
+c stop the motor if sensor crosses zero the second time
+ if (nrCrossings.eq.2) then
+ ampValueNew = zero
+ lFlagsDefine(iConcludeStep)=1
+ end if
+```
+
+
+
+```txt
+c store sensor value
+ svars(1) = vTrans_CU1
+ end if
+ end if
+ return
+ end
+```
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_014.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_014.md
new file mode 100644
index 0000000..545b362
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_014.md
@@ -0,0 +1,400 @@
+
+
+# 1.1.20 UANISOHYPER\_INV: User subroutine to define anisotropic hyperelastic material behavior using the invariant formulation.
+
+Product: Abaqus/Standard
+
+# References
+
+• “Anisotropic hyperelastic behavior,” Section 22.5.3 of the Abaqus Analysis User’s Guide
+• \*ANISOTROPIC HYPERELASTIC
+• “UANISOHYPER\_INV and VUANISOHYPER\_INV,” Section 4.1.13 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine UANISOHYPER\_INV:
+
+• can be used to define the strain energy potential of anisotropic hyperelastic materials as a function of an irreducible set of scalar invariants;
+• is called at all material calculation points of elements for which the material definition contains user-defined anisotropic hyperelastic behavior with an invariant-based formulation (“Anisotropic hyperelastic behavior,” Section 22.5.3 of the Abaqus Analysis User’s Guide);
+• can include material behavior dependent on field variables or state variables;
+• requires that the values of the derivatives of the strain energy density function of the anisotropic hyperelastic material be defined with respect to the scalar invariants; and
+• is called twice per material point in each iteration.
+
+# Enumeration of invariants
+
+To facilitate coding and provide easy access to the array of invariants passed to user subroutine UANISOHYPER\_INV, an enumerated representation of each invariant is introduced. Any scalar invariant can, therefore, be represented uniquely by an enumerated invariant, $I _ { n } ^ { * }$ , where the subscript n denotes the order of the invariant according to the enumeration scheme in the following table:
+
+| Component | 2D Case | 3D Case |
| 1 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{11}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{11}^{G}$ |
| 2 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{22}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{22}^{G}$ |
| 3 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{22}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{22}^{G}$ |
| 4 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{33}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{33}^{G}$ |
| 5 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{33}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{33}^{G}$ |
| 6 | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{33}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{33}^{G}$ |
| 7 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 8 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 9 | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 10 | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 11 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial J$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{13}^{G}$ |
| 12 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial J$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{13}^{G}$ |
| 13 | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial J$ | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{13}^{G}$ |
| 14 | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial J$ | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{13}^{G}$ |
| 15 | $\partial^{2}U/\partial J^{2}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{13}^{G}\partial\overline{\varepsilon}_{13}^{G}$ |
| 16 | | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 17 | | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 18 | | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 19 | | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 20 | | $\partial^{2}U/\partial\overline{\varepsilon}_{13}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 21 | | $\partial^{2}U/\partial\overline{\varepsilon}_{23}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 22 | | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial J$ |
| 23 | | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial J$ |
| 24 | | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial J$ |
+
+
+
+| LFLAGS (1) | Defines the procedure type. See “Results file output format,” Section 5.1.2 of the Abaqus Analysis User’s Guide, for the key used for each procedure. |
| LFLAGS (2) =0 | Small-displacement analysis. |
| LFLAGS (2) =1 | Large-displacement analysis (nonlinear geometric effects included in the step; see “General and linear perturbation procedures,” Section 6.1.3 of the Abaqus Analysis User’s Guide). |
| LFLAGS (3) =1 | Normal implicit time incrementation procedure. User subroutine UEL must define the residual vector in RHS and the Jacobian matrix in AMATRX. |
| LFLAGS (3) =2 | Define the current stiffness matrix (AMATRX = $K^{NM} = -\partial F^{N}/\partial u^{M}$ or $-\partial G^{N}/\partial u^{M}$ ) only. |
| LFLAGS (3) =3 | Define the current damping matrix (AMATRX = $C^{NM} = -\partial F^{N}/\partial \dot{u}^{M}$ or $-\partial G^{N}/\partial \dot{u}^{M}$ ) only. |
+
+
+
+| LFLAGS (3) =4 | Define the current mass matrix ( $\mathbf{AMATRX} = M^{NM} = -\partial F^{N} / \partial \ddot{u}^{M}$ ) only. Abaqus/Standard always requests an initial mass matrix at the start of the analysis. |
| LFLAGS (3) =5 | Define the current residual or load vector ( $\mathbf{RHS} = F^{N}$ ) only. |
| LFLAGS (3) =6 | Define the current mass matrix and the residual vector for the initial acceleration calculation (or the calculation of accelerations after impact). |
| LFLAGS (3) =100 | Define perturbation quantities for output. |
| LFLAGS (4) =0 | The step is a general step. |
| LFLAGS (4) =1 | The step is a linear perturbation step. |
| LFLAGS (5) =0 | The current approximations to $u^{M}$ , etc. were based on Newton corrections. |
| LFLAGS (5) =1 | The current approximations were found by extrapolation from the previous increment. |
+
+# TIME(1)
+
+Current value of step time or frequency.
+
+# TIME(2)
+
+Current value of total time.
+
+# Scalar parameters:
+
+# DTIME
+
+Time increment.
+
+# PERIOD
+
+Time period of the current step.
+
+# NDOFEL
+
+Number of degrees of freedom in the element.
+
+# MLVARX
+
+Dimensioning parameter used when several displacement or right-hand-side vectors are used.
+
+# NRHS
+
+Number of load vectors. NRHS is 1 in most nonlinear problems: it is 2 for the modified Riks static procedure (“Static stress analysis,” Section 6.2.2 of the Abaqus Analysis User’s Guide), and it is greater than 1 in some linear analysis procedures and during substructure generation.
+
+
+
+# NSVARS
+
+User-defined number of solution-dependent state variables associated with the element (“Defining the number of solution-dependent variables that must be stored within the element” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# NPROPS
+
+User-defined number of real property values associated with the element (“Defining the element properties” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# NJPROP
+
+User-defined number of integer property values associated with the element (“Defining the element properties” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# MCRD
+
+MCRD is defined as the maximum of the user-defined maximum number of coordinates needed at any node point (“Defining the maximum number of coordinates needed at any nodal point” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide) and the value of the largest active degree of freedom of the user element that is less than or equal to 3. For example, if you specify that the maximum number of coordinates is 1 and the active degrees of freedom of the user element are 2, 3, and 6, MCRD will be 3. If you specify that the maximum number of coordinates is 2 and the active degrees of freedom of the user element are 11 and 12, MCRD will be 2.
+
+# NNODE
+
+User-defined number of nodes on the element (“Defining the number of nodes associated with the element” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# JTYPE
+
+Integer defining the element type. This is the user-defined integer value n in element type Un (“Assigning an element type key to a user-defined element” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# KSTEP
+
+Current step number.
+
+# KINC
+
+Current increment number.
+
+# JELEM
+
+User-assigned element number.
+
+# NDLOAD
+
+Identification number of the distributed load or flux currently active on this element.
+
+# MDLOAD
+
+Total number of distributed loads and/or fluxes defined on this element.
+
+
+
+# NPREDF
+
+Number of predefined field variables, including temperature. For user elements Abaqus/Standard uses one value for each field variable per node.
+
+# UEL conventions
+
+The solution variables (displacement, velocity, etc.) are arranged on a node/degree of freedom basis. The degrees of freedom of the first node are first, followed by the degrees of freedom of the second node, etc.
+
+# Usage with general nonlinear procedures
+
+The values of $u ^ { N }$ (and, in direct-integration dynamic steps, $\dot { u } ^ { N }$ and $\ddot { u } ^ { N } )$ enter user subroutine UEL as their latest approximations at the end of the time increment; that is, at time $t + \Delta t$ .
+
+The values of $H ^ { \alpha }$ enter the subroutine as their values at the beginning of the time increment; that is, at time t. It is your responsibility to define suitable time integration schemes to update $H ^ { \alpha }$ . To ensure accurate, stable integration of internal state variables, you can control the time incrementation via PNEWDT.
+
+The values of $\cdot _ { p ^ { \beta } }$ enter the subroutine as the values of the total load magnitude for the th distributed load at the end of the increment. Increments in the load magnitudes are also available.
+
+In the following descriptions of the user element’s requirements, it will be assumed that LFLAGS(3)=1 unless otherwise stated.
+
+# Static analysis (LFLAGS(1)=1,2)
+
+• $F ^ { N } = F ^ { N } ( u ^ { M } , H ^ { \alpha } , p ^ { \beta } , t )$ .
+• Automatic convergence checks are applied to the force residuals corresponding to degrees of freedom 1–7.
+• You must define AMATRX $= K ^ { N M } = - \partial F ^ { N } / \partial u ^ { M }$ and $\mathbf { R H S } = F ^ { N }$ and update the state variables, $H ^ { \alpha }$ .
+
+# Modified Riks static analysis (LFLAGS(1)=1) and (NRHS=2)
+
+• $F ^ { N } = F ^ { N } ( u ^ { M } , H ^ { \alpha } , p ^ { \beta } )$ , where $p ^ { \beta } = p _ { 0 } ^ { \beta } + \lambda q ^ { \beta } , p _ { 0 } ^ { \beta }$ and $q ^ { \beta }$ are fixed load parameters, and is the Riks (scalar) load parameter.
+• Automatic convergence checks are applied to the force residuals corresponding to degrees of freedom 1–7.
+• You must define $\begin{array} { l l l } { { \tt A M A T R X } } & { = } & { K ^ { N M } } & { = } & { - \partial F ^ { N } / \partial u ^ { M } } \end{array}$ , RHS $\mathbf { \Psi } ( \mathbf { 1 } ) \ = \ \mathbf { \Psi } F ^ { N }$ , and RHS(2) $= \Delta \lambda ( \partial { \cal F } ^ { N } / \partial \lambda )$ and update the state variables, $H ^ { \alpha }$ . RHS(2) is the incremental load vector.
+
+# Direct-integration dynamic analysis (LFLAGS(1)=11, 12)
+
+• Automatic convergence checks are applied to the force residuals corresponding to degrees of freedom 1–7.
+
+
+
+• LFLAGS(3)=1: Normal time increment. Either the Hilber-Hughes-Taylor or the backward Euler time integration scheme will be used. With set to zero for the backward Euler, both schemes imply
+
+$$
+F ^ {N} = - M ^ {N M} \ddot {u} _ {t + \Delta t} + (1 + \alpha) G _ {t + \Delta t} ^ {N} - \alpha G _ {t} ^ {N},
+$$
+
+where $M ^ { N M } = M ^ { N M } ( u ^ { M } , \dot { u } ^ { M } , H ^ { \alpha } , p ^ { \beta } , t , . . . )$ and $G ^ { N } = G ^ { N } ( u ^ { M } , \dot { u } ^ { M } , H ^ { \alpha } , p ^ { \beta } , t , . . . )$ ; that is, the highest time derivative of $u ^ { M }$ in $M ^ { N M }$ and $G ^ { N }$ is $\dot { u } ^ { M }$ , so that
+
+$$
+- \frac {\partial F ^ {N}}{\partial \ddot {u} ^ {M} {} _ {t + \Delta t}} = M ^ {N M}.
+$$
+
+Therefore, you must store $G _ { t } ^ { N }$ as an internal state vector. If half-increment residual calculations are required, you must also store $G _ { t ^ { - } } ^ { N }$ as an internal state vector, where indicates the time at the beginning of the previous increment. For $\alpha = 0 , F ^ { N } = - M ^ { N M } \ddot { u } _ { t + \Delta t } + G ^ { N } { } _ { t + \Delta t }$ and $G _ { t } ^ { N }$ is not needed. You must define $\mathtt { A M A T R X } = M ^ { N M } \left( d \ddot { u } / d u \right) + \left( 1 + \alpha \right) C ^ { N M } \left( d \dot { u } / d u \right) + \left( 1 + \alpha \right) K ^ { N M }$ where $C ^ { N M } = - \partial G ^ { N } { } _ { t + \Delta t } / \partial \dot { u } ^ { M }$ and $K ^ { N M } = - \partial G ^ { N } { } _ { t + \Delta t } / \partial u ^ { M }$ . $\mathbf { R H S } = \mathbf { \nabla } F ^ { N }$ must also be defined and the state variables, $H ^ { \alpha }$ , updated. Although the value of given in the dynamic step definition is passed into UEL, the value of can vary from element to element. For example, can be set to zero for some elements in the model where numerical dissipation is not desired.
+
+• LFLAGS(3)=5: Half-increment residual $( F _ { 1 / 2 } ^ { N } )$ calculation. Abaqus/Standard will adjust the time increment so that $| F _ { 1 / 2 } ^ { N } | ~ < ~ t o l e r a n c e$ (where is specified in the dynamic step definition). The half-increment residual is defined as
+
+$$
+F _ {1 / 2} ^ {N} = - M ^ {N M} \ddot {u} _ {t + \Delta t / 2} + (1 + \alpha) G _ {t + \Delta t / 2} ^ {N} - \frac {\alpha}{2} (G _ {t} ^ {N} + G _ {t -} ^ {N}),
+$$
+
+where $t ^ { - }$ indicates the time at the beginning of the previous increment ( is a parameter of the Hilber-Hughes-Taylor time integration operator and will be set to zero if the backward Euler time integration operator is used). You must define $\mathbf { R H S } = F _ { 1 / 2 } ^ { N }$ . To evaluate $M ^ { N M }$ and ${ G ^ { N } } _ { t + \Delta t / 2 } , \mathrm { y o u }$ must calculate ${ H ^ { \alpha } } _ { t + \Delta t / 2 }$ . These half-increment values will not be saved. DTIME will still contain $\Delta t \left( { \mathrm { n o t } } \Delta t / 2 \right)$ . The values contained in U, V, A, and DU are half-increment values.
+
+• LFLAGS(3)=4: Velocity jump calculation. Abaqus/Standard solves $- M ^ { N M } \Delta \dot { u } ^ { M } = 0 \mathrm { f o r } \Delta \dot { u } ^ { M }$ , so you must define AMATRX . $\tt A M A T R X = M ^ { N M }$
+• LFLAGS(3)=6: Initial acceleration calculation. Abaqus/Standard solves $- M ^ { N M } \ddot { u } ^ { M } + G ^ { N } = 0$ for $\ddot { u } ^ { M }$ , so you must define $\mathtt { A M A T R X } = M ^ { N M }$ and $\mathbf { R } \mathbf { \bar { H } } \mathbf { S } = G ^ { N }$ .
+
+# Subspace-based dynamic analysis (LFLAGS(1)=13)
+
+• The requirements are identical to those of static analysis, except that the Jacobian (stiffness), AMATRX, is not needed. No convergence checks are performed in this case.
+
+# Quasi-static analysis (LFLAGS(1)=21)
+
+• The requirements are identical to those of static analysis.
+
+
+
+# Steady-state heat transfer analysis (LFLAGS(1)=31)
+
+• The requirements are identical to those of static analysis, except that the automatic convergence checks are applied to the heat flux residuals corresponding to degrees of freedom 11, 12, …
+
+# Transient heat transfer analysis $\left( \phantom { - } \theta _ { m a x } \right) \left( \tt L F L A G S \left( 1 \right) = 3 2 , \phantom { - } 3 3 \right)$
+
+• Automatic convergence checks are applied to the heat flux residuals corresponding to degrees of freedom 11, 12, …
+• The backward difference scheme is always used for time integration; that is, Abaqus/Standard assumes that $\dot { u } _ { t + \Delta t } = \Delta u / \Delta t .$ , where $\Delta u = u _ { t + \Delta t } - u _ { t }$ and so $d \dot { u } / d u = 1 / \Delta t$ always. For degrees of freedom 11, 12, …, $\lvert \Delta u \rvert$ will be compared against the user-prescribed maximum allowable nodal temperature change in an increment, $\Delta \theta _ { m a x }$ , for controlling the time integration accuracy.
+• You need to define $\mathtt { A M A T R X } = K ^ { N M } + ( 1 / \Delta t ) \ C ^ { N M }$ , where $C ^ { N M }$ is the heat capacity matrix and $\mathbf { R } \mathbf { \bar { H } } \mathbf { S } = F ^ { N }$ , and must update the state variables, $H ^ { \alpha }$ .
+
+# Geostatic analysis (LFLAGS(1)=61)
+
+• Identical to static analysis, except that the automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–8.
+
+# Steady-state coupled pore fluid diffusion/stress analysis (LFLAGS(1)=62, 63)
+
+• Identical to static analysis, except that the automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–8.
+
+# Transient coupled pore fluid diffusion/stress (consolidation) analysis $( ~ { u } _ { w } ^ { m a x } )$ (LFLAGS(1)=64, 65)
+
+• Automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–8.
+• The backward difference scheme is used for time integration; that is, $\dot { u } _ { t + \Delta t } ^ { M } = \Delta u ^ { M } / \Delta t$ , where $\Delta u ^ { M } = u _ { t + \Delta t } ^ { M } - u _ { t } ^ { M }$ .
+• For degree of freedom 8, $| \Delta u ^ { M } |$ will be compared against the user-prescribed maximum wetting liquid pore pressure change, $\Delta u _ { w } ^ { m a x }$ , for automatic control of the time integration accuracy.
+• You must define $\mathtt { A M A T R X } = K ^ { N M } + ( 1 / \Delta t ) \ C ^ { N M }$ , where $C ^ { N M }$ is the pore fluid capacity matrix and $\mathbf { R H S } = F ^ { N }$ , and must update the state variables, $H ^ { \alpha }$ .
+
+# Steady-state fully coupled thermal-stress analysis (LFLAGS(1)=71)
+
+• Identical to static analysis, except that the automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–7 and 11, 12, …
+
+
+
+Transient fully coupled thermal-stress analysis $( \theta _ { m a x } )$ (LFLAGS(1)=72,73)
+
+• Automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–7 and 11, 12, …
+• The backward difference scheme is used for time integration; that is, $\dot { u } _ { t + \Delta t } ^ { M } = \Delta u ^ { M } / \Delta t .$ where $\Delta u ^ { M } = u _ { t + \Delta t } ^ { M } - u _ { t } ^ { M }$ .
+• For degrees of freedom 11, 12, …, $| \Delta u ^ { M } |$ will be compared against the user-prescribed maximum allowable nodal temperature change in an increment, $\Delta \theta _ { m a x }$ , for automatic control of the time integration accuracy.
+• You must define $\mathtt { A M A T R X } = K ^ { N M } + ( 1 / \Delta t ) \ C ^ { N M }$ , where $C ^ { N M }$ is the heat capacity matrix and $\mathbf { R } \mathbf { \bar { H } } \mathbf { S } = F ^ { N }$ , and must update the state variables, $H ^ { \alpha }$ .
+
+# Steady-state coupled thermal-electrical analysis (LFLAGS(1)=75)
+
+• The requirements are identical to those of static analysis, except that the automatic convergence checks are applied to the current density residuals corresponding to degree of freedom 9, in addition to the heat flux residuals.
+
+Transient coupled thermal-electrical analysis $\left( \phantom { - } \theta _ { m a x } \right) \left( \tt L F L A G S \left( 1 \right) = 7 6 , \eta 7 7 \right)$
+
+• Automatic convergence checks are applied to the current density residuals corresponding to degree of freedom 9 and to the heat flux residuals corresponding to degree of freedom 11.
+• The backward difference scheme is always used for time integration; that is, Abaqus/Standard assumes that $\dot { u } _ { t + \Delta t } = \Delta u / \Delta t .$ , where $\Delta u = u _ { t + \Delta t } - u _ { t }$ . Therefore, $d \dot { u } / d u = 1 / \Delta t$ always. For degree of freedom 11 $| \Delta u |$ will be compared against the user-prescribed maximum allowable nodal temperature change in an increment, $\Delta \theta _ { m a x }$ , for controlling the time integration accuracy.
+• You must define $\mathtt { A M A T R X } = K ^ { N M } + ( 1 / \Delta t ) \ C ^ { N M }$ , where $C ^ { N M }$ is the heat capacity matrix and $\mathbf { R } \mathbf { \bar { H } } \mathbf { S } = F ^ { N }$ , and must update the state variables, $H ^ { \alpha }$ .
+
+Steady-state coupled thermal-electrical-structural analysis (LFLAGS(1)=102)
+
+• Identical to static analysis, except that the automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–7, 9, and 11.
+
+Transient coupled thermal-electrical-structural analysis $( \theta _ { m a x } )$ (LFLAGS(1)=103,104)
+
+• Automatic convergence checks are applied to the residuals corresponding to degrees of freedom 1–7, 9, and 11.
+• The backward difference scheme is always used for time integration; that is, Abaqus/Standard assumes that $\dot { u } _ { t + \Delta t } = \Delta u / \Delta t .$ , where $\Delta u = u _ { t + \Delta t } - u _ { t }$ . Therefore, $d \dot { u } / d u = 1 / \Delta t$ always. For degree of freedom 11 $| \Delta u |$ will be compared against the user-prescribed maximum allowable nodal temperature change in an increment, $\Delta \theta _ { m a x }$ , for controlling the time integration accuracy.
+• You must define $\mathtt { A M A T R X } = K ^ { N M } + ( 1 / \Delta t ) \ C ^ { N M }$ , where $C ^ { N M }$ is the heat capacity matrix and $\mathbf { R } \mathbf { \bar { H } } \mathbf { S } = F ^ { N }$ ; and you must update the state variables, $H ^ { \alpha }$ .
+
+
+
+“General and linear perturbation procedures,” Section 6.1.3 of the Abaqus Analysis User’s Guide, describes the linear perturbation capabilities in Abaqus/Standard. Here, base state values of variables will be denoted by $u ^ { M } , H ^ { \alpha }$ , etc. Perturbation values will be denoted by $\tilde { u } ^ { M } , \tilde { H } ^ { \alpha }$ , etc.
+
+Abaqus/Standard will not call user subroutine UEL for the eigenvalue buckling prediction procedure.
+
+For response spectrum, random response, transient modal dynamic, and mode-based steady-state dynamic procedures, user subroutine UEL is called only in a prior natural frequency extraction analysis, and the mass and stiffness contributions are taken into account during modal superposition.
+
+For direct-solution and mode-based steady-state dynamic, complex eigenvalue extraction, matrix generation, and substructure generation procedures, Abaqus/Standard will call user subroutine UEL, but only mass and stiffness contributions will be taken into account. The damping contributions will be neglected.
+
+# Static analysis (LFLAGS(1)=1, 2)
+
+• Abaqus/Standard will solve $K ^ { N M } \tilde { u } ^ { M } = \tilde { P } ^ { N }$ for $\tilde { u } ^ { M }$ , where $K ^ { N M }$ is the base state stiffness matrix and the perturbation load vector, $\tilde { P } ^ { N }$ , is a linear function of the perturbation loads, $\tilde { p } ;$ that is, $\tilde { P } ^ { N } =$ $\left( { \partial F } / { \partial \tilde { p } } \right) \tilde { p } .$ .
+• $\mathtt { L F L A G S } \left( 3 \right) = 1$ : You must define $\mathbf { A M A T R X } = K ^ { N M }$ and $\mathbf { R } \mathbf { \bar { H } } \mathbf { S } = \mathbf { \tilde { \Lambda } } \tilde { P } ^ { N }$ .
+• $\mathtt { L F L A G S } \left( 3 \right) = 1 0 0 $ : You must compute perturbations of the internal variables, $\tilde { H } ^ { \alpha }$ , and define RHS $= \tilde { P } ^ { N } - K ^ { N M } \tilde { u } ^ { M } $ for output purposes.
+
+# Eigenfrequency extraction analysis (LFLAGS(1)=41)
+
+$\bullet F ^ { N } = - M ^ { N M } \ddot { \tilde { u } } + G ^ { N } ( u ^ { M } + \tilde { u } ^ { M } , \ldots ) = - M ^ { N M } \ddot { \tilde { u } } + \left( \partial G ^ { N } / \partial u ^ { M } \right) \tilde { u } ^ { M } .$
+• Abaqus/Standard will solve $\begin{array} { r l r } { K ^ { N M } \phi _ { i } ^ { M } } & { { } = } & { \omega _ { i } ^ { 2 } M ^ { N M } \phi _ { i } ^ { M } } \end{array}$ for $\phi _ { i } ^ { N }$ and $\omega _ { i }$ , where $\begin{array} { r l } { K ^ { N M } } & { { } = } \end{array}$ $- \partial F ^ { N } / \partial u ^ { M }$ is the base state stiffness matrix and $M ^ { N M } = - \partial F ^ { N M } / \partial \ddot { u } ^ { M }$ is the base state mass matrix.
+• LFLAGS(3)=2: Define $\mathbf { A M A T R X } = K ^ { N M }$ .
+• LFLAGS(3)=4: Define $\mathtt { A M A T R X } = M ^ { N M }$ .
+
+# Example: Structural and heat transfer user element
+
+Both a structural and a heat transfer user element have been created to demonstrate the usage of subroutine UEL. These user-defined elements are applied in a number of analyses. The following excerpt is from the verification problem that invokes the structural user element in an implicit dynamics procedure:
+
+```python
+*USER ELEMENT, NODES=2, TYPE=U1, PROPERTIES=4, COORDINATES=3, VARIABLES=12
+1, 2, 3
+*ELEMENT, TYPE=U1
+101, 101, 102
+```
+
+
+
+```csv
+*ELGEN, ELSET=UTRUSS
+101, 5
+*UEL PROPERTY, ELSET=UTRUSS
+0.002, 2.1E11, 0.3, 7200.
+```
+
+The user element consists of two nodes that are assumed to lie parallel to the x-axis. The element behaves like a linear truss element. The supplied element properties are the cross-sectional area, Young’s modulus, Poisson’s ratio, and density, respectively.
+
+The next excerpt shows the listing of the subroutine. The user subroutine has been coded for use in a perturbation static analysis; general static analysis, including Riks analysis with load incrementation defined by the subroutine; eigenfrequency extraction analysis; and direct-integration dynamic analysis. The names of the verification input files associated with the subroutine and these procedures can be found in “UEL,” Section 4.1.14 of the Abaqus Verification Guide. The subroutine performs all calculations required for the relevant procedures as described earlier in this section. The flags passed in through the LFLAGS array are used to associate particular calculations with solution procedures.
+
+During a modified Riks analysis all force loads must be passed into UEL by means of distributed load definitions such that they are available for the definition of incremental load vectors; the load keys Un and UnNU must be used properly, as discussed in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide. The coding in subroutine UEL must distribute the loads into consistent equivalent nodal forces and account for them in the calculation of the RHS and ENERGY arrays.
+
+```csv
+SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
+1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,
+2 DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
+3 PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,
+4 JPROPS,NJPROP,PERIOD)
+C
+INCLUDE 'ABA_PARAM.INC'
+PARAMETER ( ZERO = 0.D0, HALF = 0.5D0, ONE = 1.D0 )
+C
+DIMENSION RHS(MLVARX, * ),AMATRX(NDOFEL,NDOFEL),
+1 SVARS(NSVARS),ENERGY(8),PROPS(*),COORDS(MCRD,NNODE),
+2 U(NDOFEL),DU(MLVARX, * ),V(NDOFEL),A(NDOFEL),TIME(2),
+3 PARAMS(3),JDLTYP(MDLOAD, * ),ADLMAG(MDLOAD, * ),
+4 DDLMAG(MDLOAD, * ),PREDEF(2,NPREDF,NNODE),LFLAGS(*),
+5 JPROPS(*)
+DIMENSION SRESID(6)
+C
+C UEL SUBROUTINE FOR A HORIZONTAL TRUSS ELEMENT
+C
+C SRESID - stores the static residual at time t+dt
+C SVARS - In 1-6, contains the static residual at time t
+C upon entering the routine. SRESID is copied to
+```
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_021.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_021.md
new file mode 100644
index 0000000..170585e
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_021.md
@@ -0,0 +1,332 @@
+
+
+```txt
+C SVARS(1-6) after the dynamic residual has been
+C calculated.
+C - For half-increment residual calculations: In 7-12,
+C contains the static residual at the beginning
+C of the previous increment. SVARS(1-6) are copied
+C into SVARS(7-12) after the dynamic residual has
+C been calculated.
+C
+AREA = PROPS(1)
+E = PROPS(2)
+ANU = PROPS(3)
+RHO = PROPS(4)
+C
+ALEN = ABS(COORDS(1,2)-COORDS(1,1))
+AK = AREA*E/ALEN
+AM = HALF*AREA*RHO*ALEN
+C
+DO K1 = 1, NDOFEL
+ SRESID(K1) = ZERO
+ DO KRHS = 1, NRHS
+ RHS(K1,KRHS) = ZERO
+ END DO
+ DO K2 = 1, NDOFEL
+ AMATRX(K2,K1) = ZERO
+ END DO
+ END DO
+C
+IF (LFLAGS(3).EQ.1) THEN
+C Normal incrementation
+ IF (LFLAGS(1).EQ.1 .OR. LFLAGS(1).EQ.2) THEN
+C *STATIC
+ AMATRX(1,1) = AK
+ AMATRX(4,4) = AK
+ AMATRX(1,4) = -AK
+ AMATRX(4,1) = -AK
+ IF (LFLAGS(4).NE.0) THEN
+ FORCE = AK*(U(4)-U(1))
+ DFORCE = AK*(DU(4,1)-DU(1,1))
+ SRESID(1) = -DFORCE
+ SRESID(4) = DFORCE
+ RHS(1,1) = RHS(1,1)-SRESID(1)
+ RHS(4,1) = RHS(4,1)-SRESID(4)
+```
+
+
+
+```matlab
+ENERGY(2) = HALF*FORCE*(DU(4,1)-DU(1,1))
+* + HALF*DFORCE*(U(4)-U(1))
+* + HALF*DFORCE*(DU(4,1)-DU(1,1))
+ELSE
+FORCE = AK*(U(4)-U(1))
+SRESID(1) = -FORCE
+SRESID(4) = FORCE
+RHS(1,1) = RHS(1,1)-SRESID(1)
+RHS(4,1) = RHS(4,1)-SRESID(4)
+DO KDLOAD = 1, NDLOAD
+IF (JDLTYP(KDLOAD,1).EQ.1001) THEN
+RHS(4,1) = RHS(4,1)+ADLMAG(KDLOAD,1)
+ENERGY(8) = ENERGY(8)+(ADLMAG(KDLOAD,1)
+* - HALF*DDLMAG(KDLOAD,1))*DU(4,1)
+IF (NRHS.EQ.2) THEN
+C
+Riks
+RHS(4,2) = RHS(4,2)+DDLMAG(KDLOAD,1)
+END IF
+END IF
+END DO
+ENERGY(2) = HALF*FORCE*(U(4)-U(1))
+END IF
+ELSE IF (LFLAGS(1).EQ.11 .OR. LFLAGS(1).EQ.12) THEN
+C
+*DYNAMIC
+ALPHA = PARAMS(1)
+BETA = PARAMS(2)
+GAMMA = PARAMS(3)
+C
+DADU = ONE/(BETA*DTIME**2)
+DVDU = GAMMA/(BETA*DTIME)
+C
+DO K1 = 1, NDOFEL
+AMATRX(K1,K1) = AM*DADU
+RHS(K1,1) = RHS(K1,1)-AM*A(K1)
+END DO
+AMATRX(1,1) = AMATRX(1,1)+(ONE+ALPHA)*AK
+AMATRX(4,4) = AMATRX(4,4)+(ONE+ALPHA)*AK
+AMATRX(1,4) = AMATRX(1,4)-(ONE+ALPHA)*AK
+AMATRX(4,1) = AMATRX(4,1)-(ONE+ALPHA)*AK
+FORCE = AK*(U(4)-U(1))
+SRESID(1) = -FORCE
+SRESID(4) = FORCE
+```
+
+
+
+```vba
+RHS(1,1) = RHS(1,1) -
+* ((ONE+ALPHA)*SRESID(1)-ALPHA*SVARS(1))
+RHS(4,1) = RHS(4,1) -
+* ((ONE+ALPHA)*SRESID(4)-ALPHA*SVARS(4))
+ENERGY(1) = ZERO
+DO K1 = 1, NDOFEL
+SVARS(K1+6) = SVARS(k1)
+SVARS(K1) = SRESID(K1)
+ENERGY(1) = ENERGY(1)+HALF*V(K1)*AM*V(K1)
+END DO
+ENERGY(2) = HALF*FORCE*(U(4)-U(1))
+END IF
+ELSE IF (LFLAGS(3).EQ.2) THEN
+C Stiffness matrix
+AMATRX(1,1) = AK
+AMATRX(4,4) = AK
+AMATRX(1,4) = -AK
+AMATRX(4,1) = -AK
+ELSE IF (LFLAGS(3).EQ.4) THEN
+C Mass matrix
+DO K1 = 1, NDOFEL
+AMATRX(K1,K1) = AM
+END DO
+ELSE IF (LFLAGS(3).EQ.5) THEN
+C Half-increment residual calculation
+ALPHA = PARAMS(1)
+FORCE = AK*(U(4)-U(1))
+SRESID(1) = -FORCE
+SRESID(4) = FORCE
+RHS(1,1) = RHS(1,1)-AM*A(1)-(ONE+ALPHA)*SRESID(1)
+* + HALF*ALPHA*(SVARS(1)+SVARS(7))
+RHS(4,1) = RHS(4,1)-AM*A(4)-(ONE+ALPHA)*SRESID(4)
+* + HALF*ALPHA*(SVARS(4)+SVARS(10))
+ELSE IF (LFLAGS(3).EQ.6) THEN
+C Initial acceleration calculation
+DO K1 = 1, NDOFEL
+AMATRX(K1,K1) = AM
+END DO
+FORCE = AK*(U(4)-U(1))
+SRESID(1) = -FORCE
+SRESID(4) = FORCE
+RHS(1,1) = RHS(1,1)-SRESID(1)
+```
+
+
+
+```matlab
+RHS(4,1) = RHS(4,1) - SRESID(4)
+ENERGY(1) = ZERO
+DO K1 = 1, NDOFEL
+SVARS(K1) = SRESID(K1)
+ENERGY(1) = ENERGY(1) + HALF*V(K1) * AM*V(K1)
+END DO
+ENERGY(2) = HALF*FORCE*(U(4) - U(1))
+ELSE IF (LFLAGS(3).EQ.100) THEN
+C Output for perturbations
+IF (LFLAGS(1).EQ.1 .OR. LFLAGS(1).EQ.2) THEN
+C *STATIC
+FORCE = AK*(U(4) - U(1))
+DFORCE = AK*(DU(4,1) - DU(1,1))
+SRESID(1) = -DFORCE
+SRESID(4) = DFORCE
+RHS(1,1) = RHS(1,1) - SRESID(1)
+RHS(4,1) = RHS(4,1) - SRESID(4)
+ENERGY(2) = HALF*FORCE*(DU(4,1) - DU(1,1))
+* + HALF*DFORCE*(U(4) - U(1))
+* + HALF*DFORCE*(DU(4,1) - DU(1,1))
+DO KVAR = 1, NSVARS
+SVARS(KVAR) = ZERO
+END DO
+SVARS(1) = RHS(1,1)
+SVARS(4) = RHS(4,1)
+ELSE IF (LFLAGS(1).EQ.41) THEN
+C *FREQUENCY
+DO KRHS = 1, NRHS
+DFORCE = AK*(DU(4, KRHS) - DU(1, KRHS))
+SRESID(1) = -DFORCE
+SRESID(4) = DFORCE
+RHS(1, KRHS) = RHS(1, KRHS) - SRESID(1)
+RHS(4, KRHS) = RHS(4, KRHS) - SRESID(4)
+END DO
+DO KVAR = 1, NSVARS
+SVARS(KVAR) = ZERO
+END DO
+SVARS(1) = RHS(1,1)
+SVARS(4) = RHS(4,1)
+END IF
+END IF
+C
+```
+
+
+
+RETURN
+END
+
+
+
+
+
+# 1.1.29 UELMAT: User subroutine to define an element with access to Abaqus materials.
+
+# Product: Abaqus/Standard
+
+WARNING: This feature is intended for advanced users only. Its use in all but the simplest test examples will require considerable coding by the user/developer. “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide, should be read before proceeding.
+
+# References
+
+• “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide
+• \*UEL PROPERTY
+• \*USER ELEMENT
+• “Accessing Abaqus materials,” Section 2.1.17
+• “Accessing Abaqus thermal materials,” Section 2.1.18
+
+# Overview
+
+# User subroutine UELMAT:
+
+• will be called for each element that is of a general user-defined element type (i.e., not defined by a linear stiffness or mass matrix read either directly or from results file data) each time element calculations are required;
+• (or subroutines called by user subroutine UELMAT) must perform all of the calculations for the element, appropriate to the current activity in the analysis;
+• can access some of the Abaqus materials through utility routines MATERIAL\_LIB\_MECH and MATERIAL\_LIB\_HT;
+• is available for a subset of the procedures supported for user subroutine UEL (see “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide); and
+• is available for plane stress and three-dimensional element types in a stress/displacement analysis and for two-dimensional and three-dimensional element types in a heat transfer analysis (see “Userdefined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# User subroutine interface
+
+SUBROUTINE UELMAT(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
+1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
+2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,
+3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD,
+4 MATERIALLIB)
+
+
+
+C
+
+INCLUDE 'ABA\_PARAM.INC'
+
+C
+
+DIMENSION RHS(MLVARX,\*),AMATRX(NDOFEL,NDOFEL),PROPS(\*),
+
+1 SVARS(\*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),
+
+2 DU(MLVARX,\*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(\*),
+
+3 JDLTYP(MDLOAD,\*),ADLMAG(MDLOAD,\*),DDLMAG(MDLOAD,\*),
+
+4 PREDEF(2,NPREDF,NNODE),LFLAGS(\*),JPROPS(\*)
+
+user coding to define RHS, AMATRX, SVARS, ENERGY, and PNEWDT
+
+RETURN
+
+END
+
+# Variables to be defined
+
+These arrays depend on the value of the LFLAGS array.
+
+# RHS
+
+An array containing the contributions of this element to the right-hand-side vectors of the overall system of equations. For most nonlinear analysis procedures, NRHS=1 and RHS should contain the residual vector. The exception is the modified Riks static procedure (“Static stress analysis,” Section 6.2.2 of the Abaqus Analysis User’s Guide), for which NRHS=2 and the first column in RHS should contain the residual vector and the second column should contain the increments of external load on the element. RHS(K1,K2) is the entry for the K1th degree of freedom of the element in the K2th right-hand-side vector.
+
+# AMATRX
+
+An array containing the contribution of this element to the Jacobian (stiffness) or other matrix of the overall system of equations. The particular matrix required at any time depends on the entries in the LFLAGS array (see below).
+
+All nonzero entries in AMATRX should be defined, even if the matrix is symmetric. If you do not specify that the matrix is unsymmetric when you define the user element, Abaqus/Standard will use the symmetric matrix defined by ${ \mathsf { \Omega } } _ { 2 } ^ { 1 } ( [ A ] + [ A ] ^ { T } )$ , where is the matrix defined as AMATRX in this subroutine. If you specify that the matrix is unsymmetric when you define the user element, Abaqus/Standard will use AMATRX directly.
+
+# SVARS
+
+An array containing the values of the solution-dependent state variables associated with this element. The number of such variables is NSVARS (see below). You define the meaning of these variables.
+
+For general nonlinear steps this array is passed into UELMAT containing the values of these variables at the start of the current increment. They should be updated to be the values at the end
+
+
+
+of the increment, unless the procedure during which UELMAT is being called does not require such an update; this requirement depends on the entries in the LFLAGS array (see below). For linear perturbation steps this array is passed into UELMAT containing the values of these variables in the base state. They should be returned containing perturbation values if you wish to output such quantities.
+
+When KINC is equal to zero, the call to UELMAT is made for zero increment output (see “Output,” Section 4.1.1 of the Abaqus Analysis User’s Guide). In this case the values returned will be used only for output purposes and are not updated permanently.
+
+# ENERGY
+
+For general nonlinear steps array ENERGY contains the values of the energy quantities associated with the element. The values in this array when UELMAT is called are the element energy quantities at the start of the current increment. They should be updated to the values at the end of the current increment. For linear perturbation steps the array is passed into UELMAT containing the energy in the base state. They should be returned containing perturbation values if you wish to output such quantities. The entries in the array are as follows:
+
+| Element type | Force and moment components |
| Three-dimensional shells (S4R, S8R, S8R5, etc.) and axisymmetric shells with asymmetric deformation (SAXA1N, SAXA2N) | $N_{11}, N_{22}, N_{12}, M_{11}, M_{22}, M_{12}$ |
| Axisymmetric shells (SAX1, SAX2, etc) | $N_{11}, N_{22}, M_{11}, M_{22}$ |
+
+There are NDI direct membrane and NSHR shear membrane components and NDI direct bending and NSHR shear bending components: a total of NSECV components. The order of the components is defined in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide.
+
+Engineering measures of shear membrane strain $( \gamma _ { 1 2 } )$ and twist $\left( \mathrm { K } _ { 1 2 } \right)$ are used.
+
+
+
+# Increments for which only the section stiffness can be defined
+
+Abaqus/Standard passes zero strain increments into user subroutine UGENS to start the first increment of all the steps and all increments of steps for which you have suppressed extrapolation in time from the previous incremental solution (“Defining an analysis,” Section 6.1.2 of the Abaqus Analysis User’s Guide). In this case you can define only the section stiffness (DDNDDE).
+
+# Stability
+
+You should ensure that the integration scheme coded in this routine is stable—no direct provision is made to include a stability limit in the time stepping scheme based on the calculations in UGENS.
+
+# Convergence rate
+
+DDNDDE must be defined accurately if rapid convergence of the overall Newton scheme is to be achieved. In most cases the accuracy of this definition is the most important factor governing the convergence rate. Unsymmetric equation solution is as much as four times as expensive as the corresponding symmetric system. Therefore, if the section stiffness matrix (DDNDDE) is only slightly unsymmetric, it may be computationally less expensive to use a symmetric approximation and accept a slightly slower rate of convergence.
+
+# Use with shells that have transverse shear and/or hourglass stiffness
+
+If user subroutine UGENS is used to describe the section behavior of shells with transverse shear, you must define the transverse shear stiffness (see “Defining the transverse shear stiffness” in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide).
+
+If user subroutine UGENS is used to describe the section behavior of shells with hourglass stiffness, you must define the hourglass stiffness parameter for hourglass control based on total stiffness (see “Specifying nondefault hourglass control parameters for reduced-integration shell elements” in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide). The hourglass stiffness parameter is not required for enhanced hourglass control, but you can define a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surface normal).
+
+# Use with continuum shell elements
+
+User subroutine UGENS cannot be used to describe the section behavior of continuum shell elements.
+
+# User subroutine interface
+
+SUBROUTINE UGENS(DDNDDE,FORCE,STATEV,SSE,SPD,PNEWDT,STRAN,
+1 DSTRAN,TSS,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CENAME,NDI,
+2 NSHR,NSECV,NSTATV,PROPS,JPROPS,NPROPS,NJPROP,COORDS,CELENT,
+3 THICK,DFGRD,CURV,BASIS,NOEL,NPT,KSTEP,KINC,NIT,LINPER)
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_027.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_027.md
new file mode 100644
index 0000000..ab3f9e7
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_027.md
@@ -0,0 +1,337 @@
+
+
+```prolog
+C
+INCLUDE 'ABA_PARAM.INC'
+CHARACTER*80 CENAME
+DIMENSION DDNDDE(NSECV,NSECV),FORCE(NSECV),STATEV(NSTATV),
+1 STRAN(NSECV),DSTRAN(NSECV),TSS(2),TIME(2),PREDEF(*),
+2 DPRED(*),PROPS(*),JPROPS(*),COORDS(3),DFGRD(3,3),
+3 CURV(2,2),BASIS(3,3)
+user coding to define DDNDDE, FORCE, STATEV, SSE, PNEWDT
+RETURN
+END
+```
+
+# Variables to be defined
+
+# DDNDDE(NSECV,NSECV)
+
+Section stiffness matrix of the shell section, , where are the section forces and moments on the shell section and are the generalized section strains in the shell. DDNDDE(I,J) defines the change in the Ith force component at the end of the time increment caused by an infinitesimal perturbation of the Jth component of the section strain increment array. The size of this matrix depends on the values of NSECV (see below for details).
+
+Unless you invoke the unsymmetric equation solution capability in the general shell section definition (“Defining whether or not the section stiffness matrices are symmetric” in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide), Abaqus/Standard will use only the symmetric part of DDNDDE. The symmetric part of the matrix is calculated by taking one half the sum of the matrix and its transpose.
+
+# FORCE(NSECV)
+
+This array is passed in as the forces and moments per unit length on the shell surface at the beginning of the increment and must be updated in this routine to be the forces and moments at the end of the increment.
+
+# STATEV(NSTATV)
+
+An array containing the solution-dependent state variables. These are passed in as the values at the beginning of the increment and must be returned as the values at the end of the increment.
+
+# SSE, SPD
+
+Elastic strain energy and plastic dissipation, respectively. These are passed in as the values at the beginning of the increment and should be updated to the corresponding energy values at the end of the increment. These values have no effect on the solution; they are used for the energy output.
+
+
+
+# PNEWDT
+
+Ratio of suggested new time increment to the time increment being used (DTIME, see below). This variable allows you to provide input to the automatic time incrementation algorithms in Abaqus/Standard (if automatic time incrementation is chosen).
+
+PNEWDT is set to a large value before each call to UGENS.
+
+If PNEWDT is redefined to be less than 1.0, Abaqus/Standard must abandon the time increment and attempt it again with a smaller time increment. The suggested new time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum value for all calls to user subroutines that allow redefinition of PNEWDT for this iteration.
+
+If PNEWDT is given a value that is greater than 1.0 for all calls to user subroutines for this iteration and the increment converges in this iteration, Abaqus/Standard may increase the time increment. The suggested new time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum value for all calls to user subroutines for this iteration.
+
+If automatic time incrementation is not selected in the analysis procedure, values of PNEWDT that are greater than 1.0 will be ignored and values of PNEWDT that are less than 1.0 will cause the job to terminate.
+
+# Variables passed in for information
+
+# STRAN(NSECV)
+
+An array containing the generalized section strains (membrane strains and curvature changes) at the beginning of the increment. The size of this array depends on the value of NSECV (see below for details).
+
+# DSTRAN(NSECV)
+
+Array of generalized section strain increments.
+
+# TSS(2)
+
+Array containing the transverse shear strains.
+
+# TIME(1)
+
+Value of step time at the beginning of the current increment.
+
+# TIME(2)
+
+Value of total time at the beginning of the current increment.
+
+# DTIME
+
+Time increment.
+
+# TEMP
+
+Temperature at the start of the increment.
+
+
+
+# DTEMP
+
+Increment of temperature.
+
+# PREDEF
+
+Array of interpolated values of predefined field variables at this point at the start of the increment, based on the values read in at the nodes.
+
+# DPRED
+
+Array of increments of predefined field variables.
+
+# CENAME
+
+User-specified element set name associated with this section, left justified.
+
+# NDI
+
+Number of direct force components at this point.
+
+# NSHR
+
+Number of shear force components at this point.
+
+# NSECV
+
+Size of the force and strain component arrays.
+
+# NSTATV
+
+User-defined number of solution-dependent state variables associated with this section (“Defining the number of solution-dependent variables that must be stored for the section” in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide).
+
+# PROPS(NPROPS)
+
+A floating point array containing the NPROPS real property values defined for use with this section.
+
+# JPROPS(NJPROP)
+
+An integer array containing the NJPROP integer property values defined for use with this section.
+
+# NPROPS
+
+User-defined number of real property values associated with this section (“Defining the section properties” in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide).
+
+# NJPROP
+
+User-defined number of integer property values associated with the element (“Defining the section properties” in “Using a general shell section to define the section behavior,” Section 29.6.6 of the Abaqus Analysis User’s Guide).
+
+# COORDS
+
+An array containing the current coordinates of this integration point.
+
+
+
+# CELENT
+
+Characteristic element length in the reference surface.
+
+# THICK
+
+Original section thickness.
+
+# DFGRD(3,3)
+
+An array containing the components of the midsurface deformation gradient, $\bar { f } _ { i j }$ . The deformation gradient curvature tensor is available for finite-strain shells (S3/S3R, S4, S4R, SAXs, and SAXAs); it is not available for small-strain shells.
+
+The deformation gradient is stored as a 3 × 3 matrix with component equivalence DFGRD $( \pmb { \tau } , \pmb { \ J } ) \Leftrightarrow$ $\bar { f } _ { i j } . \bar { f } _ { \alpha \beta }$ (Greek subscripts range from 1 to 2) are the in-plane components of the deformation gradient, and $\bar { f } _ { 3 3 }$ is the thickness change component. The components, $\bar { \bar { f } } _ { \alpha 3 }$ , are the transverse shear strains scaled by $\bar { f } _ { 3 3 }$ . The remaining components, $\bar { f } _ { 3 \beta }$ , are all zero.
+
+The tensor is provided in the local shell coordinate system.
+
+# CURV(2,2)
+
+An array containing the midsurface curvature tensor, $b _ { \alpha \beta }$ . The curvature tensor is available for finitestrain shells (S3/S3R, S4, S4R, SAXs, and SAXAs); it is not available for small-strain shells.
+
+The curvature tensor is stored as a 2 × 2 matrix with component equivalence CURV $( \pmb { \tau } , \pmb { \sigma } ) \Leftrightarrow b _ { \alpha \beta }$
+
+The tensor is provided in the local shell coordinate system.
+
+# BASIS(3,3)
+
+An array containing the direction cosines of the shell local surface coordinate system. BASIS(1,1), BASIS(2,1), and BASIS(3,1) give the (1, 2, 3) components of the first local direction, etc. The first two directions are in the plane of the element surface, and the third direction is the normal. The conventions for local directions on shell surfaces are defined in “Conventions,” Section 1.2.2 of the Abaqus Analysis User’s Guide. You can redefine the local system; see “Orientations,” Section 2.2.5 of the Abaqus Analysis User’s Guide.
+
+# NOEL
+
+Element number.
+
+# NPT
+
+Integration point number.
+
+# KSTEP
+
+Step number.
+
+# KINC
+
+Increment number.
+
+# NIT
+
+Iteration number. NIT=0 during the first assembly of the system matrix in any increment.
+
+
+
+# LINPER
+
+Linear perturbation flag. LINPER=1 if the step is a linear perturbation step. LINPER=0 if the step is a general step.
+
+
+
+
+
+# 1.1.39 UHARD: User subroutine to define the yield surface size and hardening parameters for isotropic plasticity or combined hardening models.
+
+# Product: Abaqus/Standard
+
+# References
+
+• “Classical metal plasticity,” Section 23.2.1 of the Abaqus Analysis User’s Guide
+• “Models for metals subjected to cyclic loading,” Section 23.2.2 of the Abaqus Analysis User’s Guide
+• \*CYCLIC HARDENING
+• \*PLASTIC
+
+# Overview
+
+User subroutine UHARD:
+
+• is called at all material calculation points of elements for which the material definition includes user-defined isotropic hardening or cyclic hardening for metal plasticity;
+• can be used to define a material’s isotropic yield behavior;
+• can be used to define the size of the yield surface in a combined hardening model;
+• can include material behavior dependent on field variables or state variables; and
+• requires, when appropriate, that the values of the derivatives of the yield stress (or yield surface size in combined hardening models) be defined with respect to the strain, strain rate, and temperature.
+
+# User subroutine interface
+
+```prolog
+SUBROUTINE UHARD (SYIELD, HARD, EQPLAS, EQPLASRT, TIME, DTIME, TEMP, 1 DTEMP, NOEL, NPT, LAYER, KSPT, KSTEP, KINC, CMNAME, NSTATV, 2 STATEV, NUMFIELDV, PREDEF, DPRED, NUMPROPS, PROPS)
+C
+INCLUDE 'ABA_PARAM.INC'
+C
+CHARACTER*80 CMNAME
+DIMENSION HARD(3), STATEV(NSTATV), TIME(*), $ PREDEF(NUMFIELDV), DPRED(*), PROPS(*)
+user coding to define SYIELD, HARD(1), HARD(2), HARD(3)
+RETURN
+END
+```
+
+
+
+# Variables to be defined
+
+# SYIELD
+
+. Yield stress for isotropic plasticity. Yield surface size for combined hardening.
+
+# HARD(1)
+
+Variation of SYIELD with respect to the equivalent plastic strain,
+
+# HARD(2)
+
+Variation of SYIELD with respect to the equivalent plastic strain rate,
+
+# HARD(3)
+
+Variation of SYIELD with respect to temperature, This quantity is required only in adiabatic, fully coupled temperature-displacement, and thermal-electrical-structural analyses.
+
+# STATEV(NSTATV)
+
+Array containing the user-defined solution-dependent state variables at this point. These are supplied as values at the beginning of the increment or as values updated by other user subroutines (see “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide) and must be returned as values at the end of the increment.
+
+# Variables passed in for information
+
+# EQPLAS
+
+Equivalent plastic strain,
+
+# EQPLASRT
+
+Equivalent plastic strain rate,
+
+# TIME(1)
+
+Value of step time at the beginning of the current increment.
+
+# TIME(2)
+
+Value of total time at the beginning of the current increment.
+
+# DTIME
+
+Time increment.
+
+# TEMP
+
+Temperature at the beginning of the increment.
+
+# DTEMP
+
+Increment of temperature.
+
+# NOEL
+
+Element number.
+
+
+
+# NPT
+
+Integration point number.
+
+# LAYER
+
+Layer number (for composite shells and layered solids).
+
+# KSPT
+
+Section point number within the current layer.
+
+# KSTEP
+
+Step number.
+
+# KINC
+
+Increment number.
+
+# CMNAME
+
+User-specified material name, left justified.
+
+# NSTATV
+
+User-specified number of solution-dependent state variables associated with this material (“Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
+
+# NUMFIELDV
+
+Number of field variables.
+
+# PREDEF(NUMFIELDV)
+
+Array of interpolated values of predefined field variables at this material point at the start of the increment based on the values read in at the nodes (initial values at the beginning of the analysis and current values during the analysis).
+
+# DPRED(NUMFIELDV)
+
+Array of increments of predefined field variables at this material point for this increment; this includes any values updated by user subroutine USDFLD.
+
+# NPROPS
+
+Number of hardening properties entered for this user-defined hardening definition.
+
+# PROPS(NPROPS)
+
+Array of hardening properties entered for this user-defined hardening definition.
+
+
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_028.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_028.md
new file mode 100644
index 0000000..11aef06
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_028.md
@@ -0,0 +1,302 @@
+
+
+# 1.1.40 UHYPEL: User subroutine to define a hypoelastic stress-strain relation.
+
+Product: Abaqus/Standard
+
+# References
+
+• “Hypoelastic behavior,” Section 22.4.1 of the Abaqus Analysis User’s Guide
+• \*HYPOELASTIC
+
+# Overview
+
+User subroutine UHYPEL:
+
+• can be used to define isotropic hypoelastic material behavior, thus requiring the definition of Young’s modulus, E, and Poisson’s ratio, ;
+• is called at all material calculation points of elements for which the material definition contains user-defined hypoelastic behavior;
+• can be used in conjunction with user subroutine USDFLD to redefine any field variables that are passed in (see “USDFLD,” Section 1.1.53); and
+• ignores any data specified outside the user subroutine for the associated hypoelastic material definition.
+
+# Special considerations for various element types
+
+There are several special considerations that need to be noted.
+
+# Beams and shells that calculate transverse shear energy
+
+When UHYPEL is used to define the material response of shell or beam elements that calculate transverse shear energy, Abaqus/Standard cannot calculate a default value for the transverse shear stiffness of the element. Hence, you must define the element’s transverse shear stiffness. See “Shell section behavior,” Section 29.6.4 of the Abaqus Analysis User’s Guide, and “Choosing a beam element,” Section 29.3.3 of the Abaqus Analysis User’s Guide, for guidelines on choosing this stiffness.
+
+# Elements with hourglassing modes
+
+If this capability is used to describe the material of elements with hourglassing modes, you must define the hourglass stiffness for hourglass control based on the total stiffness approach. The hourglass stiffness is not required for enhanced hourglass control, but you can define a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surface normal). See “Section controls,” Section 27.1.4 of the Abaqus Analysis User’s Guide.
+
+
+
+User subroutine interface
+```txt
+SUBROUTINE UHYPEL(E, GNU, STRAIN, NDI, NSHR, EINV1, EINV2, EINV3, 1 COORDS, NOEL, TEMP, PREDEF, CMNAME)
+C
+INCLUDE 'ABA_PARAM.INC'
+C
+CHARACTER*80 CMNAME
+C
+DIMENSION STRAIN(*), COORDS(3), PREDEF(*)
+user coding to define E and GNU
+RETURN
+END
+```
+
+Variables to be defined
+```txt
+E Young's modulus.
+GNU Poisson's ratio.
+```
+
+Variables passed in for information
+STRAIN
+Array containing the total (elastic) strains, ( $\varepsilon$ ).
+NDI
+Number of direct strain components at this point.
+NSHR
+Number of shear strain components at this point.
+
+EINV1 $I_{1} = \mathrm{trace}(\varepsilon)$ , the first strain invariant.
+
+EINV2 $I_{2} = 1 / 2(\varepsilon :\varepsilon -I_{1}^{2})$ , the second strain invariant.
+
+EINV3 $I_{3} = \operatorname{det}(\varepsilon)$ , the third strain invariant.
+
+
+
+# COORDS
+
+An array containing the coordinates of the material point. These are the current coordinates if geometric nonlinearity is accounted for during the step (see “Defining an analysis,” Section 6.1.2 of the Abaqus Analysis User’s Guide); otherwise, the array contains the original coordinates of the point.
+
+# NOEL
+
+Element number.
+
+# TEMP
+
+Current temperature at this point.
+
+# PREDEF
+
+An array containing current values of the predefined field variables at this point (initial values at the beginning of the analysis and current values during the analysis).
+
+# CMNAME
+
+User-specified material name, left justified.
+
+
+
+
+
+# 1.1.41 UHYPER: User subroutine to define a hyperelastic material.
+
+# Product: Abaqus/Standard
+
+# References
+
+• “Hyperelastic behavior of rubberlike materials,” Section 22.5.1 of the Abaqus Analysis User’s Guide
+• \*HYPERELASTIC
+• “UMAT and UHYPER,” Section 4.1.21 of the Abaqus Verification Guide
+
+# Overview
+
+# User subroutine UHYPER:
+
+• can be used to define the strain energy potential for isotropic hyperelastic material behavior;
+• is called at all material calculation points of elements for which the material definition contains user-defined hyperelastic behavior;
+• can include material behavior dependent on field variables or state variables; and
+• requires that the values of the derivatives of the strain energy density function of the hyperelastic material be defined with respect to the strain invariants.
+
+# Special considerations for various element types
+
+There are several special considerations that need to be noted.
+
+# Shells that calculate transverse shear energy
+
+When UHYPER is used to define the material response of shell elements that calculate transverse shear energy, Abaqus/Standard cannot calculate a default value for the transverse shear stiffness of the element. Hence, you must define the element’s transverse shear stiffness. See “Shell section behavior,” Section 29.6.4 of the Abaqus Analysis User’s Guide, for guidelines on choosing this stiffness.
+
+# Elements with hourglassing modes
+
+If this capability is used to describe the material of elements with hourglassing modes, you must define the hourglass stiffness for hourglass control based on the total stiffness approach. The hourglass stiffness is not required for enhanced hourglass control, but you can define a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surface normal). See “Section controls,” Section 27.1.4 of the Abaqus Analysis User’s Guide.
+
+# User subroutine interface
+
+SUBROUTINE UHYPER(BI1,BI2,AJ,U,UI1,UI2,UI3,TEMP,NOEL, 1 CMNAME,INCMPFLAG,NUMSTATEV,STATEV,NUMFIELDV,FIELDV,
+
+
+
+UHYPER
+```txt
+2 FIELDVINC, NUMPROPS, PROPS)
+C
+INCLUDE 'ABA_PARAM.INC'
+C
+CHARACTER*80 CMNAME
+DIMENSION U(2), UI1(3), UI2(6), UI3(6), STATEV(*), FIELDV(*), 2 FIELDVINC(*), PROPS(*)
+user coding to define U, UI1, UI2, UI3, STATEV
+```
+RETURN END
+
+Variables to be defined
+U(1)
+U, strain energy density function. For a compressible material, at least one derivative involving J should be nonzero. For an incompressible material, all derivatives involving J will be ignored. The strain invariants— $\overline{I}_{1}$ , $\overline{I}_{2}$ , and J—are defined in “Hyperelastic behavior of rubberlike materials,” Section 22.5.1 of the Abaqus Analysis User’s Guide.
+
+U(2) $\tilde{U}_{dev}$ , the deviatoric part of the strain energy density of the primary material response. This quantity is needed only if the current material definition also includes Mullins effect (see “Mullins effect,” Section 22.6.1 of the Abaqus Analysis User’s Guide).
+
+UI1(1) $\partial U/\partial\overline{I}_{1}$ .
+
+UI1(2) $\partial U/\partial\overline{I}_{2}$ .
+
+UI1(3) $\partial U/\partial J$ .
+
+UI2(1) $\partial^{2}U/\partial\overline{I}_{1}^{2}$ .
+
+UI2(2) $\partial^{2}U/\partial\overline{I}_{2}^{2}$ .
+
+UI2(3) $\partial^{2}U/\partial J^{2}$ .
+
+
+
+UI2 (4) $\partial^{2}U/\partial\overline{I}_{1}\partial\overline{I}_{2}.$
+
+UI2 (5) $\partial^{2}U/\partial\overline{I}_{1}\partial J.$
+
+UI2 (6) $\partial^{2}U/\partial\overline{I}_{2}\partial J.$
+
+UI3 (1) $\partial^{3}U/\partial\overline{I}_{1}^{2}\partial J.$
+
+UI3 (2) $\partial^{3}U/\partial\overline{I}_{2}^{2}\partial J.$
+
+UI3 (3) $\partial^{3}U/\partial\overline{I}_{1}\partial\overline{I}_{2}\partial J.$
+
+UI3 (4) $\partial^{3}U/\partial\overline{I}_{1}\partial J^{2}.$
+
+UI3 (5) $\partial^{3}U/\partial\overline{I}_{2}\partial J^{2}.$
+
+UI3 (6) $\partial^{3}U/\partial J^{3}.$
+
+# STATEV
+
+Array containing the user-defined solution-dependent state variables at this point. These are supplied as values at the start of the increment or as values updated by other user subroutines (see “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide) and must be returned as values at the end of the increment.
+
+# Variables passed in for information
+
+BI1 $\overline{I}_1$
+
+BI2 $\overline{I}_2$
+
+```txt
+AJ J.
+```
+
+# TEMP
+
+Current temperature at this point.
+
+
+
+# NOEL
+
+Element number.
+
+# CMNAME
+
+User-specified material name, left justified.
+
+# INCMPFLAG
+
+Incompressibility flag defined to be 1 if the material is specified as incompressible or 0 if the material is specified as compressible.
+
+# NUMSTATEV
+
+User-defined number of solution-dependent state variables associated with this material (see “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
+
+# NUMFIELDV
+
+Number of field variables.
+
+# FIELDV
+
+Array of interpolated values of predefined field variables at this material point at the end of the increment based on the values read in at the nodes (initial values at the beginning of the analysis and current values during the analysis).
+
+# FIELDVINC
+
+Array of increments of predefined field variables at this material point for this increment; this includes any values updated by the user subroutine USDFLD.
+
+# NUMPROPS
+
+Number of material properties entered for this user-defined hyperelastic material.
+
+# PROPS
+
+Array of material properties entered for this user-defined hyperelastic material.
+
+
+
+# 1.1.42 UINTER: User subroutine to define surface interaction behavior for contact surfaces.
+
+# Product: Abaqus/Standard
+
+# References
+
+• “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide
+• \*SURFACE INTERACTION
+• “UINTER,” Section 4.1.20 of the Abaqus Verification Guide
+
+# Overview
+
+# User subroutine UINTER:
+
+• is called at points on the slave surface of a contact pair with a user-defined constitutive model defining the interaction between the surfaces;
+• can be used to define the mechanical (normal and shear) and thermal (heat flux) interactions between surfaces;
+• can be used when the normal surface behavior (contact pressure versus overclosure) models (“Contact pressure-overclosure relationships,” Section 37.1.2 of the Abaqus Analysis User’s Guide) or the extended versions of the classical Coulomb friction model (“Frictional behavior,” Section 37.1.5 of the Abaqus Analysis User’s Guide) are too restrictive and a more complex definition of normal and shear transmission between contacting surfaces, including damping properties, are required;
+• must provide the entire definition of the mechanical and the thermal interaction between the contacting surfaces (hence, no additional surface behaviors can be specified in conjunction with this capability);
+• can provide the entire definition of viscous and structural damping for interactions between the contacting surfaces for direct and mode-based steady-state dynamic analysis (including the subspace projection method), transient mode-based analysis, complex eigenvalue extraction, matrix generation, and substructure generation;
+• only accounts for element damping in mode-based procedures if the SIM architecture is used;
+• can use and update solution-dependent state variables; and
+• is not available for contact elements.
+
+# User subroutine interface
+
+SUBROUTINE UINTER(STRESS,DDSDDR,DVISCOUS,DSTRUCTURAL,FLUX,DDFDDT, 1 DDSDDT,DDFDDR,STATEV,SED,SFD,SPD,SVD,SCD,PNEWDT,RDISP,
+
+
+
+```prolog
+2 DRDISP,
+3 TEMP, DTEMP, PREDEF, DPRED, TIME, DTIME, FREQR, CINAME, SLNAME,
+4 MSNAME,
+5 PROPS, COORDS, ALOCALDIR, DROT, AREA, CHRLNGTH, NODE, NDIR, NSTATV,
+6 NPRED, NPROPS, MCRD, KSTEP, KINC, KIT, LINPER, LOPENCLOSE, LSTATE,
+7 LSDI, LPRINT)
+C
+INCLUDE 'ABA_PARAM.INC'
+C
+CHARACTER*80 CINAME, SLNAME, MSNAME
+DIMENSION STRESS (NDIR), DDSDDR (NDIR, NDIR), FLUX (2), DDFDDT (2, 2),
+1 DDSDDT (NDIR, 2), DDFDDR (2, NDIR), STATEV (NSTATV),
+2 RDISP (NDIR), DRDISP (NDIR), TEMP (2), DTEMP (2), PREDEF (2, NPRED),
+3 DPRED (2, NPRED), TIME (2), PROPS (NPROPS), COORDS (MCRD),
+4 ALOCALDIR (3, 3), DROT (2, 2), DVISCOUS (NDIR, NDIR),
+5 DSTRUCTURAL (NDIR, NDIR)
+user coding to define STRESS, DDSDDR, FLUX, DDFDDT,
+DDSDDT, DDFDDR,
+and, optionally, STATEV, SED, SFD, SPD, SVD, SCD, PNEWDT,
+LOPENCLOSE, LSTATE, LSDI, DVISCOUS, DSTRUCTURAL
+RETURN
+END
+```
+
+# Variables to be defined
+
+# STRESS(NDIR)
+
+This array is passed in as the stress between the slave and master surfaces at the beginning of the increment and must be updated in this routine to be the stress at the end of the increment. The stress must be defined in a local coordinate system (see ALOCDIR). This variable must be defined for a stress/displacement, a fully coupled temperature-displacement, or a coupled thermal-electrical-structural analysis. The sign convention for stresses is that a positive stress indicates compression across contact surfaces, while a negative stress indicates tension.
+
+# DDSDDR(NDIR,NDIR)
+
+Interface stiffness matrix. DDSDDR(I,J) defines the change in the Ith stress component at the end of the time increment caused by an infinitesimal perturbation of the Jth component of the relative displacement increment array. Unless you invoke the unsymmetric equation solution capability in the contact property model definition (“Use with the unsymmetric equation solver in Abaqus/Standard” in
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_029.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_029.md
new file mode 100644
index 0000000..bc35c1d
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_029.md
@@ -0,0 +1,299 @@
+
+
+“User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide), Abaqus/Standard will use only the symmetric part of DDSDDR. For a particular off-diagonal (I,J) entry, the symmetrization is done by halving the sum of (I,J) and (J,I) components. DDSDDR must be defined for a stress/displacement, a fully coupled temperature-displacement, or a coupled thermal-electrical-structural analysis to ensure proper convergence characteristics.
+
+# FLUX(2)
+
+Magnitude of the heat flux flowing into the slave and master surfaces, respectively. This array is passed in as the value at the beginning of the increment and must be updated to the flux at the end of the increment. The convention for defining the flux is that a positive flux indicates heat flowing into a surface, while a negative flux indicates heat flowing out of the surface. This variable must be defined for a heat transfer, a fully coupled temperature-displacement, or a coupled thermal-electrical-structural analysis. The sum of these two flux terms represents the heat generated in the interface, and the difference in these flux terms represents the heat conducted through the interface.
+
+# DDFDDT(2,2)
+
+The negative of the variation of the flux at the two surfaces with respect to their respective temperatures, for a fixed relative displacement. This variable must be defined for a heat transfer, a fully coupled temperature-displacement, or a coupled thermal-electrical-structural analysis to ensure proper convergence characteristics. The entries in the first row contain the negatives of the derivatives of FLUX(1) with respect to TEMP(1) and TEMP(2), respectively. The entries in the second row contain the negatives of the corresponding derivatives of FLUX(2).
+
+# DDSDDT(NDIR,2)
+
+Variation of the stress with respect to the temperatures of the two surfaces for a fixed relative displacement. This variable is required only for thermally coupled elements (in a fully coupled temperature-displacement or a coupled thermal-electrical-structural analysis), in which the stress is a function of the surface temperatures. DDSDDT(NDIR,1) corresponds to the slave surface, and DDSDDT(NDIR,2) corresponds to the master surface.
+
+# DDFDDR(2,NDIR)
+
+Variation of the flux with respect to the relative displacement between the two surfaces. This variable is required only for thermally coupled elements (in a fully coupled temperature-displacement or a coupled thermal-electrical-structural analysis), in which the flux is a function of the relative displacement. DDFDDR(1,NDIR) corresponds to the slave surface, and DDFDDR(2,NDIR) corresponds to the master surface.
+
+# Variables that can be updated
+
+# DVISCOUS(NDIR,NDIR)
+
+Interface viscous damping matrix that can be used in direct steady-state dynamic analysis and transient and steady-state mode-based dynamic analysis (including the subspace projection method), as well as in complex eigenvalue extraction, matrix generation, and substructure generation. DVISCOUS(I,J)
+
+
+
+defines an element in the material viscous damping matrix at the current frequency. Abaqus/Standard requires that this element is defined as a damping value for each (I,J) entry in the damping matrix.
+
+Unless you invoke the unsymmetric equation solution capability in the contact property model definition (“Use with the unsymmetric equation solver in Abaqus/Standard” in “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide), Abaqus/Standard uses only the symmetric part of DVISCOUS. For a particular off-diagonal (I,J) entry the symmetrization is done by halving the sum of the (I,J) and (J,I) components.
+
+# DSTRUCTURAL(NDIR,NDIR)
+
+Interface structural damping matrix that can be used in direct steady-state dynamic analysis and steadystate mode-based dynamic analysis (including the subspace projection method), as well as in complex eigenvalue extraction, matrix generation, and substructure generation. DSTRUCTURAL(I,J) defines an element in the material structural damping matrix.
+
+Unless you invoke the unsymmetric equation solution capability in the contact property model definition (“Use with the unsymmetric equation solver in Abaqus/Standard” in “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide), Abaqus/Standard uses only the symmetric part of DSTRUCTURAL. For a particular off-diagonal (I,J) entry the symmetrization is done by halving the sum of the (I,J) and (J,I) components.
+
+# STATEV(NSTATV)
+
+An array containing the solution-dependent state variables. These are passed in as values at the beginning of the increment and must be returned as values at the end of the increment. You define the number of available state variables as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide.
+
+# SED
+
+This variable is passed in as the value of the elastic energy density at the start of the increment and should be updated to the elastic energy density at the end of the increment. This variable is used for output only and has no effect on other solution variables. It contributes to the output variable ALLSE.
+
+# SFD
+
+This variable should be defined as the incremental frictional dissipation. The units are energy per unit area. This variable is used for output only and has no effect on other solution variables. It contributes to the output variables ALLFD and SFDR (and related variables). For computing its contribution to SFDR, SFD is divided by the time increment.
+
+# SPD
+
+This variable should be defined as the incremental dissipation due to plasticity effects in the interfacial constitutive behavior. The units are energy per unit area. This variable is used for output only and has no effect on other solution variables. It contributes to the output variable ALLPD.
+
+
+
+# SVD
+
+This variable should be defined as the incremental dissipation due to viscous effects in the interfacial constitutive behavior. The units are energy per unit area. This variable is used for output only and has no effect on other solution variables. It contributes to the output variable ALLVD.
+
+# SCD
+
+This variable should be defined as the incremental dissipation due to creep effects in the interfacial constitutive behavior. The units are energy per unit area. This variable is used for output only and has no effect on other solution variables. It contributes to the output variable ALLCD.
+
+# PNEWDT
+
+Ratio of suggested new time increment to the time increment currently being used (DTIME, see below). This variable allows you to provide input to the automatic time incrementation algorithms in Abaqus/Standard (if automatic time incrementation is chosen).
+
+PNEWDT is set to a large value before each call to UINTER.
+
+If PNEWDT is redefined to be less than 1.0, Abaqus/Standard must abandon the time increment and attempt it again with a smaller time increment. The suggested new time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum value for all calls to user subroutines that allow redefinition of PNEWDT for this iteration.
+
+If PNEWDT is given a value that is greater than 1.0 for all calls to user subroutines for this iteration and the increment converges in this iteration, Abaqus/Standard may increase the time increment. The suggested new time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum value for all calls to user subroutines for this iteration.
+
+If automatic time incrementation is not selected in the analysis procedure, values of PNEWDT greater than 1.0 will be ignored and values of PNEWDT less than 1.0 will cause the job to terminate.
+
+# LOPENCLOSE
+
+An integer flag that is used to track the contact status in situations where user subroutine UINTER is used to model standard contact between two surfaces, like the default hard contact model in Abaqus/Standard. It comes in as the value at the beginning of the current iteration and should be set to the value at the end of the current iteration. It is set to −1 at the beginning of the analysis before UINTER is called. You should set it to 0 to indicate an open status and to 1 to indicate a closed status. A change in this flag from one iteration to the next will have two effects. It will result in output related to a change in contact status if you request a detailed contact printout in the message file (“The Abaqus/Standard message file” in “Output,” Section 4.1.1 of the Abaqus Analysis User’s Guide). In addition, it will also trigger a severe discontinuity iteration. Any time this flag is reset to a value of −1, Abaqus/Standard assumes that the flag is not being used. A change in this flag from −1 to another value or vice versa will not have any of the above effects.
+
+# LSTATE
+
+An integer flag that should be used in non-standard contact situations where a simple open/close status is not appropriate or enough to describe the state. It comes in as the value at the beginning of the current iteration and should be set to the value at the end of the current iteration. It is set to −1 at
+
+
+
+the beginning of the analysis before UINTER is called. It can be assigned any user-defined integer value, each corresponding to a different state. You can track changes in the value of this flag and use it to output appropriate diagnostic messages to the message file (unit 7). You may choose to output diagnostic messages only when a detailed contact printout is requested (“The Abaqus/Standard message file” in “Output,” Section 4.1.1 of the Abaqus Analysis User’s Guide). In the latter case, the LPRINT parameter is useful. In conjunction with the LSTATE flag, you may also utilize the LSDI flag to trigger a severe discontinuity iteration any time the state changes from one iteration to the next. Any time this flag is reset to a value of −1, Abaqus/Standard assumes that the flag is not being used.
+
+# LSDI
+
+This flag is set to 0 before each call to UINTER and should be set to 1 if the current iteration should be treated as a severe discontinuity iteration. This would typically be done in non-standard contact situations based on a change in the value of the LSTATE flag from one iteration to the next. The use of this flag has no effect when the LOPENCLOSE flag is also used. In that case, severe discontinuity iterations are determined based on changes in the value of LOPENCLOSE alone.
+
+# Variables passed in for information
+
+# RDISP(NDIR)
+
+An array containing the current relative positions between the two surfaces at the end of the increment. The first component is the relative position of the point on the slave surface, with respect to the master surface, in the normal direction. The second and third components, if applicable, are the accumulated incremental relative tangential displacements, measured from the beginning of the analysis. For the relative position in the normal direction a negative quantity represents an open status, while a positive quantity indicates penetration into the master surface. For open points on the slave surface for which no pairing master is found, the first component is a very large negative number $( - 1 \times 1 0 ^ { 3 6 } )$ . The local directions in which the relative displacements are defined are stored in ALOCALDIR.
+
+# DRDISP(NDIR)
+
+An array containing the increments in relative positions between the two surfaces.
+
+# TEMP(2)
+
+Temperature at the end of the increment at a point on the slave surface and the opposing master surface, respectively.
+
+# DTEMP(2)
+
+Increment in temperature at the point on the slave surface and the opposing master surface, respectively.
+
+# PREDEF(2,NPRED)
+
+An array containing pairs of values of all the predefined field variables at the end of the current increment (initial values at the beginning of the analysis and current values during the analysis). The first value in a pair, PREDEF(1,NPRED), corresponds to the value at the point on the slave surface, and the second value, PFREDEF(2,NPRED), corresponds to the value of the field variable at the nearest point on the opposing surface.
+
+
+
+# DPRED(2,NPRED)
+
+Array of increments in predefined field variables.
+
+# TIME(1)
+
+Value of step time at the end of the increment.
+
+# TIME(2)
+
+Value of total time at the end of the increment.
+
+# DTIME
+
+Current increment in time.
+
+# FREQR
+
+Current frequency for steady-state dynamic analysis in rad/time.
+
+# CINAME
+
+User-specified surface interaction name, left justified.
+
+# SLNAME
+
+Slave surface name.
+
+# MSNAME
+
+Master surface name.
+
+# PROPS(NPROPS)
+
+User-specified array of property values to define the interfacial constitutive behavior between the contacting surfaces.
+
+# COORDS(MCRD)
+
+An array containing the current coordinates of this point.
+
+# ALOCALDIR(3,3)
+
+An array containing the direction cosines of the local surface coordinate system. The directions are stored in columns. For example, ALOCALDIR(1,1), ALOCALDIR(2,1), and ALOCALDIR(3,1) give the (1, 2, 3) components of the normal direction. Thus, the first direction is the normal direction to the surface, and the remaining two directions are the local tangent directions in the plane of the surface. The local system is defined by the geometry of the master surface. The convention for the local directions is the same as the convention in situations where the model uses the built-in contact capabilities in Abaqus/Standard (described in “Contact formulations in Abaqus/Standard,” Section 38.1.1 of the Abaqus Analysis User’s Guide, for the tangential directions).
+
+# DROT(2,2)
+
+Rotation increment matrix. For contact with a three-dimensional rigid surface, this matrix represents the incremental rotation of the surface directions relative to the rigid surface. It is provided so that vector- or tensor-valued state variables can be rotated appropriately in this subroutine. Relative
+
+
+
+# UINTER
+
+displacement components are already rotated by this amount before UINTER is called. This matrix is passed in as a unit matrix for two-dimensional and axisymmetric contact problems.
+
+# AREA
+
+Surface area associated with the contact point.
+
+# CHRLNGTH
+
+Characteristic contact surface face dimension.
+
+# NODE
+
+User-defined global slave node number (or internal node number for models defined in terms of an assembly of part instances) involved with this contact point. Corresponds to the predominant slave node of the constraint if the surface-to-surface contact formulation is used.
+
+# NDIR
+
+Number of force components at this point.
+
+# NSTATV
+
+Number of solution-dependent state variables.
+
+# NPRED
+
+Number of predefined field variables.
+
+# NPROPS
+
+User-defined number of property values associated with this interfacial constitutive model (“Interfacial constants” in “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide).
+
+# MCRD
+
+Number of coordinate directions at the contact point.
+
+# KSTEP
+
+Step number.
+
+# KINC
+
+Increment number.
+
+# KIT
+
+Iteration number. KIT=0 for the first assembly, KIT=1 for the first recovery/second assembly, KIT=2 for the second recovery/third assembly, and so on.
+
+# LINPER
+
+Linear perturbation flag. LINPER=1 if the step is a linear perturbation step. LINPER=0 if the step is a general step. For a linear perturbation step, the inputs to user subroutine UINTER represent perturbation quantities about the base state. The user-defined quantities in UINTER are also perturbation quantities.
+
+
+
+The Jacobian terms should be based on the base state. No change in contact status should occur during a linear perturbation step.
+
+# LPRINT
+
+This flag is equal to 1 if a detailed contact printout to the message file is requested and 0 otherwise (“The Abaqus/Standard message file” in “Output,” Section 4.1.1 of the Abaqus Analysis User’s Guide). This flag can be used to print out diagnostic messages regarding changes in contact status selectively only when a detailed contact printout is requested.
+
+
+
+
+
+# 1.1.43 UMASFL: User subroutine to specify prescribed mass flow rate conditions for a convection/diffusion heat transfer analysis.
+
+Product: Abaqus/Standard
+
+# References
+
+• “Uncoupled heat transfer analysis,” Section 6.5.2 of the Abaqus Analysis User’s Guide
+• \*MASS FLOW RATE
+• “UTEMP, UFIELD, UMASFL, and UPRESS,” Section 4.1.25 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine UMASFL:
+
+• can be used to prescribe the mass flow rate vector at the nodes of a model as a function of position and time;
+• will be called whenever a current value of mass flow rate (per unit area) is needed for a node listed in a user-subroutine-defined mass flow rate definition (the node should belong to one or more convection/diffusion elements); and
+• will overwrite any flow rate data specified for the associated mass flow rate definition outside the user subroutine.
+
+# User subroutine interface
+
+```fortran
+SUBROUTINE UMASFL (FLOW, KFLOW, KSTEP, KINC, TIME, NODE, COORDS)
+C
+INCLUDE 'ABA_PARAM.INC'
+C
+DIMENSION FLOW (KFLOW), TIME(2), COORDS(3)
+C
+user coding to define FLOW
+RETURN
+END
+```
+
+
+
+# Variable to be defined
+
+# FLOW
+
+Total value of the mass flow rate vector at this point. The number of components in this vector is KFLOW. If KFLOW=1, give the total mass flow rate through the cross-section (for one-dimensional elements). If KFLOW=2, give the x-component and y-component of the flow rate vector as FLOW(1) and FLOW(2). If KFLOW=3, give the x-component, y-component, and z-component as FLOW(1), FLOW(2), and FLOW(3).
+
+# Variables passed in for information
+
+# KFLOW
+
+The number of components in the mass flow rate vector. UMASFL will be called only once per node.
+
+# KSTEP
+
+Step number.
+
+# KINC
+
+Increment number.
+
+# TIME(1)
+
+Current value of step time.
+
+# TIME(2)
+
+Current value of total time.
+
+# NODE
+
+Node number.
+
+# COORDS
+
+An array containing the coordinates of this node.
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_030.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_030.md
new file mode 100644
index 0000000..7d974be
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_030.md
@@ -0,0 +1,348 @@
+
+
+# 1.1.44 UMAT: User subroutine to define a material’s mechanical behavior.
+
+# Product: Abaqus/Standard
+
+WARNING: The use of this subroutine generally requires considerable expertise. You are cautioned that the implementation of any realistic constitutive model requires extensive development and testing. Initial testing on a single-element model with prescribed traction loading is strongly recommended.
+
+# References
+
+• “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide
+• “User-defined thermal material behavior,” Section 26.7.2 of the Abaqus Analysis User’s Guide
+• \*USER MATERIAL
+• “SDVINI,” Section 4.1.11 of the Abaqus Verification Guide
+• “UMAT and UHYPER,” Section 4.1.21 of the Abaqus Verification Guide
+
+# Overview
+
+# User subroutine UMAT:
+
+• can be used to define the mechanical constitutive behavior of a material;
+• will be called at all material calculation points of elements for which the material definition includes a user-defined material behavior;
+• can be used with any procedure that includes mechanical behavior;
+• can use solution-dependent state variables;
+• must update the stresses and solution-dependent state variables to their values at the end of the increment for which it is called;
+• must provide the material Jacobian matrix, $\partial \triangle \sigma / \partial \triangle \varepsilon$ , for the mechanical constitutive model;
+• can be used in conjunction with user subroutine USDFLD to redefine any field variables before they are passed in; and
+• is described further in “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide.
+
+# Storage of stress and strain components
+
+In the stress and strain arrays and in the matrices DDSDDE, DDSDDT, and DRPLDE, direct components are stored first, followed by shear components. There are NDI direct and NSHR engineering shear components. The order of the components is defined in “Conventions,” Section 1.2.2 of the Abaqus Analysis User’s Guide. Since the number of active stress and strain components varies between element types, the routine must be coded to provide for all element types with which it will be used.
+
+
+
+# Defining local orientations
+
+If a local orientation (“Orientations,” Section 2.2.5 of the Abaqus Analysis User’s Guide) is used at the same point as user subroutine UMAT, the stress and strain components will be in the local orientation; and, in the case of finite-strain analysis, the basis system in which stress and strain components are stored rotates with the material.
+
+# Stability
+
+You should ensure that the integration scheme coded in this routine is stable—no direct provision is made to include a stability limit in the time stepping scheme based on the calculations in UMAT.
+
+# Convergence rate
+
+DDSDDE and—for coupled temperature-displacement and coupled thermal-electrical-structural analyses—DDSDDT, DRPLDE, and DRPLDT must be defined accurately if rapid convergence of the overall Newton scheme is to be achieved. In most cases the accuracy of this definition is the most important factor governing the convergence rate. Since nonsymmetric equation solution is as much as four times as expensive as the corresponding symmetric system, if the constitutive Jacobian (DDSDDE) is only slightly nonsymmetric (for example, a frictional material with a small friction angle), it may be less expensive computationally to use a symmetric approximation and accept a slower convergence rate.
+
+An incorrect definition of the material Jacobian affects only the convergence rate; the results (if obtained) are unaffected.
+
+# Viscoelastic behavior in frequency domain
+
+The constitutive Jacobian (DDSDDE) must provide both the stiffness (storage modulus) and damping (loss modulus) for modeling frequency domain viscoelastic behavior.
+
+# Special considerations for various element types
+
+There are several special considerations that need to be noted.
+
+# Deformation gradient
+
+The deformation gradient is available for solid (continuum) elements, membranes, and finite-strain shells (S3/S3R, S4, S4R, SAXs, and SAXAs). It is not available for beams or small-strain shells. It is stored as a 3 × 3 matrix with component equivalence DFGRD0 $( \pmb { \tau } , \pmb { \sigma } ) \Leftrightarrow F _ { I J }$ . For fully integrated firstorder isoparametric elements (4-node quadrilaterals in two dimensions and 8-node hexahedra in three dimensions) the selectively reduced integration technique is used (also known as the technique). Thus, a modified deformation gradient
+
+$$
+\overline {{\mathbf {F}}} = \mathbf {F} \left(\frac {\overline {{J}}}{\overline {{J}}}\right) ^ {\frac {1}{n}}
+$$
+
+
+
+is passed into user subroutine UMAT. For more details, see “Solid isoparametric quadrilaterals and hexahedra,” Section 3.2.4 of the Abaqus Theory Guide.
+
+The deformation gradient, which is passed to the user subroutine, is computed with respect to the initial configuration. If a local orientation is not specified, the components of the deformation gradient are expressed in the global coordinate system. If a local orientation is used, the components of the same deformation gradient are expressed in the local coordinate system; in the case of finite-strain analysis, the basis system rotates with the material.
+
+# Beams and shells that calculate transverse shear energy
+
+If user subroutine UMAT is used to describe the material of beams or shells that calculate transverse shear energy, you must specify the transverse shear stiffness as part of the beam or shell section definition to define the transverse shear behavior. See “Shell section behavior,” Section 29.6.4 of the Abaqus Analysis User’s Guide, and “Choosing a beam element,” Section 29.3.3 of the Abaqus Analysis User’s Guide, for information on specifying this stiffness.
+
+# Open-section beam elements
+
+When user subroutine UMAT is used to describe the material response of beams with open sections (for example, an I-section), the torsional stiffness is obtained as
+
+$$
+G J = \frac {(K _ {1 3} + K _ {2 3}) J}{2 k A},
+$$
+
+where J is the torsional constant, A is the section area, k is a shear factor, and $K _ { \alpha 3 }$ is the user-specified transverse shear stiffness (see “Transverse shear stiffness definition” in “Choosing a beam element,” Section 29.3.3 of the Abaqus Analysis User’s Guide).
+
+# Elements with hourglassing modes
+
+If this capability is used to describe the material of elements with hourglassing modes, you must define the hourglass stiffness factor for hourglass control based on the total stiffness approach as part of the element section definition. The hourglass stiffness factor is not required for enhanced hourglass control, but you can define a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surface normal). See “Section controls,” Section 27.1.4 of the Abaqus Analysis User’s Guide, for information on specifying the stiffness factor.
+
+# Pipe-soil interaction elements
+
+The constitutive behavior of the pipe-soil interaction elements (see “Pipe-soil interaction elements,” Section 32.12.1 of the Abaqus Analysis User’s Guide) is defined by the force per unit length caused by relative displacement between two edges of the element. The relative-displacements are available as “strains” (STRAN and DSTRAN). The corresponding forces per unit length must be defined in the STRESS array. The Jacobian matrix defines the variation of force per unit length with respect to relative displacement.
+
+
+
+For two-dimensional elements two in-plane components of “stress” and “strain” exist (NTENS=NDI=2, and NSHR=0). For three-dimensional elements three components of “stress” and “strain” exist (NTENS=NDI=3, and NSHR=0).
+
+# Large volume changes with geometric nonlinearity
+
+If the material model allows large volume changes and geometric nonlinearity is considered, the exact definition of the consistent Jacobian should be used to ensure rapid convergence. These conditions are most commonly encountered when considering either large elastic strains or pressure-dependent plasticity. In the former case, total-form constitutive equations relating the Cauchy stress to the deformation gradient are commonly used; in the latter case, rate-form constitutive laws are generally used.
+
+For total-form constitutive laws, the exact consistent Jacobian is defined through the variation in Kirchhoff stress:
+
+$$
+\delta (J \pmb {\sigma}) = J (\mathbf {C}: \delta \mathbf {D} + \delta \mathbf {W} \cdot \pmb {\sigma} - \pmb {\sigma} \cdot \delta \mathbf {W})
+$$
+
+Here, J is the determinant of the deformation gradient, is the Cauchy stress, is the virtual rate of deformation, and is the virtual spin tensor, defined as
+
+$$
+\delta \mathbf {D} \stackrel {\mathrm{def}} {=} \mathrm{sym} (\delta \mathbf {F} \cdot \mathbf {F} ^ {- 1})
+$$
+
+and
+
+$$
+\delta \mathbf {W} \stackrel {\mathrm{def}} {=} \operatorname{asym} (\delta \mathbf {F} \cdot \mathbf {F} ^ {- 1}).
+$$
+
+For rate-form constitutive laws, the exact consistent Jacobian is given by
+
+$$
+\mathbf {C} = \frac {1}{J} \frac {\partial \Delta (J \pmb {\sigma})}{\partial \Delta \pmb {\varepsilon}}.
+$$
+
+# Use with almost incompressible or fully incompressible elastic materials
+
+For user-defined almost incompressible or incompressible elastic materials, a few different options are available depending on whether hybrid or nonhybrid elements are used. For all cases the first option should be to use user subroutine UHYPER instead of user subroutine UMAT when it is possible to do so. In user subroutine UMAT incompressible materials can be modeled via a penalty method; that is, you ensure that a finite bulk modulus is used. The bulk modulus should be large enough to model incompressibility sufficiently but small enough to avoid loss of precision. As a general guideline, the bulk modulus should be about $1 0 ^ { 4 } – 1 0 ^ { 6 }$ times the shear modulus. The tangent bulk modulus $K ^ { t }$ can be calculated from
+
+$$
+K ^ {t} = \frac {1}{9} \sum_ {\mathrm{I} = 1} ^ {3} \sum_ {\mathrm{J} = 1} ^ {3} \mathrm{DDSDDE} (\mathrm{I}, \mathrm{J}).
+$$
+
+
+
+If a hybrid element is used with user subroutine UMAT, Abaqus/Standard, by default, replaces the pressure stress calculated from your definition of STRESS with that derived from the Lagrange multiplier and modifies the Jacobian appropriately (“Hybrid incompressible solid element formulation,” Section 3.2.3 of the Abaqus Theory Guide). This approach is suitable for material models that use an incremental formulation (for example, metal plasticity) but is not consistent with a total formulation that is commonly used for hyperelastic materials. In the latter situation, the default formulation may lead to convergence problems. Such convergence problems may be observed, for example, when an almost incompressible nonlinear elastic user material is subjected to large deformations. Abaqus/Standard provides an alternate total formulation when user materials are used with hybrid elements (see “User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User’s Guide). This formulation is consistent with the native almost incompressible formulation used by Abaqus for hyperelastic materials (“Hyperelastic material behavior,” Section 4.6.1 of the Abaqus Theory Guide) and works better than the default formulation for such cases.
+
+The total hybrid formulation assumes that the response of the material can be written as the sum of its deviatoric and volumetric parts and that these parts are decoupled from each other. In particular, the volumetric response is assumed to be defined in terms of a strain energy potential, $U ( \hat { J } )$ , which is a function of an alternate variable, $\hat { J }$ in place of the actual volume change . The alternate variable is made available inside user subroutine UMAT by extending the STRESS array beyond NTENS, with the NTENS+1 entry providing read access to ${ \hat { J } } .$ You must define the hydrostatic part of the stress tensor as $\begin{array} { r } { \hat { p } = - \frac { \partial U } { \partial \hat { J } } } \end{array}$ a . The formulation also requires the additional derivatives e these additional derivatives inside user subroutine UMAT $\begin{array} { r } { \hat { K } = J \frac { \partial ^ { 2 } \hat { \boldsymbol { U } } } { \partial \hat { J } ^ { 2 } } } \end{array}$ andk and S+1 $\begin{array} { r } { \frac { \partial \hat { K } } { \partial \hat { J } } = J \frac { \partial ^ { 3 } U } { \partial \hat { J } ^ { 3 } } } \end{array}$ a3 . YouNS+2 entry, respectively, of the STRESS array. In addition, the bulk modulus of the material (contributes toward the material Jacobian matrix, DDSDDE) must be defined as $\hat { K }$ .
+
+Abaqus/Standard also provides a fully incompressible user material formulation for use with hybrid elements to define a fully incompressible user material response. This formulation is consistent with the native formulation used by Abaqus for incompressible hyperelastic materials and assumes that the deviatoric stress can be derived from a strain energy potential function. You need define only the deviatoric stress and Jacobian to define a fully incompressible material response through user subroutine UMAT.
+
+For incompressible pressure-sensitive materials the element choice is particularly important when using user subroutine UMAT. In particular, first-order wedge elements should be avoided. For these elements the $\bar { B }$ technique is not used to alter the deformation gradient that is passed into user subroutine UMAT, which increases the risk of volumetric locking.
+
+# Increments for which only the Jacobian can be defined
+
+Abaqus/Standard passes zero strain increments into user subroutine UMAT to start the first increment of all the steps and all increments of steps for which you have suppressed extrapolation (see “Defining an analysis,” Section 6.1.2 of the Abaqus Analysis User’s Guide). In this case you can define only the Jacobian (DDSDDE).
+
+
+
+# Utility routines
+
+Several utility routines may help in coding user subroutine UMAT. Their functions include determining stress invariants for a stress tensor and calculating principal values and directions for stress or strain tensors. These utility routines are discussed in detail in “Obtaining stress invariants, principal stress/strain values and directions, and rotating tensors in an Abaqus/Standard analysis,” Section 2.1.11.
+
+User subroutine interface
+SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
+1 RPL,DDSDDT,DRPLDE,DRPLDT,
+2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
+3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
+4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
+C
+INCLUDE 'ABA_PARAM.INC'
+C
+CHARACTER*80 CMNAME
+DIMENSION STRESS(NTENS),STATEV(NSTATV),
+1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
+2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
+3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
+4 JSTEP(4)
+user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
+and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
+RETURN
+END
+
+# Variables to be defined
+
+# In all situations
+
+DDSDDE(NTENS,NTENS)
+
+Jacobian matrix of the constitutive model, $\partial \Delta \sigma / \partial \Delta \varepsilon$ , where $\Delta \sigma$ are the stress increments and $\Delta \varepsilon$ are the strain increments. DDSDDE(I,J) defines the change in the Ith stress component at the end of the time increment caused by an infinitesimal perturbation of the Jth component of the strain increment array. Unless you invoke the unsymmetric equation solution capability for the user-defined material, Abaqus/Standard will use only the symmetric part of DDSDDE. The symmetric part of the matrix is calculated by taking one half the sum of the matrix and its transpose.
+
+For viscoelastic behavior in the frequency domain, the Jacobian matrix must be dimensioned as DDSDDE(NTENS,NTENS,2). The stiffness contribution (storage modulus) must be provided in
+
+
+
+DDSDDE(NTENS,NTENS,1), while the damping contribution (loss modulus) must be provided in DDSDDE(NTENS,NTENS,2).
+
+# STRESS(NTENS)
+
+This array is passed in as the stress tensor at the beginning of the increment and must be updated in this routine to be the stress tensor at the end of the increment. If you specified initial stresses (“Initial conditions in Abaqus/Standard and Abaqus/Explicit,” Section 34.2.1 of the Abaqus Analysis User’s Guide), this array will contain the initial stresses at the start of the analysis. The size of this array depends on the value of NTENS as defined below. In finite-strain problems the stress tensor has already been rotated to account for rigid body motion in the increment before UMAT is called, so that only the corotational part of the stress integration should be done in UMAT. The measure of stress used is “true” (Cauchy) stress.
+
+If the UMAT utilizes a hybrid formulation that is total (as opposed to the default incremental behavior), the stress array is extended beyond NTENS. The first NTENS entries of the array contain the stresses, as described above. The additional quantities are as follows:
+
+STRESS(NTENS+1)
+
+Read only: ,
+
+STRESS(NTENS+2)
+
+Write only: , $\begin{array} { r } { \hat { K } = J \frac { \partial ^ { 2 } U } { \partial \hat { J } ^ { 2 } } } \end{array}$ and
+
+STRESS(NTENS+3)
+
+Write only: $\begin{array} { r } { \frac { \partial \hat { K } } { \partial \hat { J } } = J \frac { \partial ^ { 3 } U } { \partial \hat { J } ^ { 3 } } } \end{array}$ =Jou, , where is the volumetric part of the strain energy density potential.
+
+# STATEV(NSTATV)
+
+An array containing the solution-dependent state variables. These are passed in as the values at the beginning of the increment unless they are updated in user subroutines USDFLD or UEXPAN, in which case the updated values are passed in. In all cases STATEV must be returned as the values at the end of the increment. The size of the array is defined as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide.
+
+In finite-strain problems any vector-valued or tensor-valued state variables must be rotated to account for rigid body motion of the material, in addition to any update in the values associated with constitutive behavior. The rotation increment matrix, DROT, is provided for this purpose.
+
+# SSE, SPD, SCD
+
+Specific elastic strain energy, plastic dissipation, and “creep” dissipation, respectively. These are passed in as the values at the start of the increment and should be updated to the corresponding specific energy values at the end of the increment. They have no effect on the solution, except that they are used for energy output.
+
+# Only in a fully coupled thermal-stress or a coupled thermal-electrical-structural analysis
+
+# RPL
+
+Volumetric heat generation per unit time at the end of the increment caused by mechanical working of the material.
+
+
+
+# DDSDDT(NTENS)
+
+Variation of the stress increments with respect to the temperature.
+
+# DRPLDE(NTENS)
+
+Variation of RPL with respect to the strain increments.
+
+# DRPLDT
+
+Variation of RPL with respect to the temperature.
+
+# Only in a geostatic stress procedure or a coupled pore fluid diffusion/stress analysis for pore pressure cohesive elements
+
+# RPL
+
+RPL is used to indicate whether or not a cohesive element is open to the tangential flow of pore fluid. Set RPL equal to 0 if there is no tangential flow; otherwise, assign a nonzero value to RPL if an element is open. Once opened, a cohesive element will remain open to the fluid flow.
+
+# Variable that can be updated
+
+# PNEWDT
+
+Ratio of suggested new time increment to the time increment being used (DTIME, see discussion later in this section). This variable allows you to provide input to the automatic time incrementation algorithms in Abaqus/Standard (if automatic time incrementation is chosen). For a quasi-static procedure the automatic time stepping that Abaqus/Standard uses, which is based on techniques for integrating standard creep laws (see “Quasi-static analysis,” Section 6.2.5 of the Abaqus Analysis User’s Guide), cannot be controlled from within the UMAT subroutine.
+
+PNEWDT is set to a large value before each call to UMAT.
+
+If PNEWDT is redefined to be less than 1.0, Abaqus/Standard must abandon the time increment and attempt it again with a smaller time increment. The suggested new time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum value for all calls to user subroutines that allow redefinition of PNEWDT for this iteration.
+
+If PNEWDT is given a value that is greater than 1.0 for all calls to user subroutines for this iteration and the increment converges in this iteration, Abaqus/Standard may increase the time increment. The suggested new time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum value for all calls to user subroutines for this iteration.
+
+If automatic time incrementation is not selected in the analysis procedure, values of PNEWDT that are greater than 1.0 will be ignored and values of PNEWDT that are less than 1.0 will cause the job to terminate.
+
+# Variables passed in for information
+
+# STRAN(NTENS)
+
+An array containing the total strains at the beginning of the increment. If thermal expansion is included in the same material definition, the strains passed into UMAT are the mechanical strains only (that is,
+
+
+
+the thermal strains computed based upon the thermal expansion coefficient have been subtracted from the total strains). These strains are available for output as the “elastic” strains.
+
+In finite-strain problems the strain components have been rotated to account for rigid body motion in the increment before UMAT is called and are approximations to logarithmic strain.
+
+# DSTRAN(NTENS)
+
+Array of strain increments. If thermal expansion is included in the same material definition, these are the mechanical strain increments (the total strain increments minus the thermal strain increments).
+
+# TIME(1)
+
+Value of step time at the beginning of the current increment or frequency.
+
+# TIME(2)
+
+Value of total time at the beginning of the current increment.
+
+# DTIME
+
+Time increment.
+
+# TEMP
+
+Temperature at the start of the increment.
+
+# DTEMP
+
+Increment of temperature.
+
+# PREDEF
+
+Array of interpolated values of predefined field variables at this point at the start of the increment, based on the values read in at the nodes.
+
+# DPRED
+
+Array of increments of predefined field variables.
+
+# CMNAME
+
+User-defined material name, left justified. Some internal material models are given names starting with the “ABQ\_” character string. To avoid conflict, you should not use “ABQ\_” as the leading string for CMNAME.
+
+# NDI
+
+Number of direct stress components at this point.
+
+# NSHR
+
+Number of engineering shear stress components at this point.
+
+# NTENS
+
+Size of the stress or strain component array (NDI + NSHR).
+
+
+
+# NSTATV
+
+Number of solution-dependent state variables that are associated with this material type (defined as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
+
+# PROPS(NPROPS)
+
+User-specified array of material constants associated with this user material.
+
+# NPROPS
+
+User-defined number of material constants associated with this user material.
+
+# COORDS
+
+An array containing the coordinates of this point. These are the current coordinates if geometric nonlinearity is accounted for during the step (see “Defining an analysis,” Section 6.1.2 of the Abaqus Analysis User’s Guide); otherwise, the array contains the original coordinates of the point.
+
+# DROT(3,3)
+
+Rotation increment matrix. This matrix represents the increment of rigid body rotation of the basis system in which the components of stress (STRESS) and strain (STRAN) are stored. It is provided so that vector- or tensor-valued state variables can be rotated appropriately in this subroutine: stress and strain components are already rotated by this amount before UMAT is called. This matrix is passed in as a unit matrix for small-displacement analysis and for large-displacement analysis if the basis system for the material point rotates with the material (as in a shell element or when a local orientation is used).
+
+# CELENT
+
+Characteristic element length, which is a typical length of a line across an element for a first-order element; it is half of the same typical length for a second-order element. For beams and trusses it is a characteristic length along the element axis. For membranes and shells it is a characteristic length in the reference surface. For axisymmetric elements it is a characteristic length in the plane only. For cohesive elements it is equal to the constitutive thickness.
+
+# DFGRD0(3,3)
+
+Array containing the deformation gradient at the beginning of the increment. If a local orientation is defined at the material point, the deformation gradient components are expressed in the local coordinate system defined by the orientation at the beginning of the increment. For a discussion regarding the availability of the deformation gradient for various element types, see “Deformation gradient.”
+
+# DFGRD1(3,3)
+
+Array containing the deformation gradient at the end of the increment. If a local orientation is defined at the material point, the deformation gradient components are expressed in the local coordinate system defined by the orientation. This array is set to the identity matrix if nonlinear geometric effects are not included in the step definition associated with this increment. For a discussion regarding the availability of the deformation gradient for various element types, see “Deformation gradient.”
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_031.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_031.md
new file mode 100644
index 0000000..c42179e
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_031.md
@@ -0,0 +1,381 @@
+
+
+# NOEL
+
+Element number.
+
+# NPT
+
+Integration point number.
+
+# LAYER
+
+Layer number (for composite shells and layered solids).
+
+# KSPT
+
+Section point number within the current layer.
+
+# JSTEP(1)
+
+Step number.
+
+# JSTEP(2)
+
+Procedure type key (see “Results file output format,” Section 5.1.2 of the Abaqus Analysis User’s Guide).
+
+# JSTEP(3)
+
+1 if NLGEOM=YES for the current step; 0 otherwise.
+
+# JSTEP(4)
+
+1 if current step is a linear perturbation procedure; 0 otherwise.
+
+# KINC
+
+Increment number.
+
+# Example: Using more than one user-defined mechanical material model
+
+To use more than one user-defined mechanical material model, the variable CMNAME can be tested for different material names inside user subroutine UMAT as illustrated below:
+
+```txt
+IF (CMNAME(1:4) .EQ. 'MAT1') THEN
+ CALL UMAT_MAT1(argument_list)
+ELSE IF(CMNAME(1:4) .EQ. 'MAT2') THEN
+ CALL UMAT_MAT2(argument_list)
+END IF
+```
+
+UMAT\_MAT1 and UMAT\_MAT2 are the actual user material subroutines containing the constitutive material models for each material MAT1 and MAT2, respectively. Subroutine UMAT merely acts as a directory here. The argument list may be the same as that used in subroutine UMAT.
+
+
+
+As a simple example of the coding of user subroutine UMAT, consider the linear, viscoelastic model shown in Figure 1.1.44–1. Although this is not a very useful model for real materials, it serves to illustrate how to code the routine.
+
+The behavior of the one-dimensional model shown in the figure is
+
+$$
+\sigma + \frac {\mu_ {1}}{(E _ {1} + E _ {2})} \dot {\sigma} = \frac {\mu_ {1}}{(1 + E _ {1} / E _ {2})} \dot {\varepsilon} + \frac {1}{(1 / E _ {1} + 1 / E _ {2})} \varepsilon ,
+$$
+
+where $\dot { \sigma }$ and $\dot { \varepsilon }$ are the time rates of change of stress and strain. This can be generalized for small straining of an isotropic solid as
+
+$$
+\sigma_ {x x} + \tilde {\nu} \dot {\sigma} _ {x x} = \lambda \varepsilon_ {V} + 2 \mu \varepsilon_ {x x} + \tilde {\lambda} \dot {\varepsilon} _ {V} + 2 \tilde {\mu} \dot {\varepsilon} _ {x x}, \quad \mathrm{etc.},
+$$
+
+and
+
+$$
+\sigma_ {x y} + \tilde {\nu} \dot {\sigma} _ {x y} = \mu \gamma_ {x y} + \tilde {\mu} \dot {\gamma} _ {x y}, \quad \mathrm{etc.},
+$$
+
+where
+
+$$
+\varepsilon_ {V} = \varepsilon_ {x x} + \varepsilon_ {y y} + \varepsilon_ {z z},
+$$
+
+and $\tilde { \nu } , \lambda , \mu , \tilde { \lambda } ,$ , and $\tilde { \mu }$ are material constants ( and $\mu$ are the Lamé constants).
+
+A simple, stable integration operator for this equation is the central difference operator:
+
+$$
+\dot {f} _ {t + \frac {1}{2} \Delta t} = \frac {\Delta f}{\Delta t},
+$$
+
+$$
+f _ {t + \frac {1}{2} \Delta t} = f _ {t} + \frac {\Delta f}{2},
+$$
+
+where $f \mathrm { i }$ is some function, $f _ { t }$ is its value at the beginning of the increment, $\Delta f$ is the change in the function over the increment, and $\Delta t$ is the time increment.
+
+Applying this to the rate constitutive equations above gives
+
+$$
+(\frac {\Delta t}{2} + \tilde {\nu}) \Delta \sigma_ {x x} = (\Delta t \frac {\lambda}{2} + \tilde {\lambda}) \Delta \varepsilon_ {V} + (\Delta t \mu + 2 \tilde {\mu}) \Delta \varepsilon_ {x x} + \Delta t (\lambda \varepsilon_ {V} + 2 \mu \varepsilon_ {x x} - \sigma_ {x x}) _ {t}, \quad \mathrm{etc.,}
+$$
+
+and
+
+$$
+(\frac {\Delta t}{2} + \tilde {\nu}) \Delta \sigma_ {x y} = (\Delta t \frac {\mu}{2} + \tilde {\mu}) \Delta \gamma_ {x y} + \Delta t (\mu \gamma_ {x y} - \sigma_ {x y}) _ {t}, \quad \mathrm{etc.},
+$$
+
+so that the Jacobian matrix has the terms
+
+
+
+
+
+| lFlagsDefine (iComputeSecDeriv) | If set to 1, you must provide the computation of the amplitude second derivative. The default is 0, which means that Abaqus computes the second derivative automatically. |
| lFlagsDefine (iComputeInteg) | If set to 1, you must provide the computation of the amplitude incremental integral. The default is 0, which means that Abaqus computes the incremental integral automatically. |
| lFlagsDefine (iStopAnalysis) | If set to 1, the analysis will be stopped and an error message will be issued. The default is 0, which means that Abaqus will not stop the analysis. |
| lFlagsDefine (iConcludeStep) | If set to 1, Abaqus will conclude the step execution and advance to the next step (if a next step is available). The default is 0. |
+
+# svars
+
+An array containing the values of the solution-dependent state variables associated with this amplitude definition. The number of such variables is nsvars (see above). You define the meaning of these variables.
+
+This array is passed into VUAMP containing the values of these variables at the start of the current increment. In most cases they should be updated to be the values at the end of the increment.
+
+# AmpDerivative
+
+Current value of the amplitude derivative.
+
+# AmpSecDerivative
+
+Current value of the amplitude second derivative.
+
+# AmpIncIntegral
+
+Current value of the amplitude incremental integral.
+
+# Variables passed in for information
+
+# ampName
+
+User-specified amplitude name, left justified.
+
+# time(iStepTime)
+
+Current value of step time.
+
+
+
+time(iTotalTime)
+
+Current value of total time.
+
+ampValueOld
+
+Old value of the amplitude from the previous increment.
+
+dt
+
+Current stable time increment.
+
+nprops
+
+User-defined number of properties associated with this amplitude definition.
+
+props(nprops)
+
+User-supplied amplitude properties.
+
+nSvars
+
+User-defined number of solution-dependent state variables associated with this amplitude definition.
+
+lFlagsInfo
+
+Integer flag array with information regrading the current call to VUAMP:
+
+lFlagsInfo(iInitialization)
+
+This flag is equal to 1 if VUAMP is called from the initialization phase of each step and is set to 0 otherwise.
+
+lFlagsInfo(iRegularInc)
+
+This flag is equal to 1 if VUAMP is called from a regular increment and is set to 0 otherwise.
+
+lFlagsInfo(ikStep)
+
+Step number.
+
+nSensor
+
+Total number of sensors in the model.
+
+sensorValues
+
+Array with sensor values at the end of the previous increment. Each sensor value corresponds to a history output variable associated with the output database request defining the sensor.
+
+sensorNames
+
+Array with user-defined sensor names in the entire model, left justified. Each sensor name corresponds to a sensor value provided with the output database request. All names will be converted to uppercase characters if lowercase or mixed-case characters were used in their definition.
+
+jSensorLookUpTable
+
+Variable that must be passed into the utility functions IVGETSENSORID and VGETSENSORVALUE.
+
+
+
+Example: Amplitude definition using sensor and state variables
+```txt
+c user amplitude subroutine
+Subroutine VUAMP(
+C passed in for information and state variables
+* ampName, time, ampValueOld, dt, nprops, props, nSvars,
+* svars, lFlagsInfo, nSensor, sensorValues, sensorNames,
+* jSensorLookUpTable,
+C to be defined
+* ampValueNew,
+* lFlagsDefine,
+* AmpDerivative, AmpSecDerivative, AmpIncIntegral)
+include 'vaba_param.inc'
+C svars - additional state variables, similar to (V)UEL
+ dimension sensorValues(nSensor), props(nprops),
+* svars(nSvars)
+ character*80 sensorNames(nSensor)
+ character*80 ampName
+C time indices
+ parameter( iStepTime = 1,
+* iTotalTime = 2,
+* nTime = 2)
+C flags passed in for information
+ parameter( iInitialization = 1,
+* iRegularInc = 2,
+* ikStep = 3,
+* nFlagsInfo = 3)
+C optional flags to be defined
+ parameter( iComputeDeriv = 1,
+* iComputeSecDeriv = 2,
+* iComputeInteg = 3,
+* iStopAnalysis = 4,
+* iConcludeStep = 5,
+* nFlagsDefine = 5)
+ parameter( tStep=0.18, tAccelerateMotor = .00375,
+* omegaFinal=23.26)
+```
+c Alternatively, assign the user-defined amplitude
+
+
+
+```txt
+c properties on the data lines rather than using a parameter
+c definition above.
+c tStep = props(1)
+c tAccelerateMotor = props(2)
+c omegaFinal = props(3)
+
+ dimension time(nTime), lFlagsInfo(nFlagsInfo),
+ * lFlagsDefine(nFlagsDefine)
+ dimension jSensorLookUpTable(*)
+
+ lFlagsDefine(iComputeDeriv) = 1
+ lFlagsDefine(iComputeSecDeriv) = 1
+
+c get sensor value
+ vTrans_CU1 = vGetSensorValue('HORIZ_TRANSL_MOTION',
+ * jSensorLookUpTable,
+ * sensorValues)
+
+ if (ampName(1:22) .eq. 'MOTOR_WITH_STOP_SENSOR') then
+ if (lFlagsInfo(iInitialization).eq.1) then
+ ampValueNew = ampValueOld
+
+ svars(1) = 0.0
+ svars(2) = 0.0
+ else
+ tim = time(iStepTime)
+
+c ramp up the angular rot velocity of the electric
+c motor after which hold constant
+ if (tim .le. tAccelerateMotor) then
+ ampValueNew = omegaFinal*tim/tAccelerateMotor
+
+ else
+ ampValueNew = omegaFinal
+ end if
+
+c retrieve old sensor value
+ vTrans_CU1_old = svars(1)
+
+c detect a zero crossing and count the number of
+c crossings
+ if (vTrans_CU1_old*vTrans_CU1 .le. 0.0 .and.
+```
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_045.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_045.md
new file mode 100644
index 0000000..2bfb596
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_045.md
@@ -0,0 +1,307 @@
+
+
+```txt
+* tim .gt. tAccelerateMotor ) then
+ svars(2) = svars(2) + 1.0
+ end if
+ nrCrossings = int(svars(2))
+
+c stop the motor if sensor crosses zero the second
+c time
+ if (nrCrossings.eq.2) then
+ ampValueNew = 0.0
+ lFlagsDefine(iConcludeStep)=1
+ end if
+
+c store sensor value
+ svars(1) = vTrans_CU1
+ end if
+ end if
+
+return
+end
+```
+
+
+
+
+
+# 1.2.10 VUANISOHYPER\_INV: User subroutine to define anisotropic hyperelastic material behavior using the invariant formulation.
+
+Product: Abaqus/Explicit
+
+# References
+
+• “Anisotropic hyperelastic behavior,” Section 22.5.3 of the Abaqus Analysis User’s Guide
+• \*ANISOTROPIC HYPERELASTIC
+• “UANISOHYPER\_INV and VUANISOHYPER\_INV,” Section 4.1.13 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine VUANISOHYPER\_INV:
+
+• can be used to define the strain energy potential of anisotropic hyperelastic materials as a function of an irreducible set of scalar invariants;
+• will be called for blocks of material calculation points for which the material definition contains user-defined anisotropic hyperelastic behavior with invariant-based formulation (“Anisotropic hyperelastic behavior,” Section 22.5.3 of the Abaqus Analysis User’s Guide);
+• can use and update solution-dependent state variables;
+• can use any field variables that are passed in; and
+• requires that the values of the derivatives of the strain energy density function be defined with respect to the scalar invariants.
+
+# Enumeration of invariants
+
+To facilitate coding and provide easy access to the array of invariants passed to user subroutine VUANISOHYPER\_INV, an enumerated representation of each invariant is introduced. Any scalar invariant can, therefore, be represented uniquely by an enumerated invariant, $I _ { n } ^ { * }$ , where the subscript n denotes the order of the invariant according to the enumeration scheme in the following table:
+
+| Component | 2D Case | 3D Case |
| 1 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{11}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{11}^{G}$ |
| 2 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{22}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{22}^{G}$ |
| 3 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{22}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{22}^{G}$ |
| 4 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{33}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{33}^{G}$ |
| 5 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{33}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{33}^{G}$ |
| 6 | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{33}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{33}^{G}$ |
| 7 | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 8 | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 9 | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 10 | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{12}^{G}$ | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{12}^{G}$ |
| 11 | | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 12 | | $\partial^{2}U/\partial\overline{\varepsilon}_{22}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 13 | | $\partial^{2}U/\partial\overline{\varepsilon}_{33}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 14 | | $\partial^{2}U/\partial\overline{\varepsilon}_{12}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 15 | | $\partial^{2}U/\partial\overline{\varepsilon}_{23}^{G}\partial\overline{\varepsilon}_{23}^{G}$ |
| 16 | | $\partial^{2}U/\partial\overline{\varepsilon}_{11}^{G}\partial\overline{\varepsilon}_{31}^{G}$ |
+
+
+
+| predef (kblock, k2, 1, k4) | Temperature. |
| predef (kblock, k2, 2, k4) | First predefined field variable. |
| predef (kblock, k2, 3, k4) | Second predefined field variable. |
| Etc. | Any other predefined field variable. |
+
+
+
+predef(kblock,k2,k3,k4)
+
+Value of the (k3–1)th predefined field variable at the k2th node of the element at the beginning or the end of the increment.
+
+predef(kblock,k2,k3,1)
+
+Values of the variables at the end of the current increment.
+
+predef(kblock,k2,k3,2)
+
+Values of the variables at the beginning of the current increment.
+
+# lflags
+
+An array containing the flags that define the current solution procedure and requirements for element calculations.
+
+lflags(iProcedure)
+
+Defines the procedure type. See
+
+“Results file output format,”
+
+Section 5.1.2 of the Abaqus Analysis
+
+User’s Guide, for the key used for each procedure.
+
+lflags(iNlgeom)=0
+
+Small-displacement analysis.
+
+lflags(iNlgeom)=1
+
+Large-displacement analysis (nonlinear geometric effects included in the step; see “General and linear perturbation procedures,” Section 6.1.3 of the Abaqus Analysis User’s Guide).
+
+lflags(iOpCode)=jMassCalc
+
+Define the mass matrix amass in the beginning of the analysis.
+
+lflags(iOpCode)=jIntForceAnd-DtStable
+
+Define the element internal force.
+
+Define the stable time increment as well.
+
+lflags(iOpCode)=jExternForce
+
+Define the distributed load effect on the external force associated with the element.
+
+# dMassScaleFactor
+
+An array containing the mass scale factors for each element.
+
+# time(iStepTime)
+
+Current value of step time.
+
+
+
+time(iTotalTime)
+
+Current value of total time.
+
+# Scalar parameters:
+
+nblock
+
+Number of user elements to be processed in this call to VUEL.
+
+dtimeCur
+
+Current time increment.
+
+dtimePrev
+
+Previous time increment.
+
+period
+
+Time period of the current step.
+
+ndofel
+
+Number of degrees of freedom in the elements processed.
+
+nsvars
+
+User-defined number of solution-dependent state variables associated with the element (“Defining the number of solution-dependent variables that must be stored within the element” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+nprops
+
+User-defined number of real property values associated with the elements processed (“Defining the element properties” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+njprops
+
+User-defined number of integer property values associated with the elements processed (“Defining the element properties” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+mcrd
+
+mcrd is defined as the maximum of the user-defined maximum number of coordinates needed at any node point (“Defining the maximum number of coordinates needed at any nodal point” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide) and the value of the largest active degree of freedom of the user element that is less than or equal to 3. For example, if you specify that the maximum number of coordinates is 1 and the active degrees of freedom of the user element are 2, 3, and 6 mcrd will be 3. If you specify that the maximum number of coordinates is 2 and the active degree of freedom of the user element is 11, mcrd will be 2.
+
+nnode
+
+User-defined number of nodes on the elements (“Defining the number of nodes associated with the element” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+
+
+# jtype
+
+Integer defining the element type. This is the user-defined integer value n in element type VUn (“Assigning an element type key to a user-defined element” in “User-defined elements,” Section 32.17.1 of the Abaqus Analysis User’s Guide).
+
+# kstep
+
+Current step number.
+
+# kinc
+
+Current increment number.
+
+# npredef
+
+Number of predefined field variables, including temperature. For user elements Abaqus/Explicit uses one value for each field variable per node.
+
+# VUEL conventions
+
+The solution variables (displacement, velocity, etc.) are arranged on a node/degree of freedom basis. The degrees of freedom of the first node are first, followed by the degrees of freedom of the second node, etc. The degrees of freedom that will be updated automatically in Abaqus/Explicit are: 1–3 (displacements), 4–6 (rotations), 8 (pressure), and 11 (temperature). Depending on the procedure type (see below), only some of the degrees of freedom listed above will be updated. Other degrees of freedom will not be updated by the time integration procedure in Abaqus/Explicit and, hence, should not be used.
+
+The mass matrix defined in user subroutine VUEL must be symmetric. In addition, the following requirements apply:
+
+• The mass matrix entries associated with the translational degrees of freedom for a particular node must be diagonal. Moreover, these diagonal entries must be equal to each other.
+• There must be no coupling (off-diagonal) entries specified between degrees of freedom of different kinds. For example, you cannot specify nonzero mass matrix entries to couple the translational degrees of freedom to the rotational degrees of freedom.
+• There must be no coupling (off-diagonal) entries specified between degrees of freedom belonging to different nodes.
+
+You must be using appropriate lumping techniques to provide a mass matrix that follows these requirements. For the rotational degrees of freedom at a particular node in three-dimensional analyses, you can specify a fully populated symmetric 3 × 3 inertia tensor.
+
+# Usage with general nonlinear procedures
+
+The following illustrates the use in explicit dynamic procedures:
+
+# Direct-integration explicit dynamic analysis (lflags(iProcedure)=17)
+
+• Automatic updates for degrees of freedom 1–6, 8, and 11.
+
+
+
+• The governing equations are as described in “Explicit dynamic analysis,” Section 6.3.3 of the Abaqus Analysis User’s Guide.
+• Coding for the operational code lflags(iOpCode)=jExternForceis optional.
+
+# Transient fully coupled thermal-stress analysis (lflags(iProcedure)=74)
+
+• Automatic updates for degrees of freedom 1–6 and 11.
+• The governing equations are as described in “Fully coupled thermal-stress analysis in Abaqus/Explicit” in “Fully coupled thermal-stress analysis,” Section 6.5.3 of the Abaqus Analysis User’s Guide.
+• Coding for the operational code lflags(iOpCode)=jExternForce is optional.
+
+# Example: Structural user element
+
+A structural user element has been created to demonstrate the usage of subroutine VUEL. These userdefined elements are applied in a number of analyses. The following excerpt is from the verification problem that invokes the structural user element in an explicit dynamic procedure:
+
+```txt
+*USER ELEMENT, NODES=2, TYPE=VU1, PROPERTIES=4, COORDINATES=3, VARIABLES=12
+1, 2, 3
+*ELEMENT, TYPE=VU1
+101, 101, 102
+*ELGEN, ELSET=VUTRUSS
+101, 5
+*UEL PROPERTY, ELSET=VUTRUSS
+0.002, 2.1E11, 0.3, 7200.
+```
+
+The user element consists of two nodes that are assumed to lie parallel to the x-axis. The element behaves similarly to a linear truss element. The supplied element properties are the cross-sectional area, Young’s modulus, Poisson’s ratio, and density, respectively.
+
+The next excerpt shows the listing of the subroutine. The user subroutine has been coded for use in an explicit dynamic analysis. The names of the verification input files associated with the subroutine and these procedures can be found in “VUEL,” Section 4.1.33 of the Abaqus Verification Guide.
+
+```csv
+subroutine vuel(
+* nblock,
+* rhs, amass, dtimeStable,
+* svars, nsvars,
+* energy,
+* nnode, ndofel,
+* props, nprops,
+* jprops, njprops,
+* coords, ncrd,
+* u, du, v, a,
+```
+
+
+
+```txt
+* jtype, jElem,
+* time, period, dtimeCur, dtimePrev, kstep, kinc, lflags,
+* dMassScaleFactor,
+* predef, npredef,
+* ndload, adImag)
+
+include 'vaba_param.inc'
+
+c operation code
+parameter ( jMassCalc = 1,
+* jIntForceAndDtStable = 4)
+
+c flags
+parameter (iProcedure = 1,
+* iNlgeom = 2,
+* iOpCode = 3,
+* nFlags = 3)
+
+c procedure flags
+parameter ( jDynExplicit = 17 )
+
+c time
+parameter (iStepTime = 1,
+* iTotalTime = 2,
+* nTime = 2)
+
+c energies
+parameter ( iElPd = 1,
+* iElCd = 2,
+* iElIe = 3,
+* iElTs = 4,
+* iElDd = 5,
+* iElBv = 6,
+* iElDe = 7,
+* iElHe = 8,
+* iElKe = 9,
+* iElTh = 10,
+* iElDmd = 11,
+* iElDc = 12,
+* nElEnergy = 12)
+
+parameter (factorStable = 0.99d0)
+```
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_050.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_050.md
new file mode 100644
index 0000000..09141b0
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_050.md
@@ -0,0 +1,374 @@
+
+
+```python
+parameter (zero = 0.d0, half = 0.5d0, one = 1.d0, two=2.d0)
+c
+ dimension rhs(nblock,ndofel), amass(nblock,ndofel,ndofel),
+ * dtimeStable(nblock),
+ * svars(nblock,nsvars), energy(nblock,nElEnergy),
+ * props(nprops), jprops(njprops),
+ * jElem(nblock), time(nTime), l(nFlags),
+ * coords(nblock,nnode,ncrd), u(nblock,ndofel),
+ * du(nblock,ndofel), v(nblock,ndofel), a(nblock, ndofel),
+ * predef(nblock, nnode, npred, nPred), adlmag(nblock),
+ * dMassScaleFactor(nblock)
+
+c Notes:
+c Define only nonzero entries; the arrays to be defined
+c have been zeroed out just before this call
+
+if (jtype .eq. 1001 .and.
+ * lflags(iProcedure).eq.jDynExplicit) then
+
+ area0 = props(1)
+ eMod = props(2)
+ anu = props(3)
+ rho = props(4)
+
+ eDampTra = zero
+ amassFact0 = half*area0*rho
+
+ if ( lflags(iOpCode).eq.jMassCalc ) then
+ do kblock = 1, nblock
+
+c use original distance to compute mass
+ alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
+ alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
+ alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
+ alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
+ * alenZ0*alenZ0)
+ am0 = amassFact0*alen0
+ amass(kblock,1,1) = am0
+ amass(kblock,2,2) = am0
+ amass(kblock,3,3) = am0
+ amass(kblock,4,4) = am0
+ amass(kblock,5,5) = am0
+```
+
+
+
+```txt
+amass(kblock,6,6) = am0
+end do
+else if ( lflags(iOpCode) .eq.
+* jIntForceAndDtStable) then
+do kblock = 1, nblock
+ alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
+ alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
+ alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
+ alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
+* alenZ0*alenZ0)
+ vol0 = area0*alen0
+ amElem0 = two*amassFact0*alen0
+ alenX = alenX0
+* + (u(kblock,4) - u(kblock,1))
+ alenY = alenY0
+* + (u(kblock,5) - u(kblock,2))
+ alenZ = alenZ0
+* + (u(kblock,6) - u(kblock,3))
+ alen = sqrt(alenX*alenX + alenY*alenY + alenZ*alenZ)
+ area = vol0/alen
+ ak = area*eMod/alen
+c stable time increment for translations
+ dtimeStable(kblock) = factorStable*sqrt(amElem0/ak)
+c force = E * logarithmic strain *current area
+ strainLog = log(alen/alen0)
+ fElasTra = eMod*strainLog*area
+ forceTra = fElasTra
+c assemble internal load in RHS
+ rhs(kblock,1) = -forceTra
+ rhs(kblock,4) = forceTra
+c internal energy calculation
+ alenOld = svars(kblock,1)
+ fElasTraOld = svars(kblock,2)
+ energy(kblock, iElIe) = energy(kblock, iElIe) +
+```
+
+
+
+```fortran
+* half*(fElasTra+fElasTraOld)*(alen - alenOld)
+c update state variables
+ svars(kblock,1) = alen
+ svars(kblock,2) = fElasTra
+ end do
+ end if
+ end if
+c
+ return
+ end
+```
+
+
+
+
+
+# 1.2.15 VUEOS: User subroutine to define equation of state material model.
+
+# Product: Abaqus/Explicit
+
+# References
+
+• “Equation of state,” Section 25.2.1 of the Abaqus Analysis User’s Guide
+• \*EOS
+“Equation of state material,” Section 2.2.20 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine VUEOS:
+
+• can be used to define the hydrodynamic material model in which the material’s volumetric response is determined by the user-defined equation of state;
+• will be called for blocks of material calculation points for which the material definition contains a user-defined equation of state;
+• can use and update solution-dependent state variables; and
+• can use any field variables that are passed in.
+
+# User subroutine interface
+
+```python
+subroutine vueos (
+C Read only (unmodifiable) variables -
+ 1 nblock,
+ 2 jElem, kIntPt, kLayer, kSecPt,
+ 3 steTime, totalTime, dt, cmname,
+ 4 nstatev, nfieldv, nprops,
+ 5 props, tempOld, tempNew, fieldOld, fieldNew,
+ 6 stateOld, charLength, coordMp,
+ 7 densityMean, refDensity, densityNew,
+ 8 dkk, Em,
+C Write only (modifiable) variables -
+ 8 press, dPdRho, dPdEm,
+ 9 stateNew )
+C
+ include 'vaba_param.inc'
+C
+ dimension props(nprops),
+ 1 tempOld(nblock),
+ 2 fieldOld(nblock, nfieldv),
+```
+
+
+
+```matlab
+3 stateOld(nblock, nstatev),
+4 tempNew(nblock),
+5 fieldNew(nblock, nfieldv),
+6 charLength(nblock), coordMp(nblock, *),
+7 densityMean(nblock), refDensity(nblock),
+8 densityNew(nblock),
+9 dkk(nblock), Em(nblock),
+1 press(nblock), dPdRho(nblock), dPdEm(nblock),
+2 stateNew(nblock)
+C
+ character*80 cmname
+C
+ do 100 km = 1, nblock
+ user coding to define/update press, dPdRho, dPdEm
+100 continue
+ return
+ end
+```
+
+# Variables to be defined
+
+press(nblock)
+
+The material point pressure stress, p.
+
+dPdRho(nblock)
+
+The derivative of the pressure with respect to the density, $\partial p / \partial \rho .$ This quantity is needed for the evaluation of the effective moduli of the material, which enters the stable time increment calculation.
+
+dPdEm(nblock)
+
+The derivative of the pressure with respect to the internal energy, $\partial p / \partial E _ { m }$ . This quantity is needed for the iterative Newton loop used outside of the user subroutine to solve for pressure.
+
+# Variable that can be updated
+
+stateNew(nblock,nstatev)
+
+State variables at each material point at the end of the increment. You define the size of this array by allocating space for it (see “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide, for more information).
+
+
+
+# Variables passed in for information
+
+nblock
+
+Number of material points to be processed in this call to VUEOS.
+
+jElem(nblock)
+
+Array of element numbers.
+
+kIntPt
+
+Integration point number.
+
+kLayer
+
+Layer number.
+
+kSecPt
+
+Section point number within the current layer.
+
+stepTime
+
+Value of time since the step began.
+
+totalTime
+
+Value of total time. The time at the beginning of the step is given by totalTime - stepTime.
+
+dt
+
+Time increment size.
+
+cmname
+
+User-specified material name, left justified. It is passed in as an uppercase character string. Some internal material models are given names starting with the “ABQ\_” character string. To avoid conflict, you should not use “ABQ\_” as the leading string for cmname.
+
+nstatev
+
+Number of user-defined state variables that are associated with this material type (you define the number as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
+
+nfieldv
+
+Number of user-defined external field variables.
+
+nprops
+
+User-specified number of user-defined material properties.
+
+props(nprops)
+
+User-supplied material properties.
+
+
+
+tempOld(nblock)
+
+Temperatures at each material point at the beginning of the increment.
+
+tempNew(nblock)
+
+Temperatures at each material point at the end of the increment.
+
+fieldOld(nblock,nfieldv)
+
+Values of the user-defined field variables at each material point at the beginning of the increment.
+
+fieldNew(nblock,nfieldv)
+
+Values of the user-defined field variables at each material point at the end of the increment.
+
+stateOld(nblock,nstatev)
+
+State variables at each material point at the beginning of the increment.
+
+charLength(nblock)
+
+Characteristic element length, which is either the default value based on the geometric mean or the user-defined characteristic element length defined in user subroutine VUCHARLENGTH. The default value is a typical length of a line across an element for a first-order element; it is half the same typical length for a second-order element. For beams, pipes, and trusses, the default value is a characteristic length along the element axis. For membranes and shells it is a characteristic length in the reference surface. For axisymmetric elements it is a characteristic length in the r–z plane only. For cohesive elements it is equal to the constitutive thickness.
+
+coordMp(nblock,\*)
+
+Material point coordinates.
+
+densityMean(nblock)
+
+The mean density.
+
+refDensity(nblock)
+
+The reference density.
+
+densityNew(nblock)
+
+The current density for this increment.
+
+dkk(nblock)
+
+The volumetric strain increment.
+
+Em(nblock)
+
+The element specific internal energy (per unit mass)
+
+Example: User subroutine VUEOS to reproduce results obtained with \*EOS, TYPE=USUP
+
+As a simple example of coding of user subroutine VUEOS, consider the following form of the Mie-Grüneisen equation of state with 0.0 and a linear dependency between pressure and internal energy:
+
+
+
+$$
+p = \rho_ {0} c _ {0} ^ {2} \eta (1 - \frac {\Gamma_ {0} \eta}{2}) + \Gamma_ {0} \rho_ {0} E _ {m},
+$$
+
+where $\eta = 1 - \rho _ { 0 } / \rho .$ . Therefore, the results obtained with user subroutine VUEOS should be the same as the results obtained with the linear $U _ { s } - U _ { p }$ type of EOS already available in “Equation of state material,” Section 2.2.20 of the Abaqus Verification Guide.
+
+The code in user subroutine VUEOS must return the pressure, , as in the above equation; the derivative of the pressure with respect to the density, $\partial p / \partial \rho ;$ and the derivative of the pressure with respect to the energy, $\partial p / \partial E _ { m }$ . For the case considered here, these values are
+
+$$
+\frac {\partial p}{\partial \rho} = \frac {\rho_ {0} ^ {2} c _ {0} ^ {2}}{\rho^ {2}} (1 - \Gamma_ {0} \eta),
+$$
+
+$$
+\frac {\partial p}{\partial E _ {m}} = \Gamma_ {0} \rho_ {0}.
+$$```fortran
+subroutine vueos (
+C Read only (unmodifiable) variables -
+ 1 nblock,
+ 2 jElem, kIntPt, kLayer, kSecPt,
+ 3 steTime, totalTime, dt, cmname,
+ 4 nstatev, nfieldv, nprops,
+ 5 props, tempOld, tempNew, fieldOld, fieldNew,
+ 6 stateOld, charLength, coordMp,
+ 7 densityMean, refDensity, densityNew,
+ 8 dkk, Em,
+C Write only (modifiable) variables -
+ 8 press, dPdRho, dPdEm,
+ 9 stateNew)
+C
+ include 'vaba_param.inc'
+C
+ dimension props(nprops),
+ 1 tempOld(nblock),
+ 2 fieldOld(nblock,nfieldv),
+ 3 stateOld(nblock,nstatev),
+ 4 tempNew(nblock),
+ 5 fieldNew(nblock,nfieldv),
+ 6 charLength(nblock), coordMp(nblock,*),
+ 7 densityMean(nblock), refDensity(nblock),
+ 8 densityNew(nblock),
+ 9 dkk(nblock), Em(nblock),
+ 1 press(nblock), dPdRho(nblock), dPdEm(nblock),
+```
+
+
+
+```txt
+2 stateNew(nblock)
+C
+ character*80 cmname
+C
+ parameter ( zero = 0.d0, one = 1.d0, half = 0.5d0 )
+C
+ c0 = props(1)
+ gamma0 = props(2)
+ c02 = c0*c0
+C
+ do k=1, nblock
+ rho0 = refDensity(k)
+ eta = one - rho0/densityNew(k)
+ f1 = rho0*c02*eta*(one-half*gamma0*eta)
+ f2 = gamma0*rho0
+ press(k) = f1 + f2*Em(k)
+C dP/dEm
+ dPdEm(k) = f2
+C dP/dRho
+ dPdRho(k) = c02*(rho0/densityNew(k))**2*(one-gamma0*eta)
+ end do
+C
+ return
+ end
+```
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_051.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_051.md
new file mode 100644
index 0000000..d852eaa
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_051.md
@@ -0,0 +1,379 @@
+
+
+# 1.2.16 VUFIELD: User subroutine to specify predefined field variables.
+
+# Product: Abaqus/Explicit
+
+# References
+
+• “Predefined fields,” Section 34.6.1 of the Abaqus Analysis User’s Guide
+• \*FIELD
+
+# Overview
+
+User subroutine VUFIELD:
+
+• allows you to prescribe predefined field variables at the nodes of a model—the predefined field variables at a node can be updated individually, or a number of field variables at the nodes can be updated simultaneously;
+• can be called for blocks of nodes for which the field variable values are defined in the subroutine;
+• ignores any field variable values specified directly;
+• can be used to modify field variable values read from a results file; and
+• can be used in conjunction with user subroutine VUSDFLD such that the field variables that are passed in from VUFIELD and interpolated to the material points can be modified (such changes are local to material point values, and nodal field variable values remain unaffected).
+
+# Updating field variables
+
+Two different methods are provided for updating field variables.
+
+# Individual variable updates
+
+By default, only one field variable is updated at a time for given nodes or a given node set in user subroutine VUFIELD. The user subroutine is called whenever a current value of a field variable is needed for the nodes that are listed in the field variable definition. This method is ideal for cases in which the field variables are independent of each other.
+
+# Simultaneous variable updates
+
+User subroutine VUFIELD can also be used to update multiple field variables simultaneously for given nodes or a given node set. This method is well-suited for cases in which there are dependencies between some of the field variables. In this case you must specify the number of field variables to be updated simultaneously, and the user subroutine will be called each time the field variable values are needed.
+
+
+
+User subroutine interface
+```txt
+SUBROUTINE VUFIELD(FIELD, NBLOCK, NFIELD, KFIELD, NCOMP, 1 KSTEP, JFLAGS, JNODEID, TIME, 2 COORDS, U, V, A)
+
+C
+INCLUDE 'VABA_PARAM.INC'
+
+C indices for the time array TIME
+PARAMETER( i_ufld_Current = 1,
+* i_ufld_Increment = 2,
+* i_ufld_Period = 3,
+* i_ufld_Total = 4 )
+
+C indices for the coordinate array COORDS
+PARAMETER( i_ufld_CoordX = 1,
+* i_ufld_CoordY = 2,
+* i_ufld_CoordZ = 3 )
+
+C indices for the displacement array U
+PARAMETER( i_ufld_SpaDisplX = 1,
+* i_ufld_SpaDisplY = 2,
+* i_ufld_SpaDisplZ = 3,
+* i_ufld_RotDisplX = 4,
+* i_ufld_RotDisplY = 5,
+* i_ufld_RotDisplZ = 6,
+* i_ufld_AcoPress = 7,
+* i_ufld_Temp = 8 )
+
+C indices for the velocity array V
+PARAMETER( i_ufld_SpaVelX = 1,
+* i_ufld_SpaVelY = 2,
+* i_ufld_SpaVelZ = 3,
+* i_ufld_RotVelX = 4,
+* i_ufld_RotVelY = 5,
+* i_ufld_RotVelZ = 6,
+* i_ufld_DAcoPress = 7,
+* i_ufld_DTemp = 8 )
+
+C indices for the acceleration array A
+PARAMETER( i_ufld_SpaAccelX = 1,
+```
+
+
+
+```c
+* i_ufld_SpaAccelY = 2,
+* i_ufld_SpaAccelZ = 3,
+* i_ufld_RotAccelX = 4,
+* i_ufld_RotAccelY = 5,
+* i_ufld_RotAccelZ = 6,
+* i_ufld_DDAcoPress = 7,
+* i_ufld_DDTemp = 8)
+
+C indices for JFLAGS
+PARAMETER( i_ufld_kInc = 1,
+* i_ufld_kPass = 2 )
+
+C
+DIMENSION FIELD(NBLOCK, NCOMP, NFIELD)
+DIMENSION JFLAGS(2), JNODEID(NBLOCK), TIME(4),
+* COORDS(3, NBLOCK)
+DIMENSION U(8, NBLOCK), V(8, NBLOCK), A(8, NBLOCK)
+
+C
+user coding to define FIELD
+
+RETURN
+END
+```
+
+# Variable to be defined
+
+# FIELD(NBLOCK,NCOMP,NFIELD)
+
+Array of field variable values at a collective number of nodes NBLOCK (see NBLOCK below). When updating one field variable at a time, only the value of the specified field variable KFIELD must be returned. In this case NFIELD is passed into user subroutine VUFIELD with a value of 1, and FIELD is thus dimensioned as FIELD(NBLOCK,NCOMP,1). When updating all field variables simultaneously, the values of the specified number of field variables must be returned. In this case FIELD is dimensioned as FIELD(NBLOCK,NCOMP,NFIELD), where NFIELD is the number of field variables specified and KFIELD, which is set to −1, has no meaning.
+
+If fields are applied to nodes that are not part of pipe, beam, or shell elements, only one value of each field variable is required (NCOMP=1), and the user subroutine is invoked in a single pass. For nodes that are part of pipe, beam, or shell elements, VUFIELD is invoked in two passes per increment for these elements, and the number of values to be returned depends on the mode of temperature and field variable input selected for the beam or shell section. The following cases are possible:
+
+1. Field variables are given as values at the points on the shell or beam section. For a beam section the number of values required is determined by the particular section type specified, as described in “Beam cross-section library,” Section 29.3.9 of the Abaqus Analysis User’s Guide. For a shell section temperatures and field variables are given as values at n equally spaced points through
+
+
+
+each layer of a shell section. In the first pass NCOMP is passed in with the value of 1 to define the field variable values at the first point. The second pass is used to define field variables at the remaining points.
+
+2. Field variables for the shell or beam section are given as values at the origin of the cross-section together with gradients along the cross-section. The number of gradient values required is 2 for three-dimensional beams, 1 for two-dimensional beams, and 1 for shells. In the first pass NCOMP is passed in with the value of 1 to define the field variable values at the origin of the cross-section. The gradients are defined in the second pass.
+
+Because field variables can also be defined directly, it is important to understand the hierarchy used in situations with conflicting information (see “Predefined fields,” Section 34.6.1 of the Abaqus Analysis User’s Guide).
+
+When the array FIELD is passed into user subroutine VUFIELD, it will contain either the field variable values from the previous increment or those values obtained from the results file if this method was used. You can then modify these values within this subroutine.
+
+# Variables passed in for information
+
+# NBLOCK
+
+User-specified number of nodes to be processed as a block in this call to VUFIELD. The value is equal to the total number of nodes given in a node set when blocking is disabled. When blocking is enabled, NBLOCK is equal to a predefined number set in Abaqus/Explicit. You can also modify NBLOCK by specifying a blocking size in the Abaqus/Explicit analysis.
+
+# NFIELD
+
+User-specified number of field variables to be updated. The default value is 1.
+
+# KFIELD
+
+User-specified field variable number. This variable is meaningful only when updating individual field variables at a time; otherwise, the value is set to −1.
+
+# NCOMP
+
+Maximum number of section values to be defined for any node in the model in the current pass. The first pass to user subroutine VUFIELD has NCOMP passed in with the value of 1.
+
+# KSTEP
+
+Current step number.
+
+# JFLAGS(i\_ufld\_kInc)
+
+Increment number for step KSTEP.
+
+# JFLAGS(i\_ufld\_kPass)
+
+This flag is equal to 1 for the first pass to user subroutine VUFIELD and is equal to 2 for the second pass.
+
+
+
+# JNODEUID(NBLOCK)
+
+Array for user-defined node numbers. This array is dimensioned based on the size of NBLOCK, and the contained node numbers are identical to those defined in the input file. You can perform additional interdependent field variable operations by using nodal indices stored in this array.
+
+# TIME(4)
+
+Array for information of analysis time. You can retrieve any time information from this array by using the parameters given above. TIME(i\_ufld\_Current) stores the current analysis time, TIME(i\_ufld\_Increment) gives the time increment at this instance, TIME(i\_ufld\_Period) is the time period of the current step, and TIME(i\_ufld\_Total) is the total analysis time up to this point. You can use this time information to perform possible time-dependent field variable operations.
+
+# COORDS(3,NBLOCK)
+
+Coordinates for nodes in the array JNODEUID. This array stores current coordinates of nodes in which the order of coordinates stored corresponds to the order of nodes listed in the array JNODEUID. The coordinates can be retrieved by using the parameters given above. You can make use of COORDS to define possible position-dependent field variable operations.
+
+# U(8,NBLOCK), V(8,NBLOCK), and A(8,NBLOCK)
+
+Arrays containing solution variables of displacements, rotations, temperatures, and pressures and their corresponding temporal derivatives. The order in which these solutions are stored follows the order defined in the array JNODEUID. For a specific node its solution variables can be retrieved by using the parameter indices given above. Depending on the degrees of freedom, some solution variables are not valid for a given node. The displacement values correspond to the current increment. However, the acceleration is from a configuration that is one increment behind, and the velocity is such that it is consistent with the displacement increment and the time increment between the two successive configurations.
+
+
+
+
+
+# 1.2.17 VUFLUIDEXCH: User subroutine to define the mass flow rate/heat energy flow rate for fluid exchange.
+
+# Product: Abaqus/Explicit
+
+# References
+
+• “Fluid exchange definition,” Section 11.5.3 of the Abaqus Analysis User’s Guide
+• \*FLUID EXCHANGE
+• \*FLUID EXCHANGE ACTIVATION
+• \*FLUID EXCHANGE PROPERTY
+
+# Overview
+
+User subroutine VUFLUIDEXCH:
+
+• can be used to define mass flow rate and/or heat energy flow rate for fluid exchange;
+• can be used when built-in fluid exchange property types cannot satisfactorily model the mass/heat energy flow;
+• can use and update solution-dependent state variables;
+• can use any field variables that are passed in; and
+• requires that the derivatives of mass/heat energy flow rates be defined with respect to pressure and temperature in the primary and secondary fluid cavities.
+
+# Conventions for defining mass flow/heat energy flow rate
+
+A positive mass/heat energy flow rate indicates flow from the primary fluid cavity to the secondary fluid cavity. A negative value for mass flow rate will be ignored if the fluid exchange is between a cavity and its environment.
+
+# User subroutine interface
+
+```txt
+subroutine vufluidexch(
+C Read only (unmodifiable)variables -
+1 nstatev, nfieldv, nprops,
+2 steppTime, totalTime, dt,
+3 jCavType, fluExchName, effArea, amplitude,
+4 props, lExchEnv, pcavNew, pcavOld,
+5 ctempNew, ctempOld, cvol, cmass,
+6 rMix, CpMix, DCpDtemp,
+7 field, stateOld,
+```
+
+
+
+```python
+C Write only (modifiable) variables
+ 8 stateNew, rMassRate, rEneRate,
+ 9 DMassRateDPcav, DMassRateDTemp,
+ * DEneRateDPcav, DEneRateDTemp)
+c
+ include 'vaba_param.inc'
+c
+ dimension props(nprops),
+ 1 pcavNew(2), pcavOld(2),
+ 2 ctempNew(2), ctempOld(2), cvol(2), cmass(2),
+ 3 rMix(2), CpMix(2), dCpDtemp(2),
+ 4 field(nfieldv),
+ 5 stateOld(nstatev), stateNew(nstatev),
+ 6 DMassRateDPcav(2), DMassRateDTemp(2),
+ 7 DEneRateDPcav(2), DEneRateDTemp(2)
+
+c Fluid cavity type
+ parameter( iHydraulic = 1,
+ * iAdiabaticGas = 2,
+ * iIsothermalGas = 3)
+
+ character*80 fluExchName
+
+c User coding to calculate mass flow rate,
+c heat energy flow rate and its derivatives with respect
+c to fluid cavity pressure and temperature.
+
+ return
+ end
+```
+
+# Variables to be defined
+
+# rMassRate
+
+Mass flow rate. The mass flow rate is negative if the flow is into the primary cavity.
+
+# DMassRateDPcav(2)
+
+Derivative of mass flow rate with respect to pressure in primary and secondary fluid cavities.
+
+# DMassRateDTemp(2)
+
+Derivative of mass flow rate with respect to temperature in primary and secondary fluid cavities.
+
+# rEneRate
+
+Heat energy flow rate. The energy flow rate is negative if the flow is into the primary cavity.
+
+
+
+# DEneRateDPcav(2)
+
+Derivative of heat energy flow rate with respect to pressure in primary and secondary fluid cavities.
+
+# DEneRateDTemp(2)
+
+Derivative of heat energy flow rate with respect to temperature in primary and secondary fluid cavities.
+
+# Variable that can be updated
+
+# stateNew(nstatev)
+
+State variable for fluid exchange at the end of the increment. You define the size of this array by allocating space for it (see “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide, for more information).
+
+# Variables passed in for information
+
+# nstatev
+
+Number of user-defined state variables that are associated with this fluid exchange (you define this as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
+
+# nfieldv
+
+Number of user-defined external field variables.
+
+# nprops
+
+User-specified number of user-defined fluid exchange properties required to define mass/heat energy flow rate.
+
+# stepTime
+
+Value of time since the step began.
+
+# totalTime
+
+Value of total time. The time at the beginning of the step is given by totalTime−stepTime.
+
+# dt
+
+Time increment size.
+
+# jCavType
+
+Indicator of fluid cavity type: 1 for fluid cavity with hydraulic fluids, 2 for fluid cavity with adiabatic gases, and 3 for fluid cavity with isothermal gases.
+
+# fluExchName
+
+User-specified fluid exchange name.
+
+# effArea
+
+Effective area for fluid exchange.
+
+
+
+amplitude
+
+Current value of the amplitude referenced for this fluid exchange. You must multiply the flow rates by the current amplitude value within the user subroutine if the amplitude is required.
+
+props(nprop)
+
+User-defined fluid exchange properties.
+
+lExchEnv
+
+The fluid exchange is to the environment if lExchEnv=1 and to another fluid cavity if lExchEnv=0.
+
+pcavNew(2)
+
+Pressure in primary and secondary fluid cavities at the end of the increment.
+
+pcavOld(2)
+
+Pressure in primary and secondary fluid cavities at the beginning of the increment.
+
+ctempNew(2)
+
+Temperature in primary and secondary fluid cavities at the end of the increment.
+
+ctempOld(2)
+
+Temperature in primary and secondary fluid cavities at the beginning of the increment.
+
+cvol(2)
+
+Volume of primary and secondary fluid cavities.
+
+cmass(2)
+
+Mass of fluid in primary and secondary fluid cavities.
+
+rMix(2)
+
+Gas constant of mixture in primary and secondary fluid cavities.
+
+CpMix(2)
+
+Specific heat of mixture in primary and secondary fluid cavities.
+
+DCpDtemp(2)
+
+Derivative of specific heat with respect to temperature for primary and secondary fluid cavities.
+
+field(nfieldv)
+
+Field variables at orifice.
+
+stateOld(nstatev)
+
+State variables for fluid exchange at the beginning of the increment.
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_052.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_052.md
new file mode 100644
index 0000000..c14268f
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_052.md
@@ -0,0 +1,371 @@
+
+
+# 1.2.18 VUFLUIDEXCHEFFAREA: User subroutine to define the effective area for fluid exchange.
+
+# Product: Abaqus/Explicit
+
+# References
+
+• “Fluid exchange definition,” Section 11.5.3 of the Abaqus Analysis User’s Guide
+• \*FLUID EXCHANGE
+
+# Overview
+
+User subroutine VUFLUIDEXCHEFFAREA:
+
+• can be used to define an effective area for fluid exchange that depends on the material state in the underlying elements on the fluid exchange surface;
+• will be called for blocks of material calculation points on the fluid exchange surface;
+• can be used only if the specified surface over which the fluid exchange occurs is a surface defined over membrane elements; and
+• can be used with any fluid exchange property type.
+
+# Defining effective area
+
+The contribution of each material point can be defined as a function of:
+
+• the original area associated with the material point;
+• the current material state in the underlying elements; and
+• the temperature and pressure in the primary fluid cavity and the secondary fluid cavity or environment.
+
+The effective area for fabric materials can depend on the nominal strain in the yarn directions and the change in angle between the two yarn directions, as well as the current angle between the two yarn directions. For nonfabric materials the effective area can depend on the material point strain.
+
+# User subroutine interface
+
+```csv
+subroutine vufluidexcheffarea(
+C Read only (unmodifiable) variables -
+1 nblock, nprop, props,
+2 steppTime, totalTime, fluExchName,
+3 cMatName, lFabric, braidAngle,
+4 strain, origArea,
+```
+
+
+
+```txt
+5 pcav, ctemp,
+C Write only (modifiable) variables
+6 effArea)
+c
+ include 'vaba_param.inc'
+c
+ parameter (ndir = 3, nshr=1)
+c
+c pointers for retrieving fabric constitutive strains
+ parameter( iFiberStrain1 = 1,
+ * iFiberStrain2 = 2,
+ * iFiberChangeAng = 4)
+c
+ dimension props(nprop),
+ 1 braidAngle(nblock),
+ 2 strain(nblock, ndir+nshr),
+ 3 origArea(nblock),
+ 4 pcav(2),ctemp(2),
+ 5 effArea(nblock)
+
+ character*80 fluExchName, cMatName
+
+c do k = 1, nblock
+c User coding to update effArea(k) = area associated with
+c material point contributing to area for fluid exchange
+c (leakage).
+c end do
+
+ return
+ end
+```
+
+# Variable to be defined
+
+# effArea(nblock)
+
+Area associated with the material point contributing to the total effective area for fluid exchange. The subroutine is called with effArea set to the current area associated with the material point and should be updated to reflect the area that contributes to fluid exchange.
+
+# Variables passed in for information
+
+# nBlock
+
+Number of material points to be processed in this call to VUFLUIDEXCHEFFAREA.
+
+
+
+# nprop
+
+User-specified number of user-defined fluid exchange properties required to define the effective area.
+
+# props(nprop)
+
+User-defined fluid exchange properties.
+
+# stepTime
+
+Value of time since the step began.
+
+# totalTime
+
+Value of total time. The time at the beginning of the step is given by totalTime−stepTime.
+
+# fluExchName
+
+User-specified fluid exchange name.
+
+# cMatName
+
+User-specified material name associated with material points processed in this call.
+
+# lFabric
+
+Flag indicating whether the subroutine is called for material points on the fluid exchange surface with a fabric material (lFabric=1 if fabric material, lFabric=0 otherwise).
+
+# braidAngle(nblock)
+
+Angle in radians between the two yarn directions for fabric materials.
+
+# strain(nblock, ndir+nshr)
+
+Fabric constitutive strains (nominal strain in the yarn directions and change in angle between the two yarn directions) or strains for nonfabric materials at current location.
+
+# origArea(nblock)
+
+Original area associated with current material point.
+
+# pcav(2)
+
+Absolute pressure in primary and secondary (or ambient) fluid cavities at the start of the increment.
+
+# ctemp(2)
+
+Temperature in primary and secondary (or ambient) fluid cavities at the start of the increment.
+
+
+
+
+
+# 1.2.19 VUHARD: User subroutine to define the yield surface size and hardening parameters for isotropic plasticity or combined hardening models.
+
+# Product: Abaqus/Explicit
+
+# References
+
+• “Classical metal plasticity,” Section 23.2.1 of the Abaqus Analysis User’s Guide
+• “Models for metals subjected to cyclic loading,” Section 23.2.2 of the Abaqus Analysis User’s Guide
+• \*CYCLIC HARDENING
+• \*PLASTIC
+• “Deformation of a sandwich plate under CONWEP blast loading,” Section 9.1.9 of the Abaqus Example Problems Guide
+• “VUHARD,” Section 4.1.35 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine VUHARD:
+
+• is called at all material points of elements for which the material definition includes user-defined isotropic hardening or cyclic hardening for metal plasticity;
+• can be used to define a material’s isotropic yield behavior;
+• can be used to define the size of the yield surface in a combined hardening model;
+• can include material behavior dependent on field variables or state variables; and
+• requires that the derivatives of the yield stress (or yield surface size in combined hardening models) be defined with respect to the appropriate independent variables, such as strain, strain rate, and temperature.
+
+# User subroutine interface
+
+```c
+subroutine vuhard(
+C Read only -
+* nblock,
+* jElem, kIntPt, kLayer, kSecPt,
+* lAnneal, stepTime, totalTime, dt, cmname,
+* nstatev, nfieldv, nprops,
+* props, tempOld, tempNew, fieldOld, fieldNew,
+* stateOld,
+* eqps, eqpsRate,
+C Write only -
+```
+
+
+
+```txt
+* yield, dyieldDtemp, dyieldDeqps,
+* stateNew)
+C
+ include 'vaba_param.inc'
+C
+ dimension props(nprops), tempOld(nblock), tempNew(nblock),
+ 1 fieldOld(nblock, nfieldv), fieldNew(nblock, nfieldv),
+ 2 stateOld(nblock, nstatev), eqps(nblock), eqpsRate(nblock),
+ 3 yield(nblock), dyieldDtemp(nblock), dyieldDeqps(nblock, 2),
+ 4 stateNew(nblock, nstatev), jElem(nblock)
+C
+ character*80 cmname
+C
+ do 100 km = 1, nblock
+ user coding
+100 continue
+C
+ return
+ end
+```
+
+# Variables to be defined
+
+# yield(nblock)
+
+Array containing the yield stress (for isotropic plasticity) or yield surface size (for combined hardening) at the material points.
+
+# dyieldDeqps(nblock,1)
+
+Array containing the derivative of the yield stress or yield surface size with respect to the equivalent plastic strain at the material points.
+
+# dyieldDeqps(nblock,2)
+
+Array containing the derivative of the yield stress with respect to the equivalent plastic strain rate at the material points.
+
+# dyieldDtemp(nblock)
+
+Array containing the derivative of the yield stress or yield surface size with respect to temperature at the material points. This quantity is required only in adiabatic and fully coupled temperature-displacement analyses.
+
+# stateNew(nblock,nstatev)
+
+Array containing the state variables at the material points at the end of the increment. The allocation of this array is described in “Solution-dependent state variables” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide.
+
+
+
+nblock
+
+Number of material points to be processed in this call to VUHARD.
+
+jElem(nblock)
+
+Array of element numbers.
+
+kIntPt
+
+Integration point number.
+
+kLayer
+
+Layer number (for composite shells).
+
+kSecPt
+
+Section point number within the current layer.
+
+lanneal
+
+Flag indicating whether the routine is being called during an annealing process. lanneal=0 indicates that the routine is being called during a normal mechanics increment. lanneal=1 indicates that this is an annealing process and the internal state variables, stateNew, should be reinitialized if necessary. Abaqus/Explicit will automatically set the stresses, stretches, and state to a value of zero during the annealing process.
+
+stepTime
+
+Value of time since the step began.
+
+totalTime
+
+Value of total time. The time at the beginning of the step is given by totalTime-stepTime.
+
+dt
+
+Time increment size.
+
+cmname
+
+Material name, left justified. It is passed in as an uppercase character string. Some internal material models are given names starting with the “ABQ\_” character string. To avoid conflict, “ABQ\_” should not be used as the leading string for cmname.
+
+nstatev
+
+Number of user-defined state variables that are associated with this material type (see “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
+
+nfieldv
+
+Number of user-defined external field variables.
+
+
+
+nprops
+
+User-specified number of user-defined material properties.
+
+tempOld(nblock)
+
+Temperatures at the material points at the beginning of the increment.
+
+tempNew(nblock)
+
+Temperatures at the material points at the end of the increment.
+
+fieldOld(nblock,nfieldv)
+
+Values of the user-defined field variables at the material points at the beginning of the increment.
+
+fieldNew(nblock,nfieldv)
+
+Values of the user-defined field variables at the material points at the end of the increment.
+
+stateOld(nblock,nstatev)
+
+State variables at the material points at the beginning of the increment.
+
+eqps(nblock)
+
+Equivalent plastic strain at the material points.
+
+eqpsRate(nblock)
+
+Equivalent plastic strain rate at the material points.
+
+
+
+# 1.2.20 VUINTER: User subroutine to define the interaction between contact surfaces.
+
+# Product: Abaqus/Explicit
+
+# References
+
+• “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide
+• \*SURFACE INTERACTION
+• “VUINTER,” Section 4.1.36 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine VUINTER:
+
+• can be used to define the mechanical and thermal interaction between contacting surfaces;
+• must provide the entire definition of the interaction between the contacting surfaces;
+• can use and update solution-dependent state variables; and
+• must be used with the penalty contact constraint algorithm.
+
+# Terminology
+
+The use of user subroutine VUINTER requires familiarity with the following terminology.
+
+# Surface node numbers
+
+The “surface node number” refers to the position of a particular node in the list of nodes on the surface. For example, there are nSlvNod nodes on the slave surface. Number nSlvNod, is the surface node number of the nth node in this list; jSlvUid is the user-defined global number of this node. An Abaqus/Explicit model can be defined in terms of an assembly of part instances (see “Defining an assembly,” Section 2.10.1 of the Abaqus Analysis User’s Guide). In such models a node number in jSlvUid is an internally generated node number. If the original node number and part instance name are required, call the utility routine VGETPARTINFO (see “Obtaining part information,” Section 2.1.5).
+
+# Local coordinate system
+
+The array alocaldir defines the direction cosines of a local coordinate system for each slave node. The first local direction corresponds to the contact normal direction from the perspective of the slave node. For a two-dimensional VUINTER model the second local direction is the tangent direction defined by the cross product of the vector into the plane of the model (0., 0., −1.0) and the slave normal. For a three-dimensional VUINTER model the second and third local directions correspond to two orthogonal tangent directions $\mathbf { t } _ { 1 }$ and $\mathbf { t } _ { 2 }$ , which are set as follows:
+
+
+
+• If the master surface is a cylindrical analytical surface, the second local direction corresponds to the generator direction (see “Analytical rigid surface definition,” Section 2.3.4 of the Abaqus Analysis User’s Guide), and the third local direction is the cross product of the first and second local directions.
+• If the master surface is an analytical surface of revolution, the third local direction corresponds to the hoop direction, and the second local direction is the cross product of the third and first local directions.
+• If the master surface is a three-dimensional, element-based surface, the tangent directions are based on the slave normal, using the standard convention for calculating surface tangents (see “Conventions,” Section 1.2.2 of the Abaqus Analysis User’s Guide).
+
+For the two cases listed above involving three-dimensional analytical surfaces, the local tangent directions will reflect a rotation of the master surface. For the last case (three-dimensional, element-based master surface) the tangent directions may not follow the rotation of either the master or slave surfaces; for example, the local system would remain fixed with respect to the global system if a slave node and its surrounding facets rotate about an axes parallel to the slave normal.
+
+The 2 × 2 array stored in drot for each slave node represents the incremental rotation of the tangent directions within the tangent plane corresponding to the tracked point of a three-dimensional master surface. (This incremental rotation array is equal to a unit matrix if nDir is equal to 2.) This incremental rotation matrix is provided so that vector- or tensor-valued state variables defined within the tangent plane can be rotated in this subroutine. For example, the second and third components of the rdisp array (i.e., the relative slip components) are rotated by this amount before VUINTER is called. However, as already mentioned, the rotation of the tangent directions may not reflect a physical rotation of the master or slave surface.
+
+# Conventions for heat flux and stress
+
+A positive flux indicates heat flowing into a surface, and a negative flux denotes heat leaving the surface. Flux must be specified for both surfaces, and they need not be equal and opposite so that effects such as frictional dissipation and differential surface heating can be modeled.
+
+A positive normal stress denotes a pressure directed into the surface (opposite the local normal direction). Positive shear stresses denote shear tractions in the direction of the local surface tangents.
+
+# User subroutine interface
+
+```txt
+subroutine vuinter(
+C Write only
+ 1 sfd, scd, spd, svd,
+C Read/Write -
+ 2 stress, fluxSlv, fluxMst, sed, statev,
+C Read only -
+ 3 kStep, kInc, nFacNod, nSlvNod, nMstNod, nSurfDir,
+ 4 nDir, nStateVar, nProps, nTemp, nPred, numDefTfv,
+```
diff --git a/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_053.md b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_053.md
new file mode 100644
index 0000000..6d26256
--- /dev/null
+++ b/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_053.md
@@ -0,0 +1,390 @@
+
+
+```julia
+5 jSlvUid, jMstUid, jConMstid, timStep, timGlb,
+6 dTimCur, surfInt, surfSlv, surfMst,
+7 rdisp, drdisp, drot, stiffDflt, condDflt,
+8 shape, coordSlv, coordMst, alocaldir, props,
+9 areaSlv, tempSlv, dtempSlv, preDefSlv, dpreDefSlv,
+1 tempMst, dtempMst, preDefMst, dpreDefMst)
+C
+ include `vaba_param.inc'
+C
+ character*80 surfInt, surfSlv, surfMst
+C
+ dimension props(nProps), statev(nStateVar,nSlvNod),
+ 1 drot(2,2,nSlvNod), sed(nSlvNod), sfd(nSlvNod),
+ 2 scd(nSlvNod), spd(nSlvNod), svd(nSlvNod),
+ 3 rdisp(nDir,nSlvNod), drdisp(nDir,nSlvNod),
+ 4 stress(nDir,nSlvNod), fluxSlv(nSlvNod),
+ 5 fluxMst(nSlvNod), areaSlv(nSlvNod),
+ 6 stiffDflt(nSlvNod), condDflt(nSlvNod),
+ 7 alocaldir(nDir,nDir,nSlvNod), shape(nFacNod,nSlvNod),
+ 8 coordSlv(nDir,nSlvNod), coordMst(nDir,nMstNod),
+ 9 jSlvUid(nSlvNod), jMstUid(nMstNod),
+ 1 jConMstid(nFacNod,nSlvNod), tempSlv(nSlvNod),
+ 2 dtempSlv(nSlvNod), preDefSlv(nPred,nSlvNod),
+ 3 dpreDefSlv(nPred,nSlvNod), tempMst(numDefTfv),
+ 4 dtempMst(numDefTfv), preDefMst(nPred,numDefTfv),
+ 5 dpreDefMst(nPred,numDefTfv)
+
+ user coding to define stress,
+ and, optionally, fluxSlv, fluxMst, statev, sed, sfd, scd, spd,
+ and svd
+
+ return
+ end
+```
+
+# Variable to be defined
+
+stress(nDir, nSlvNod)
+
+On entry this array contains the stress at the interface during the previous time increment. It must be updated to the stress at the interface in the current time increment.
+
+
+
+# fluxSlv(nSlvNod)
+
+On entry this array contains the flux entering the slave surface during the previous time increment. It must be updated to the flux entering the slave surface during the current increment.
+
+# fluxMst(nSlvNod)
+
+On entry this array contains the flux entering the master surface during the previous time increment. It must be updated to the flux entering the master surface during the current time increment.
+
+# sfd(nSlvNod)
+
+This array can be updated to contain the increment in frictional dissipation at each node (units of energy per unit area). These values contribute to the output variables SFDR and ALLFD and have no effect on other solution variables.
+
+# scd(nSlvNod)
+
+This array can be updated to contain the increment in creep dissipation at each node (units of energy per unit area). These values contribute to the output variables SFDR and ALLCD and have no effect on other solution variables.
+
+# spd(nSlvNod)
+
+This array can be updated to contain the increment in plastic dissipation at each node (units of energy per unit area). These values contribute to the output variables SFDR and ALLPD and have no effect on other solution variables.
+
+# svd(nSlvNod)
+
+This array can be updated to contain the increment in viscous dissipation at each node (units of energy per unit area). These values contribute to the output variables SFDR and ALLVD and have no effect on other solution variables.
+
+# sed(nSlvNod)
+
+On entry this array contains the elastic energy density at the slave nodes at the beginning of the increment. It can be updated to contain the elastic energy density at the end of the current time increment. These values contribute to the output variable ALLSE and have no effect on other solution variables.
+
+# statev(nstateVar, nSlvNod)
+
+This array contains the user-defined solution-dependent state variables for all the nodes on the slave surface. You define the size of this array (see “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide, for more information). This array will be passed in containing the values of these variables prior to the call to user subroutine VUINTER. If any of the solution-dependent state variables is being used in conjunction with the surface interaction, it must be updated in this subroutine.
+
+
+
+kStep
+
+Step number.
+
+kInc
+
+Increment number.
+
+nFacNod
+
+Number of nodes on each master surface facet. nFacNod is 2 for two-dimensional surfaces, and nFacNod is 4 for three-dimensional surfaces (the first and last nodes are the same for triangular facets). If the master surface is an analytical rigid surface, this variable is passed in as 0.
+
+nSlvNod
+
+Number of slave nodes.
+
+nMstNod
+
+Number of master surface nodes, if the master surface is made up of facets. If the master surface is an analytical rigid surface, this variable is passed in as 0.
+
+nSurfDir
+
+Number of tangent directions at the contact points (nSurfDir = nDir - 1).
+
+nDir
+
+Number of coordinate directions at the contact points. (In a three-dimensional model nDir will be 2 if the surfaces in the contact pair are two-dimensional analytical rigid surfaces or are formed by two-dimensional elements.)
+
+nStateVar
+
+Number of user-defined state variables.
+
+nProps
+
+User-specified number of property values associated with this surface interaction model.
+
+nTemp
+
+1 if the temperature is defined and 0 if the temperature is not defined.
+
+nPred
+
+Number of predefined field variables.
+
+numDefTfv
+
+Equal to nSlvNod if the master surface is made up of facets. If the master surface is an analytical rigid surface, this variable is passed in as 1.
+
+
+
+jSlvUid(nSlvNod)
+
+This array lists the user-defined global node numbers (or internal node numbers for models defined in terms of an assembly of part instances) of the nodes on the slave surface.
+
+jMstUid(nMstNod)
+
+This array lists the user-defined global node numbers (or internal node numbers for models defined in terms of an assembly of part instances) of the nodes on the master surface. If the master surface is an analytical rigid surface, this array is passed in as a dummy array.
+
+jConMstid(nFacNod, nSlvNod)
+
+This array lists the surface node numbers of the master surface nodes that make up the facet onto which each slave node projects. If the master surface is an analytical rigid surface, this array is passed in as a dummy array.
+
+timStep
+
+Value of step time.
+
+timGlb
+
+Value of total time.
+
+dtimCur
+
+Current increment in time from $t = t _ { c u r r } - \Delta t { \mathrm { t } } 0 t = t _ { c u r r } .$
+
+surfInt
+
+User-specified surface interaction name, left justified.
+
+surfSlv
+
+Slave surface name.
+
+surfMst
+
+Master surface name.
+
+rdisp(nDir, nSlvNod)
+
+An array containing the relative positions between the two surfaces. The first component is the relative position of the slave node, with respect to the master surface, in the normal direction (a positive value indicates a penetration, and a negative value indicates a gap). The second and third components, if applicable, are the accumulated incremental relative tangential displacements of the slave node, measured from the beginning of the step in which the contact pair is defined. The local directions in which the relative displacements are defined are stored in alocaldir. If the master surface is an analytical surface, the elements in rdisp are set to r\_MaxVal for the slave nodes that are far from the master surface.
+
+
+
+drdisp(nDir, nSlvNod)
+
+An array containing the increments in relative positions between the two surfaces during the current time increment. If the master surface is an analytical surface, the elements in drdisp are set to r\_MaxVal for the slave nodes that are far from the master surface.
+
+drot(2, 2, nSlvNod)
+
+Rotation increment matrix. This matrix represents the incremental rotation of the local surface tangent directions for a three-dimensional surface. This rotation matrix for each slave node is defined as a unit matrix for two-dimensional surfaces. If the master surface is an analytical surface, the elements in drot are set to r\_MaxVal for the slave nodes that are far from the master surface.
+
+stiffDflt(nSlvNod)
+
+Values of the default penalty stiffnesses for each slave node (units of FL ).
+
+condDflt(nSlvNod)
+
+Values of the default penalty conductances for each slave node (units of $J \theta ^ { - 1 } \mathrm { T } ^ { - 1 } )$ .
+
+shape(nFacNod, nSlvNod)
+
+For each contact point this array contains the shape functions of the nodes of its master surface facet, evaluated at the location of the contact point. If the master surface is an analytical rigid surface, this array is passed in as a dummy array.
+
+coordSlv(nDir, nSlvNod)
+
+Array containing the nDir components of the current coordinates of the slave nodes.
+
+coordMst(nDir, nMstNod)
+
+Array containing the nDir components of the current coordinates of the master nodes. If the master surface is an analytical rigid surface, this array is passed in as the coordinates of the contact points on the master surface.
+
+alocaldir(nDir, nDir, nSlvNod)
+
+Direction cosines of the local surface coordinate system. The first array index corresponds to the components of the local directions, and the second array index corresponds to the local direction number. The first direction (alocaldir(1..nDir,1,...)) is the normal to the surface. The second direction (alocaldir(1..nDir,2,...)) is the first surface tangent. For a three-dimensional surface, the third direction (alocaldir(1..3,3,...)) is the second surface tangent. If the master surface is an analytical rigid surface, the numbers in alocaldir are valid only if the corresponding parts in rdisp are valid (i.e., not equal to r\_MaxVal).
+
+props(nProps)
+
+User-specified vector of property values to define the behavior between the contacting surfaces.
+
+areaSlv(nSlvNod)
+
+Area associated with the slave nodes (equal to 1 for node-based surface nodes).
+
+
+
+tempSlv(nSlvNod)
+
+Current temperature at the slave nodes.
+
+dtempSlv(nSlvNod)
+
+Increment in temperature during the previous time increment at the slave nodes.
+
+preDefSlv(nPred, nSlvNod)
+
+Current user-specified predefined field variables at the slave nodes (initial values at the beginning of the analysis and current values during the analysis).
+
+dpreDefSlv(nPred, nSlvNod)
+
+Increment in the predefined field variables at the slave nodes during the previous time increment.
+
+tempMst(numDefTfv)
+
+Current temperature at the nearest points on the master surface.
+
+dtempMst(numDefTfv)
+
+Increment in temperature during the previous time increment at the nearest points on the master surface.
+
+preDefMst(nPred, numDefTfv)
+
+Current user-specified predefined field variables at the nearest points on the master surface (initial values at the beginning of the analysis and current values during the analysis).
+
+dpreDefMst(nPred, numDefTfv)
+
+Increment in the predefined field variables during the previous time increment at the nearest points on the master surface.
+
+
+
+# 1.2.21 VUINTERACTION: User subroutine to define the contact interaction between surfaces with the general contact algorithm.
+
+Product: Abaqus/Explicit
+
+# References
+
+• “User-defined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide
+• \*SURFACE INTERACTION
+• “VUINTERACTION,” Section 4.1.37 of the Abaqus Verification Guide
+
+# Overview
+
+User subroutine VUINTERACTION:
+
+• can be used to define the mechanical and thermal interaction between contact surfaces;
+• must provide the entire definition of the interaction between the contact surfaces;
+• can utilize a user-specified tracking thickness to determine potential points of interaction on a surface (and thus which nodes should be passed into the subroutine);
+• can use and update solution-dependent state variables for node-to-face contact and node-toanalytical rigid surface contact; and
+• must be used with the general contact algorithm.
+
+# Terminology
+
+The use of user subroutine VUINTERACTION requires familiarity with the following terminology.
+
+# Tracking thickness
+
+For efficiency, user subroutine VUINTERACTION considers only regions of two surfaces that are likely to be in contact or come into contact in a given increment. This likelihood is defined by a tracking thickness: only portions of surfaces separated by less than the tracking thickness in a given increment are passed into the subroutine; portions of the surfaces with a separation larger than the tracking thickness are ignored for the current increment. Surface thicknesses are accounted for in the separation calculations.
+
+Abaqus/Explicit provides an internal default value for the tracking thickness, but a nondefault value can be specified; see “Tracking thickness when VUINTER or VUINTERACTION is used” in “Userdefined interfacial constitutive behavior,” Section 37.1.6 of the Abaqus Analysis User’s Guide. The tracking thickness is passed into VUINTERACTION using the variable rData(4).
+
+# Proximity points
+
+A proximity point is a potential point of interaction for you to consider in user subroutine VUINTERACTION. Each proximity point is primarily associated with a slave node or a point along
+
+
+
+a slave edge; the proximity point also references a corresponding, locally nearest point on the master surface within the tracking thickness. A proximity point exists for each pairing of slave node and proximal master surface point. Therefore, more than one proximity point may reference the same node on the slave surface but different points on the master surface if multiple local minimum distances to the slave node exist on the master surface; this phenomenon commonly occurs near the corners of a master surface. No proximity points exist for a slave node that is separated from the master surface by more than the tracking thickness. A two-dimensional representation for contact between portions of shell surfaces is shown in Figure 1.2.21–1.
+
+
+
+