c$Id: elmt40.f,v 1.1 2008/12/12 23:58:47 rlt Exp $
      subroutine elmt40(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

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

c....  Copyright (c) 1984-2011: Sanjay Govindjee
c                               All rights reserved

c-----[--.----+----.----+----.-----------------------------------------]
c     Modification log                                Date (dd/mm/year)
c       Original version                                    02/04/2011
c-----[--.----+----.----+----.-----------------------------------------]
c     Purpose: 
c              Radiation boundary damper
c              2-D Lysmer-Kuhlemeyer theory for Cartessian and axis-symm
c

c     Inputs:
c            Cp -- density times P-wave speed
c            Cs -- density times S-wave speed
c            ax -- axis-symmetry flag (= 0 -- off, = 1 -- on)

c-----[--.----+----.----+----.-----------------------------------------]
      implicit none

      include 'eldata.h'
      include 'iofile.h'

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

      save

      if(isw.eq.0)  then
        if(ior.lt.0) then
         write(*,*) '   Elmt 40: 2-d Lysmer-Kuhlemeyer damper'
        endif

      elseif(isw.eq.1) then

        call inpt40(d)

      elseif(isw.eq.3 .or. isw.eq.6) then

        call stif40(d,ul,xl,s,p,ndf,ndm,nst,isw)

      endif

      end


c-----[--.----+----.----+----.-----------------------------------------]
c
c   Stiffness and RHS routine
c
c-----[--.----+----.----+----.-----------------------------------------]
      subroutine stif40(d,ul,xl,s,p,ndf,ndm,nst,isw)

      implicit none

      include 'cdata.h'
      include 'eldata.h'
      include 'eltran.h'

      integer isw,ndf,ndm,nst
      integer i,j,l,lint,j1,i1

      real*8  d(*),ul(ndf,nen,*),xl(ndm,nen),s(nst,*),p(*)
      real*8  L2,LL,cp,cs,xsj,r,v(2),t(2),mn(2,2)
      real*8  no(2),m(2,2),sg(2,4),shp(2,nel)

      save
 
      data lint /0/

c     Get quadrature points
      if(lint.ne.2) then
        lint = 2
        call int1d(lint,sg)
      endif

c     Extract damper values = velocities*rho 
      cp = d(1)
      cs = d(2)

c     Compute element length
      L2 = 0.d0
      do i = 1,ndm
        L2 = L2 + (xl(i,2) - xl(i,1))**2
      end do
      LL = sqrt(L2)

c     Compute outward normal
      no(1) =  (xl(2,2) - xl(2,1)) / LL
      no(2) = -(xl(1,2) - xl(1,1)) / LL

c     Damping matrix [ cp*n*n' + cs*(iden - n*n') ]
      m(1,1) = cp*no(1)*no(1) + cs*no(2)*no(2)
      m(1,2) = (cp-cs)*no(1)*no(2) 
      m(2,1) = (cp-cs)*no(1)*no(2) 
      m(2,2) = cs*no(1)*no(1) + cp*no(2)*no(2)


      do l = 1,lint
        call shp1d(sg(1,l),xl,shp,ndm,nel,xsj)

        xsj = xsj*sg(2,l)

c       Adjust if axis-symmetric
        if (d(3).eq.1) then
          r = 0.d0
          do j=1,nel
           r = r + xl(1,j)*shp(2,j)
          end do
          xsj = xsj * r
        endif

c       Compute velocities
        v(1) = 0.d0
        v(2) = 0.d0
        do j = 1,nel
          v(1) = v(1) + ul(1,j,4)*shp(2,j)
          v(2) = v(2) + ul(2,j,4)*shp(2,j)
        end do

c       Compute damping tractions
        t(1) = m(1,1)*v(1) + m(1,2)*v(2)
        t(2) = m(2,1)*v(1) + m(2,2)*v(2)


c       Compute RHS
        j1 = 0
        do j = 1,nel
          p(1+j1)     = p(1+j1)     - shp(2,j)*t(1)*xsj
          p(2+j1)     = p(2+j1)     - shp(2,j)*t(2)*xsj
          j1 = j1 +ndf
        end do

c       Compute Damping matrix for implicit
        if(isw.eq.3) then
          i1 = 0
	  xsj = xsj * ctan(2)
          do i = 1,nel

            mn(1,1) = m(1,1)*shp(2,i)*xsj
            mn(1,2) = m(1,2)*shp(2,i)*xsj
            mn(2,1) = m(2,1)*shp(2,i)*xsj
            mn(2,2) = m(2,2)*shp(2,i)*xsj

            j1 = 0
            do j = 1,nel
              s(i+i1,j+j1)     = s(i+i1,j+j1)     + shp(2,j)*mn(1,1)
              s(i+i1,j+j1+1)   = s(i+i1,j+j1+1)   + shp(2,j)*mn(1,2)
              s(i+i1+1,j+j1)   = s(i+i1+1,j+j1)   + shp(2,j)*mn(2,1)
              s(i+i1+1,j+j1+1) = s(i+i1+1,j+j1+1) + shp(2,j)*mn(2,2)
              j1 = j1 + ndf
            end do

            i1 = i1 + ndf
          end do
        end if 

      end do


      end

      subroutine inpt40(d)
c-------------------------------------------------------------------
c
c   Material Properties:
c
c          d(1) = rho*Vp
c          d(2) = rho*Vs
c          d(3) = axis-symmetric flag
c
c-------------------------------------------------------------------

      implicit none

      logical pcomp, errck, tinput
      character*15 text
      real*8 td,d(*)

      include 'iofile.h'

      save

      d(3) = 0
      text = 'start'
      do while (.not.pcomp(text,'    ',4))

          errck = tinput(text,1,td,1)
          
          if (pcomp(text,'cp',2)) then
            d(1) = td
          elseif (pcomp(text,'cs',2)) then
            d(2) = td
          elseif (pcomp(text,'ax',2)) then
            d(3) = 1
          endif

      end do

      if( d(3) .eq. 1 ) then
        if(ior.lt.0) then
         write(*,3010)
        endif
        write(iow,3010)
      endif

      if(ior.lt.0) then
          write(*,3000) d(1),d(2)
      endif
      write(iow,3000) d(1),d(2)


3000  format(/5x,'2 - D   Lysmer-Kuhlemeyer Damper'//
     &       10x,'Cp = ',e12.5/
     &       10x,'Cs = ',e12.5/x)

3010  format(10x,'AXIS-SYMMETRIC'//)

      end