subroutine orient(mu,nu,kcu,gg)
c
c CONTACT: ap@marc.de
c      
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 cylindrical for 3D problems
c  x1(ijk)   = coordinates of 1st point of cylinder axis
c  x2(ijk)   = coordinates of 2nd point of cylinder axis
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 axial
c
c USAGE
c  Two points that define the cylinder axis must be
c  defined below in x1 and x2.
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 ijk=1,ncrd
       xip(ijk) = vars(la1+ijk-1)
      enddo
c
c 3rd pref. direction is along cylinder axis
      anorm = 0d0
      do ijk=1,ncrd
       gg(3,ijk) = x2(ijk) - x1(ijk)
       anorm = anorm + gg(3,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       gg(3,ijk) = gg(3,ijk)/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 ijk=1,ncrd
       anorm = anorm + gg(2,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       gg(2,ijk) = gg(2,ijk)/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