VUMAT

Component Ordering in Tensors

The component ordering depends on whether the tensor is symmetric or nonsymmetric.

Symmetric Tensors

For symmetric tensors such as the stress and strain tensors, there are ndir+nshr components, and the component order is given as a natural permutation of the indices of the tensor. The direct components are first and then the indirect components, beginning with the 12-component. For example, a stress tensor contains ndir direct stress components and nshr shear stress components, which are passed in as


Component	2D Case	3D Case
1	$σ_{11}$	$σ_{11}$
2	$σ_{22}$	$σ_{22}$
3	$σ_{33}$	$σ_{33}$
4	$σ_{12}$	$σ_{12}$
5		$σ_{23}$
6		$σ_{31}$

The shear strain components in user subroutine VUMAT are stored as tensor components and not as engineering components; this is different from user subroutine UMAT in Abaqus/Standard, which uses engineering components.

Nonsymmetric Tensors

For nonsymmetric tensors there are ndir+2*nshr components, and the component order is given as a natural permutation of the indices of the tensor. The direct components are first and then the indirect components, beginning with the 12-component. For example, the deformation gradient is passed as


Component	2D Case	3D Case
1	$F_{11}$	$F_{11}$
2	$F_{22}$	$F_{22}$
3	$F_{33}$	$F_{33}$
4	$F_{12}$	$F_{12}$
5	$F_{21}$	$F_{23}$
6		$F_{31}$
7		$F_{21}$
8		$F_{32}$
9		$F_{13}$

Initial Calculations and Checks

In the data check phase of the analysis Abaqus/Explicit calls user subroutine VUMAT with a set of fictitious strains and a totalTime and stepTime both equal to 0.0. This is done as a check on your constitutive relation and to calculate the equivalent initial material properties, based on which the transverse shear stiffness (if VUMAT is used to define the material response of shell elements and material transverse shear modulus is defined) and the initial elastic wave speeds are computed.

Defining Local Orientations

All stresses, strains, stretches, and state variables are in the orientation of the local material axes. These local material axes form a basis system in which stress and strain components are stored. This represents a corotational coordinate system in which the basis system rotates with the material. If a user-specified coordinate system (Orientations) is used, it defines the local material axes in the undeformed configuration.

Special Considerations for Various Element Types

In addition to defining the stresses and internal state variables, the use of user subroutine VUMAT requires special consideration for various element types.

Shell and Plane Stress Elements

In the case of shell or plane stress elements, NDIR=3 and NSHR=1; you must define the thickness strain increment, strainInc(*,3), inside VUMAT (typically by enforcing the plane stress constraint or the incompressibility constraint for fully incompressible materials). When the constitutive model employs the right stretch tensor $U$ (such as in hyperelastic formulations), you should also define the thickness components of the stretch and deformation gradient tensors, stretchNew (*, 3) and defgradNew (*,3). Optionally, you can define the internal energies. If you do not define them, Abaqus/Explicit cannot provide a meaningful energy balance.

Shell Elements

When VUMAT is used to define the material response of shell elements, Abaqus/Explicit cannot calculate a default value for the transverse shear stiffness of the element. Hence, you must define the material transverse shear modulus (see Defining the Elastic Transverse Shear Modulus) or the element's transverse shear stiffness (see Shell Section Behavior).

Beam Elements

For beam elements the stretch tensor and the deformation gradient tensor are not available. For beams in space you must define the thickness strains, strainInc(*,2) and strainInc(*,3). strainInc(*,4) is the shear strain associated with twist. Thickness stresses, stressNew(*,2) and stressNew(*,3), are assumed to be zero, and any values you assign are ignored.

Pipe Elements

For pipe elements the stretch tensor and the deformation gradient tensor are not available. The axial strain, strainInc(*,1), and the shear strain, strainInc(*,4), associated with twist are provided along with the hoop stress, stressNew(*,2). The hoop stress is predefined based on your pipe internal and external pressure load definitions (PE, PI, HPE, HPI, PENU, and PINU), and it should not be modified here. The thickness stress, stressNew(*,3), is assumed to be zero and any value you assign is ignored. You must define the axial stress, stressNew(*,1), and the shear stress, stressNew(*,4). You must also define hoop strain, strainInc(*,2), and the pipe thickness strain, strainInc(*,3).

Deformation Gradient

The polar decomposition of the deformation gradient is written as $F = R \cdot U = V \cdot R$ , where $U$ and $V$ are the right and left symmetric stretch tensors, respectively. The constitutive model is defined in a corotational coordinate system in which the basis system rotates with the material. All stress and strain tensor quantities are defined with respect to the corotational basis system. The right stretch tensor, $U$ , is used. The relative spin tensor $W - Ω$ represents the spin (the antisymmetric part of the velocity gradient) defined with respect to the corotational basis system.

Special Considerations for Hyperelasticity

Hyperelastic constitutive models in VUMAT should be defined in a corotational coordinate system in which the basis system rotates with the material. This is most effectively accomplished by formulating the hyperelastic constitutive model in terms of the stretch tensor, $U$ , instead of in terms of the deformation gradient, $F = R \cdot U$ . Using the deformation gradient can present some difficulties because the deformation gradient includes the rotation tensor and the resulting stresses would need to be rotated back to the corotational basis.

Objective Stress Rates

The Green-Naghdi stress rate is used when the mechanical behavior of the material is defined using user subroutine VUMAT. The stress rate obtained with user subroutine VUMAT might differ from that obtained with a built-in Abaqus material model. For example, most material models used with solid (continuum) elements in Abaqus/Explicit employ the Jaumann stress rate. This difference in the formulation will cause significant differences in the results only if finite rotation of a material point is accompanied by finite shear. For a discussion of the objective stress rates used in Abaqus, see Stress rates.

Defining Effective Modulus to Control Time Incrementation in Abaqus/Explicit

The stable time increment in Abaqus/Explicit is a function of the dilatational wave speed of the material and the characteristic length of the element. By default, the dilatational wave speed is determined automatically by Abaqus/Explicit using a numerical calculation of the effective hypoelastic bulk and shear moduli from the user material's constitutive response. In situations when the user material is highly nonlinear, this calculation might not yield a conservative value of the time increment, leading to unstable behavior. To avoid these situations, you can define the values of the effective bulk and shear moduli directly inside user subroutine VUMAT when the EFFMOD parameter is specified in the USER MATERIAL option. These values are then used by Abaqus/Explicit to compute the dilatational wave speed and the stable time increment. For more information, see Stability.

Material Point Deletion

Material points that satisfy a user-defined failure criterion can be deleted from the model (see User-Defined Mechanical Material Behavior). You must specify the state variable number controlling the element deletion flag when you allocate space for the solution-dependent state variables, as explained in User-Defined Mechanical Material Behavior. The deletion state variable should be set to a value of one or zero in VUMAT. A value of one indicates that the material point is active, while a value of zero indicates that Abaqus/Explicit should delete the material point from the model by setting the stresses to zero. The structure of the block of material points passed to user subroutine VUMAT remains unchanged during the analysis; deleted material points are not removed from the block. Abaqus/Explicit will pass zero stresses and strain increments for all deleted material points. Once a material point has been flagged as deleted, it cannot be reactivated.

User Subroutine Interface

      subroutine vumat(
C Read only (unmodifiable)variables -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, jInfoArray,
     2  stepTime, totalTime, dtArray, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
      parameter (i_info_AnnealFlag = 1, 
     *     i_info_Intpt    = 2, ! Integration station number
     *     i_info_layer  = 3, ! Layer number
     *     i_info_kspt   = 4, ! Section point number in current layer
     *     i_info_effModDefn = 5, ! =1 if Bulk/ShearMod need to be defined
     *     i_info_ElemNumStartLoc   = 6) ! Start loc of user element number
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(nblock), dtArray(2*(nblock)+1), strainInc(nblock,ndir+nshr),
     2  relSpinInc(nblock,nshr), tempOld(nblock), 
     3  stretchOld(nblock,ndir+nshr),
     4  defgradOld(nblock,ndir+nshr+nshr),
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(nblock),
     8  stretchNew(nblock,ndir+nshr),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock), jInfoArray(*)
C
      character*80 cmname
C
      pointer (ptrjElemNum, jElemNum)
      dimension jElemNum(nblock)
C
      lAnneal = jInfoArray(i_info_AnnealFlag) 
      iLayer = jInfoArray(i_info_layer)
      kspt   = jInfoArray(i_info_kspt)
      intPt  = jInfoArray(i_info_Intpt)
      iUpdateEffMod = jInfoArray(i_info_effModDefn)
      iElemNumStartLoc = jInfoArray(i_info_ElemNumStartLoc)
      ptrjElemNum = loc(jInfoArray(iElemNumStartLoc))

      do 100 km = 1,nblock
        user coding
  100 continue

      return
      end

Variables to Be Defined

stressNew (nblock, ndir+nshr): Stress tensor at each material point at the end of the increment.
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 About User Subroutines and Utilities for more information).

Variables That Can Be Updated

enerInternNew (nblock): Internal energy per unit mass at each material point at the end of the increment.
enerInelasNew (nblock): Dissipated inelastic energy per unit mass at each material point at the end of the increment.
dtArray(2 : nblock+1): Shear modulus at each material point for nblock elements in dtArray(2 : nblock+1). The shear modulus returned by the user subroutine is used to compute the time increment only when the EFFMOD parameter is specified in the USER MATERIAL option.
dtArray(nblock+2 : 2*(nblock)+1): Bulk modulus at each material point for nblock elements in dtArray(nblock+2 : 2*(nblock)+1). The bulk modulus returned by the subroutine is used to compute the time increment only when the EFFMOD parameter is specified in the USER MATERIAL option.

Variables Passed in for Information

nblock

Number of material points to be processed in this call to VUMAT.

ndir

Number of direct components in a symmetric tensor.

nshr

Number of indirect components in a symmetric tensor.

nstatev

Number of user-defined state variables that are associated with this material type (you define this as described in Allocating Space for Solution-Dependent State Variables).

nfieldv

Number of user-defined external field variables.

nprops

User-specified number of user-defined material properties.

jInfoArray


jInfoArray(i_info_AnnealFlag)	Flag indicating whether the routine is being called during an annealing process. jInfoArray(i_info_AnnealFlag)=0 indicates that the routine is being called during a normal mechanics increment. jInfoArray(i_info_AnnealFlag)=1 indicates that this is an annealing process and you should reinitialize the internal state variables, stateNew, if necessary. Abaqus/Explicit automatically sets the stresses, stretches, and state to a value of zero during the annealing process.
jInfoArray(i_info_Intpt)	Integration location number.
jInfoArray(i_info_layer)	Layer number.
jInfoArray(i_info_kspt)	Section point number in the current layer.
jInfoArray(i_info_effModDefn)	Flag indicating whether the bulk modulus and shear modulus must be updated. jInfoArray(i_info_effModDefn) = 1 indicates that you must update the values for the bulk and shear moduli.
jInfoArray(i_info_ElemNumStartLoc)	Start location in jInfoArray for user-defined element numbers of all elements in nblock.

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.

dtArray(1)

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.

coordMp(nblock,*)

Material point coordinates. It is the midplane material point for shell elements and the centroid for beam and pipe elements.

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 of 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.

props(nprops)

User-supplied material properties.

density(nblock)

Current density at the material points in the midstep configuration. This value might be inaccurate in problems where the volumetric strain increment is very small. If an accurate value of the density is required in such cases, the analysis should be run in double precision. This value of the density is not affected by mass scaling.

strainInc (nblock, ndir+nshr)

Strain increment tensor at each material point.

relSpinInc (nblock, nshr)

Incremental relative rotation vector at each material point defined in the corotational system. Defined as $Δ t (W - Ω)$ , where $W$ is the antisymmetric part of the velocity gradient, $L$ , and $Ω = \dot{R} \cdot R^{T}$ . Stored in 3D as $(32, 13, 21)$ and in 2D as $(21)$ .

tempOld(nblock)

Temperatures at each material point at the beginning of the increment.

stretchOld (nblock, ndir+nshr)

Stretch tensor, $U$ , at each material point at the beginning of the increment defined from the polar decomposition of the deformation gradient by $F = R \cdot U$ .

defgradOld (nblock,ndir+2*nshr)

Deformation gradient tensor at each material point at the beginning of the increment. Stored in 3D as ( $F_{11}$ , $F_{22}$ , $F_{33}$ , $F_{12}$ , $F_{23}$ , $F_{31}$ , $F_{21}$ , $F_{32}$ , $F_{13}$ ) and in 2D as ( $F_{11}$ , $F_{22}$ , $F_{33}$ , $F_{12}$ , $F_{21}$ ).

fieldOld (nblock, nfieldv)

Values of the user-defined field variables at each material point at the beginning of the increment.

stressOld (nblock, ndir+nshr)

Stress tensor at each material point at the beginning of the increment.

stateOld (nblock, nstatev)

State variables at each material point at the beginning of the increment.

enerInternOld (nblock)

Internal energy per unit mass at each material point at the beginning of the increment.

enerInelasOld (nblock)

Dissipated inelastic energy per unit mass at each material point at the beginning of the increment.

tempNew(nblock)

Temperatures at each material point at the end of the increment.

stretchNew (nblock, ndir+nshr)

Stretch tensor, $U$ , at each material point at the end of the increment defined from the polar decomposition of the deformation gradient by $F = R \cdot U$ .

defgradNew (nblock,ndir+2*nshr)

Deformation gradient tensor at each material point at the end of the increment. Stored in 3D as ( $F_{11}$ , $F_{22}$ , $F_{33}$ , $F_{12}$ , $F_{23}$ , $F_{31}$ , $F_{21}$ , $F_{32}$ , $F_{13}$ ) and in 2D as ( $F_{11}$ , $F_{22}$ , $F_{33}$ , $F_{12}$ , $F_{21}$ ).

fieldNew (nblock, nfieldv)

Values of the user-defined field variables at each material point at the end of the increment.

Example: Using More than One User-Defined Material Model

To use more than one user-defined material model, the variable cmname can be tested for different material names inside user subroutine VUMAT, as illustrated below:

if (cmname(1:4) .eq. 'MAT1') then
   call VUMAT_MAT1(argument_list)
else if (cmname(1:4) .eq. 'MAT2') then
   call VUMAT_MAT2(argument_list)
end if

VUMAT_MAT1 and VUMAT_MAT2 are the actual user material subroutines containing the constitutive material models for each material MAT1 and MAT2, respectively. Subroutine VUMAT merely acts as a directory here. The argument list can be the same as that used in subroutine VUMAT. The material names must be in uppercase characters since cmname is passed in as an uppercase character string.

Example: Elastic/Plastic Material with Kinematic Hardening

As a simple example of the coding of subroutine VUMAT, consider the generalized plane strain case for an elastic/plastic material with kinematic hardening. The basic assumptions and definitions of the model are as follows.

Let $σ$ be the current value of the stress, and define $S$ to be the deviatoric part of the stress. The center of the yield surface in deviatoric stress space is given by the tensor $α$ , which has initial values of zero. The stress difference, $ξ$ , is the stress measured from the center of the yield surface and is given by

ξ = S - α .

The von Mises yield surface is defined as

f (σ) = \frac{1}{2} ξ : ξ - \frac{1}{3} σ_{0}^{2},

where $σ_{0}$ is the uniaxial equivalent yield stress. The von Mises yield surface is a cylinder in deviatoric stress space with a radius of

R = \sqrt{\frac{2}{3}} σ_{0} .

For the kinematic hardening model, R is a constant. The normal to the Mises yield surface can be written as

Q = \sqrt{\frac{3}{2}} \frac{ξ}{σ_{0}} .

We decompose the strain rate into an elastic and plastic part using an additive decomposition:

\dot{ϵ} = {\dot{ϵ}}^{e l} + {\dot{ϵ}}^{p l} .

The plastic part of the strain rate is given by a normality condition

{\dot{ϵ}}^{p l} = \dot{γ} Q,

where the scalar multiplier $\dot{γ}$ must be determined. A scalar measure of equivalent plastic strain rate is defined by

{\dot{\bar{ϵ}}}^{p l} = \sqrt{\frac{2}{3} {\dot{ϵ}}^{p l} : {\dot{ϵ}}^{p l}} .

The stress rate is assumed to be purely due to the elastic part of the strain rate and is expressed in terms of Hooke's law by

\dot{σ} = λ trace ({\dot{ϵ}}^{e l}) I + 2 μ {\dot{ϵ}}^{e l},

where $λ$ and $2 μ$ are the Lamés constants for the material.

The evolution law for $α$ is given as

\dot{α} = \frac{2}{3} \dot{γ} H Q,

where H is the slope of the uniaxial yield stress versus plastic strain curve.

During active plastic loading the stress must remain on the yield surface, so that

\sqrt{Q : Q} = 1.

The equivalent plastic strain rate is related to $\dot{γ}$ by

{\dot{\bar{ϵ}}}^{p l} = \sqrt{\frac{2}{3}} \dot{γ} .

The kinematic hardening constitutive model is integrated in a rate form as follows. A trial elastic stress is computed as

σ_{n e w}^{t r i a l} = σ_{o l d} + λ trace (Δ ϵ) I + 2 μ Δ ϵ,

where the subscripts $o l d$ and $n e w$ refer to the beginning and end of the increment, respectively. If the trial stress does not exceed the yield stress, the new stress is set equal to the trial stress. If the yield stress is exceeded, plasticity occurs in the increment. We then write the incremental analogs of the rate equations as

σ_{n e w} = σ_{n e w}^{t r i a l} - 2 μ Δ ϵ^{p l} = σ_{n e w}^{t r i a l} - 2 μ Δ γ Q,

α_{n e w} = α_{o l d} + \frac{2}{3} H Δ γ Q,

{\bar{ϵ}}_{n e w}^{p l} = {\bar{ϵ}}_{o l d}^{p l} + \sqrt{\frac{2}{3}} Δ γ,

where

Δ γ = \dot{γ} Δ t .

From the definition of the normal to the yield surface at the end of the increment, $Q$ ,

α_{n e w} + \sqrt{\frac{2}{3}} σ_{0} Q = S_{n e w} .

This can be expanded using the incremental equations as

α_{o l d} + \frac{2}{3} H Δ γ Q + \sqrt{\frac{2}{3}} σ_{0} Q = S_{n e w}^{t r i a l} - Δ γ 2 μ Q .

Taking the tensor product of this equation with $Q$ , using the yield condition at the end of the increment, and solving for $Δ γ$ :

Δ γ = \frac{1}{2 μ (1 + H / 3 μ)} ({(ξ_{n e w}^{t r i a l} : ξ_{n e w}^{t r i a l})}^{1 / 2} - \sqrt{\frac{2}{3}} σ_{0}) .

The value for $Δ γ$ is used in the incremental equations to determine $σ_{n e w}$ , $α_{n e w}$ , and ${\bar{ϵ}}_{n e w}^{p l}$ .

This algorithm is often referred to as an elastic predictor, radial return algorithm because the correction to the trial stress under the active plastic loading condition returns the stress state to the yield surface along the direction defined by the vector from the center of the yield surface to the elastic trial stress. The subroutine would be coded as follows:

       subroutine vumat(
C Read only -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, jInfoArray,
     2  stepTime, totalTime, dtArray, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     3  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
     5  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
C J2 Mises Plasticity with kinematic hardening for plane 
C strain case.
C Elastic predictor, radial corrector algorithm.
C
C The state variables are stored as:
C      STATE(*,1) = back stress component 11
C      STATE(*,2) = back stress component 22
C      STATE(*,3) = back stress component 33
C      STATE(*,4) = back stress component 12
C      STATE(*,5) = equivalent plastic strain
C
C
C All arrays dimensioned by (*) are not used in this algorithm
      dimension props(nprops), density(nblock),
     1  coordMp(nblock,*),
     2  charLength(*), dtArray(*), strainInc(nblock,ndir+nshr),
     3  relSpinInc(*), tempOld(*),
     4  stretchOld(*), defgradOld(*),
     5  fieldOld(*),  stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(*),
     8  stretchNew(*), defgradNew(*), fieldNew(*),
     9  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock), jInfoArray(*)
C
      parameter (i_info_AnnealFlag = 1, 
     *     i_info_Intpt    = 2, ! Integration station number
     *     i_info_layer  = 3, ! Layer number
     *     i_info_kspt   = 4, ! Section point number in current layer
     *     i_info_effModDefn = 5, ! =1 if Bulk/ShearMod need to be defined
     *     i_info_ElemNumStartLoc   = 6) ! Start loc of user element number
C
      parameter( zero = 0., one = 1., two = 2., three = 3.,
     1  third = one/three, half = .5, twoThirds = two/three,
     2  threeHalfs = 1.5 )

C
      character*80 cmname
C
      pointer (ptrjElemNum, jElemNum)
      dimension jElemNum(nblock)
C
      lAnneal = jInfoArray(i_info_AnnealFlag) 
      iLayer = jInfoArray(i_info_layer)
      kspt   = jInfoArray(i_info_kspt)
      intPt  = jInfoArray(i_info_Intpt)
      iUpdateEffMod = jInfoArray(i_info_effModDefn)
      iElemNumStartLoc = jInfoArray(i_info_ElemNumStartLoc)
      ptrjElemNum = loc(jInfoArray(iElemNumStartLoc))
c      
      e     = props(1)
      xnu   = props(2)
      yield = props(3)
      hard  = props(4)
C
      twomu  = e / ( one + xnu )
      thremu = threeHalfs * twomu
      sixmu  = three * twomu
      alamda = twomu * ( e - twomu ) / ( sixmu - two * e )
      term   = one / ( twomu * ( one + hard/thremu ) )
      con1   = sqrt( twoThirds )
C
      do 100 i = 1,nblock
C
C Trial stress
        trace  = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
        sig1 = stressOld(i,1) + alamda*trace + twomu*strainInc(i,1)
        sig2 = stressOld(i,2) + alamda*trace + twomu*strainInc(i,2)
        sig3 = stressOld(i,3) + alamda*trace + twomu*strainInc(i,3)
        sig4 = stressOld(i,4)                + twomu*strainInc(i,4)
C
C Trial stress measured from the back stress
        s1 = sig1 - stateOld(i,1)
        s2 = sig2 - stateOld(i,2)
        s3 = sig3 - stateOld(i,3)
        s4 = sig4 - stateOld(i,4)
C
C Deviatoric part of trial stress measured from the back stress
        smean = third * ( s1 + s2 + s3 )
        ds1 = s1 - smean
        ds2 = s2 - smean
        ds3 = s3 - smean
C
C Magnitude of the deviatoric trial stress difference
        dsmag = sqrt( ds1**2 + ds2**2 + ds3**2 + 2.*s4**2 )
C
C Check for yield by determining the factor for plasticity,
C zero for elastic, one for yield
        radius = con1 * yield
        facyld = zero
        if(  dsmag - radius .ge. zero ) facyld = one
C
C Add a protective addition factor to prevent a divide by zero 
C when dsmag is zero. If dsmag is zero, we will not have exceeded
C the yield stress and facyld will be zero.
        dsmag  = dsmag + ( one - facyld )
C
C Calculated increment in gamma (this explicitly includes the 
C time step)
        diff   = dsmag - radius
        dgamma = facyld * term * diff
C
C Update equivalent plastic strain
        deqps  = con1 * dgamma
        stateNew(i,5) = stateOld(i,5) + deqps
C
C Divide dgamma by dsmag so that the deviatoric stresses are
C explicitly converted to tensors of unit magnitude in the
C following calculations
        dgamma = dgamma / dsmag
C
C Update back stress
        factor  = hard * dgamma * twoThirds
        stateNew(i,1) = stateOld(i,1) + factor * ds1
        stateNew(i,2) = stateOld(i,2) + factor * ds2
        stateNew(i,3) = stateOld(i,3) + factor * ds3
        stateNew(i,4) = stateOld(i,4) + factor *  s4
C
C Update the stress
        factor   = twomu * dgamma
        stressNew(i,1) = sig1 - factor * ds1
        stressNew(i,2) = sig2 - factor * ds2
        stressNew(i,3) = sig3 - factor * ds3
        stressNew(i,4) = sig4 - factor *  s4
C
C Update the specific internal energy -
        stressPower = half * (
     1    ( stressOld(i,1)+stressNew(i,1) )*strainInc(i,1)
     1    +     ( stressOld(i,2)+stressNew(i,2) )*strainInc(i,2)
     1    +     ( stressOld(i,3)+stressNew(i,3) )*strainInc(i,3)
     1    + two*( stressOld(i,4)+stressNew(i,4) )*strainInc(i,4) )
C
        enerInternNew(i) = enerInternOld(i)
     1    + stressPower / density(i)
C
C Update the dissipated inelastic specific energy -
        plasticWorkInc = dgamma * half * (
     1    ( stressOld(i,1)+stressNew(i,1) )*ds1
     1    +     ( stressOld(i,2)+stressNew(i,2) )*ds2
     1    +     ( stressOld(i,3)+stressNew(i,3) )*ds3
     1    + two*( stressOld(i,4)+stressNew(i,4) )*s4 )
        enerInelasNew(i) = enerInelasOld(i)
      1   + plasticWorkInc / density(i)
  100 continue
C
      return
      end