c$Id: cubic.f,v 1.1 2002/08/14 23:32:48 sanjay Exp $
        subroutine cubic(isw,f,detf,ud,sig,kk)
c      * * F E A P * * A Finite Element Analysis Program

c....  Copyright (c) 1999: Sanjay Govindjee

c-----[--.----+----.----+----.-----------------------------------------]
c     Purpose: User Constitutive Model Routine for a single
c              crystal cubic material St. Venant Kirchhoff formulation

c     Input:
c          isw = task switch (3,6,4,6) = (stress+tangent,stress,stress,stress)
c            f = deformation gradient
c         detf = Jacobian determinent
c          ud  = user properties [ alpha, beta, gamma, a(3), b(3), c(3) ]

c     Output:
c          sig = Cauchy stresses
c          kk  = spatial tangent

c-----[--.----+----.----+----.-----------------------------------------]



	implicit none

        integer isw, i, j, k, l, m, n
        integer vtot(2,6), ttov(3,3)

        real*8 ud(*)
	real*8 f(3,3),sig(6),kk(6,6),detf

        real*8 alpha, beta, gamma, tmpx
        real*8 aaxis(3), baxis(3), caxis(3)
        real*8 Fa(3),    Fb(3),    Fc(3)
        real*8 bb(6), c(6), e(6), pe(6)
        real*8 eaa,ebb,ecc,tre

        save

        data vtot /1,1, 2,2, 3,3, 1,2, 2,3, 3,1/
        data ttov /1,4,6,4,2,5,6,5,3/

c       Basic Tolerance


        if(isw.eq.3 .or. isw.eq.6 .or. isw.eq.4 .or. 
     &          isw.eq.8 ) then

c       Set up the moduli

        alpha = ud(1)
        beta  = ud(2)
        gamma = ud(3)
        tmpx = alpha - beta - 2.d0*gamma

        aaxis(1) = ud(4)
        aaxis(2) = ud(5)
        aaxis(3) = ud(6)
        baxis(1) = ud(7)
        baxis(2) = ud(8)
        baxis(3) = ud(9)
        caxis(1) = ud(10)
        caxis(2) = ud(11)
        caxis(3) = ud(12)


c       Compute Left-Cauchy Green deformation tensor

        bb(1) = f(1,1)*f(1,1) + f(1,2)*f(1,2) + f(1,3)*f(1,3)
        bb(2) = f(2,1)*f(2,1) + f(2,2)*f(2,2) + f(2,3)*f(2,3)
        bb(3) = f(3,1)*f(3,1) + f(3,2)*f(3,2) + f(3,3)*f(3,3)
        bb(4) = f(1,1)*f(2,1) + f(1,2)*f(2,2) + f(1,3)*f(2,3)
        bb(5) = f(2,1)*f(3,1) + f(2,2)*f(3,2) + f(2,3)*f(3,3)
        bb(6) = f(1,1)*f(3,1) + f(1,2)*f(3,2) + f(1,3)*f(3,3)

c       Compute Right-Cauchy Green deformation tensor

        c(1) = f(1,1)*f(1,1) + f(2,1)*f(2,1) + f(3,1)*f(3,1)
        c(2) = f(1,2)*f(1,2) + f(2,2)*f(2,2) + f(3,2)*f(3,2)
        c(3) = f(1,3)*f(1,3) + f(2,3)*f(2,3) + f(3,3)*f(3,3)
        c(4) = f(1,1)*f(1,2) + f(2,1)*f(2,2) + f(3,1)*f(3,2)
        c(5) = f(1,2)*f(1,3) + f(2,2)*f(2,3) + f(3,2)*f(3,3)
        c(6) = f(1,1)*f(1,3) + f(2,1)*f(2,3) + f(3,1)*f(3,3)

c       Compute Green-Lagrange strain tensor and its push

        e(1) = 0.5d0*(c(1) - 1.d0)
        e(2) = 0.5d0*(c(2) - 1.d0)
        e(3) = 0.5d0*(c(3) - 1.d0)
        e(4) = 0.5d0*c(4)
        e(5) = 0.5d0*c(5)
        e(6) = 0.5d0*c(6)

        tre = e(1)+e(2)+e(3)

        pe(1) = 0.5d0*(bb(1)*bb(1)+bb(4)*bb(4)+bb(6)*bb(6)- bb(1))
        pe(2) = 0.5d0*(bb(2)*bb(2)+bb(4)*bb(4)+bb(5)*bb(5)- bb(2))
        pe(3) = 0.5d0*(bb(3)*bb(3)+bb(5)*bb(5)+bb(6)*bb(6)- bb(3))
        pe(4) = 0.5d0*(bb(1)*bb(4)+bb(4)*bb(2)+bb(6)*bb(5)- bb(4))
        pe(5) = 0.5d0*(bb(6)*bb(4)+bb(5)*bb(2)+bb(3)*bb(5)- bb(5))
        pe(6) = 0.5d0*(bb(1)*bb(6)+bb(4)*bb(5)+bb(6)*bb(3)- bb(6))

c       Cubic axes strains

        eaa = aaxis(1)*aaxis(1)*e(1) + aaxis(2)*aaxis(2)*e(2) +
     &        aaxis(3)*aaxis(3)*e(3) + 2.d0*aaxis(1)*aaxis(2)*e(4) +
     &        2.d0*aaxis(2)*aaxis(3)*e(5) + 2.d0*aaxis(1)*aaxis(3)*e(6)
        ebb = baxis(1)*baxis(1)*e(1) + baxis(2)*baxis(2)*e(2) +
     &        baxis(3)*baxis(3)*e(3) + 2.d0*baxis(1)*baxis(2)*e(4) +
     &        2.d0*baxis(2)*baxis(3)*e(5) + 2.d0*baxis(1)*baxis(3)*e(6)
        ecc = caxis(1)*caxis(1)*e(1) + caxis(2)*caxis(2)*e(2) +
     &        caxis(3)*caxis(3)*e(3) + 2.d0*caxis(1)*caxis(2)*e(4) +
     &        2.d0*caxis(2)*caxis(3)*e(5) + 2.d0*caxis(1)*caxis(3)*e(6)

c       Push cubic axes

        do i = 1,3
          Fa(i) = f(i,1)*aaxis(1) + f(i,2)*aaxis(2) + f(i,3)*aaxis(3)
          Fb(i) = f(i,1)*baxis(1) + f(i,2)*baxis(2) + f(i,3)*baxis(3)
          Fc(i) = f(i,1)*caxis(1) + f(i,2)*caxis(2) + f(i,3)*caxis(3)
        end do

        do i = 1,6
          sig(i) = (beta*tre*bb(i)+2.d0*gamma*pe(i)
     &           + tmpx*(eaa*Fa(vtot(1,i))*Fa(vtot(2,i))
     &           +       ebb*Fb(vtot(1,i))*Fb(vtot(2,i))
     &           +       ecc*Fc(vtot(1,i))*Fc(vtot(2,i))))/detf
        end do

c       Moduli computation

        call pzero(kk,36)

        do m = 1,6
         i = vtot(1,m)
         j = vtot(2,m)
         do n = 1,6
           l = vtot(1,n)
           k = vtot(2,n)
           kk(m,n) = kk(m,n) + (beta*bb(m)*bb(n)
     &             + gamma*(bb(ttov(i,k))*bb(ttov(j,l))
     &             +        bb(ttov(i,l))*bb(ttov(j,k)))
     &             + tmpx*(Fa(i)*Fa(j)*Fa(k)*Fa(l)  
     &             +       Fb(i)*Fb(j)*Fb(k)*Fb(l)  
     &             +       Fc(i)*Fc(j)*Fc(k)*Fc(l)))/detf
         end do
        end do

        else
c        write(*,*) 'ISW ',isw,' not implemented'
        endif


	end ! subroutine