Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  ertwertewr
Language: FORTRAN
Code:
            subroutine vumat( 
c Read only - 
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, 
     2  stepTime, totalTime, dt, 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 - 
     7  stressNew, stateNew, enerInternNew, enerInelasNew ) 

      include 'vaba_param.inc' 

c 3D Orthotropic Elasticity with Hashin 3d Failure criterion 

c The state variables are stored as: 
c    state(*,1)   = material point status 
c    state(*,2:7) = damping stresses 

c User defined material properties are stored as 
c  * First line: 
c     props(1) --> Young's modulus in 1-direction, E1 
c     props(2) --> Young's modulus in 2-direction, E2 
c     props(3) --> Young's modulus in 3-direction, E3 
c     props(4) --> Poisson's ratio, nu12 
c     props(5) --> Poisson's ratio, nu13 
c     props(6) --> Poisson's ratio, nu23 
c     props(7) --> Shear modulus, G12 
c     props(8) --> Shear modulus, G13 

c  * Second line: 
c     props(9)  --> Shear modulus, G23 
c     props(10) --> beta damping parameter 
c     props(11) --> "not used" 
c     props(12) --> "not used" 
c     props(13) --> "not used" 
c     props(14) --> "not used" 
c     props(15) --> "not used" 
c     props(16) --> "not used" 

c  * Third line: 
c     props(17) --> Ultimate tens stress in 1-direction, sigu1t 
c     props(18) --> Ultimate comp stress in 1-direction, sigu1c 
c     props(19) --> Ultimate tens stress in 2-direction, sigu2t 
c     props(20) --> Ultimate comp stress in 2-direction, sigu2c 
c     props(21) --> Ultimate tens stress in 2-direction, sigu3t 
c     props(22) --> Ultimate comp stress in 2-direction, sigu3c 
c     props(23) --> "not used" 
c     props(24) --> "not used" 

c  * Fourth line: 
c     props(25) --> Ultimate shear stress, sigu12 
c     props(26) --> Ultimate shear stress, sigu13 
c     props(27) --> Ultimate shear stress, sigu23 
c     props(28) --> "not used" 
c     props(29) --> "not used" 
c     props(30) --> "not used" 
c     props(31) --> "not used" 
c     props(32) --> "not used" 


      dimension props(nprops), density(nblock), 
     1  coordMp(nblock,*), 
     2  charLength(*), strainInc(nblock,ndir+nshr), 
     3  relSpinInc(nblock,nshr), tempOld(nblock), 
     4  stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr), 
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), 
     6  stateOld(nblock,nstatev), enerInternOld(nblock), 
     7  enerInelasOld(nblock), tempNew(*), 
     8  stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), 
     9  fieldNew(nblock,nfieldv), stressNew(nblock,ndir+nshr), 
     1  stateNew(nblock,nstatev), 
     2  enerInternNew(nblock), enerInelasNew(nblock) 

      character*80 cmname 

      parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0 ) 

      parameter( 
     *     i_svd_DmgFiberT   = 1, 
     *     i_svd_DmgFiberC   = 2, 
     *     i_svd_DmgMatrixT  = 3, 
     *     i_svd_DmgMatrixC  = 4, 
     *     i_svd_statusMp   = 5, 
     *     i_svd_dampStress = 6, 
c     *    i_svd_dampStressXx = 6, 
c     *    i_svd_dampStressYy = 7, 
c     *    i_svd_dampStressZz = 8, 
c     *    i_svd_dampStressXy = 9, 
c     *    i_svd_dampStressYz = 10, 
c     *    i_svd_dampStressZx = 11, 
     *     i_svd_Strain   = 12, 
c     *    i_svd_StrainXx = 12, 
c     *    i_svd_StrainYy = 13, 
c     *    i_svd_StrainZz = 14, 
c     *    i_svd_StrainXy = 15, 
c     *    i_svd_StrainYz = 16, 
c     *    i_svd_StrainZx = 17, 
     *     n_svd_required = 17 ) 

      parameter( 
     *     i_s33_Xx = 1, 
     *     i_s33_Yy = 2, 
     *     i_s33_Zz = 3, 
     *     i_s33_Xy = 4, 
     *     i_s33_Yz = 5, 
     *     i_s33_Zx = 6 ) 

* Structure of property array 
      parameter ( 
     *     i_pro_E1    = 1, 
     *     i_pro_E2    = 2, 
     *     i_pro_E3    = 3, 
     *     i_pro_nu12  = 4, 
     *     i_pro_nu13  = 5, 
     *     i_pro_nu23  = 6, 
     *     i_pro_G12   = 7, 
     *     i_pro_G13   = 8, 
     *     i_pro_G23   = 9, 

     *     i_pro_beta  = 10, 

     *     i_pro_sigu1t = 17, 
     *     i_pro_sigu1c = 18, 
     *     i_pro_sigu2t = 19, 
     *     i_pro_sigu2c = 20, 
     *     i_pro_sigu3t = 21, 
     *     i_pro_sigu3c = 22, 
     *     i_pro_sigu12 = 25, 
     *     i_pro_sigu13 = 26, 
     *     i_pro_sigu23 = 27 ) 
* Temporary arrays 
      dimension eigen(maxblk*3) 

* Read material properties 

      E1 = props(i_pro_E1) 
      E2 = props(i_pro_E2) 
      E3 = props(i_pro_E3) 
      xnu12 = props(i_pro_nu12) 
      xnu13 = props(i_pro_nu13) 
      xnu23 = props(i_pro_nu23) 
      G12 = props(i_pro_G12) 
      G13 = props(i_pro_G13) 
      G23 = props(i_pro_G23) 

      xnu21 = xnu12 * E2 / E1 
      xnu31 = xnu13 * E3 / E1 
      xnu32 = xnu23 * E3 / E2 


* Compute terms of stiffness matrix 
      gg = one / ( one - xnu12*xnu21 - xnu23*xnu32 - xnu31*xnu13 
     *     - two*xnu21*xnu32*xnu13 ) 
      C11  = E1 * ( one - xnu23*xnu32 ) * gg 
      C22  = E2 * ( one - xnu13*xnu31 ) * gg 
      C33  = E3 * ( one - xnu12*xnu21 ) * gg 
      C12  = E1 * ( xnu21 + xnu31*xnu23 ) * gg 
      C13  = E1 * ( xnu31 + xnu21*xnu32 ) * gg 
      C23  = E2 * ( xnu32 + xnu12*xnu31 ) * gg 

      f1t = props(i_pro_sigu1t) 
      f1c = props(i_pro_sigu1c) 
      f2t = props(i_pro_sigu2t) 
      f2c = props(i_pro_sigu2c) 
      f3t = props(i_pro_sigu3t) 
      f3c = props(i_pro_sigu3c) 
      f12 = props(i_pro_sigu12) 
      f13 = props(i_pro_sigu13) 
      f23 = props(i_pro_sigu23) 

      beta = props(i_pro_beta) 

* Assume purely elastic material at the beginning of the analysis 
*       
      if ( totalTime .eq. zero ) then 
         if (nstatev .lt. n_svd_Required) then 
            call xplb_abqerr(-2,'Subroutine VUMAT requires the '// 
     *           'specification of %I state variables. Check the '// 
     *           'definition of *DEPVAR in the input file.', 
     *           n_svd_Required,zero,' ') 
            call xplb_exit 
         end if 
         call OrthoEla3dExp ( nblock, 
     *        stateOld(1,i_svd_DmgFiberT), 
     *        stateOld(1,i_svd_DmgFiberC), 
     *        stateOld(1,i_svd_DmgMatrixT), 
     *        stateOld(1,i_svd_DmgMatrixC), 
     *        C11, C22, C33, C12, C23, C13, G12, G23, G13, 
     *        strainInc, 
     *        stressNew ) 
         return 
      end if 

*  Update total elastic strain 
      call strainUpdate ( nblock, strainInc, 
     *     stateOld(1,i_svd_strain), stateNew(1,i_svd_strain) ) 

* Stress update 
      call OrthoEla3dExp ( nblock, 
     *     stateOld(1,i_svd_DmgFiberT), 
     *     stateOld(1,i_svd_DmgFiberC), 
     *     stateOld(1,i_svd_DmgMatrixT), 
     *     stateOld(1,i_svd_DmgMatrixC), 
     *     C11, C22, C33, C12, C23, C13, G12, G23, G13, 
     *     stateNew(1,i_svd_strain), 
     *     stressNew ) 

* Failure evaluation 

      call copyr ( nblock, 
     *     stateOld(1,i_svd_DmgFiberT), stateNew(1,i_svd_DmgFiberT) ) 
      call copyr ( nblock, 
     *     stateOld(1,i_svd_DmgFiberC), stateNew(1,i_svd_DmgFiberC) ) 
      call copyr ( nblock, 
     *     stateOld(1,i_svd_DmgMatrixT), stateNew(1,i_svd_DmgMatrixT) ) 
      call copyr ( nblock, 
     *     stateOld(1,i_svd_DmgMatrixC), stateNew(1,i_svd_DmgMatrixC) ) 
      nDmg = 0 
      call eig33Anal ( nblock, stretchNew, eigen ) 
      call Hashin3d  ( nblock, nDmg, 
     *     f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13, 
     *     stateNew(1,i_svd_DmgFiberT), 
     *     stateNew(1,i_svd_DmgFiberC), 
     *     stateNew(1,i_svd_DmgMatrixT), 
     *     stateNew(1,i_svd_DmgMatrixC), 
     *     stateNew(1,i_svd_statusMp), 
     *     stressNew, eigen ) 
*     -- Recompute stresses if new Damage is occurring 
      if ( nDmg .gt. 0 ) then 
         call OrthoEla3dExp ( nblock, 
     *        stateNew(1,i_svd_DmgFiberT), 
     *        stateNew(1,i_svd_DmgFiberC), 
     *        stateNew(1,i_svd_DmgMatrixT), 
     *        stateNew(1,i_svd_DmgMatrixC), 
     *        C11, C22, C33, C12, C23, C13, G12, G23, G13, 
     *        stateNew(1,i_svd_strain), 
     *        stressNew ) 
      end if 

* Beta damping 
      if ( beta .gt. zero ) then 
         call betaDamping3d ( nblock, 
     *        beta, dt, strainInc, 
     *        stressOld, stressNew, 
     *        stateNew(1,i_svd_statusMp), 
     *        stateOld(1,i_svd_dampStress), 
     *        stateNew(1,i_svd_dampStress) ) 
      end if 

* Integrate the internal specific energy (per unit mass) 

      call EnergyInternal3d ( nblock, stressOld, stressNew, 
     *   strainInc, density, enerInternOld, enerInternNew ) 

      return 
      end 


************************************************************ 
*   OrthoEla3dExp: Orthotropic elasticity - 3d             * 
************************************************************ 
      subroutine OrthoEla3dExp ( nblock, 
     *     dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC, 
     *     C11, C22, C33, C12, C23, C13, G12, G23, G13, 
     *     strain, stress ) 

      include 'vaba_param.inc' 

*  Orthotropic elasticity, 3D case - 

      parameter( zero = 0.d0, one = 1.d0, two = 2.d0) 
      parameter( 
     *     i_s33_Xx = 1, 
     *     i_s33_Yy = 2, 
     *     i_s33_Zz = 3, 
     *     i_s33_Xy = 4, 
     *     i_s33_Yz = 5, 
     *     i_s33_Zx = 6, 
     *     n_s33_Car = 6 ) 

      dimension  strain(nblock,n_s33_Car), 
     *     dmgFiberT(nblock), dmgFiberC(nblock), 
     *     dmgMatrixT(nblock), dmgMatrixC(nblock), 
     *     stress(nblock,n_s33_Car) 
*     -- shear fraction in matrix tension and compression mode 
      parameter ( smt = 0.9d0, smc = 0.5d0 ) 

      do k = 1, nblock 
*     -- Compute damaged stiffness 
         dft = dmgFiberT(k) 
         dfc = dmgFiberC(k) 
         dmt = dmgMatrixT(k) 
         dmc = dmgMatrixC(k) 
         df = one - ( one - dft ) * ( one - dfc ) 

         dC11 = ( one - df ) * C11 
         dC22 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C22 
         dC33 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C33 
         dC12 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C12 
         dC23 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C23 
         dC13 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C13 
         dG12 = ( one - df ) 
     *        * ( one - smt*dmt ) * ( one - smc*dmc ) * G12 
         dG23 = ( one - df ) 
     *        * ( one - smt*dmt ) * ( one - smc*dmc ) * G23 
         dG13 = ( one - df ) 
     *        * ( one - smt*dmt ) * ( one - smc*dmc ) * G13 
*     -- Stress update 
         stress(k,i_s33_Xx) = dC11 * strain(k,i_s33_Xx) 
     *        + dC12 * strain(k,i_s33_Yy) 
     *        + dC13 * strain(k,i_s33_Zz) 
         stress(k,i_s33_Yy) = dC12 * strain(k,i_s33_Xx) 
     *        + dC22 * strain(k,i_s33_Yy) 
     *        + dC23 * strain(k,i_s33_Zz) 
         stress(k,i_s33_Zz) = dC13 * strain(k,i_s33_Xx) 
     *        + dC23 * strain(k,i_s33_Yy) 
     *        + dC33 * strain(k,i_s33_Zz) 
         stress(k,i_s33_Xy) = two * dG12 * strain(k,i_s33_Xy) 
         stress(k,i_s33_Yz) = two * dG23 * strain(k,i_s33_Yz) 
         stress(k,i_s33_Zx) = two * dG13 * strain(k,i_s33_Zx) 
      end do 
*     
      return 
      end 

************************************************************ 
*   strainUpdate: Update total strain                      * 
************************************************************ 
      subroutine strainUpdate ( nblock, 
     *     strainInc, strainOld, strainNew ) 

      include 'vaba_param.inc' 

      parameter( 
     *     i_s33_Xx = 1, 
     *     i_s33_Yy = 2, 
     *     i_s33_Zz = 3, 
     *     i_s33_Xy = 4, 
     *     i_s33_Yz = 5, 
     *     i_s33_Zx = 6, 
     *     n_s33_Car = 6 ) 

      dimension strainInc(nblock,n_s33_Car), 
     *     strainOld(nblock,n_s33_Car), 
     *     strainNew(nblock,n_s33_Car) 

      do k = 1, nblock 
         strainNew(k,i_s33_Xx)= strainOld(k,i_s33_Xx) 
     *                        + strainInc(k,i_s33_Xx) 
         strainNew(k,i_s33_Yy)= strainOld(k,i_s33_Yy) 
     *                        + strainInc(k,i_s33_Yy) 
         strainNew(k,i_s33_Zz)= strainOld(k,i_s33_Zz) 
     *                        + strainInc(k,i_s33_Zz) 
         strainNew(k,i_s33_Xy)= strainOld(k,i_s33_Xy) 
     *                        + strainInc(k,i_s33_Xy) 
         strainNew(k,i_s33_Yz)= strainOld(k,i_s33_Yz) 
     *                        + strainInc(k,i_s33_Yz) 
         strainNew(k,i_s33_Zx)= strainOld(k,i_s33_Zx) 
     *                        + strainInc(k,i_s33_Zx) 
      end do 

      return 
      end 


************************************************************ 
*   Hashin3d w/ Modified Puck: Evaluate Hashin 3d failure  * 
*   criterion for fiber, Puck for matrix                   * 
************************************************************ 
      subroutine Hashin3d ( nblock, nDmg, 
     *     f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13, 
     *     dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC, 
     *     statusMp, stress, eigen ) 

      include 'vaba_param.inc' 

      parameter( zero = 0.d0, one = 1.d0, half = 0.5d0, three =3.d0 ) 
      parameter( 
     *     i_s33_Xx = 1, 
     *     i_s33_Yy = 2, 
     *     i_s33_Zz = 3, 
     *     i_s33_Xy = 4, 
     *     i_s33_Yz = 5, 
     *     i_s33_Zx = 6, 
     *     n_s33_Car = 6 ) 

      parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 ) 
      parameter(n_v3d_Car=3 ) 

      parameter ( eMax = 1.00d0, eMin = -0.8d0 ) 

      dimension  dmgFiberT(nblock), dmgFiberC(nblock), 
     *     dmgMatrixT(nblock), dmgMatrixC(nblock), 
     *     stress(nblock,n_s33_Car), 
     *     eigen(nblock,n_v3d_Car), 
     *     statusMp(nblock) 

      f1tInv = zero 
      f2tInv = zero 
      f3tInv = zero 
      f1cInv = zero 
      f2cInv = zero 
      f3cInv = zero 
      f12Inv = zero 
      f23Inv = zero 
      f13Inv = zero 

      if ( f1t .gt. zero ) f1tInv = one / f1t 
      if ( f2t .gt. zero ) f2tInv = one / f2t 
      if ( f3t .gt. zero ) f3tInv = one / f3t 
      if ( f1c .gt. zero ) f1cInv = one / f1c 
      if ( f2c .gt. zero ) f2cInv = one / f2c 
      if ( f3c .gt. zero ) f3cInv = one / f3c 
      if ( f12 .gt. zero ) f12Inv = one / f12 
      if ( f23 .gt. zero ) f23Inv = one / f23 
      if ( f13 .gt. zero ) f13Inv = one / f13 

      do k = 1, nblock 
         if ( statusMp(k) .eq. one ) then 
*     
         lFail = 0 

         s11 = stress(k,i_s33_Xx) 
         s22 = stress(k,i_s33_Yy) 
         s33 = stress(k,i_s33_Zz) 
         s12 = stress(k,i_s33_Xy) 
         s23 = stress(k,i_s33_Yz) 
         s13 = stress(k,i_s33_Zx) 

*     Evaluate Fiber modes 
         if ( s11 .gt. zero ) then 
*     -- Tensile Fiber Mode 
            rft = (s11*f1tInv )**2 + (s12*f12Inv )**2 + (s13*f13Inv )**2 
            if ( rft .ge. one ) then 
               lDmg = 1 
               dmgFiberT(k) = one 
            end if 
         else if ( s11 .lt. zero ) then 
*     -- Compressive Fiber Mode 
            rfc = abs(s11) * f1cInv 
            if ( rfc .ge. one ) then 
               lDmg = 1 
               dmgFiberC(k) = one 
            end if 
         end if 

*     Evaluate Matrix Modes 
         if ( ( s22 + s33 ) .gt. zero ) then 
*     -- Tensile Matrix mode 
            rmt = ( s11 * half * f1tInv )**2 
     *           + ( s22**2 * abs(f2tInv * f2cInv) ) 
     *           + ( s12 * f12Inv )**2 
     *           + ( s22 * (f2tInv + f2cInv) ) 
            if ( rmt .ge. one ) then 
               lDmg = 1 
               dmgMatrixT(k) = one 
            end if 
         else if ( ( s22 + s33 ) .lt. zero ) then 
*     -- Compressive Matrix Mode 
            rmc = ( s11 * half * f1tInv )**2 
     *           + ( s22**2 * abs(f2tInv * f2cInv) ) 
     *           + ( s12 * f12Inv )**2 
     *           + ( s22 * (f2tInv + f2cInv) ) 
            if ( rmc .ge. one ) then 
               lDmg = 1 
               dmgMatrixC(k) = one 
            end if 
         end if 

         eigMax=max(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z)) 
         eigMin=min(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z)) 
         enomMax = eigMax - one 
         enomMin = eigMin - one 

         if ( enomMax .gt. eMax .or. 
     *        enomMin .lt. eMin .or. 
     *        dmgFiberT(k) .eq. one ) then 
            statusMp(k) = zero 
         end if 

         nDmg = nDmk + lDmg 

         end if 

      end do 

      return 
      end 


************************************************************ 
*   betaDamping: Add beta damping                          * 
************************************************************ 
      subroutine betaDamping3d ( nblock, 
     *     beta, dt, strainInc, sigOld, sigNew, 
     *     statusMp, sigDampOld, sigDampNew ) 

      include 'vaba_param.inc' 

      parameter( 
     *     i_s33_Xx = 1, 
     *     i_s33_Yy = 2, 
     *     i_s33_Zz = 3, 
     *     i_s33_Xy = 4, 
     *     i_s33_Yz = 5, 
     *     i_s33_Zx = 6, 
     *     n_s33_Car = 6 ) 

      dimension  sigOld(nblock,n_s33_Car), 
     *     sigNew(nblock,n_s33_Car), 
     *     strainInc(nblock,n_s33_Car), 
     *     statusMp(nblock), 
     *     sigDampOld(nblock,n_s33_Car), 
     *     sigDampNew(nblock,n_s33_Car)       

      parameter ( zero = 0.d0, one = 1.d0, two=2.0d0, 
     *     half = 0.5d0, third = 1.d0/3.d0 ) 
      parameter ( asmall = 1.d-16 ) 
*     
      betaddt =  beta / dt 

      do k =1 , nblock 
         sigDampNew(k,i_s33_Xx) = betaddt * statusMp(k) * 
     *        ( sigNew(k,i_s33_Xx) 
     *        - ( sigOld(k,i_s33_Xx) - sigDampOld(k,i_s33_Xx) ) ) 
         sigDampNew(k,i_s33_Yy) = betaddt * statusMp(k) * 
     *        ( sigNew(k,i_s33_Yy) 
     *        - ( sigOld(k,i_s33_Yy) - sigDampOld(k,i_s33_Yy) ) ) 
         sigDampNew(k,i_s33_Zz) = betaddt * statusMp(k) * 
     *        ( sigNew(k,i_s33_Zz) 
     *        - ( sigOld(k,i_s33_Zz) - sigDampOld(k,i_s33_Zz) ) ) 
         sigDampNew(k,i_s33_Xy) = betaddt * statusMp(k) * 
     *        ( sigNew(k,i_s33_Xy) 
     *        - ( sigOld(k,i_s33_Xy) - sigDampOld(k,i_s33_Xy) ) ) 
         sigDampNew(k,i_s33_Yz) = betaddt * statusMp(k) * 
     *        ( sigNew(k,i_s33_Yz) 
     *        - ( sigOld(k,i_s33_Yz) - sigDampOld(k,i_s33_Yz) ) ) 
         sigDampNew(k,i_s33_Zx) = betaddt * statusMp(k) * 
     *        ( sigNew(k,i_s33_Zx) 
     *        - ( sigOld(k,i_s33_Zx) - sigDampOld(k,i_s33_Zx) ) ) 

         sigNew(k,i_s33_Xx) = sigNew(k,i_s33_Xx)+sigDampNew(k,i_s33_Xx) 
         sigNew(k,i_s33_Yy) = sigNew(k,i_s33_Yy)+sigDampNew(k,i_s33_Yy) 
         sigNew(k,i_s33_Zz) = sigNew(k,i_s33_Zz)+sigDampNew(k,i_s33_Zz) 
         sigNew(k,i_s33_Xy) = sigNew(k,i_s33_Xy)+sigDampNew(k,i_s33_Xy) 
         sigNew(k,i_s33_Yz) = sigNew(k,i_s33_Yz)+sigDampNew(k,i_s33_Yz) 
         sigNew(k,i_s33_Zx) = sigNew(k,i_s33_Zx)+sigDampNew(k,i_s33_Zx) 

      end do 
*     
      return 
      end 


************************************************************ 
*   EnergyInternal3d: Compute internal energy for 3d case  * 
************************************************************ 
      subroutine EnergyInternal3d(nblock, sigOld, sigNew , 
     *   strainInc, curDensity, enerInternOld, enerInternNew) 

      include 'vaba_param.inc' 

      parameter( 
     *     i_s33_Xx = 1, 
     *     i_s33_Yy = 2, 
     *     i_s33_Zz = 3, 
     *     i_s33_Xy = 4, 
     *     i_s33_Yz = 5, 
     *     i_s33_Zx = 6, 
     *     n_s33_Car = 6 ) 

      parameter( two = 2.d0, half = .5d0 ) 

      dimension sigOld (nblock,n_s33_Car), sigNew (nblock,n_s33_Car), 
     *     strainInc (nblock,n_s33_Car), curDensity (nblock), 
     *     enerInternOld(nblock), enerInternNew(nblock) 

      do k = 1, nblock 
         stressPower  = half * ( 
     *        ( sigOld(k,i_s33_Xx) + sigNew(k,i_s33_Xx) ) 
     *        * ( strainInc(k,i_s33_Xx) ) 
     *        +       ( sigOld(k,i_s33_Yy) + sigNew(k,i_s33_Yy) ) 
     *        * ( strainInc(k,i_s33_Yy)) 
     *        +       ( sigOld(k,i_s33_Zz) + sigNew(k,i_s33_Zz) ) 
     *        * ( strainInc(k,i_s33_Zz)) 
     *        + two * ( sigOld(k,i_s33_Xy) + sigNew(k,i_s33_Xy) ) 
     *        * strainInc(k,i_s33_Xy) 
     *        + two * ( sigOld(k,i_s33_Yz) + sigNew(k,i_s33_Yz) ) 
     *        * strainInc(k,i_s33_Yz) 
     *        + two * ( sigOld(k,i_s33_Zx) + sigNew(k,i_s33_Zx) ) 
     *        * strainInc(k,i_s33_Zx) ) 
*     
         enerInternNew(k) = enerInternOld(k) + stressPower/curDensity(k) 
      end do 
*     
      return   
      end   

************************************************************ 
*   CopyR: Copy from one array to another                  * 
************************************************************ 
      subroutine CopyR(nCopy, from, to ) 

      include 'vaba_param.inc' 

      dimension from(nCopy), to(nCopy) 

      do k = 1, nCopy 
         to(k) = from(k) 
      end do 

      return 
      end 

********************************************************************* 
***** 
* eig33Anal: Compute eigen values of a 3x3 symmetric matrix analytically * 
********************************************************************* 
***** 
      subroutine eig33Anal( nblock, sMat, eigVal ) 

      include 'vaba_param.inc' 

      parameter(i_s33_Xx=1,i_s33_Yy=2,i_s33_Zz=3 ) 
      parameter(i_s33_Xy=4,i_s33_Yz=5,i_s33_Zx=6 ) 
      parameter(i_s33_Yx=i_s33_Xy ) 
      parameter(i_s33_Zy=i_s33_Yz ) 
      parameter(i_s33_Xz=i_s33_Zx,n_s33_Car=6 ) 

      parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 ) 
      parameter(n_v3d_Car=3 ) 

      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, 
     *     three = 3.d0, half = 0.5d0, third = one / three, 
     *     pi23 = 2.094395102393195d0, 
     *     fuzz = 1.d-8, 
     *     preciz = fuzz * 1.d4 ) 

      dimension eigVal(nblock,n_v3d_Car), sMat(nblock,n_s33_Car) 

      do k = 1, nblock 
        sh  = third*(sMat(k,i_s33_Xx)+sMat(k,i_s33_Yy)+sMat(k,i_s33_Zz)) 
        s11 = sMat(k,i_s33_Xx) - sh 
        s22 = sMat(k,i_s33_Yy) - sh 
        s33 = sMat(k,i_s33_Zz) - sh 
        s12 = sMat(k,i_s33_Xy) 
        s13 = sMat(k,i_s33_Xz) 
        s23 = sMat(k,i_s33_Yz) 

        fac  = max(abs(s11), abs(s22), abs(s33)) 
        facs = max(abs(s12), abs(s13), abs(s23)) 
        if( facs .lt. (preciz*fac) ) then 
          eigVal(k,i_v3d_X) = sMat(k,i_s33_Xx) 
          eigVal(k,i_v3d_Y) = sMat(k,i_s33_Yy) 
          eigVal(k,i_v3d_Z) = sMat(k,i_s33_Zz) 
        else 
          q = third*((s12**2+s13**2+s23**2)+half*(s11**2+s22**2+s33**2)) 
          fac = two * sqrt(q) 
          if( fac .gt. fuzz ) then 
            ofac = two/fac 
          else 
            ofac = zero 
          end if 
          s11 = ofac*s11 
          s22 = ofac*s22 
          s33 = ofac*s33 
          s12 = ofac*s12 
          s13 = ofac*s13 
          s23 = ofac*s23 
          r = s12*s13*s23 
     *         + half*(s11*s22*s33-s11*s23**2-s22*s13**2-s33*s12**2) 
          if( r .ge. one-fuzz ) then 
            cos1 = -half 
            cos2 = -half 
            cos3 = one 
          else if( r .le. fuzz-one ) then 
            cos1 = -one 
            cos2 = half 
            cos3 = half 
          else 
            ang = third * acos(r) 
            cos1 = cos(ang) 
            cos2 = cos(ang+pi23) 
            cos3 =-cos1-cos2 
          end if 
          eigVal(k,i_v3d_X) = sh + fac*cos1 
          eigVal(k,i_v3d_Y) = sh + fac*cos2 
          eigVal(k,i_v3d_Z) = sh + fac*cos3 
        end if 
      end do 

      return 
      end          
Comments: