c********************************************************************************
c
c           Transformation of the 2nd PK stresses to Cauchy Stresses 
c
c********************************************************************************
c
      subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,idum0,layer,idum1,
     *                 idum2,jpltcd)
c* * * * * *
c
c     select a variable contour plotting (user subroutine).
c
c     v          variable
c     s (idss)   stress array
c     sp         stresses in preferred direction
c     etot       total strain (generalized)
c     eplas      total plastic strain
c     ecreep     total creep strain
c     t          current temperature
c     m          element number
c     nn         integration point number
c     layer      layer number
c     ndi (3)    number of direct stress components
c     nshear (3) number of shear stress components
c     lsub       control variable for the part of element assembly function
c     jpltcd     plot variable for post processing
c
c* * * * * *
c
      implicit real*8 (a-h,o-z)                                               dp
      dimension s(1),etot(1),eplas(1),ecreep(1),sp(1)
      dimension m(2),ff(3,3),dispt(24),pk2(3,3),cauchy(3,3)
      include '../common/blnk'
      include '../common/deveub'
      include '../common/dimen'
      include '../common/elmcom'
      include '../common/pival'
      include '../common/space'
      include '../common/parpts'
      include '../common/concom'
      include '../common/heat'
      include '../common/arrays'
      include '../common/array5'
      include '../common/lass'
      include '../common/array2'
c
c
c Input file must include LARGE DISP option (but no Update flag !!!!)
c Can be used with/ without the processor card (ie. processor,1,1,1) in the input
c
c
      if(jtype.eq.7) then
      lofr=(n-1)*nelstr
      itel=3
      if(nn.eq.1) then
      lovl=4
c
c Wt. functions
c
      lsub=5
c
      call ellib
      if(jparel.eq.1) then
        call scla(vars(ibp),0.d0,ngenel,n1*8,0)
        if(ismall.eq.0) call scla(vars(idumpp),0.d0,nldsrc,n1*8,0)
        if(ismall.eq.0) call scla(vars(idump1p),0.d0,nldsrc,n1*8,0)
      endif
      endif
c
c Linear beta matrix
c
      lsub=2
      if(jparel.eq.0) then
        call scla(vars(idump),0.d0,n1,1,0)
        lsub=4
      endif
      call ellib
      call elvec(dispt,nnode,vars(idsxt),ndegmx,ndegel,numnp,5,lm)
      if(jparel.eq.1) then
        call dispdp(dispt,ff,vars(idump1p),itel,1.d0,ndi,ityp,n1,nn)
      else
        call dispdr(dispt,ff,vars(idump1),itel,1.d0,ndi,ityp,n1,nn)
      endif
c
c
      if((m(1).eq.1).and.(nn.eq.1)) then
      write(50,*) 'increment=',inc,' iteration=',ncycle
      write(50,*) 'element  =',m(1),' int. point=',nn
      write(50,*) 'total disp.'
      write(50,*)  dispt
      write(50,*) 'deformation gradient'
      write(50,'(3(e15.8,1x))') ((ff(i,j),j=1,3),i=1,3)
      endif
c
c
      call pushs(s,pk2,ff,cauchy,detf)
c
c
      if((m(1).eq.1).and.(nn.eq.1)) then
      write(50,*) 'Jacobian=',detf
      write(50,*) '2nd PK stresses'
      write(50,'(3(e15.8,1x))') ((pk2(i,j),j=1,3),i=1,3)
      write(50,*) 'Cauchy stresses'
      write(50,'(3(e15.8,1x))') ((cauchy(i,j),j=1,3),i=1,3)
      endif
      endif
        if(jpltcd.eq.1) then
         v=cauchy(1,1)
        else if(jpltcd.eq.2) then
         v=cauchy(2,2)
        endif
c
      return
      end
c
c
c********************************************************************************
c
c     Push forward stresses           
c
c********************************************************************************
c
      subroutine pushs(s,pk2,ff,cauchy,detf) 
c
      implicit real*8 (a-h,o-z)
c
      dimension s(1),pk2(3,3),ff(3,3),temp(3,3),store(3,3),cauchy(3,3)
c
      pk2(1,1)=s(1)
      pk2(1,2)=s(4)
      pk2(1,3)=s(6)
      pk2(2,1)=s(4)
      pk2(2,2)=s(2)
      pk2(2,3)=s(5)
      pk2(3,1)=s(6)
      pk2(3,2)=s(5)
      pk2(3,3)=s(3)
c
        do 10 i=1,3
         do 10 j=1,3
         temp(i,j)=0.
          cauchy(i,j)=0.
 10     continue

      do 300 i=1,3
       do 200 j=1,3
        do 100 k=1,3
         temp(i,j)=temp(i,j)+pk2(i,k)*ff(j,k)
100     continue
200    continue
300   continue 
c
      call mcpy(ff,store,3,3,0)
      call inv3x3(store,store,detf,0)
c
      do 600 i=1,3
       do 500 j=1,3
        do 400 k=1,3
         cauchy(i,j)=cauchy(i,j)+ff(i,k)*temp(k,j)/detf
400     continue
500    continue
600   continue 
c
      return
      end
c
c
c********************************************************************************
c
c     Example user sub. : HYPELA 
c
c********************************************************************************
c
      SUBROUTINE HYPELA(D,G,E,DE,S,T,DT,NGENS,N,NNN,LAYER,MAT,
     * NDI,NSHEAR)
      implicit real*8 (a-h,o-z)                                               dp
      DIMENSION E(1),DE(1),T(1),DT(1),G(1),D(NGENS,NGENS),S(1),
     C ET(6)
      DO 100 I=1,NGENS
      ET(I)= E(I)+0.5*DE(I)
100   CONTINUE
      POIS=.45
      CONS1=1.0d0/SQRT(1.0d0 + 2.0d0*POIS*pois)
      CONS2=1.0d0 - POIS - 2.0d0*POIS*pois
      EVM=CONS1*SQRT(ET(1)*ET(1)+ET(2)*ET(2)+ET(3)*ET(3)+
     c (.5d0*(ET(4)*ET(4)+ET(5)*ET(5)+ET(6)*ET(6))))
      IF(EVM.GT..105) GO TO 200
      EMOD=.97352d0
      GO TO 500
200   IF(EVM.GT..22) GO TO 300
      EMOD=.464d0
      GO TO 500
300   IF(EVM.GT..345) GO TO 400
      EMOD=.296d0
      GO TO 500
400   IF(EVM.GT..48) GO TO 410
      EMOD=.179d0
      GO TO 500
410   EMOD=.047d0
500   CONTINUE
      SHMOD=EMOD/(2.0d0*(1.0d0+POIS))
      DO 510 I=1,NGENS
      DO 510 J=1,NGENS
510   D(I,J)=0.0d0
      D(1,1)=EMOD*(1.0d0 - POIS)/CONS2
      D(2,2)=D(1,1)
      D(3,3)=D(1,1)
      D(1,2)=EMOD*POIS/CONS2
      D(1,3)=D(1,2)
      D(2,1)=D(1,2)
      D(2,3)=D(1,2)
      D(3,1)=D(1,2)
      D(3,2)=D(1,2)
      D(4,4)=SHMOD
      D(5,5)=SHMOD
      D(6,6)=SHMOD
      DO 520 I=1,NGENS
      DO 515 J=1,NGENS
      S(I)=S(I)+D(I,J)*DE(J)
515   CONTINUE
520   CONTINUE
      RETURN
      END