subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,
     * nshear,jpltcd)
c
c CONTACT: ap@marc.de
c      
c* * * * * *
c
c PURPOSE:
c   transform the global stress tensor components to
c   cylindrical components.
c
c NOTE:
c   -1- use postcodes -1 through -6 for the user variables.
c   -2- define the cylinder axis in subroutine "orient_plot".
c       check the description of "orient_plotv" below for the
c       definition of the cylindrical coordinate system.
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
c* * * * * *
      implicit real*8 (a-h,o-z)                                               dp
      dimension s(1),etot(1),eplas(1),ecreep(1),sp(1)
c
      dimension q1(3,3),q2(3,3),s1(3,3),st(3,3),t1(6)
      save t1
c
      if (jpltcd.eq.1) then
c
c compute cylindrical transformation matrix in q1
       call orient_plotv(m,nn,layer,q1)
c
       do i1=1,3
        do i2=1,3
c
c intialize stress tensor
         s1(i1,i2) = 0d0
         st(i1,i2) = 0d0
        enddo
       enddo
c
c store stress components in stress tensor
       do i1=1,ndi
        s1(i1,i1) = s(i1)
       enddo
       if (nshear.eq.1) then
        s1(1,2) = s(ndi+1)
        s1(2,1) = s1(1,2)
       endif
       if (nshear.eq.2) then
        s1(1,2) = s(ndi+1)
        s1(2,1) = s1(1,2)
        s1(2,3) = s(ndi+2)
        s1(3,2) = s1(2,3)
       endif
       if (nshear.eq.3) then
        s1(1,2) = s(ndi+1)
        s1(2,1) = s1(1,2)
        s1(2,3) = s(ndi+2)
        s1(3,2) = s1(2,3)
        s1(3,1) = s(ndi+3)
        s1(1,3) = s1(3,1)
       endif
c
c compute q x stress x q_transposed: s1 = q1xs1xq1_transposed
       do i1=1,3
        do i2=1,3
         st(i1,i2) = 0d0
         do i3=1,3
          st(i1,i2) = st(i1,i2) + q1(i1,i3)*s1(i3,i2)
         enddo
        enddo
       enddo
       do i1=1,3
        do i2=1,3
         s1(i1,i2) = 0d0
         do i3=1,3
          s1(i1,i2) = s1(i1,i2) + st(i1,i3)*q1(i2,i3)
         enddo
        enddo
       enddo
c
c store transformed stresses into vector t1
c t1(1) = s_11; t1(2) = s_22; t1(3) = s_33
c t1(4) = s_12; t1(5) = s_23; t1(6) = s_31
       t1(1) = s1(1,1)
       t1(2) = s1(2,2)
       t1(3) = s1(3,3)
       t1(4) = s1(1,2)
       t1(5) = s1(2,3)
       t1(6) = s1(3,1)
      endif
c
c assign the desired transformed stress component to v
      v = t1(jpltcd)
c
      return
      end
      subroutine orient_plotv(mu,nu,kcu,gg)
c* * * * * *
c
c     mu           element number
c     nu           integration point number
c     kcu          layer number
c     gg           orientation matrix
c
c* * * * * *
c
c  preferred system is cylindrical for 2D and 3D problems
c  used to transform stress tensor to cylindrical system
c  i1 runs from 1 to 3 in definitions of x1,x2,xip and gg
c  x1(i1)   = coordinates of 1st point of cylinder axis
c  x2(i1)   = coordinates of 2nd point of cylinder axis
c  xip(i1)  = integration point coordinates
c  gg(1,i1) = 1st preferred direction is radial
c  gg(2,i1) = 2nd preferred direction is circumferential
c  gg(3,i1) = 3rd preferred direction is axial
c
c* * * * * *
c
      implicit real*8 (a-h,o-z)                                               dp
      dimension gg(3,3),m(2),xip(12),x1(3),x2(3)
c
      include '../common/space'
      include '../common/heat'
      include '../common/lass'
      include '../common/dimen'
      include '../common/arrays'
      include '../common/array4'
c
c define two points on cylinder axis
      x1(1) = 0d0
      x1(2) = 0d0
      x1(3) = 0d0
      x2(1) = 0d0
      x2(2) = 0d0
      x2(3) = 1d0
c
c use el. integr. point coordinates to define the transformations
      la1 = icrxpt + (n-1)*nelstr + (nn-1)*ncrd
      do i1=1,ncrd
       xip(i1) = vars(la1+i1-1)
      enddo
      if (ncrd.eq.2) xip(3) = 0d0
c
c 3rd pref. direction is along cylinder axis
      anorm = 0d0
      do i1=1,3
       gg(3,i1) = x2(i1) - x1(i1)
       anorm = anorm + gg(3,i1)**2
      enddo
      anorm = dsqrt(anorm)
      do i1=1,3
       gg(3,i1) = gg(3,i1)/anorm
      enddo
c
c 2nd pref. direction is circumferential
      gg(2,1) = gg(3,2)*(xip(3)-x1(3)) - gg(3,3)*(xip(2)-x1(2))
      gg(2,2) = gg(3,3)*(xip(1)-x1(1)) - gg(3,1)*(xip(3)-x1(3))
      gg(2,3) = gg(3,1)*(xip(2)-x1(2)) - gg(3,2)*(xip(1)-x1(1))
      anorm = 0d0
      do i1=1,3
       anorm = anorm + gg(2,i1)**2
      enddo
      anorm = dsqrt(anorm)
      do i1=1,3
       gg(2,i1) = gg(2,i1)/anorm
      enddo
c
c 1st pref. direction is radial
      gg(1,1) = gg(2,2)*gg(3,3) - gg(2,3)*gg(3,2)
      gg(1,2) = gg(2,3)*gg(3,1) - gg(2,1)*gg(3,3)
      gg(1,3) = gg(2,1)*gg(3,2) - gg(2,2)*gg(3,1)
c
      return
      end