c$Id:$
      subroutine frame3d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw)

c      * * F E A P * * A Finite Element Analysis Program

c....  Copyright (c) 1984-2005: Robert L. Taylor
c                               All rights reserved

c     Three dimensional frame element

c     Control data:
c         ndm - 3 (x,y,z)
c         ndf - 6 (u,v,w, theta_x,theta_y,theta_z)
c         nen - 3 or more (see below)

c      Beam end nodes 1 and 2
c      Plane defined by nodes 1, 2, 3 contains z-axis
c                            (perpendicular to x-axis)
c      Vector products: e_1 = (x_2 - x_1)/|x_2 - x_1|
c                       e_2 = (x_2 - x_1) x ( x_3 - x_1)
c                       e_3 =  e_2 x e_1


c                  z (e_3)         x (e_1)
c                3 o - - - - - - /
c                  |           o 2
c                  |         /
c                  |       /
c                  |     /
c                  |   /
c                  | /
c                  / ------------- y (e_2)
c              1 o

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

      implicit  none

      include  'cdata.h'
      include  'eldata.h'
      include  'hdata.h'
      include  'iofile.h'
      include  'pmod2d.h'

      include  'comblk.h'

      integer   ndf,ndm,nst,isw, tdof, i

      integer   ix(*)
      real*8    d(*),xl(ndm,*),ul(ndf,nen,*),s(nst,*),p(nst)

      save

c     Output descriptor

      if(isw.eq.0 .and. ior.lt.0) then

        write(*,*) '   Frame3d: 3-d Frame Elastic'
        return

      else

        if(isw.eq.1) then
          write(iow,2000)
          if(ior.lt.0) write(*,2000)
          call inmate(d,tdof,12,3)
          call pltln2(iel)

          if(ndf.lt.6.or.ndm.ne.3) then
            write(iow,2001) ndf,ndm
            call plstop()
          end if

c         Deactivate dof in element for dof > 6

          do i = 7,ndf
            ix(i) = 0
          end do

c         Check for reference node definition

          if(int(d(96)).lt.1 .or. int(d(96)).gt.2) then
            write(iow,2002) int(d(96))
            call plstop()
          endif

        endif

        if(d(18).gt.0.0d0) then
          call frams3d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw)
        else
          write(*,*) ' *ERROR* No 3-d finite deformation frame'
          call plstop()
        endif

      endif

c     Format statements

2000  format(5x,'T h r e e   D i m e n s i o n a l   F r a m e'/)

2001  format(/5x,'Error in FRAME3D. Any of following values is wrong:'
     &       /6x,'ndf(should be 6) =',i3/6x,'ndm (should be 3) =',i3)

2002  format(/5x,'Error in FRAME3D. No reference vector definition for'
     &          ,' axes: lref =',i3)

      end

      subroutine frams3d(d,ul,xl,ix,s,r,ndf,ndm,nst,isw)

c-----[--.----+----.----+----.-----------------------------------------]
c     Small Deformation Three dimensional frame element

c     Control data:
c         ndm - 3 (x,y,z)
c         ndf - 6 (u,v,w, theta_x,theta_y,theta_z)
c         nen - 3 or more (see below)

c      Beam end nodes 1 and 2
c      Plane defined by nodes 1, 2, 3 contains z-axis
c                            (perpendicular to x-axis)
c      Vector products: e_1 = (x_2 - x_1)/|x_2 - x_1|
c                       e_2 = (x_2 - x_1) x ( x_3 - x_1)
c                       e_3 =  e_2 x e_1


c                  z (e_3)         x (e_1)
c                3 o - - - - - - /
c                  |           o 2
c                  |         /
c                  |       /
c                  |     /
c                  |   /
c                  | /
c                  / ------------- y (e_2)
c              1 o
c-----[--.----+----.----+----.-----------------------------------------]
      implicit  none

      include  'bdata.h'
      include  'cdata.h'
      include  'eldata.h'
      include  'eltran.h'
      include  'iofile.h'
      include  'prld1.h'
      include  'prstrs.h'
      include  'refnd.h'

      include  'comblk.h'

      integer   ndf,ndm,nst,isw, i,i1,i2, j,j1,j2, k,l
      real*8    le,dl,dl2,dl3,dn,ctan1,ctan3

      integer   ix(*)
      real*8    d(*),xl(ndm,*),ul(ndf,nen,*),s(nst,*),r(ndf,*)
      real*8    b(6),t(3,3),fi(6,2),sm(12,12),pm(12)

      save

      data      b / 6*0.0d0 /

c     Transfer to correct processor

      go to (1,3,3,3,3,3,7,3), isw

1     return

c     Compute direction cosine terms and member length

3     lref    = d(96)
      refx(1) = d(97)
      refx(2) = d(98)
      refx(3) = d(99)

      t(1,1) = xl(1,2) - xl(1,1)
      t(1,2) = xl(2,2) - xl(2,1)
      t(1,3) = xl(3,2) - xl(3,1)
      le  = sqrt(t(1,1)*t(1,1)+t(1,2)*t(1,2)+t(1,3)*t(1,3))
      dl  = 1.0d0/le
      dl2 = dl*dl
      dl3 = dl*dl2
      t(1,1) = t(1,1)*dl
      t(1,2) = t(1,2)*dl
      t(1,3) = t(1,3)*dl
      if(nel.eq.2) then
        if    (lref.eq.1) then
          t(3,1) = refx(1) - xl(1,1)
          t(3,2) = refx(2) - xl(2,1)
          t(3,3) = refx(3) - xl(3,1)
        elseif(lref.eq.2) then
          t(3,1) = refx(1)
          t(3,2) = refx(2)
          t(3,3) = refx(3)
        else
          write(iow,3000)
          if(ior.lt.0) then
            write(*,3000)
          endif
          call plstop()
        endif
      else
        t(3,1) = xl(1,3) - xl(1,1)
        t(3,2) = xl(2,3) - xl(2,1)
        t(3,3) = xl(3,3) - xl(3,1)
      endif

      if(isw.eq.2) return

      t(2,1) = (t(3,2)*t(1,3) - t(3,3)*t(1,2))
      t(2,2) = (t(3,3)*t(1,1) - t(3,1)*t(1,3))
      t(2,3) = (t(3,1)*t(1,2) - t(3,2)*t(1,1))
      dn  = 1.0d0/sqrt(t(2,1)*t(2,1)+t(2,2)*t(2,2)+t(2,3)*t(2,3))
      t(2,1) = t(2,1)*dn
      t(2,2) = t(2,2)*dn
      t(2,3) = t(2,3)*dn
      t(3,1) = t(1,2)*t(2,3) - t(1,3)*t(2,2)
      t(3,2) = t(1,3)*t(2,1) - t(1,1)*t(2,3)
      t(3,3) = t(1,1)*t(2,2) - t(1,2)*t(2,1)

c     Branch for mass matrix

      if(isw.eq.5) go to 5

c     Compute axial stiffness terms

      i2       =  ndf + 1
      s( 1, 1) =  d(21)*d(32)*dl
      s(i2, 1) = -s(1,1)
      s( 1,i2) = -s(1,1)
      s(i2,i2) =  s(1,1)

c     Compute torsional stiffness terms

      i2       =  ndf + 4
      s( 4, 4) =  d(21)*d(36)*dl/(1.d0 + d(2))*0.5d0
      s(i2, 4) = -s(4,4)
      s( 4,i2) = -s(4,4)
      s(i2,i2) =  s(4,4)

c     Compute bending stiffness terms for z-displacements

      i1       = ndf + 3
      i2       = ndf + 5

      s( 3, 3) =  12.0d0*d(21)*d(33)*dl3
      s( 3,i1) = -s(3,3)
      s(i1, 3) = -s(3,3)
      s(i1,i1) =  s(3,3)

      s( 5,i2) = 2.d0*d(21)*d(33)*dl
      s( 5, 5) =  s(5,i2) + s(5,i2)
      s(i2, 5) =  s(5,i2)
      s(i2,i2) =  s(5,5)

      s( 3, 5) = -6.d0*d(21)*d(33)*dl2
      s( 5, 3) =  s(3,5)
      s( 3,i2) =  s(3,5)
      s(i2, 3) =  s(3,i2)
      s( 5,i1) = -s(3,5)
      s(i1, 5) =  s(5,i1)
      s(i1,i2) = -s(3,5)
      s(i2,i1) =  s(i1,i2)

c     Compute bending stiffness terms for y-displacement

      i1       = ndf + 2
      i2       = ndf + 6

      s( 2, 2) =  12.d0*d(21)*d(34)*dl3
      s( 2,i1) = -s(2,2)
      s(i1, 2) = -s(2,2)
      s(i1,i1) =  s(2,2)

      s( 6,i2) = 2.d0*d(21)*d(34)*dl
      s( 6, 6) =  s(6,i2) + s(6,i2)
      s(i2, 6) =  s(6,i2)
      s(i2,i2) =  s(6,6)

      s( 2, 6) =  6.d0*d(21)*d(34)*dl2
      s( 6, 2) =  s(2,6)
      s( 2,i2) =  s(2,6)
      s(i2, 2) =  s(2,i2)
      s( 6,i1) = -s(2,6)
      s(i1, 6) =  s(6,i1)
      s(i1,i2) = -s(2,6)
      s(i2,i1) =  s(i1,i2)

c     Compute bending stiffness terms for yz-displacement

      if(d(35).ne.0.0d0) then

        i1       = ndf + 3
        i2       = ndf + 5
        j1       = ndf + 2
        j2       = ndf + 6

        s( 2, 3) =  12.d0*d(21)*d(35)*dl3
        s( 3, 2) =  s(2,3)

        s( 2,i1) = -s(2,3)
        s(i1, 2) = -s(2,3)

        s(j1, 3) = -s(2,3)
        s( 3,j1) = -s(2,3)

        s(j1,i1) =  s(2,3)
        s(i1,j1) =  s(2,3)

        s( 6,i2) = 2.d0*d(21)*d(35)*dl
        s(i2, 6) =  s(6,i2)

        s( 6, 5) =  s(6,i2) + s(6,i2)
        s( 5, 6) =  s(6,5)

        s(j2, 5) =  s(6,i2)
        s( 5,j2) =  s(6,i2)

        s(j2,i2) =  s(6,5)
        s(i2,j2) =  s(6,5)

        s( 2, 5) =  6.d0*d(21)*d(35)*dl2
        s( 5, 2) =  s(2,5)

        s( 6, 3) =  s(2,5)
        s( 3, 6) =  s(2,5)

        s( 2,i2) =  s(2,5)
        s(i2, 2) =  s(2,5)

        s(j2, 3) =  s(2,i2)
        s( 3,j2) =  s(2,i2)

        s( 6,i1) = -s(2,5)
        s(i1, 6) = -s(2,5)

        s(j1, 5) =  s(6,i1)
        s( 5,j1) =  s(6,i1)

        s(j1,i2) = -s(2,5)
        s(i2,j1) =  s(j1,i2)

        s(i1,j2) = -s(2,5)
        s(j2,i1) =  s(j1,i2)

      endif

c     Transform to global coordinate displacements

      if(isw.eq.3 .or. isw.eq.6) then
        call tran13(s,t,nst,ndf)

        do i = 1,6

c         Do deformation state

          do j = 1,6
            r(i,1) = r(i,1)
     &             - s(i    ,j    )*(ul(j,1,1) + d(78)*ul(j,1,4))
     &             - s(i    ,j+ndf)*(ul(j,2,1) + d(78)*ul(j,2,4))
            r(i,2) = r(i,2)
     &             - s(i+ndf,j    )*(ul(j,1,1) + d(78)*ul(j,1,4))
     &             - s(i+ndf,j+ndf)*(ul(j,2,1) + d(78)*ul(j,2,4))
          end do ! j
        end do ! i

        ctan1 = ctan(1) + d(78)*ctan(2)
        do j = 1,nst
          do i = 1,nst
            s(i,j) = s(i,j)*ctan1
          end do ! i
        end do ! j

c       Inertial and body force effects

        dn = d(4)*d(32)*le*0.5d0

c       Set body loading factors

        if(int(d(74)).gt.0) then
          b(1) = 0.5*le*(d(11) + prldv(int(d(74)))*d(71))
        else
          b(1) = 0.5*le*d(11)*dm
        endif
        if(int(d(75)).gt.0) then
          b(2) = 0.5*le*(d(12) + prldv(int(d(75)))*d(72))
        else
          b(2) = 0.5*le*d(12)*dm
        endif
        if(int(d(76)).gt.0) then
          b(3) = 0.5*le*(d(13) + prldv(int(d(76)))*d(73))
        else
          b(3) = 0.5*le*d(13)*dm
        endif

c       Lumped mass computation

        ctan3 = ctan(3) + d(77)*ctan(2)
        if(d(7).eq.0.0d0) then

          do i = 1,3
            r(i,1)         = r(i,1) - dn*(ul(i,1,5) + d(77)*ul(i,1,4))
     &                     + b(i)
            r(i,2)         = r(i,2) - dn*(ul(i,2,5) + d(77)*ul(i,2,4))
     &                     + b(i)
            s(i    ,i    ) = s(i    ,i    ) + ctan3*dn
            s(i+ndf,i+ndf) = s(i+ndf,i+ndf) + ctan3*dn
          end do

c       Consistent mass computation

        else
          do i = 1,12
            do j = 1,12
              sm(j,i) = 0.0d0
            end do ! i
            pm(i) = 0.0d0
          end do ! i
          call massf3(sm,pm,d(7),d,le,12,ndm,6)
          call tran13(sm,t,12,6)
          i1 = 0
          i2 = 0
          do i = 1,2
            do k = 1,6
              r(k,i) = r(k,i) + b(k)
              j1 = 0
              j2 = 0
              do j = 1,2
                do l = 1,6
                  r(k,i)  = r(k,i)
     &                    - sm(k+i2,l+j2)*(ul(l,j,5) + d(77)*ul(l,j,4))
                  s(k+i1,l+j1) = s(k+i1,l+j1) + sm(k+i2,l+j2)*ctan3
                end do
                j1 = j1 + ndf
                j2 = j2 + 6
              end do
            end do
            i1 = i1 + ndf
            i2 = i2 + 6
          end do

        endif

      endif

c     Output member forces

      if(isw.eq.4 .or. isw.eq.8) then

c       Transform displacements

        do i = 1,ndf
          ul(i,1,2) = ul(i,1,1)
          ul(i,2,2) = ul(i,2,1)
        end do
        do i = 1,3
          ul(i  ,1,1) = t(i,1)*ul(1,1,2)
     &                + t(i,2)*ul(2,1,2)
     &                + t(i,3)*ul(3,1,2)
          ul(i+3,1,1) = t(i,1)*ul(4,1,2)
     &                + t(i,2)*ul(5,1,2)
     &                + t(i,3)*ul(6,1,2)
          ul(i  ,2,1) = t(i,1)*ul(1,2,2)
     &                + t(i,2)*ul(2,2,2)
     &                + t(i,3)*ul(3,2,2)
          ul(i+3,2,1) = t(i,1)*ul(4,2,2)
     &                + t(i,2)*ul(5,2,2)
     &                + t(i,3)*ul(6,2,2)
        end do ! i

c       Compute member force if necessary

        do i = 1,6
          r(i,1) = 0.0d0
          r(i,2) = 0.0d0
          do j = 1,6
            r(i,1) = r(i,1) - s(i    ,j    )*ul(j,1,1)
     &                      - s(i    ,j+ndf)*ul(j,2,1)
            r(i,2) = r(i,2) - s(i+ndf,j    )*ul(j,1,1)
     &                      - s(i+ndf,j+ndf)*ul(j,2,1)
          end do ! j
        end do ! i

        do i = 1,6
          fi(i,1) =  r(i,1)
          fi(i,2) = -r(i,2)
        end do

        if(isw.eq.4) then
          mct = mct - 1
          if(mct.le.0) then
            write(iow,2001) o,head
            if(ior.lt.0) then
              write(*,2001) o,head
            endif
            mct = 50
          endif

          write(iow,2002) n,ma,(xl(i,1),i=1,3),(fi(i,1),i=1,6),
     &                         (xl(i,2),i=1,3),(fi(i,2),i=1,6)
          if(ior.lt.0) then
            write(*,2002) n,ma,(xl(i,1),i=1,3),(fi(i,1),i=1,6),
     &                         (xl(i,2),i=1,3),(fi(i,2),i=1,6)
          endif

c       Project stress resultants to nodes

        else

          call frcn3d(ix,fi,hr(nph),hr(nph+numnp),numnp)

        endif

      endif

      return

c     Compute mass matrix (lumped only)

5     if(d(7).eq.0.0d0) then
        dn = d(4)*d(32)*le*0.5d0
        do i = 1,3
          r(i,1)         = dn
          r(i,2)         = dn
          s(i    ,i    ) = dn
          s(i+ndf,i+ndf) = dn
        end do

      else
        call massf3(s,r,d(7),d,le,nst,ndm,ndf)
        call tran13(s,t,nst,ndf)
      endif

7     return

c     Format statements

2001  format(a1,20a4//5x,'3-D Frame Element Forces'//
     &   '    Elmt  Mat     x-Coor     y-Coor     z-Coor'/
     & 7x,'I-end:      Force    1-Shear    2-Shear   1-Torque',
     &     '   1-Moment   2-Moment'/
     &   '                  x-Coor     y-Coor     z-Coor'/
     & 7x,'J-end:      Force    1-Shear    2-Shear   1-Torque',
     &     '   1-Moment   2-Moment'/1x,78('-'))

2002  format(i8,i5,1p,3e11.3/13x,1p,6e11.3/
     +       13x,  1p,3e11.3/13x,1p,6e11.3/1x)

3000  format('    *ERROR*  No Reference point for element')

      end

      subroutine tran13(s,t,nst,ndf)

      implicit  none

      integer   nst,ndf, i,j,j1,j2,k, nsiz

      real*8    s(nst,nst),ss(6),t(3,3)

      save

c     Transform if not identity

      if(t(1,1)+t(2,2)+t(3,3).lt.2.999999d+00) then

c     Postmultiply local stiffness by transformation array

        j1 = 0
        nsiz = ndf + ndf
        do k = 1,2
          do i = 1,nsiz
            do j = 1,6
              ss(j) = s(i,j+j1)
            end do
            j2 = j1 + 3
            do j = 1,3
              s(i,j+j1) = ss(1)*t(1,j) + ss(2)*t(2,j) + ss(3)*t(3,j)
              s(i,j+j2) = ss(4)*t(1,j) + ss(5)*t(2,j) + ss(6)*t(3,j)
            end do
          end do
          j1 = j1 + ndf
        end do

c     Premultiply result by transpose of transformation array

        j1 = 0
        do k = 1,2
          do i = 1,nsiz
            do j = 1,6
              ss(j) = s(j+j1,i)
            end do
            j2 = j1 + 3
            do j = 1,3
              s(j+j1,i) = t(1,j)*ss(1) + t(2,j)*ss(2) + t(3,j)*ss(3)
              s(j+j2,i) = t(1,j)*ss(4) + t(2,j)*ss(5) + t(3,j)*ss(6)
            end do
          end do
          j1 = j1 + ndf
        end do

      endif

      end

      subroutine frcn3d(ix,fi,dt,st,numnp)

      implicit  none

      integer   numnp
      integer   i,j,ll

      integer   ix(*)
      real*8    dt(numnp),st(numnp,*),fi(6,*)

      save

      do i = 1,2

        ll = ix(i)
        if(ll.gt.0) then

          dt(ll) = dt(ll) + 1.d0

c         Stress projections

          do j = 1,6
            st(ll,j) = st(ll,j) + fi(j,i)
          end do
        endif
      end do

      end

      subroutine massf3(s,r,cfac,d,le,nst,ndm,ndf)

c     Frame mass matrix

      implicit  none

      integer   nst,ndm,ndf,i,j,k,l, i1,j1, ll
      real*8    cfac,lfac,ra,le,lr,t,dv,s1,s2,s3
      real*8    d(*), r(ndf,*),s(nst,nst),sg(4),ag(4),bb(2,2),db(2,2)
      real*8    nt(2,6,2),nr(3,6,2),ntak(2,6),nrik(3,6),in(3,3)

      save


c     Set inertial properties

      ra       = d(4)*d(32)
      in(1,2)  = 0.0d0
      in(1,3)  = 0.0d0
      in(2,1)  = 0.0d0
      in(3,1)  = 0.0d0

c     Lumped mass matrix

      lr       = 1.d0/le
      t        = 0.5d0*ra*le
      do i = 1,ndm
        r(i,1) = t
        r(i,2) = t
      end do
      do i = 1,6
        do j = 1,2
          nt(1,i,j) = 0.0d0
          nt(2,i,j) = 0.0d0
          nr(1,i,j) = 0.0d0
          nr(2,i,j) = 0.0d0
          nr(3,i,j) = 0.0d0
        end do
      end do

c     Consistent mass matrix

      s(1    ,1    ) = 0.6666666666666667d0*t
      s(1    ,ndf+1) = 0.3333333333333333d0*t
      s(ndf+1,1    ) = s(1,ndf+1)
      s(ndf+1,ndf+1) = s(1,1)

      call int1d(4,sg,ag)

      do ll = 1,4

        dv     = t*ag(ll)
        in(1,1)  = d(4)*d(36)*0.5d0*le*ag(ll)
        in(2,2)  = d(4)*d(33)*0.5d0*le*ag(ll)
        in(2,3)  = d(4)*d(35)*0.5d0*le*ag(ll)
        in(3,2)  = d(4)*d(35)*0.5d0*le*ag(ll)
        in(3,3)  = d(4)*d(34)*0.5d0*le*ag(ll)

        s1     = 0.5d0 + 0.5d0*sg(ll)
        s2     = s1*s1
        s3     = s1*s2

        bb(1,2)  = 3.d0*s2 - s3 - s3
        bb(2,2)  = le*(s3 - s2)
        bb(1,1)  = 1.d0 - bb(1,2)
        bb(2,1)  = le*(s1 - s2) + bb(2,2)

        db(1,2)  = 6.d0*(s1 - s2)*lr
        db(2,2)  = 3.d0*s2 - 2.d0*s1
        db(1,1)  = -db(1,2)
        db(2,1)  = 1.d0 -2.d0*s1 + db(2,2)

        nr(1,4,1) = 1.d0 - s1
        nr(1,4,2) = s1
        do k = 1,2

          nt(1,2,k) =  bb(1,k)
          nt(1,6,k) =  bb(2,k)
          nt(2,3,k) =  bb(1,k)
          nt(2,5,k) = -bb(2,k)

          nr(2,3,k) =  db(1,k)
          nr(2,5,k) = -db(2,k)
          nr(3,2,k) = -db(1,k)
          nr(3,6,k) = -db(2,k)

        end do
        i1 = 0
        do k = 1,2
          ntak(1,2) = dv*nt(1,2,k)
          ntak(1,6) = dv*nt(1,6,k)
          ntak(2,3) = dv*nt(2,3,k)
          ntak(2,5) = dv*nt(2,5,k)
          do i = 1,6
            do j = 1,3
              nrik(j,i) = nr(1,i,k)*in(1,j)
     &                  + nr(2,i,k)*in(2,j)
     &                  + nr(3,i,k)*in(3,j)
            end do ! j
          end do ! i

          j1 = 0
          do l = 1,2

            do i = 1,6
              do j = 1,6
                s(i+i1,j+j1) = s(i+i1,j+j1)
     &                       + ntak(1,i)*nt(1,j,l)
     &                       + ntak(2,i)*nt(2,j,l)
     &                       + nrik(1,i)*nr(1,j,l)
     &                       + nrik(2,i)*nr(2,j,l)
     &                       + nrik(3,i)*nr(3,j,l)
              end do ! j
            end do ! i
            j1 = j1 + ndf
          end do ! l
          i1 = i1 + ndf
        end do ! k
      end do ! ll

c     Interpolate mass between lumped and consistent
c     Consistent Mass: cfac = 1.0 ; Lumped Mass: cfac = 0.0

      lfac = 1.d0 - cfac
      do i = 1,nst
        do j = 1,nst
          s(i,j) = cfac*s(i,j)
        end do
      end do
      do i = 1,ndm
        s(i    ,i    ) = s(i    ,i    ) + lfac*r(i,1)
        s(i+ndf,i+ndf) = s(i+ndf,i+ndf) + lfac*r(i,2)
      end do

      end