SUBROUTINE VUMAT(NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV,
1 NPROPS,LANNEAL,STEPTIME, TOTALTIME, DT, CMNAME, COORDMP,
2 CHARLENGTH,PROPS, DENSITY, STRAININC, RELSPININC,
3 TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
4 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
5 TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,
6 STRESSNEW, STATENEW, ENERINTERNNEW,ENERINELASNEW)
C -------------------VUMAT Interface Variable Description-----------------------
C NBLOCK——Number of material integration points processed in this VUMAT call
C NDIR——Number of diagonal tensor components for stress/strain; 3 for plane/3D problems
C NSHR——Number of off-diagonal shear tensor components; 1 for plane, 3 for 3D
C NSTATEV——Number of user-defined state dependent variables (SDVs)
C NFIELDV——Number of external user field variables
C NPROPS——Number of material constitutive parameters
C LANNEAL——Annealing flag to reinitialize internal state variables
C STEPTIME——Time value at the start of current increment
C TOTALTIME——Total accumulated analysis time
C DT——Time increment size of current step
C CMNAME——Material name character string for distinguishing multiple materials
C COORDMP——Spatial coordinates of material integration points
C CHARLENGTH——Characteristic element length at material point
C PROPS——Array storing user-defined material constants (same as UMAT material table)
C DENSITY——Material mass density defined in material input block
C STRAININC——Incremental strain tensor at each material point
C RELSPININC——Incremental relative rotation tensor under reference rotating coordinate system
C TEMPOLD——Temperature at the beginning of current increment
C STRETCHOLD——Left stretch tensor U at increment start; suffix OLD = start of increment, NEW = end of increment
C DEFGRADOLD——Deformation gradient tensor at increment start
C FIELDOLD——User external field variables at increment start
C STRESSOLD——Cauchy stress tensor at increment start
C STATEOLD——State dependent variables at increment start
C ENERINTERNOLD——Internal energy density at increment start
C ENERINELASOLD——Elastic strain energy density at increment start
C TEMPNEW——Temperature at the end of current increment
C STRETCHNEW——Left stretch tensor U at increment end
C DEFGRADNEW——Deformation gradient tensor at increment end
C FIELDNEW——User external field variables at increment end
C STRESSNEW——Updated Cauchy stress tensor at increment end (output variable)
C STATENEW——Updated state dependent variables at increment end (output variable)
C ENERINTERNNEW——Internal energy density at increment end
C ENERINELASNEW——Elastic strain energy density at increment end
C -------------------------------------------------------------------------------
INCLUDE 'VABA_PARAM.INC'
C VARIABLE DECLARATION
DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK,*),
1 CHARLENGTH(NBLOCK), 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),
9 DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
1 FIELDNEW(NBLOCK,NFIELDV),
2 STRESSNEW(NBLOCK,NDIR+NSHR), STATENEW(NBLOCK,NSTATEV),
3 ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)
CHARACTER*80 CMNAME
!---------------------------------------------------------------------
! LOCAL VARIABLE DECLARATION
!---------------------------------------------------------------------
INTEGER ZERO, ONE, TWO, NTENS,
& NDI, NSTATV,
& I, NSHRUMAT, NPROPSUMAT
DOUBLE PRECISION DTIME
DOUBLE PRECISION STRESS(NDIR+NSHR), STATEV(NSTATEV),
& STRAN(NDIR+NSHR), DSTRAN(NDIR+NSHR),
& TIME(2),
& DFGRD0(3,3), DFGRD1(3,3)
DOUBLE PRECISION PROPSUMAT(NPROPS)
PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0)
!*************************************************************************
C Initialize specific SDVs at zero total time (initial analysis step)
IF (TOTALTIME .EQ. ZERO) THEN
DO KM = 1, NBLOCK
STATENEW(KM,2) = ZERO
STATENEW(KM,7) = STATEOLD(KM,7)
ENDDO
ENDIF
!*************************************************************************
C Loop over all material integration points in current block
DO 100 KM = 1,NBLOCK
C Map old stress and strain increment to local UMAT-style array
DO I = 1, NDIR
STRESS(I) = STRESSOLD(KM,I)
DSTRAN(I) = STRAININC(KM,I)
ENDDO
STRESS(4) = STRESSOLD(KM,4)
DSTRAN(4) = TWO * STRAININC(KM,4)
IF (NSHR .GT. 1) THEN
STRESS(5) = STRESSOLD(KM,6)
DSTRAN(5) = TWO * STRAININC(KM,6)
STRESS(6) = STRESSOLD(KM,5)
DSTRAN(6) = TWO * STRAININC(KM,5)
ENDIF
C Assign deformation gradient at increment end DFGRD1 (3x3 full tensor)
DFGRD1(1,1) = DEFGRADNEW(KM,1)
DFGRD1(2,2) = DEFGRADNEW(KM,2)
DFGRD1(3,3) = DEFGRADNEW(KM,3)
DFGRD1(1,2) = DEFGRADNEW(KM,4)
DFGRD1(2,3) = DEFGRADNEW(KM,5)
DFGRD1(3,1) = DEFGRADNEW(KM,6)
DFGRD1(2,1) = DEFGRADNEW(KM,7)
DFGRD1(3,2) = DEFGRADNEW(KM,8)
DFGRD1(1,3) = DEFGRADNEW(KM,9)
C Assign deformation gradient at increment start DFGRD0 (3x3 full tensor)
DFGRD0(1,1) = DEFGRADOLD(KM,1)
DFGRD0(2,2) = DEFGRADOLD(KM,2)
DFGRD0(3,3) = DEFGRADOLD(KM,3)
DFGRD0(1,2) = DEFGRADOLD(KM,4)
DFGRD0(2,3) = DEFGRADOLD(KM,5)
DFGRD0(3,1) = DEFGRADOLD(KM,6)
DFGRD0(2,1) = DEFGRADOLD(KM,7)
DFGRD0(3,2) = DEFGRADOLD(KM,8)
DFGRD0(1,3) = DEFGRADOLD(KM,9)
C Pass global VUMAT variables to local UMAT-style variables for crystal plasticity subroutine
NSTATV = NSTATEV
STATEV(:) = STATEOLD(KM,:)
DTIME = DT
NDI = NDIR
NSHRUMAT = NSHR
NTENS = NDIR + NSHR
PROPSUMAT = PROPS
NPROPSUMAT = NPROPS
TIME(1)=TOTALTIME
TIME(2)=STEPTIME
C Call Huang Yonggang single crystal plasticity constitutive subroutine to update stress and state variables
CALL CRYSTALPLASTICITY(STRESS,STATEV,STRAN,DSTRAN,
& TIME,DTIME,CMNAME,NDI,NSHRUMAT,NTENS,NSTATV,PROPSUMAT,
& NPROPSUMAT,DFGRD0,DFGRD1)
C Write updated stress back to VUMAT output array STRESSNEW
DO I = 1, NDIR
STRESSNEW(KM,I) = STRESS(I)
ENDDO
STRESSNEW(KM,4) = STRESS(4)
IF( NSHRUMAT .GT. 1 ) THEN
STRESSNEW(KM,5) = STRESS(6)
STRESSNEW(KM,6) = STRESS(5)
ENDIF
C Write updated state dependent variables back to VUMAT output array STATENEW
STATENEW(KM,:) = STATEV(:)
100 CONTINUE ! End of material point loop
RETURN
END
C Huang Yonggang Single Crystal Plasticity Constitutive Subroutine
C Trimmed version retaining only stress and state variable update logic
SUBROUTINE CRYSTALPLASTICITY(STRESS,STATEV,STRAN,DSTRAN,
2 TIME,DTIME,CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,
3 DFGRD0,DFGRD1)
C----- Single precision compilation note for Cray machines:
C (1) Delete statement "IMPLICIT*8 (A-H,O-Z)";
C (2) Change "REAL*8 FUNCTION" to "FUNCTION";
C (3) Replace double precision intrinsic DSIGN with SIGN.
C
C----- Internal Auxiliary Subroutines:
C
C ROTATION -- Construct crystal orientation rotation matrix;
C computes direction cosines of cubic crystal [100], [010], [001]
C axes under global coordinate system at initial state
C
C SLIPSYS -- Generate independent slip systems, unit slip direction vectors
C and unit slip plane normal vectors for cubic crystals at initial state
C
C GSLPINIT -- Assign initial slip system critical shear strength values
C
C STRAINRATE -- Evaluate slip shear strain rate based on resolved shear stress
C and current slip system hardening strength following power-law viscoplasticity
C
C LATENTHARDEN -- Assemble self-hardening and latent hardening interaction matrix
C
C ITERATION -- Construct Jacobian arrays for Newton-Rhapson implicit iteration
C
C LUDCMP -- Perform LU matrix decomposition for linear system solving
C
C LUBKSB -- Solve linear system via precomputed LU decomposition
C
C
C----- Internal Auxiliary Function:
C
C F -- Slip system shear strain rate function (power-law viscoplasticity)
C
C----- Subroutine Input/Output Variables:
C
C STRESS -- Cauchy stress tensor (INPUT & OUTPUT)
C Finite deformation framework adopts true Cauchy stress
C STATEV -- State dependent solution variables (INPUT & OUTPUT)
C
C----- Passed auxiliary variables for reference:
C
C STRAN -- Integral logarithmic strain tensor for finite deformation
C Equivalent to time integral of symmetric velocity gradient
C DSTRAN -- Incremental strain tensor
C CMNAME -- Material name string defined in *MATERIAL keyword block
C NDI -- Count of direct normal stress tensor components
C NSHR -- Count of engineering shear stress tensor components
C NTENS -- Total tensor component count = NDI + NSHR
C NSTATV -- Total number of state dependent variables defined via *DEPVAR
C PROPS -- User material constants defined under *USER MATERIAL keyword
C NPROPS -- Total number of user material constants
C
C----- Constitutive Theory Overview:
C This subroutine implements finite deformation single crystal plasticity for ABAQUS.
C Crystal plastic slip follows Schmid's resolved shear stress criterion.
C Total strain increment decomposes additively into elastic lattice stretch strain and plastic slip strain.
C Elastic strain increment corresponds to pure lattice stretching; plastic strain is the superposition
C of shear slip over all activated slip systems.
C Slip shear strain increment is a power-law function of resolved shear stress normalized by slip system strength.
C Slip system hardening strength increment couples with accumulated shear strain via self/latent hardening interaction.
C
C----- Time Integration Scheme:
C Implicit backward integration algorithm proposed by Peirce, Shih & Needleman (1984) is adopted.
C Optional nested Newton-Rhapson iteration is available to converge stress and internal state variables per increment.
C
C----- Crystal System Restriction:
C Core implementation for single cubic crystals (FCC/BCC). Extensions to HCP, tetragonal, orthotropic
C lattices only require modifications to ROTATION and SLIPSYS subroutines to incorporate lattice aspect ratios.
C
C----- Critical User Setup Requirements:
C
C (1) Minimum required state variable count NSTATV:
C NSTATV >= 10 * NSLPTL + 5
C NSLPTL = total independent slip systems across all slip families
C Slip systems (s,-m), (-s,m), (-s,-m) are treated as dependent duplicates of (s,m) and excluded
C Cubic slip family examples: {110}<111> contains 12 independent slip systems
C If additional constitutive state parameters are required (e.g. Zarka model), extend to:
C NSTATV >= NPARMT + 10 * NSLPTL + 5
C
C (2) Tangent stiffness matrix asymmetry:
C Latent hardening introduces non-symmetric consistent tangent stiffness.
C Must declare keyword "UNSYMM" under *USER MATERIAL in ABAQUS input deck.
C
PARAMETER (ND=150)
C----- ND defines maximum array dimension for slip system storage
C Default value 150 supports up to three cubic slip families fully activated.
C Reduce ND to NSLPTL if fewer slip families are used (e.g. ND=12 for single {110}<111> family).
C
include 'aba_param.inc'
C
CHARACTER*8 CMNAME
EXTERNAL F
DIMENSION STRESS(NTENS),STATEV(NSTATV),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),
3 PROPS(NPROPS),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),
2 SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),
3 SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),
4 DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),
5 FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),
6 H(ND,ND), DDGDDE(ND,6),
7 DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),
8 DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),
9 WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)
DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),
2 GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),
3 DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),
4 DSPNRO(3,ND), DSPDRO(3,ND),
5 DHDGDG(ND,ND)
DOUBLE PRECISION LV(3,3),
4 LVT(3,3),WV(3,3),ROTA(3,3),ROTAT(3,3),DETADG(3,3),
5 DGINV(3,3),TERMIDW(3,3),TERMIPW(3,3),TERMIDWINV(3,3)
C----- NSLIP -- Independent slip system count per slip family
C----- SLPDIR -- Unit slip direction vectors in initial crystal local coordinate system
C----- SLPNOR -- Unit slip plane normal vectors in initial crystal local coordinate system
C----- SLPDEF -- Slip deformation tensor (Schmid factor matrix, Voigt 6-component format)
C SLPDEF(1,i) = SLPDIR(1,i)*SLPNOR(1,i)
C SLPDEF(2,i) = SLPDIR(2,i)*SLPNOR(2,i)
C SLPDEF(3,i) = SLPDIR(3,i)*SLPNOR(3,i)
C SLPDEF(4,i) = SLPDIR(1,i)*SLPNOR(2,i)+SLPDIR(2,i)*SLPNOR(1,i)
C SLPDEF(5,i) = SLPDIR(1,i)*SLPNOR(3,i)+SLPDIR(3,i)*SLPNOR(1,i)
C SLPDEF(6,i) = SLPDIR(2,i)*SLPNOR(3,i)+SLPDIR(3,i)*SLPNOR(2,i)
C Index i denotes ith independent slip system
C----- SLPSPN -- Slip spin tensor components (only required for finite rotation formulation)
C SLPSPN(1,i) = 0.5*(SLPDIR(1,i)*SLPNOR(2,i)-SLPDIR(2,i)*SLPNOR(1,i))
C SLPSPN(2,i) = 0.5*(SLPDIR(3,i)*SLPNOR(1,i)-SLPDIR(1,i)*SLPNOR(3,i))
C SLPSPN(3,i) = 0.5*(SLPDIR(2,i)*SLPNOR(3,i)-SLPDIR(3,i)*SLPNOR(2,i))
C----- DSPDIR -- Incremental update of slip direction vectors under finite rotation
C----- DSPNOR -- Incremental update of slip plane normal vectors under finite rotation
C
C----- DLOCAL -- Anisotropic elastic stiffness matrix defined in crystal local coordinate system
C----- D -- Anisotropic elastic stiffness matrix rotated to global Cartesian coordinate system
C----- ROTD -- Voigt rotation transformation matrix mapping DLOCAL to global D
C
C----- ROTATE -- Crystal orientation rotation matrix; stores direction cosines of crystal [100], [010], [001] axes
C relative to global X/Y/Z axes at initial state
C
C----- FSLIP -- Current shear strain rate magnitude for each slip system
C----- DFDXSP -- Derivative dF/dX where X = resolved shear stress / slip system hardening strength
C
C----- DDEMSD -- Double dot product of elastic stiffness tensor with Schmid tensor, plus
C spin-stress coupling term exclusively for finite rotation kinematics
C
C----- H -- Slip system hardening interaction matrix
C H(i,i) = Self-hardening modulus of ith slip system
C H(i,j) = Latent hardening modulus on system i induced by slip on system j (i≠j)
C
C----- DDGDDE -- Derivative of slip shear strain increment with respect to macroscopic strain increment
C
C----- DSTRES -- Jaumann corotational stress increment tensor co-rotated with material spin
C----- DELATS -- Lattice elastic stretch strain increment (macro strain minus plastic slip strain)
C DELATS(1-3) = Normal elastic strain increments
C DELATS(4-6) = Engineering elastic shear strain increments
C----- DSPIN -- Material element spin increment tensor components
C DSPIN(1) = 12 component of spin tensor
C DSPIN(2) = 31 component of spin tensor
C DSPIN(3) = 23 component of spin tensor
C
C----- DVGRAD -- Incremental velocity gradient tensor (velocity gradient × time increment)
C
C----- DGAMMA -- Incremental plastic shear strain magnitude for each slip system
C----- DTAUSP -- Incremental resolved shear stress on each slip system
C----- DGSLIP -- Incremental hardening strength increase for each slip system
C
C----- Iteration temporary storage arrays:
C FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,
C DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO,
C DHDGDG
C
C----- STATEV State Variable Storage Layout Definition:
C NSLPTL = total independent slip systems across all slip families
C
C STATEV(1 : NSLPTL) : Current slip system hardening strength g_α
C STATEV(NSLPTL+1 : 2*NSLPTL) : Accumulated shear strain γ_α on each slip system
C STATEV(2*NSLPTL+1 : 3*NSLPTL) : Current resolved shear stress τ_α on each slip system
C
C STATEV(3*NSLPTL+1 : 6*NSLPTL) : Current updated slip plane normal vectors m_α (3 components per system)
C STATEV(6*NSLPTL+1 : 9*NSLPTL) : Current updated slip direction vectors s_α (3 components per system)
C
C STATEV(9*NSLPTL+1 : 10*NSLPTL) : Cumulative absolute shear strain |γ_α| for each individual slip system
C
C STATEV(10*NSLPTL+1) : Global total cumulative absolute shear strain sum(|γ_α|) over all slip systems
C
C STATEV(10*NSLPTL+2 : NSTATV-4) : User-extended auxiliary constitutive state parameters (optional)
C
C STATEV(NSTATV-3) : Slip system count of first slip family
C STATEV(NSTATV-2) : Slip system count of second slip family
C STATEV(NSTATV-1) : Slip system count of third slip family
C STATEV(NSTATV) : Total independent slip system count NSLPTL
C
C----- PROPS Material Constant Array Layout Definition:
C
C PROPS(1) - PROPS(21) -- Anisotropic elastic stiffness constants
C
C Isotropic elastic model: PROPS(i)=0 for i>2
C PROPS(1) = Young's Modulus E
C PROPS(2) = Poisson's Ratio ν
C
C Cubic elastic model: PROPS(i)=0 for i>3
C PROPS(1) = C11
C PROPS(2) = C12
C PROPS(3) = C44
C
C Orthotropic elastic model: PROPS(1)-PROPS(9) match ABAQUS orthotropic elastic input order
C D1111, D1122, D2222, D1133, D2233, D3333, D1212, D1313, D2323
C
C General fully anisotropic elastic model: PROPS(1)-PROPS(21) full 4th-order stiffness tensor Voigt components
C
C
C PROPS(25) - PROPS(56) -- Slip family definition parameters for cubic crystal
C
C PROPS(25) -- Number of distinct slip families (maximum 3, input as floating point value e.g. 3.0)
C
C PROPS(33) - PROPS(35) -- Miller indices of reference slip plane normal for slip family 1 e.g. (1,1,0)
C PROPS(36) - PROPS(38) -- Miller indices of reference slip direction for slip family 1 e.g. [1,-1,1]
C
C PROPS(41) - PROPS(43) -- Reference slip plane normal Miller indices for slip family 2
C PROPS(44) - PROPS(46) -- Reference slip direction Miller indices for slip family 2
C
C PROPS(49) - PROPS(51) -- Reference slip plane normal Miller indices for slip family 3
C PROPS(52) - PROPS(54) -- Reference slip direction Miller indices for slip family 3
C
C
C PROPS(57) - PROPS(72) -- Crystal orientation definition parameters
C Two non-parallel vectors required to construct orientation rotation matrix
C
C PROPS(57) - PROPS(59) -- Local crystal coordinate vector 1 Miller indices e.g. [1,1,0]
C PROPS(60) - PROPS(62) -- Global Cartesian coordinate vector 1 components (non-unit vector allowed)
C
C PROPS(65) - PROPS(67) -- Local crystal coordinate vector 2 Miller indices
C PROPS(68) - PROPS(70) -- Global Cartesian coordinate vector 2 components
C
C
C PROPS(73) - PROPS(96) -- Viscoplastic power-law slip rate parameters for each slip family
C
C PROPS(73) - PROPS(80) -- Power-law parameters for slip family 1
C PROPS(81) - PROPS(88) -- Power-law parameters for slip family 2
C PROPS(89) - PROPS(96) -- Power-law parameters for slip family 3
C
C
C PROPS(97) - PROPS(144) -- Self & latent hardening law parameters for each slip family
C
C PROPS(97) - PROPS(104)-- Self-hardening parameters slip family 1
C PROPS(105)- PROPS(112)-- Latent hardening interaction parameters slip family 1
C
C PROPS(113)- PROPS(120)-- Self-hardening parameters slip family 2
C PROPS(121)- PROPS(128)-- Latent hardening interaction parameters slip family 2
C
C PROPS(129)- PROPS(136)-- Self-hardening parameters slip family 3
C PROPS(137)- PROPS(144)-- Latent hardening interaction parameters slip family 3
C
C
C PROPS(145)- PROPS(152)-- Time integration and finite deformation switch parameters
C
C PROPS(145) -- Implicit integration weighting factor θ (0 ≤ θ ≤ 1)
C θ=0: Explicit forward Euler integration
C θ=0.5: Recommended midpoint integration
C θ=1.0: Fully implicit backward Euler integration
C
C PROPS(146) -- Finite geometry activation flag NLGEOM
C 0.0 = Small deformation infinitesimal theory
C Non-zero = Finite rotation & finite strain kinematics (requires *NLGEOM step keyword)
C
C
C PROPS(153)- PROPS(160)-- Newton iteration control parameters
C
C PROPS(153) -- Iteration enable flag ITRATN
C 0.0 = No nested iteration, single step solve
C Non-zero = Activate Newton-Rhapson iteration
C
C PROPS(154) -- Maximum allowed iteration count ITRMAX
C
C PROPS(155) -- Convergence tolerance GAMERR for slip shear strain residual
C
C ***************************************
C Supplementary routine to compute incremental rotation tensor DROT
C ***************************************
C Compute deformation gradient increment DFGRAD1 - DFGRD0
DO I = 1,3
DO J = 1,3
DETADG(I,J) = DFGRD1(I,J)-DFGRD0(I,J)
END DO
END DO
C Copy end-of-step deformation gradient for matrix inversion
DO I = 1,3
DO J = 1,3
DGINV(I,J) = DFGRD1(I,J)
END DO
END DO
C Call subroutine to calculate inverse of deformation gradient
CALL GET_INV_DET(DGINV,GARB,0)
C Compute incremental velocity gradient LV = ΔF · F⁻¹
DO I = 1,3
DO J = 1,3
LV(I,J) = 0.0D0
DO K = 1,3
LV(I,J) = LV(I,J)+DETADG(I,K)*DGINV(K,J)
END DO
END DO
END DO
C Extract skew-symmetric spin tensor component WV = Ω·dt
DO I =1,3
WV(I,I) = 0.0D0
END DO
WV(1,2) =0.5D0*( LV(1,2)-LV(2,1))
WV(2,1) = -WV(1,2)
WV(1,3) =0.5D0*( LV(1,3)-LV(3,1))
WV(3,1) = -WV(1,3)
WV(2,3) =0.5D0*( LV(2,3)-LV(3,2))
WV(3,2) = -WV(2,3)
C Zero temporary matrix storage
CALL CLEAR(TERMIDW,9)
CALL CLEAR(TERMIPW,9)
C Construct I ± 0.5Ωdt matrix for midpoint rotation integration
DO I = 1,3
TERMIDW(I,I) = 1.0D0
TERMIPW(I,I) = 1.0D0
DO J = 1,3
TERMIDW(I,J) = TERMIDW(I,J)-0.5D0*WV(I,J)
TERMIPW(I,J) = TERMIPW(I,J)+0.5D0*WV(I,J)
END DO
END DO
C Copy matrix for inversion
DO I=1,3
DO J = 1,3
TERMIDWINV(I,J) = TERMIDW(I,J)
END DO
END DO
C Compute inverse of (I + 0.5Ωdt)
CALL GET_INV_DET(TERMIDWINV,GARB,0)
C Calculate incremental rotation tensor DROT = (I+0.5Ωdt)⁻¹ · (I-0.5Ωdt)
DO I = 1,3
DO J = 1,3
ROTA(I,J) = 0.0D0
DO K =1,3
ROTA(I,J) = ROTA(I,J)+TERMIDWINV(I,K)*TERMIPW(K,J)
END DO
END DO
END DO
C Assign final incremental rotation matrix to DROT
DO I = 1,3
DO J = 1,3
DROT(I,J) = ROTA(I,J)
END DO
END DO
C----- Construct anisotropic elastic stiffness matrix DLOCAL in crystal local coordinate system
DO J=1,6
DO I=1,6
DLOCAL(I,J)=0.
END DO
END DO
C Check if fully anisotropic elastic constants are provided
CHECK=0.
DO J=10,21
CHECK=CHECK+ABS(PROPS(J))
END DO
IF (CHECK.EQ.0.) THEN
C Check orthotropic elastic input
DO J=4,9
CHECK=CHECK+ABS(PROPS(J))
END DO
IF (CHECK.EQ.0.) THEN
C Isotropic or cubic elastic model
IF (PROPS(3).EQ.0.) THEN
C----- Isotropic elastic constitutive model
GSHEAR=PROPS(1)/2./(1.+PROPS(2))
E11=2.*GSHEAR*(1.-PROPS(2))/(1.-2.*PROPS(2))
E12=2.*GSHEAR*PROPS(2)/(1.-2.*PROPS(2))
DO J=1,3
DLOCAL(J,J)=E11
DO I=1,3
IF (I.NE.J) DLOCAL(I,J)=E12
END DO
DLOCAL(J+3,J+3)=GSHEAR
END DO
ELSE
C----- Cubic elastic stiffness matrix C11, C12, C44
DO J=1,3
DLOCAL(J,J)=PROPS(1)
DO I=1,3
IF (I.NE.J) DLOCAL(I,J)=PROPS(2)
END DO
DLOCAL(J+3,J+3)=PROPS(3)
END DO
END IF
ELSE
C----- Orthotropic elastic stiffness matrix
DLOCAL(1,1)=PROPS(1)
DLOCAL(1,2)=PROPS(2)
DLOCAL(2,1)=PROPS(2)
DLOCAL(2,2)=PROPS(3)
DLOCAL(1,3)=PROPS(4)
DLOCAL(3,1)=PROPS(4)
DLOCAL(2,3)=PROPS(5)
DLOCAL(3,2)=PROPS(5)
DLOCAL(3,3)=PROPS(6)
DLOCAL(4,4)=PROPS(7)
DLOCAL(5,5)=PROPS(8)
DLOCAL(6,6)=PROPS(9)
END IF
ELSE
C----- General fully anisotropic elastic stiffness matrix (Voigt symmetric storage)
ID=0
DO J=1,6
DO I=1,J
ID=ID+1
DLOCAL(I,J)=PROPS(ID)
DLOCAL(J,I)=DLOCAL(I,J)
END DO
END DO
END IF
C----- Call orientation subroutine to build crystal rotation matrix ROTATE
CALL ROTATION (PROPS(57), ROTATE)
C----- Construct Voigt rotation transformation matrix ROTD for 6-component stiffness rotation
DO J=1,3
J1=1+J/3
J2=2+J/2
DO I=1,3
I1=1+I/3
I2=2+I/2
ROTD(I,J)=ROTATE(I,J)**2
ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTATE(I,J2)
ROTD(I+3,J)=ROTATE(I1,J)*ROTATE(I2,J)
ROTD(I+3,J+3)=ROTATE(I1,J1)*ROTATE(I2,J2)+
2 ROTATE(I1,J2)*ROTATE(I2,J1)
END DO
END DO
C----- Rotate local crystal stiffness DLOCAL to global Cartesian stiffness D
C D = ROTD · DLOCAL · ROTDᵀ
DO J=1,6
DO I=1,6
D(I,J)=0.
END DO
END DO
DO J=1,6
DO I=1,J
DO K=1,6
DO L=1,6
D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L)
END DO
END DO
D(J,I)=D(I,J)
END DO
END DO
C----- Read number of independent slip families NSET
NSET=NINT(PROPS(25))
IF (NSET.LT.1) THEN
WRITE (6,*) '***ERROR - zero slip families defined'
STOP
ELSE IF (NSET.GT.3) THEN
WRITE (6,*)
2 '***ERROR - maximum supported slip family count is 3'
STOP
END IF
C----- Implicit integration weighting factor θ
THETA=PROPS(145)
C----- Finite geometry switch NLGEOM
IF (PROPS(146).EQ.0.) THEN
NLGEOM=0
ELSE
NLGEOM=1
END IF
C----- Newton iteration enable flag ITRATN
IF (PROPS(153).EQ.0.) THEN
ITRATN=0
ELSE
ITRATN=1
END IF
ITRMAX=NINT(PROPS(154))
GAMERR=PROPS(155)
C Initialize iteration residual storage
NITRTN=-1
DO I=1,NTENS
DSOLD(I)=0.
END DO
DO J=1,ND
DGAMOD(J)=0.
DTAUOD(J)=0.
DGSPOD(J)=0.
DO I=1,3
DSPNRO(I,J)=0.
DSPDRO(I,J)=0.
END DO
END DO
C----- Compute material spin increment DSPIN from incremental rotation matrix DROT (finite deformation only)
IF (NLGEOM.NE.0) THEN
DO J=1,3
DO I=1,3
TERM(I,J)=DROT(J,I)
TRM0(I,J)=DROT(J,I)
END DO
TERM(J,J)=TERM(J,J)+1.D0
TRM0(J,J)=TRM0(J,J)-1.D0
END DO
CALL LUDCMP (TERM, 3, 3, ITRM, DDCMP)
DO J=1,3
CALL LUBKSB (TERM, 3, 3, ITRM, TRM0(1,J))
END DO
DSPIN(1)=TRM0(2,1)-TRM0(1,2)
DSPIN(2)=TRM0(1,3)-TRM0(3,1)
DSPIN(3)=TRM0(3,2)-TRM0(2,3)
END IF
C----- Volumetric strain increment trace
DEV=0.D0
DO I=1,NDI
DEV=DEV+DSTRAN(I)
END DO
C----- Newton-Rhapson iteration loop entry label
1000 CONTINUE
NITRTN=NITRTN+1
C----- Initialize slip system data at first analysis increment (TOTALTIME=0)
IF (STATEV(1).EQ.0.) THEN
NSLPTL=0
DO I=1,NSET
ISPNOR(1)=NINT(PROPS(25+8*I))
ISPNOR(2)=NINT(PROPS(26+8*I))
ISPNOR(3)=NINT(PROPS(27+8*I))
ISPDIR(1)=NINT(PROPS(28+8*I))
ISPDIR(2)=NINT(PROPS(29+8*I))
ISPDIR(3)=NINT(PROPS(30+8*I))
C Generate all independent slip systems for current slip family
CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),
2 SLPNOR(1,NSLPTL+1), ROTATE)
NSLPTL=NSLPTL+NSLIP(I)
END DO
C Check array dimension ND is sufficient for total slip systems
IF (ND.LT.NSLPTL) THEN
WRITE (6,*)
2 '***ERROR - Parameter ND smaller than total independent slip systems NSLPTL'
STOP
END IF
C----- Assemble Schmid slip deformation tensor SLPDEF for all slip systems
DO J=1,NSLPTL
SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)
SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)
SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)
SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)
SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)
SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J)
END DO
C----- Write slip system metadata to state variables
STATEV(NSTATV)=FLOAT(NSLPTL)
DO I=1,NSET
STATEV(NSTATV-4+I)=FLOAT(NSLIP(I))
END DO
C Store initial slip plane normals and slip directions to SDVs
IDNOR=3*NSLPTL
IDDIR=6*NSLPTL
DO J=1,NSLPTL
DO I=1,3
IDNOR=IDNOR+1
STATEV(IDNOR)=SLPNOR(I,J)
IDDIR=IDDIR+1
STATEV(IDDIR)=SLPDIR(I,J)
END DO
END DO
C Assign initial critical slip strength g0 via GSLPINIT
CALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))
C Initialize accumulated shear strain and cumulative absolute slip
DO I=1,NSLPTL
STATEV(NSLPTL+I)=0.
STATEV(9*NSLPTL+I)=0.
END DO
STATEV(10*NSLPTL+1)=0.
C Calculate initial resolved shear stress τ_α = σ : S_α
DO I=1,NSLPTL
TERM1=0.
DO J=1,NTENS
IF (J.LE.NDI) THEN
TERM1=TERM1+SLPDEF(J,I)*STRESS(J)
ELSE
TERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J)
END IF
END DO
STATEV(2*NSLPTL+I)=TERM1
END DO
ELSE
C----- Post-initialization: Read existing slip system data from state variables
NSLPTL=NINT(STATEV(NSTATV))
DO I=1,NSET
NSLIP(I)=NINT(STATEV(NSTATV-4+I))
END DO
C Recover current slip plane normals and slip directions from SDVs
IDNOR=3*NSLPTL
IDDIR=6*NSLPTL
DO J=1,NSLPTL
DO I=1,3
IDNOR=IDNOR+1
SLPNOR(I,J)=STATEV(IDNOR)
IDDIR=IDDIR+1
SLPDIR(I,J)=STATEV(IDDIR)
END DO
END DO
C Rebuild Schmid tensor SLPDEF from updated slip vectors
DO J=1,NSLPTL
SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)
SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)
SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)
SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)
SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)
SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J)
END DO
END IF
C----- Construct slip spin tensor SLPSPN (finite rotation kinematics only)
IF (NLGEOM.NE.0) THEN
DO J=1,NSLPTL
SLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-
2 SLPDIR(2,J)*SLPNOR(1,J))
SLPSPN(2,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-
2 SLPDIR(1,J)*SLPNOR(3,J))
SLPSPN(3,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-
2 SLPDIR(3,J)*SLPNOR(2,J))
END DO
END IF
C----- Compute DDEMSD coupling matrix D:S_α + σ×W_α (finite rotation extra term)
DO J=1,NSLPTL
DO I=1,6
DDEMSD(I,J)=0.
DO K=1,6
DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)
END DO
END DO
END IF
IF (NLGEOM.NE.0) THEN
DO J=1,NSLPTL
DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*STRESS(1)
DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*STRESS(1)
IF (NDI.GT.1) THEN
DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)
DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(3,J)*STRESS(2)
END IF
IF (NDI.GT.2) THEN
DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(2,J)*STRESS(3)
DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(3,J)*STRESS(3)
END IF
IF (NSHR.GE.1) THEN
DDEMSD(1,J)=DDEMSD(1,J)+SLPSPN(1,J)*STRESS(NDI+1)
DDEMSD(2,J)=DDEMSD(2,J)-SLPSPN(1,J)*STRESS(NDI+1)
DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(3,J)*STRESS(NDI+1)
DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(2,J)*STRESS(NDI+1)
END IF
IF (NSHR.GE.2) THEN
DDEMSD(1,J)=DDEMSD(1,J)-SLPSPN(2,J)*STRESS(NDI+2)
DDEMSD(3,J)=DDEMSD(3,J)+SLPSPN(2,J)*STRESS(NDI+2)
DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(3,J)*STRESS(NDI+2)
DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(1,J)*STRESS(NDI+2)
END IF
IF (NSHR.EQ.3) THEN
DDEMSD(2,J)=DDEMSD(2,J)+SLPSPN(3,J)*STRESS(NDI+3)
DDEMSD(3,J)=DDEMSD(3,J)-SLPSPN(3,J)*STRESS(NDI+3)
DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(2,J)*STRESS(NDI+3)
DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(1,J)*STRESS(NDI+3)
END IF
END DO
END IF
C----- Evaluate slip shear strain rate FSLIP and derivative dF/dX via STRAINRATE
ID=1
DO I=1,NSET
IF (I.GT.1) ID=ID+NSLIP(I-1)
CALL STRAINRATE (STATEV(NSLPTL+ID), STATEV(2*NSLPTL+ID),
2 STATEV(ID), NSLIP(I), FSLIP(ID), DFDXSP(ID),
3 PROPS(65+8*I))
END DO
C----- Assemble self/latent hardening interaction matrix H
CALL LATENTHARDEN (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1),
2 STATEV(1), STATEV(9*NSLPTL+1),
3 STATEV(10*NSLPTL+1), NSLIP, NSLPTL,
4 NSET, H(1,1), PROPS(97), ND)
C----- Build Jacobian matrix for slip shear strain increment solve
TERM1=THETA*DTIME
DO I=1,NSLPTL
TAUSLP=STATEV(2*NSLPTL+I)
GSLIP=STATEV(I)
X=TAUSLP/GSLIP
TERM2=TERM1*DFDXSP(I)/GSLIP
TERM3=TERM1*X*DFDXSP(I)/GSLIP
DO J=1,NSLPTL
TERM4=0.
DO K=1,6
TERM4=TERM4+DDEMSD(K,I)*SLPDEF(K,J)
END DO
WORKST(I,J)=TERM2*TERM4+H(I,J)*TERM3*DSIGN(1.D0,FSLIP(J))
IF (NITRTN.GT.0) WORKST(I,J)=WORKST(I,J)+TERM3*DHDGDG(I,J)
END DO
WORKST(I,I)=WORKST(I,I)+1.
END DO
C LU decomposition of slip system Jacobian matrix
CALL LUDCMP (WORKST, NSLPTL, ND, INDX, DDCMP)
C----- Increment of shear strain in a slip system: DGAMMA
TERM1=THETA*DTIME
DO I=1,NSLPTL
IF (NITRTN.EQ.0) THEN
TAUSLP=STATEV(2*NSLPTL+I)
GSLIP=STATEV(I)
X=TAUSLP/GSLIP
TERM2=TERM1*DFDXSP(I)/GSLIP
DGAMMA(I)=0.
DO J=1,NDI
DGAMMA(I)=DGAMMA(I)+DDEMSD(J,I)*DSTRAN(J)
END DO
IF (NSHR.GT.0) THEN
DO J=1,NSHR
DGAMMA(I)=DGAMMA(I)+DDEMSD(J+3,I)*DSTRAN(J+NDI)
END DO
END IF
DGAMMA(I)=DGAMMA(I)*TERM2+FSLIP(I)*DTIME
ELSE
DGAMMA(I)=TERM1*(FSLIP(I)-FSLIP1(I))+FSLIP1(I)*DTIME
2 -DGAMOD(I)
END IF
END DO
CALL LUBKSB (WORKST, NSLPTL, ND, INDX, DGAMMA)
DO I=1,NSLPTL
DGAMMA(I)=DGAMMA(I)+DGAMOD(I)
END DO
C----- Update the shear strain in a slip system: STATEV(NSLPTL+1) -
C STATEV(2*NSLPTL)
C
DO I=1,NSLPTL
STATEV(NSLPTL+I)=STATEV(NSLPTL+I)+DGAMMA(I)-DGAMOD(I) ! 更新切应变
END DO
C----- Increment of current strength in a slip system: DGSLIP
DO I=1,NSLPTL
DGSLIP(I)=0.
DO J=1,NSLPTL
DGSLIP(I)=DGSLIP(I)+H(I,J)*ABS(DGAMMA(J)) !求滑移系强度增量
END DO
END DO
C----- Update the current strength in a slip system: STATEV(1) -
C STATEV(NSLPTL)
C
DO I=1,NSLPTL
STATEV(I)=STATEV(I)+DGSLIP(I)-DGSPOD(I) !并更新强度
END DO
C----- Increment of strain associated with lattice stretching: DELATS
DO J=1,6
DELATS(J)=0.
END DO
DO J=1,3
IF (J.LE.NDI) DELATS(J)=DSTRAN(J)
DO I=1,NSLPTL
DELATS(J)=DELATS(J)-SLPDEF(J,I)*DGAMMA(I)
END DO
END DO
DO J=1,3
IF (J.LE.NSHR) DELATS(J+3)=DSTRAN(J+NDI)
DO I=1,NSLPTL
DELATS(J+3)=DELATS(J+3)-SLPDEF(J+3,I)*DGAMMA(I)
END DO
END DO
C----- Increment of deformation gradient associated with lattice
C stretching in the current state, i.e. the velocity gradient
C (associated with lattice stretching) times the increment of time:
C DVGRAD (only needed for finite rotation)
C
IF (NLGEOM.NE.0) THEN
DO J=1,3
DO I=1,3
IF (I.EQ.J) THEN
DVGRAD(I,J)=DELATS(I)
ELSE
DVGRAD(I,J)=DELATS(I+J+1)
END IF
END DO
END DO
DO J=1,3
DO I=1,J
IF (J.GT.I) THEN
IJ2=I+J-2
IF (MOD(IJ2,2).EQ.1) THEN
TERM1=1.
ELSE
TERM1=-1.
END IF
DVGRAD(I,J)=DVGRAD(I,J)+TERM1*DSPIN(IJ2)
DVGRAD(J,I)=DVGRAD(J,I)-TERM1*DSPIN(IJ2)
DO K=1,NSLPTL
DVGRAD(I,J)=DVGRAD(I,J)-TERM1*DGAMMA(K)*
2 SLPSPN(IJ2,K)
DVGRAD(J,I)=DVGRAD(J,I)+TERM1*DGAMMA(K)*
2 SLPSPN(IJ2,K)
END DO
END IF
END DO
END DO
END IF
C----- Increment of resolved shear stress in a slip system: DTAUSP
DO I=1,NSLPTL
DTAUSP(I)=0.
DO J=1,6
DTAUSP(I)=DTAUSP(I)+DDEMSD(J,I)*DELATS(J)
END DO
END DO
C----- Update the resolved shear stress in a slip system:
C STATEV(2*NSLPTL+1) - STATEV(3*NSLPTL)
C
DO I=1,NSLPTL
STATEV(2*NSLPTL+I)=STATEV(2*NSLPTL+I)+DTAUSP(I)-DTAUOD(I)
END DO
C----- Increment of stress: DSTRES
IF (NLGEOM.EQ.0) THEN
DO I=1,NTENS
DSTRES(I)=0.
END DO
ELSE
DO I=1,NTENS
DSTRES(I)=-STRESS(I)*DEV
END DO
END IF
DO I=1,NDI
DO J=1,NDI
DSTRES(I)=DSTRES(I)+D(I,J)*DSTRAN(J)
END DO
IF (NSHR.GT.0) THEN
DO J=1,NSHR
DSTRES(I)=DSTRES(I)+D(I,J+3)*DSTRAN(J+NDI)
END DO
END IF
DO J=1,NSLPTL
DSTRES(I)=DSTRES(I)-DDEMSD(I,J)*DGAMMA(J)
END DO
END DO
IF (NSHR.GT.0) THEN
DO I=1,NSHR
DO J=1,NDI
DSTRES(I+NDI)=DSTRES(I+NDI)+D(I+3,J)*DSTRAN(J)
END DO
DO J=1,NSHR
DSTRES(I+NDI)=DSTRES(I+NDI)+D(I+3,J+3)*DSTRAN(J+NDI)
END DO
DO J=1,NSLPTL
DSTRES(I+NDI)=DSTRES(I+NDI)-DDEMSD(I+3,J)*DGAMMA(J)
END DO
END DO
END IF
C----- Update the stress: STRESS
DO I=1,NTENS
STRESS(I)=STRESS(I)+DSTRES(I)-DSOLD(I)
END DO
C----- Increment of normal to a slip plane and a slip direction (only
C needed for finite rotation)
C
IF (NLGEOM.NE.0) THEN
DO J=1,NSLPTL
DO I=1,3
DSPNOR(I,J)=0.
DSPDIR(I,J)=0.
DO K=1,3
DSPNOR(I,J)=DSPNOR(I,J)-SLPNOR(K,J)*DVGRAD(K,I)
DSPDIR(I,J)=DSPDIR(I,J)+SLPDIR(K,J)*DVGRAD(I,K)
END DO
END DO
END DO
C----- Update the normal to a slip plane and a slip direction (only
C needed for finite rotation)
C
IDNOR=3*NSLPTL
IDDIR=6*NSLPTL
DO J=1,NSLPTL
DO I=1,3
IDNOR=IDNOR+1
STATEV(IDNOR)=STATEV(IDNOR)+DSPNOR(I,J)-DSPNRO(I,J)
IDDIR=IDDIR+1
STATEV(IDDIR)=STATEV(IDDIR)+DSPDIR(I,J)-DSPDRO(I,J)
END DO
END DO
END IF
C----- Iteration ?
IF (ITRATN.NE.0) THEN
C----- Save solutions (without iteration):
C Shear strain-rate in a slip system FSLIP1
C Current strength in a slip system GSLP1
C Shear strain in a slip system GAMMA1
C Resolved shear stress in a slip system TAUSP1
C Normal to a slip plane SPNOR1
C Slip direction SPDIR1
C Stress STRES1
C Jacobian matrix DDSDE1
C
IF (NITRTN.EQ.0) THEN
IDNOR=3*NSLPTL
IDDIR=6*NSLPTL
DO J=1,NSLPTL
FSLIP1(J)=FSLIP(J)
GSLP1(J)=STATEV(J)
GAMMA1(J)=STATEV(NSLPTL+J)
TAUSP1(J)=STATEV(2*NSLPTL+J)
DO I=1,3
IDNOR=IDNOR+1
SPNOR1(I,J)=STATEV(IDNOR)
IDDIR=IDDIR+1
SPDIR1(I,J)=STATEV(IDDIR)
END DO
END DO
END IF
C----- Increments of stress DSOLD, and solution dependent state
C variables DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO (for the next
C iteration)
C
DO I=1,NTENS
DSOLD(I)=DSTRES(I)
END DO
DO J=1,NSLPTL
DGAMOD(J)=DGAMMA(J)
DTAUOD(J)=DTAUSP(J)
DGSPOD(J)=DGSLIP(J)
DO I=1,3
DSPNRO(I,J)=DSPNOR(I,J)
DSPDRO(I,J)=DSPDIR(I,J)
END DO
END DO
C----- Check if the iteration solution converges
IDBACK=0
ID=0
DO I=1,NSET
DO J=1,NSLIP(I)
ID=ID+1
X=STATEV(2*NSLPTL+ID)/STATEV(ID)
RESIDU=THETA*DTIME*F(X,PROPS(65+8*I))+DTIME*(1.0-THETA)*
2 FSLIP1(ID)-DGAMMA(ID)
IF (ABS(RESIDU).GT.GAMERR) IDBACK=1
END DO
END DO
IF (IDBACK.NE.0.AND.NITRTN.LT.ITRMAX) THEN
C----- Iteration: arrays for iteration
CFIXA
CALL ITERATION (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1),
2 STATEV(1), STATEV(9*NSLPTL+1),
3 STATEV(10*NSLPTL+1), NSLPTL,
4 NSET, NSLIP, ND, PROPS(97), DGAMOD,
5 DHDGDG)
CFIXB
GO TO 1000
ELSE IF (NITRTN.GE.ITRMAX) THEN
C----- Solution not converge within maximum number of iteration (the
C solution without iteration will be used)
IDNOR=3*NSLPTL
IDDIR=6*NSLPTL
DO J=1,NSLPTL
STATEV(J)=GSLP1(J)
STATEV(NSLPTL+J)=GAMMA1(J)
STATEV(2*NSLPTL+J)=TAUSP1(J)
DO I=1,3
IDNOR=IDNOR+1
STATEV(IDNOR)=SPNOR1(I,J)
IDDIR=IDDIR+1
STATEV(IDDIR)=SPDIR1(I,J)
END DO
END DO
END IF
END IF
C----- Total cumulative shear strains on all slip systems (sum of the
C absolute values of shear strains in all slip systems)
CFIX-- Total cumulative shear strains on each slip system (sum of the
CFIX absolute values of shear strains in each individual slip system)
C
DO I=1,NSLPTL
CFIXA
STATEV(10*NSLPTL+1)=STATEV(10*NSLPTL+1)+ABS(DGAMMA(I))
STATEV(9*NSLPTL+I)=STATEV(9*NSLPTL+I)+ABS(DGAMMA(I))
CFIXB
END DO
RETURN
END
C---------------------以下是VUMAT调用的函数--------------------------------
C
SUBROUTINE ROTATION (PROP, ROTATE)
C----- This subroutine calculates the rotation matrix, i.e. the
C direction cosines of cubic crystal [100], [010] and [001]
C directions in global system
C----- The rotation matrix is stored in the array ROTATE.
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION PROP(16), ROTATE(3,3), TERM1(3,3), TERM2(3,3), INDX(3)
C----- Subroutines:
C
C CROSS -- cross product of two vectors
C
C LUDCMP -- LU decomposition
C
C LUBKSB -- linear equation solver based on LU decomposition
C method (must call LUDCMP first)
C----- PROP -- constants characterizing the crystal orientation
C (INPUT)
C
C PROP(1) - PROP(3) -- direction of the first vector in
C local cubic crystal system
C PROP(4) - PROP(6) -- direction of the first vector in
C global system
C
C PROP(9) - PROP(11)-- direction of the second vector in
C local cubic crystal system
C PROP(12)- PROP(14)-- direction of the second vector in
C global system
C
C----- ROTATE -- rotation matrix (OUTPUT):
C
C ROTATE(i,1) -- direction cosines of direction [1 0 0] in
C local cubic crystal system
C ROTATE(i,2) -- direction cosines of direction [0 1 0] in
C local cubic crystal system
C ROTATE(i,3) -- direction cosines of direction [0 0 1] in
C local cubic crystal system
C----- local matrix: TERM1
CALL CROSS (PROP(1), PROP(9), TERM1, ANGLE1)
C----- LU decomposition of TERM1
CALL LUDCMP (TERM1, 3, 3, INDX, DCMP)
C----- inverse matrix of TERM1: TERM2
DO J=1,3
DO I=1,3
IF (I.EQ.J) THEN
TERM2(I,J)=1.
ELSE
TERM2(I,J)=0.
END IF
END DO
END DO
DO J=1,3
CALL LUBKSB (TERM1, 3, 3, INDX, TERM2(1,J))
END DO
C----- global matrix: TERM1
CALL CROSS (PROP(4), PROP(12), TERM1, ANGLE2)
C----- Check: the angle between first and second vector in local and
C global systems must be the same. The relative difference must be
C less than 0.1%.
C
IF (ABS(ANGLE1/ANGLE2-1.).GT.0.001) THEN
WRITE (6,*)
2 '***ERROR - angles between two vectors are not the same'
STOP
END IF
C----- rotation matrix: ROTATE
DO J=1,3
DO I=1,3
ROTATE(I,J)=0.
DO K=1,3
ROTATE(I,J)=ROTATE(I,J)+TERM1(I,K)*TERM2(K,J)
END DO
END DO
END DO
RETURN
END
C-----------------------------------
SUBROUTINE CROSS (A, B, C, ANGLE)
C----- (1) normalize vectors A and B to unit vectors
C (2) store A, B and A*B (cross product) in C
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(3), B(3), C(3,3)
SUM1=SQRT(A(1)**2+A(2)**2+A(3)**2)
SUM2=SQRT(B(1)**2+B(2)**2+B(3)**2)
IF (SUM1.EQ.0.) THEN
WRITE (6,*) '***ERROR - first vector is zero'
STOP
ELSE
DO I=1,3
C(I,1)=A(I)/SUM1
END DO
END IF
IF (SUM2.EQ.0.) THEN
WRITE (6,*) '***ERROR - second vector is zero'
STOP
ELSE
DO I=1,3
C(I,2)=B(I)/SUM2
END DO
END IF
ANGLE=0.
DO I=1,3
ANGLE=ANGLE+C(I,1)*C(I,2)
END DO
ANGLE=ACOS(ANGLE)
C(1,3)=C(2,1)*C(3,2)-C(3,1)*C(2,2)
C(2,3)=C(3,1)*C(1,2)-C(1,1)*C(3,2)
C(3,3)=C(1,1)*C(2,2)-C(2,1)*C(1,2)
SUM3=SQRT(C(1,3)**2+C(2,3)**2+C(3,3)**2)
IF (SUM3.LT.1.E-8) THEN
WRITE (6,*)
2 '***ERROR - first and second vectors are parallel'
STOP
END IF
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE SLIPSYS (ISPDIR, ISPNOR, NSLIP, SLPDIR, SLPNOR,
2 ROTATE)
C----- This subroutine generates all slip systems in the same set for
C a CUBIC crystal. For other crystals (e.g., HCP, Tetragonal,
C Orthotropic, ...), it has to be modified to include the effect of
C crystal aspect ratio.
C----- Denote s as a slip direction and m as normal to a slip plane.
C In a cubic crystal, (s,-m), (-s,m) and (-s,-m) are NOT considered
C independent of (s,m).
C----- Subroutines: LINE1 and LINE
C----- Variables:
C
C ISPDIR -- a typical slip direction in this set of slip systems
C (integer) (INPUT)
C ISPNOR -- a typical normal to slip plane in this set of slip
C systems (integer) (INPUT)
C NSLIP -- number of independent slip systems in this set
C (OUTPUT)
C SLPDIR -- unit vectors of all slip directions (OUTPUT)
C SLPNOR -- unit normals to all slip planes (OUTPUT)
C ROTATE -- rotation matrix (INPUT)
C ROTATE(i,1) -- direction cosines of [100] in global system
C ROTATE(i,2) -- direction cosines of [010] in global system
C ROTATE(i,3) -- direction cosines of [001] in global system
C
C NSPDIR -- number of all possible slip directions in this set
C NSPNOR -- number of all possible slip planes in this set
C IWKDIR -- all possible slip directions (integer)
C IWKNOR -- all possible slip planes (integer)
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION ISPDIR(3), ISPNOR(3), SLPDIR(3,50), SLPNOR(3,50),
* ROTATE(3,3), IWKDIR(3,24), IWKNOR(3,24), TERM(3)
NSLIP=0
NSPDIR=0
NSPNOR=0
C----- Generating all possible slip directions in this set
C
C Denote the slip direction by [lmn]. I1 is the minimum of the
C absolute value of l, m and n, I3 is the maximum and I2 is the
C mode, e.g. (1 -3 2), I1=1, I2=2 and I3=3. I1<=I2<=I3.
I1=MIN(IABS(ISPDIR(1)),IABS(ISPDIR(2)),IABS(ISPDIR(3)))
I3=MAX(IABS(ISPDIR(1)),IABS(ISPDIR(2)),IABS(ISPDIR(3)))
I2=IABS(ISPDIR(1))+IABS(ISPDIR(2))+IABS(ISPDIR(3))-I1-I3
RMODIR=SQRT(FLOAT(I1*I1+I2*I2+I3*I3))
C I1=I2=I3=0
IF (I3.EQ.0) THEN
WRITE (6,*) '***ERROR - slip direction is [000]'
STOP
C I1=I2=0, I3>0 --- [001] type
ELSE IF (I2.EQ.0) THEN
NSPDIR=3
DO J=1,3
DO I=1,3
IWKDIR(I,J)=0
IF (I.EQ.J) IWKDIR(I,J)=I3
END DO
END DO
C I1=0, I3>=I2>0
ELSE IF (I1.EQ.0) THEN
C I1=0, I3=I2>0 --- [011] type
IF (I2.EQ.I3) THEN
NSPDIR=6
DO J=1,6
DO I=1,3
IWKDIR(I,J)=I2
IF (I.EQ.J.OR.J-I.EQ.3) IWKDIR(I,J)=0
IWKDIR(1,6)=-I2
IWKDIR(2,4)=-I2
IWKDIR(3,5)=-I2
END DO
END DO
C I1=0, I3>I2>0 --- [012] type
ELSE
NSPDIR=12
CALL LINE1 (I2, I3, IWKDIR(1,1), 1)
CALL LINE1 (I3, I2, IWKDIR(1,3), 1)
CALL LINE1 (I2, I3, IWKDIR(1,5), 2)
CALL LINE1 (I3, I2, IWKDIR(1,7), 2)
CALL LINE1 (I2, I3, IWKDIR(1,9), 3)
CALL LINE1 (I3, I2, IWKDIR(1,11), 3)
END IF
C I1=I2=I3>0 --- [111] type
ELSE IF (I1.EQ.I3) THEN
NSPDIR=4
CALL LINE (I1, I1, I1, IWKDIR)
C I3>I2=I1>0 --- [112] type
ELSE IF (I1.EQ.I2) THEN
NSPDIR=12
CALL LINE (I1, I1, I3, IWKDIR(1,1))
CALL LINE (I1, I3, I1, IWKDIR(1,5))
CALL LINE (I3, I1, I1, IWKDIR(1,9))
C I3=I2>I1>0 --- [122] type
ELSE IF (I2.EQ.I3) THEN
NSPDIR=12
CALL LINE (I1, I2, I2, IWKDIR(1,1))
CALL LINE (I2, I1, I2, IWKDIR(1,5))
CALL LINE (I2, I2, I1, IWKDIR(1,9))
C I3>I2>I1>0 --- [123] type
ELSE
NSPDIR=24
CALL LINE (I1, I2, I3, IWKDIR(1,1))
CALL LINE (I3, I1, I2, IWKDIR(1,5))
CALL LINE (I2, I3, I1, IWKDIR(1,9))
CALL LINE (I1, I3, I2, IWKDIR(1,13))
CALL LINE (I2, I1, I3, IWKDIR(1,17))
CALL LINE (I3, I2, I1, IWKDIR(1,21))
END IF
C----- Generating all possible slip planes in this set
C
C Denote the normal to slip plane by (pqr). J1 is the minimum of
C the absolute value of p, q and r, J3 is the maximum and J2 is the
C mode, e.g. (1 -2 1), J1=1, J2=1 and J3=2. J1<=J2<=J3.
J1=MIN(IABS(ISPNOR(1)),IABS(ISPNOR(2)),IABS(ISPNOR(3)))
J3=MAX(IABS(ISPNOR(1)),IABS(ISPNOR(2)),IABS(ISPNOR(3)))
J2=IABS(ISPNOR(1))+IABS(ISPNOR(2))+IABS(ISPNOR(3))-J1-J3
RMONOR=SQRT(FLOAT(J1*J1+J2*J2+J3*J3))
IF (J3.EQ.0) THEN
WRITE (6,*) '***ERROR - slip plane is [000]'
STOP
C (001) type
ELSE IF (J2.EQ.0) THEN
NSPNOR=3
DO J=1,3
DO I=1,3
IWKNOR(I,J)=0
IF (I.EQ.J) IWKNOR(I,J)=J3
END DO
END DO
ELSE IF (J1.EQ.0) THEN
C (011) type
IF (J2.EQ.J3) THEN
NSPNOR=6
DO J=1,6
DO I=1,3
IWKNOR(I,J)=J2
IF (I.EQ.J.OR.J-I.EQ.3) IWKNOR(I,J)=0
IWKNOR(1,6)=-J2
IWKNOR(2,4)=-J2
IWKNOR(3,5)=-J2
END DO
END DO
C (012) type
ELSE
NSPNOR=12
CALL LINE1 (J2, J3, IWKNOR(1,1), 1)
CALL LINE1 (J3, J2, IWKNOR(1,3), 1)
CALL LINE1 (J2, J3, IWKNOR(1,5), 2)
CALL LINE1 (J3, J2, IWKNOR(1,7), 2)
CALL LINE1 (J2, J3, IWKNOR(1,9), 3)
CALL LINE1 (J3, J2, IWKNOR(1,11), 3)
END IF
C (111) type
ELSE IF (J1.EQ.J3) THEN
NSPNOR=4
CALL LINE (J1, J1, J1, IWKNOR)
C (112) type
ELSE IF (J1.EQ.J2) THEN
NSPNOR=12
CALL LINE (J1, J1, J3, IWKNOR(1,1))
CALL LINE (J1, J3, J1, IWKNOR(1,5))
CALL LINE (J3, J1, J1, IWKNOR(1,9))
C (122) type
ELSE IF (J2.EQ.J3) THEN
NSPNOR=12
CALL LINE (J1, J2, J2, IWKNOR(1,1))
CALL LINE (J2, J1, J2, IWKNOR(1,5))
CALL LINE (J2, J2, J1, IWKNOR(1,9))
C (123) type
ELSE
NSPNOR=24
CALL LINE (J1, J2, J3, IWKNOR(1,1))
CALL LINE (J3, J1, J2, IWKNOR(1,5))
CALL LINE (J2, J3, J1, IWKNOR(1,9))
CALL LINE (J1, J3, J2, IWKNOR(1,13))
CALL LINE (J2, J1, J3, IWKNOR(1,17))
CALL LINE (J3, J2, J1, IWKNOR(1,21))
END IF
C----- Generating all slip systems in this set
C
C----- Unit vectors in slip directions: SLPDIR, and unit normals to
C slip planes: SLPNOR in local cubic crystal system
C
WRITE (6,*) ' '
WRITE (6,*) ' # Slip plane Slip direction'
DO J=1,NSPNOR
DO I=1,NSPDIR
IDOT=0
DO K=1,3
IDOT=IDOT+IWKDIR(K,I)*IWKNOR(K,J)
END DO
IF (IDOT.EQ.0) THEN
NSLIP=NSLIP+1
DO K=1,3
SLPDIR(K,NSLIP)=IWKDIR(K,I)/RMODIR
SLPNOR(K,NSLIP)=IWKNOR(K,J)/RMONOR
END DO
WRITE (6,10) NSLIP,
2 (IWKNOR(K,J),K=1,3), (IWKDIR(K,I),K=1,3)
END IF
END DO
END DO
10 FORMAT(1X,I2,9X,'(',3(1X,I2),1X,')',10X,'[',3(1X,I2),1X,']')
WRITE (6,*) 'Number of slip systems in this set = ',NSLIP
WRITE (6,*) ' '
IF (NSLIP.EQ.0) THEN
WRITE (6,*)
* 'There is no slip direction normal to the slip planes!'
STOP
ELSE
C----- Unit vectors in slip directions: SLPDIR, and unit normals to
C slip planes: SLPNOR in global system
C
DO J=1,NSLIP
DO I=1,3
TERM(I)=0.
DO K=1,3
TERM(I)=TERM(I)+ROTATE(I,K)*SLPDIR(K,J)
END DO
END DO
DO I=1,3
SLPDIR(I,J)=TERM(I)
END DO
DO I=1,3
TERM(I)=0.
DO K=1,3
TERM(I)=TERM(I)+ROTATE(I,K)*SLPNOR(K,J)
END DO
END DO
DO I=1,3
SLPNOR(I,J)=TERM(I)
END DO
END DO
END IF
RETURN
END
C----------------------------------
SUBROUTINE LINE (I1, I2, I3, IARRAY)
C----- Generating all possible slip directions <lmn> (or slip planes
C {lmn}) for a cubic crystal, where l,m,n are not zeros.
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION IARRAY(3,4)
DO J=1,4
IARRAY(1,J)=I1
IARRAY(2,J)=I2
IARRAY(3,J)=I3
END DO
DO I=1,3
DO J=1,4
IF (J.EQ.I+1) IARRAY(I,J)=-IARRAY(I,J)
END DO
END DO
RETURN
END
C-----------------------------------
SUBROUTINE LINE1 (J1, J2, IARRAY, ID)
C----- Generating all possible slip directions <0mn> (or slip planes
C {0mn}) for a cubic crystal, where m,n are not zeros and m does
C not equal n.
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION IARRAY(3,2)
IARRAY(ID,1)=0
IARRAY(ID,2)=0
ID1=ID+1
IF (ID1.GT.3) ID1=ID1-3
IARRAY(ID1,1)=J1
IARRAY(ID1,2)=J1
ID2=ID+2
IF (ID2.GT.3) ID2=ID2-3
IARRAY(ID2,1)=J2
IARRAY(ID2,2)=-J2
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE GSLPINIT (GSLIP0, NSLIP, NSLPTL, NSET, PROP)
C----- This subroutine calculates the initial value of current
C strength for each slip system in a rate-dependent single crystal.
C Two sets of initial values, proposed by Asaro, Pierce et al, and
C by Bassani, respectively, are used here. Both sets assume that
C the initial values for all slip systems are the same (initially
C isotropic).
C----- These initial values are assumed the same for all slip systems
C in each set, though they could be different from set to set, e.g.
C <110>{111} and <110>{100}.
C----- Users who want to use their own initial values may change the
C function subprogram GSLP0. The parameters characterizing these
C initial values are passed into GSLP0 through array PROP.
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL GSLP0
DIMENSION GSLIP0(NSLPTL), NSLIP(NSET), PROP(16,NSET)
C----- Function subprograms:
C
C GSLP0 -- User-supplied function subprogram given the initial
C value of current strength at initial state
C----- Variables:
C
C GSLIP0 -- initial value of current strength (OUTPUT)
C
C NSLIP -- number of slip systems in each set (INPUT)
C NSLPTL -- total number of slip systems in all the sets (INPUT)
C NSET -- number of sets of slip systems (INPUT)
C
C PROP -- material constants characterizing the initial value of
C current strength (INPUT)
C
C For Asaro, Pierce et al's law
C PROP(1,i) -- initial hardening modulus H0 in the ith
C set of slip systems
C PROP(2,i) -- saturation stress TAUs in the ith set of
C slip systems
C PROP(3,i) -- initial critical resolved shear stress
C TAU0 in the ith set of slip systems
C
C For Bassani's law
C PROP(1,i) -- initial hardening modulus H0 in the ith
C set of slip systems
C PROP(2,i) -- stage I stress TAUI in the ith set of
C slip systems (or the breakthrough stress
C where large plastic flow initiates)
C PROP(3,i) -- initial critical resolved shear stress
C TAU0 in the ith set of slip systems
C
ID=0
DO I=1,NSET
ISET=I
DO J=1,NSLIP(I)
ID=ID+1
GSLIP0(ID)=GSLP0(NSLPTL,NSET,NSLIP,PROP(1,I),ID,ISET)
END DO
END DO
RETURN
END
C----------------------------------
C----- Use single precision on cray
C
REAL*8 FUNCTION GSLP0(NSLPTL,NSET,NSLIP,PROP,ISLIP,ISET)
C----- User-supplied function subprogram given the initial value of
C current strength at initial state
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION NSLIP(NSET), PROP(16)
GSLP0=PROP(3)
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE STRAINRATE (GAMMAR, TAUSLP, GSLIP, NSLIP, FSLIP,
2 DFDXSP, PROP)
C----- This subroutine calculates the shear strain-rate in each slip
C system for a rate-dependent single crystal. The POWER LAW
C relation between shear strain-rate and resolved shear stress
C proposed by Hutchinson, Pan and Rice, is used here.
C----- The power law exponents are assumed the same for all slip
C systems in each set, though they could be different from set to
C set, e.g. <110>{111} and <110>{100}. The strain-rate coefficient
C in front of the power law form are also assumed the same for all
C slip systems in each set.
C----- Users who want to use their own constitutive relation may
C change the function subprograms F and its derivative DFDX,
C where F is the strain hardening law, dGAMMA/dt = F(X),
C X=TAUSLP/GSLIP. The parameters characterizing F are passed into
C F and DFDX through array PROP.
C----- Function subprograms:
C
C F -- User-supplied function subprogram which gives shear
C strain-rate for each slip system based on current
C values of resolved shear stress and current strength
C
C DFDX -- User-supplied function subprogram dF/dX, where x is the
C ratio of resolved shear stress over current strength
C----- Variables:
C
C GAMMAR -- shear strain in each slip system at the start of time
C step (INPUT)
C TAUSLP -- resolved shear stress in each slip system (INPUT)
C GSLIP -- current strength (INPUT)
C NSLIP -- number of slip systems in this set (INPUT)
C
C FSLIP -- current value of F for each slip system (OUTPUT)
C DFDXSP -- current value of DFDX for each slip system (OUTPUT)
C
C PROP -- material constants characterizing the strain hardening
C law (INPUT)
C
C For the current power law strain hardening law
C PROP(1) -- power law hardening exponent
C PROP(1) = infinity corresponds to a rate-independent
C material
C PROP(2) -- coefficient in front of power law hardening
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL F, DFDX
DIMENSION GAMMAR(NSLIP), TAUSLP(NSLIP), GSLIP(NSLIP),
2 FSLIP(NSLIP), DFDXSP(NSLIP), PROP(8)
DO I=1,NSLIP
X=TAUSLP(I)/GSLIP(I)
FSLIP(I)=F(X,PROP)
DFDXSP(I)=DFDX(X,PROP)
END DO
RETURN
END
C-----------------------------------
C----- Use single precision on cray
C
REAL*8 FUNCTION F(X,PROP)
C----- User-supplied function subprogram which gives shear
C strain-rate for each slip system based on current values of
C resolved shear stress and current strength
C
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION PROP(8)
F=PROP(2)*(ABS(X))**PROP(1)*DSIGN(1.D0,X)
RETURN
END
C-----------------------------------
C----- Use single precision on cray
C
REAL*8 FUNCTION DFDX(X,PROP)
C----- User-supplied function subprogram dF/dX, where x is the
C ratio of resolved shear stress over current strength
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION PROP(8)
DFDX=PROP(1)*PROP(2)*(ABS(X))**(PROP(1)-1.)
RETURN
END
C----------------------------------------------------------------------
CFIXA
SUBROUTINE LATENTHARDEN (GAMMAR, TAUSLP, GSLIP, GMSLTL, GAMTOL,
2 NSLIP, NSLPTL, NSET, H, PROP, ND)
CFIXB
C----- This subroutine calculates the current self- and latent-
C hardening moduli for all slip systems in a rate-dependent single
C crystal. Two kinds of hardening law are used here. The first
C law, proposed by Asaro, and Pierce et al, assumes a HYPER SECANT
C relation between self- and latent-hardening moduli and overall
C shear strain. The Bauschinger effect has been neglected. The
C second is Bassani's hardening law, which gives an explicit
C expression of slip interactions between slip systems. The
C classical three stage hardening for FCC single crystal could be
C simulated.
C----- The hardening coefficients are assumed the same for all slip
C systems in each set, though they could be different from set to
C set, e.g. <110>{111} and <110>{100}.
C----- Users who want to use their own self- and latent-hardening law
C may change the function subprograms HSELF (self hardening) and
C HLATNT (latent hardening). The parameters characterizing these
C hardening laws are passed into HSELF and HLATNT through array
C PROP.
C----- Function subprograms:
C
C HSELF -- User-supplied self-hardening function in a slip
C system
C
C HLATNT -- User-supplied latent-hardening function
C----- Variables:
C
C GAMMAR -- shear strain in all slip systems at the start of time
C step (INPUT)
C TAUSLP -- resolved shear stress in all slip systems (INPUT)
C GSLIP -- current strength (INPUT)
CFIX GMSLTL -- total cumulative shear strains on each individual slip system
CFIX (INPUT)
C GAMTOL -- total cumulative shear strains over all slip systems
C (INPUT)
C NSLIP -- number of slip systems in each set (INPUT)
C NSLPTL -- total number of slip systems in all the sets (INPUT)
C NSET -- number of sets of slip systems (INPUT)
C
C H -- current value of self- and latent-hardening moduli
C (OUTPUT)
C H(i,i) -- self-hardening modulus of the ith slip system
C (no sum over i)
C H(i,j) -- latent-hardening molulus of the ith slip
C system due to a slip in the jth slip system
C (i not equal j)
C
C PROP -- material constants characterizing the self- and latent-
C hardening law (INPUT)
C
C For the HYPER SECANT hardening law
C PROP(1,i) -- initial hardening modulus H0 in the ith
C set of slip systems
C PROP(2,i) -- saturation stress TAUs in the ith set of
C slip systems
C PROP(3,i) -- initial critical resolved shear stress
C TAU0 in the ith set of slip systems
C PROP(9,i) -- ratio of latent to self-hardening Q in the
C ith set of slip systems
C PROP(10,i)-- ratio of latent-hardening from other sets
C of slip systems to self-hardening in the
C ith set of slip systems Q1
C
C For Bassani's hardening law
C PROP(1,i) -- initial hardening modulus H0 in the ith
C set of slip systems
C PROP(2,i) -- stage I stress TAUI in the ith set of
C slip systems (or the breakthrough stress
C where large plastic flow initiates)
C PROP(3,i) -- initial critical resolved shear stress
C TAU0 in the ith set of slip systems
C PROP(4,i) -- hardening modulus during easy glide Hs in
C the ith set of slip systems
C PROP(5,i) -- amount of slip Gamma0 after which a given
C interaction between slip systems in the
C ith set reaches peak strength
C PROP(6,i) -- amount of slip Gamma0 after which a given
C interaction between slip systems in the
C ith set and jth set (i not equal j)
C reaches peak strength
C PROP(7,i) -- representing the magnitude of the strength
C of interaction in the ith set of slip
C system
C PROP(8,i) -- representing the magnitude of the strength
C of interaction between the ith set and jth
C set of system
C PROP(9,i) -- ratio of latent to self-hardening Q in the
C ith set of slip systems
C PROP(10,i)-- ratio of latent-hardening from other sets
C of slip systems to self-hardening in the
C ith set of slip systems Q1
C
C ND -- leading dimension of arrays defined in subroutine UMAT
C (INPUT)
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL HSELF, HLATNT
CFIXA
DIMENSION GAMMAR(NSLPTL), TAUSLP(NSLPTL), GMSLTL(NSLPTL),
2 GSLIP(NSLPTL), NSLIP(NSET), PROP(16,NSET),
3 H(ND,NSLPTL)
CFIXB
CHECK=0.
DO I=1,NSET
DO J=4,8
CHECK=CHECK+ABS(PROP(J,I))
END DO
END DO
C----- CHECK=0 -- HYPER SECANT hardening law
C otherwise -- Bassani's hardening law
ISELF=0
DO I=1,NSET
ISET=I
DO J=1,NSLIP(I)
ISELF=ISELF+1
DO LATENT=1,NSLPTL
IF (LATENT.EQ.ISELF) THEN
CFIXA
H(LATENT,ISELF)=HSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,
2 NSET,NSLIP,PROP(1,I),CHECK,
3 ISELF,ISET)
CFIXB
ELSE
CFIXA
H(LATENT,ISELF)=HLATNT(GAMMAR,GMSLTL,GAMTOL,NSLPTL,
2 NSET,NSLIP,PROP(1,I),CHECK,
3 ISELF,ISET,LATENT)
CFIXB
END IF
END DO
END DO
END DO
RETURN
END
C-----------------------------------
C----- Use single precision on cray
CFIXA
REAL*8 FUNCTION HSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,
2 NSLIP,PROP,CHECK,ISELF,ISET)
CFIXB
C----- User-supplied self-hardening function in a slip system
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
CFIXA
DIMENSION GAMMAR(NSLPTL), NSLIP(NSET), PROP(16),
2 GMSLTL(NSLPTL)
CFIXB
IF (CHECK.EQ.0.) THEN
C----- HYPER SECANT hardening law by Asaro, Pierce et al
TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
HSELF=PROP(1)*TERM2**2
ELSE
C----- Bassani's hardening law
CFIXA
TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)
ID=0
G=1.
DO I=1,NSET
IF (I.EQ.ISET) THEN
GAMMA0=PROP(5)
FAB=PROP(7)
ELSE
GAMMA0=PROP(6)
FAB=PROP(8)
END IF
DO J=1,NSLIP(I)
ID=ID+1
IF (ID.NE.ISELF) THEN
CFIXA
G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
END IF
END DO
END DO
HSELF=F*G
END IF
RETURN
END
C-----------------------------------
C----- Use single precision on cray
CFIXA
REAL*8 FUNCTION HLATNT(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,
2 NSLIP,PROP,CHECK,ISELF,ISET,LATENT)
CFIXB
C----- User-supplied latent-hardening function
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
CFIXA
DIMENSION GAMMAR(NSLPTL), NSLIP(NSET), PROP(16),
2 GMSLTL(NSLPTL)
CFIXB
ILOWER=0
IUPPER=NSLIP(1)
IF (ISET.GT.1) THEN
DO K=2,ISET
ILOWER=ILOWER+NSLIP(K-1)
IUPPER=IUPPER+NSLIP(K)
END DO
END IF
IF (LATENT.GT.ILOWER.AND.LATENT.LE.IUPPER) THEN
Q=PROP(9)
ELSE
Q=PROP(10)
END IF
IF (CHECK.EQ.0.) THEN
C----- HYPER SECANT hardening law by Asaro, Pierce et al
TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
HLATNT=PROP(1)*TERM2**2*Q
ELSE
C----- Bassani's hardening law
CFIXA
TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)
ID=0
G=1.
DO I=1,NSET
IF (I.EQ.ISET) THEN
GAMMA0=PROP(5)
FAB=PROP(7)
ELSE
GAMMA0=PROP(6)
FAB=PROP(8)
END IF
DO J=1,NSLIP(I)
ID=ID+1
IF (ID.NE.ISELF) THEN
CFIXA
G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
END IF
END DO
END DO
HLATNT=F*G*Q
END IF
RETURN
END
C----------------------------------------------------------------------
CFIXA
SUBROUTINE ITERATION (GAMMAR, TAUSLP, GSLIP, GMSLTL, GAMTOL,
2 NSLPTL, NSET, NSLIP, ND, PROP, DGAMOD,
3 DHDGDG)
CFIXB
C----- This subroutine generates arrays for the Newton-Rhapson
C iteration method.
C----- Users who want to use their own self- and latent-hardening law
C may change the function subprograms DHSELF (self hardening) and
C DHLATN (latent hardening). The parameters characterizing these
C hardening laws are passed into DHSELF and DHLATN through array
C PROP.
C----- Function subprograms:
C
C DHSELF -- User-supplied function of the derivative of self-
C hardening moduli
C
C DHLATN -- User-supplied function of the derivative of latent-
C hardening moduli
C----- Variables:
C
C GAMMAR -- shear strain in all slip systems at the start of time
C step (INPUT)
C TAUSLP -- resolved shear stress in all slip systems (INPUT)
C GSLIP -- current strength (INPUT)
CFIX GMSLTL -- total cumulative shear strains on each individual slip system
CFIX (INPUT)
C GAMTOL -- total cumulative shear strains over all slip systems
C (INPUT)
C NSLPTL -- total number of slip systems in all the sets (INPUT)
C NSET -- number of sets of slip systems (INPUT)
C NSLIP -- number of slip systems in each set (INPUT)
C ND -- leading dimension of arrays defined in subroutine UMAT
C (INPUT)
C
C PROP -- material constants characterizing the self- and latent-
C hardening law (INPUT)
C
C For the HYPER SECANT hardening law
C PROP(1,i) -- initial hardening modulus H0 in the ith
C set of slip systems
C PROP(2,i) -- saturation stress TAUs in the ith set of
C slip systems
C PROP(3,i) -- initial critical resolved shear stress
C TAU0 in the ith set of slip systems
C PROP(9,i) -- ratio of latent to self-hardening Q in the
C ith set of slip systems
C PROP(10,i)-- ratio of latent-hardening from other sets
C of slip systems to self-hardening in the
C ith set of slip systems Q1
C
C For Bassani's hardening law
C PROP(1,i) -- initial hardening modulus H0 in the ith
C set of slip systems
C PROP(2,i) -- stage I stress TAUI in the ith set of
C slip systems (or the breakthrough stress
C where large plastic flow initiates)
C PROP(3,i) -- initial critical resolved shear stress
C TAU0 in the ith set of slip systems
C PROP(4,i) -- hardening modulus during easy glide Hs in
C the ith set of slip systems
C PROP(5,i) -- amount of slip Gamma0 after which a given
C interaction between slip systems in the
C ith set reaches peak strength
C PROP(6,i) -- amount of slip Gamma0 after which a given
C interaction between slip systems in the
C ith set and jth set (i not equal j)
C reaches peak strength
C PROP(7,i) -- representing the magnitude of the strength
C of interaction in the ith set of slip
C system
C PROP(8,i) -- representing the magnitude of the strength
C of interaction between the ith set and jth
C set of system
C PROP(9,i) -- ratio of latent to self-hardening Q in the
C ith set of slip systems
C PROP(10,i)-- ratio of latent-hardening from other sets
C of slip systems to self-hardening in the
C ith set of slip systems Q1
C
C----- Arrays for iteration:
C
C DGAMOD (INPUT)
C
C DHDGDG (OUTPUT)
C
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL DHSELF, DHLATN
CFIXA
DIMENSION GAMMAR(NSLPTL), TAUSLP(NSLPTL), GMSLTL(NSLPTL),
2 GSLIP(NSLPTL), NSLIP(NSET), PROP(16,NSET),
3 DGAMOD(NSLPTL), DHDGDG(ND,NSLPTL)
CFIXB
CHECK=0.
DO I=1,NSET
DO J=4,8
CHECK=CHECK+ABS(PROP(J,I))
END DO
END DO
C----- CHECK=0 -- HYPER SECANT hardening law
C otherwise -- Bassani's hardening law
ISELF=0
DO I=1,NSET
ISET=I
DO J=1,NSLIP(I)
ISELF=ISELF+1
DO KDERIV=1,NSLPTL
DHDGDG(ISELF,KDERIV)=0.
DO LATENT=1,NSLPTL
IF (LATENT.EQ.ISELF) THEN
CFIXA
DHDG=DHSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,
2 NSLIP,PROP(1,I),CHECK,ISELF,ISET,
3 KDERIV)
CFIXB
ELSE
CFIXA
DHDG=DHLATN(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,
2 NSLIP,PROP(1,I),CHECK,ISELF,ISET,
3 LATENT,KDERIV)
CFIXB
END IF
DHDGDG(ISELF,KDERIV)=DHDGDG(ISELF,KDERIV)+
2 DHDG*ABS(DGAMOD(LATENT))
END DO
END DO
END DO
END DO
RETURN
END
C-----------------------------------
C----- Use single precision on cray
CFIXA
REAL*8 FUNCTION DHSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,
2 NSLIP,PROP,CHECK,ISELF,ISET,
3 KDERIV)
CFIXB
C----- User-supplied function of the derivative of self-hardening
C moduli
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
CFIXA
DIMENSION GAMMAR(NSLPTL), GMSLTL(NSLPTL),
2 NSLIP(NSET), PROP(16)
CFIXB
IF (CHECK.EQ.0.) THEN
C----- HYPER SECANT hardening law by Asaro, Pierce et al
TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
TERM3=PROP(1)/(PROP(2)-PROP(3))*DSIGN(1.D0,GAMMAR(KDERIV))
DHSELF=-2.*PROP(1)*TERM2**2*TANH(TERM1)*TERM3
ELSE
C----- Bassani's hardening law
CFIXA
TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
TERM3=(PROP(1)-PROP(4))/(PROP(2)-PROP(3))
IF (KDERIV.EQ.ISELF) THEN
F=-2.*(PROP(1)-PROP(4))*TERM2**2*TANH(TERM1)*TERM3
ID=0
G=1.
DO I=1,NSET
IF (I.EQ.ISET) THEN
GAMMA0=PROP(5)
FAB=PROP(7)
ELSE
GAMMA0=PROP(6)
FAB=PROP(8)
END IF
DO J=1,NSLIP(I)
ID=ID+1
CFIXA
IF (ID.NE.ISELF) G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
END DO
END DO
ELSE
F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)
ILOWER=0
IUPPER=NSLIP(1)
IF (ISET.GT.1) THEN
DO K=2,ISET
ILOWER=ILOWER+NSLIP(K-1)
IUPPER=IUPPER+NSLIP(K)
END DO
END IF
IF (KDERIV.GT.ILOWER.AND.KDERIV.LE.IUPPER) THEN
GAMMA0=PROP(5)
FAB=PROP(7)
ELSE
GAMMA0=PROP(6)
FAB=PROP(8)
END IF
CFIXA
TERM4=GMSLTL(KDERIV)/GAMMA0
CFIXB
TERM5=2.*EXP(-TERM4)/(1.+EXP(-2.*TERM4))
G=FAB/GAMMA0*TERM5**2
END IF
DHSELF=F*G
END IF
RETURN
END
C-----------------------------------
C----- Use single precision on cray
CFIXA
REAL*8 FUNCTION DHLATN(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,
2 NSLIP,PROP,CHECK,ISELF,ISET,LATENT,
3 KDERIV)
CFIXB
C----- User-supplied function of the derivative of latent-hardening
C moduli
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
CFIXA
DIMENSION GAMMAR(NSLPTL), GMSLTL(NSLPTL), NSLIP(NSET),
2 PROP(16)
CFIXB
ILOWER=0
IUPPER=NSLIP(1)
IF (ISET.GT.1) THEN
DO K=2,ISET
ILOWER=ILOWER+NSLIP(K-1)
IUPPER=IUPPER+NSLIP(K)
END DO
END IF
IF (LATENT.GT.ILOWER.AND.LATENT.LE.IUPPER) THEN
Q=PROP(9)
ELSE
Q=PROP(10)
END IF
IF (CHECK.EQ.0.) THEN
C----- HYPER SECANT hardening law by Asaro, Pierce et al
TERM1=PROP(1)*GAMTOL/(PROP(2)-PROP(3))
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
TERM3=PROP(1)/(PROP(2)-PROP(3))*DSIGN(1.D0,GAMMAR(KDERIV))
DHLATN=-2.*PROP(1)*TERM2**2*TANH(TERM1)*TERM3*Q
ELSE
C----- Bassani's hardening law
CFIXA
TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)/(PROP(2)-PROP(3))
CFIXB
TERM2=2.*EXP(-TERM1)/(1.+EXP(-2.*TERM1))
TERM3=(PROP(1)-PROP(4))/(PROP(2)-PROP(3))
IF (KDERIV.EQ.ISELF) THEN
F=-2.*(PROP(1)-PROP(4))*TERM2**2*TANH(TERM1)*TERM3
ID=0
G=1.
DO I=1,NSET
IF (I.EQ.ISET) THEN
GAMMA0=PROP(5)
FAB=PROP(7)
ELSE
GAMMA0=PROP(6)
FAB=PROP(8)
END IF
DO J=1,NSLIP(I)
ID=ID+1
CFIXA
IF (ID.NE.ISELF) G=G+FAB*TANH(GMSLTL(ID)/GAMMA0)
CFIXB
END DO
END DO
ELSE
F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)
ILOWER=0
IUPPER=NSLIP(1)
IF (ISET.GT.1) THEN
DO K=2,ISET
ILOWER=ILOWER+NSLIP(K-1)
IUPPER=IUPPER+NSLIP(K)
END DO
END IF
IF (KDERIV.GT.ILOWER.AND.KDERIV.LE.IUPPER) THEN
GAMMA0=PROP(5)
FAB=PROP(7)
ELSE
GAMMA0=PROP(6)
FAB=PROP(8)
END IF
CFIXA
TERM4=GMSLTL(KDERIV)/GAMMA0
CFIXB
TERM5=2.*EXP(-TERM4)/(1.+EXP(-2.*TERM4))
G=FAB/GAMMA0*TERM5**2
END IF
DHLATN=F*G*Q
END IF
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE LUDCMP (A, N, NP, INDX, D)
C----- LU decomposition
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NMAX=200, TINY=1.0E-20)
DIMENSION A(NP,NP), INDX(N), VV(NMAX)
D=1.
DO I=1,N
AAMAX=0.
DO J=1,N
IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
END DO
IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
VV(I)=1./AAMAX
END DO
DO J=1,N
DO I=1,J-1
SUM=A(I,J)
DO K=1,I-1
SUM=SUM-A(I,K)*A(K,J)
END DO
A(I,J)=SUM
END DO
AAMAX=0.
DO I=J,N
SUM=A(I,J)
DO K=1,J-1
SUM=SUM-A(I,K)*A(K,J)
END DO
A(I,J)=SUM
DUM=VV(I)*ABS(SUM)
IF (DUM.GE.AAMAX) THEN
IMAX=I
AAMAX=DUM
END IF
END DO
IF (J.NE.IMAX) THEN
DO K=1,N
DUM=A(IMAX,K)
A(IMAX,K)=A(J,K)
A(J,K)=DUM
END DO
D=-D
VV(IMAX)=VV(J)
END IF
INDX(J)=IMAX
IF (A(J,J).EQ.0.) A(J,J)=TINY
IF (J.NE.N) THEN
DUM=1./A(J,J)
DO I=J+1,N
A(I,J)=A(I,J)*DUM
END DO
END IF
END DO
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE LUBKSB (A, N, NP, INDX, B)
C----- Linear equation solver based on LU decomposition
C----- Use single precision on cray
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(NP,NP), INDX(N), B(N)
II=0
DO I=1,N
LL=INDX(I)
SUM=B(LL)
B(LL)=B(I)
IF (II.NE.0) THEN
DO J=II,I-1
SUM=SUM-A(I,J)*B(J)
END DO
ELSE IF (SUM.NE.0.) THEN
II=I
END IF
B(I)=SUM
END DO
DO I=N,1,-1
SUM=B(I)
IF (I.LT.N) THEN
DO J=I+1,N
SUM=SUM-A(I,J)*B(J)
END DO
END IF
B(I)=SUM/A(I,I)
END DO
RETURN
END
C ********************* 25. GET_INV_DET *********************
c this subroutine inverts a given matrix and return it
SUBROUTINE GET_INV_DET(XS,DET,NFLAG)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XS(3,3),A(3,3)
DO 10 I=1,3
I1=I+1
I2=I+2
IF(I1.GT.3) I1=I1-3
IF(I2.GT.3) I2=I2-3
DO 10 J=1,3
J1=J+1
J2=J+2
IF(J1.GT.3) J1=J1-3
IF(J2.GT.3) J2=J2-3
10 A(I,J)=XS(I1,J1)*XS(I2,J2)-XS(I1,J2)*XS(I2,J1)
DET=0.D0
DO 20 I=1,3
20 DET=DET+A(1,I)*XS(1,I)
NE=1
IF (NFLAG.EQ.0.AND.DET.LE.0.0D0) THEN
DET=0.00001
!WRITE(*,*) 'ERROR IN THE JACOBAIN, DETERMINANT', DET
C PAUSE'ERROR IN THE JACOBIAN'
ENDIF
DO 30 I=1,3
DO 30 J=1,3
30 XS(I,J)=A(J,I)/DET
RETURN
END
C ********************* 29. CLEAR *********************
c this subroutine initializes a real matrix to zero
SUBROUTINE CLEAR(A,N)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(N)
DO 10 I=1,N
10 A(I)=0.D0
RETURN
END
