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   spherical components.
c
c NOTE:
c   -1- use postcodes -1 through -6 for the user variables.
c   -2- define the sphere origin and north pole in subroutine 
c       "orient_plot". check the description of "orient_plotv"
c       below for the definition of the spherical 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 spherical 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     user subroutine for input of orientation for anisotropic option.
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 spherical for 3D problems
c  x1(ijk)   = coordinates of center point of sphere
c  x2(ijk)   = coordinates of north pole of sphere
c  xip(ijk)  = integration point coordinates
c  gg(1,ijk) = 1st preferred direction is radial
c  gg(2,ijk) = 2nd preferred direction is circumferential
c  gg(3,ijk) = 3rd preferred direction is normal to 1 and 2
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 center point of sphere
      x1(1) = 0d0
      x1(2) = 0d0
      x1(3) = 0d0
c
c define the north pole of sphere
      x2(1) = 0d0
      x2(2) = 0d0
      x2(3) = 1d0
c
c define some epsilon value
      epsilon = 1d-6
c
c use el. integr. point coordinates to define the transformations
      la1 = icrxpt + (n-1)*nelstr + (nn-1)*ncrd
      do ijk=1,ncrd
       xip(ijk) = vars(la1+ijk-1)
      enddo
      if (ncrd.eq.2) xip(3) = 0d0
c
c 1st pref. direction is radially from center to integr. point
      anorm = 0d0
      do ijk=1,ncrd
       gg(1,ijk) = xip(ijk) - x1(ijk)
       anorm = anorm + gg(1,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      if (anorm.lt.epsilon) then
c if integr. point and center point coincide don't transform
       gg(1,1) = 1d0
       gg(1,2) = 0d0
       gg(1,3) = 0d0
       gg(2,1) = 0d0
       gg(2,2) = 1d0
       gg(2,3) = 0d0
       gg(3,1) = 0d0
       gg(3,2) = 0d0
       gg(3,3) = 1d0
       return
      endif
      do ijk=1,ncrd
       gg(1,ijk) = gg(1,ijk)/anorm
      enddo
c
c sphere axis is directed from center to north pole and stored in x1
      anorm = 0d0
      do ijk=1,ncrd
       x1(ijk) = x2(ijk) - x1(ijk)
       anorm = anorm + x1(ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       x1(ijk) = x1(ijk)/anorm
      enddo
c
      anorm = 0d0
      do ijk=1,ncrd
       anorm = anorm + x1(ijk)*gg(1,ijk)
      enddo
      anorm = dabs(anorm)
      if (anorm.gt.1d0-epsilon) then
c
c if 1st dir. and polar axis dir. coincide then choose some 2nd dir.
       if (dabs(gg(1,1)).lt.epsilon) then
        gg(2,1) =  0d0
        gg(2,2) = -gg(1,3)
        gg(2,3) =  gg(1,2)
       else
        gg(2,1) =  gg(1,2)
        gg(2,2) = -gg(1,1)
        gg(2,3) =  0d0
       endif
      else
c
c 2nd pref. direction is circumferential around sphere axis
       gg(2,1) = x1(2)*gg(1,3) - x1(3)*gg(1,2)
       gg(2,2) = x1(3)*gg(1,1) - x1(1)*gg(1,3)
       gg(2,3) = x1(1)*gg(1,2) - x1(2)*gg(1,1)
      endif
      anorm = 0d0
      do ijk=1,ncrd
       anorm = anorm + gg(2,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       gg(2,ijk) = gg(2,ijk)/anrom
      enddo
c
c 3rd pref. direction is right handed normal to the other two
      gg(3,1) = gg(1,2)*gg(2,3) - gg(1,3)*gg(2,2)
      gg(3,2) = gg(1,3)*gg(2,1) - gg(1,1)*gg(2,3)
      gg(3,3) = gg(1,1)*gg(2,2) - gg(1,2)*gg(2,1)
c
      return
      end