program fem1d
C-----------------------------------------------------------------------
C                                                                       
C*Program: fem1d - 1D FEM harness
C 
C*Author:  Mark A. Christon 
C
C Date  :  09/28/1999                                 
C
C Date     |    Modification
C =========|============================================================
C xx/xx/xx | 
C                  
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
#include "cntlvec.h"
      character *80 cfn(3), exenm
      character *24 dat

C ... set unit numbers
      lum = 5
      lus = 6 
      luo = 7
      lup = 8
      icntl(91) = lum
      icntl(92) = luo
      icntl(93) = lus
      icntl(94) = lup

#ifdef LINUX
      call fdate( dat )
#endif
      if( iargc() .ne. 3 ) then
         write( lus, 9000 )
         stop
      endif
      call getarg( 0, exenm )
      call getarg( 1, cfn(1) )
      call getarg( 2, cfn(2) )
      call getarg( 3, cfn(3) )
      write( lus, 1000 ) exenm, dat
      write( lus, 1010 ) cfn(1)
      write( lus, 1020 ) cfn(2)
      write( lus, 1025 ) cfn(3)

C ... Open all files & setup the parser
C     lum=mesh file, luo=output, lup=plot, lus=screen
      iunit  = 5
      ncolum = 82
      open ( lum, file=cfn(1), form='formatted', status='unknown' )
      open ( luo, file=cfn(2), form='formatted', status='unknown' )
      open ( lup, file=cfn(3), form='formatted', status='unknown' )
      call intprs( iunit, luo, -lus, ncolum )

C ... Read & echo the header
      write( luo, 1000 ) exenm, dat
      write( luo, 1010 ) cfn(1)
      write( luo, 1020 ) cfn(2)
      write( luo, 1025 ) cfn(3)
      call rdhead( lus, luo )

C ... Allocate memory
      call alloct

C ... Read & Echo nodal coordinates, connectivity and nodesets
C     ip1=x
C     ip2=f
C     ip4=np
C     ip5=ndset
      nnp    = icntl( 1)
      nel    = icntl( 2)
      ndim   = icntl( 4)
      nnpe   = icntl( 6)
      nwix   = icntl( 7)
      nndset = icntl( 9)
      nsdset = icntl(10)
      call rdcoor( ra(ip1), nnp, ndim, lum, luo )
      call rdconn( ia(ip4), nnp, nel, nnpe, nwix, lum, luo )
      if (nndset.gt.0) call rwnset( ia(ip5), nndset, lum, luo, lus )

C ... Compute the 1/2-bandwidth & allocate for coefficient matrix
      call bandwd( ia(ip4), nel, nnpe, nwix,  nbw )
      ip5    = locrnd
      locrnd = locrnd+nnp*nbw
      if (locrnd .gt. MAXMEM) then
        write(lus,1000) locrnd,MAXMEM
        stop
      endif

C ... Form K and F
      T=30000.0
      ltype=0
      write( lus, 1030 ) nbw, T, ltype
      write( luo, 1030 ) nbw, T, ltype
      call formkf( T, ltype, ra(ip1), ra(ip5), ra(ip2), ia(ip4), 
     2     nnp, nel, nwix, nbw )

C ... Set EBC's based on nodesets
      call blstbc( ra(ip5), ra(ip2), nnp )

C ... Solve Ku=F
      call udu( ra(ip5), ra(ip2), 1, nnp, nbw, nnp, nbw, lus )

C ... Output the displacements
      call writeu( ra(ip1), ra(ip2), nnp, luo, lup )

C ... Close all files
      close(lum)
      close(luo)
      close(lup)

      stop
 1000 format(//,
     2       5x,'FEM1D       : ',a60,/
     3       5x,'Date & Time : ',a24,/)
 1010 format(5x,'Mesh   file : ',a40)
 1020 format(5x,'Output file : ',a40)
 1025 format(5x,'Plot   file : ',a40)
 1030 format(//,
     2 5x,'1/2-bandwidth= ',i4,/,
     3 5x,'Tension      = ',1pe10.4,/,
     4 5x,'Load Type    = ',i4,//)
 2000 format(//,
     2  5x,'!!!!! Memory Exceeded !!!!!',/,
     3  5x,'During allocation for stiffnes ...',/,
     4  5x,'Requested: ',i10,/,
     5  5x,'Available: ',i10,/,
     6  5x,'!!!!! Memory Exceeded !!!!!',//)
 9000 format(//,5x,'Usage: fem1d mesh output plot',//)
      end

      subroutine alloct
C-----------------------------------------------------------------------
C                                                                       
C* Routine : alloct - allocate for coordinates, connectivity, nodeset
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
#include "cntlvec.h"

      nnp    = icntl( 1)
      nel    = icntl( 2)
      ndim   = icntl( 4)
      nnpe   = icntl( 6)
      nwix   = icntl( 7)
      nndset = icntl( 9)
      nsdset = icntl(10)
      lus    = icntl(93)

C ... Allocate reals from real a-array, int's from int a-array !!!

C ... pointer for x-array from ra()
      locrnd = 1
      ip1    = locrnd
      nw     = nnp*ndim
      locrnd = locrnd+nw
C ... pointer for the solution vector & rhs from ra()
      ip2    = locrnd
      locrnd = locrnd+nnp
      if (locrnd .gt. MAXMEM) then
        write(lus,1000) locrnd,MAXMEM
        stop
      endif

C ... connectivity from ia()
      locind = 1
      ip4    = locind
      nw     = nel*nwix
      locind = locind+nw

C ... nodeset data - part #1
      ip5    = locind
      locind = locind + nndset*3

      if (locind .gt. MAXMEM) then
        write(lus,1000) locind,MAXMEM
        stop
      endif

      return
 1000 format(//,
     2  5x,'!!!!! Memory Exceeded !!!!!',/,
     3  5x,'Requested: ',i10,/,
     4  5x,'Available: ',i10,/,
     5  5x,'!!!!! Memory Exceeded !!!!!',//)
      end

      subroutine rdhead(   luo,   lus)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : rdhead - read the header block
C
C Argument |    Meaning
C =========|============================================================
C luo      | output LUN
C lus      | screen LUN
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
#include "cntlvec.h"
      character *8  cmd
      
C ... Parse the header info 
      nnp  = 0
      nnpe = 0
      nel  = 0
      nmat = 0
      nndset = 0
      nsdset = 0
  100 continue
      call getsym( 16 , 'Mesh  Command : ',  cmd,     8 )
      if( cmd .eq. 'nnpe' ) then
        call getnum( 16 , 'No. Node/El.  : ',    r,  nnpe )
        goto 100
      endif
      if( cmd .eq. 'nel'  ) then
        call getnum( 16 , 'No. Elements  : ',    r,  nel  )
        goto 100
      endif
      if( cmd .eq. 'nnp'  ) then
        call getnum( 16 , 'No. Nodes     : ',    r,  nnp  )
        goto 100
      endif
      if( cmd .eq. 'nmat'  ) then
        call getnum( 16 , 'No. Materials : ',    r,  nmat )
        goto 100
      endif
      if( cmd .eq. 'nnd_sets'  ) then
        call getnum( 16 , 'No. Node sets : ',    r,  nndset )
        goto 100
      endif
      if( cmd .eq. 'nsd_sets'  ) then
        call getnum( 16 , 'No. Node sets : ',    r,  nsdset )
        goto 100
      endif
      if( cmd .ne. 'end' ) goto 100
      ndim = 3
C ... nnpe=8 for Hex8, nnpe=4 for Quad4, nnpe=2 for Oned2 element
      if( nnpe .eq. 8 ) ndim = 3
      if( nnpe .eq. 4 ) ndim = 2
      if( nnpe .eq. 2 ) ndim = 1
      nwix=nnpe+1

      write( lus, 1000 ) ndim, nnp, nel, nmat, nnpe, nndset, nsdset
      write( luo, 1000 ) ndim, nnp, nel, nmat, nnpe, nndset, nsdset

      icntl( 1) = nnp
      icntl( 2) = nel
      icntl( 3) = nmat
      icntl( 4) = ndim
      icntl( 5) = 1
      icntl( 6) = nnpe
      icntl( 7) = nwix
      icntl( 8) = nelblk
      icntl( 9) = nndset
      icntl(10) = nsdset

      return
 1000 format(5x,'Ndim        : ',i10,/,
     2       5x,'Nnp         : ',i10,/,
     3       5x,'Nel         : ',i10,/,
     4       5x,'Nmat        : ',i10,/,
     5       5x,'Nnpe        : ',i10,/,
     6       5x,'Nnd_sets    : ',i10,/,
     7       5x,'Nsd_sets    : ',i10)
      end

      subroutine rdconn(    np,   nnp,   nel,  nnpe,  nwix,  lum,
     2                     luo)
C-----------------------------------------------------------------------
C                                                                       
C* Routine: rdconn - read in the connectivity
C
C Argument |    Meaning
C =========|============================================================
C np       | connectivity array
C nnp      | # of nodes
C nel      | # of elements
C nnpe     | # of nodes per element
C nwix     | size of np
C lum      | mesh   LUN
C luo      | output LUN
C
C* Author:  Mark A. Christon 
C
C Date  :  01/21/00                                   
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      dimension np(nel,nwix)
      character *80 txt

C ... Mesh I/O
      if( nnpe .eq. 2 ) then
C ..... Read 1-D connectivity
        write( luo, 2000 )
        do i = 1, nel
          call gttxsg( txt, nl, lum )
          read ( txt, 1000 ) mat, ( np(i,j), j = 1, nnpe )
          write( luo, 2010 ) mat, ( np(i,j), j = 1, nnpe )
        end do
      else if( nnpe .eq. 4 ) then
C ..... Read 2-D connectivity
        write( luo, 2020 )
        do i = 1, nel
          call gttxsg( txt, nl, lum )
          read ( txt, 1010 ) mat, ( np(i,j), j = 1, nnpe )
          write( luo, 2030 ) mat, ( np(i,j), j = 1, nnpe )
        end do
      else if( nnpe .eq. 8 ) then
C ..... Read 3-D connectivity
        write( luo, 2040 )
        do i = 1, nel
          call gttxsg( txt, nl, lum )
          read ( txt, 1020 ) mat, ( np(i,j), j = 1, nnpe )
          write( luo, 2050 ) mat, ( np(i,j), j = 1, nnpe )
        end do
      endif

      return
 1000 format(i5,2i8)
 1010 format(i5,4i8)
 1020 format(i5,8i8)
 2000 format(//,
     2    5x,' F o u r   N o d e  E l e m e n t s'//
     3    2x,'element material  node1  node2')
 2010 format(11x,i5,2i8)
 2020 format(//,
     2    5x,' F o u r   N o d e  E l e m e n t s'//
     3    2x,'element material  node1  node2  node3  node4')
 2030 format(11x,i5,4i8)
 2040 format(//,
     2    5x,' E i g h t  N o d e  E l e m e n t s'//
     3    2x,'element material  node1  node2  node3',
     4       '  node4  node5  node6  node7  node8')
 2050 format(11x,i5,8i8)
      end

      subroutine rdcoor(     x,   nnp,  ndim,  lum,   luo)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : rdcoor - Read in the coordinates
C
C Argument |    Meaning
C =========|============================================================
C x        | nodal coordinates
C nnp      | # of nodes
C ndim     | # of dimensions
C lum      | mesh   lun
C luo      | output lun
C                                                                       
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      dimension x(nnp,ndim)
      character *80 txt

      if( ndim .eq. 1 ) then
C ..... Read 1-D coordinates
        write( luo, 2000 )
        do i = 1, nnp
          call gttxsg( txt, nl, lum )
          read ( txt, 1000 ) nd, x(nd,1)
          write( luo, 2010 ) nd, x(nd,1)
        end do
      else if( ndim .eq. 2 ) then
C ..... Read 2-D coordinates
        write( luo, 2020 )
        do i = 1, nnp
          call gttxsg( txt, nl, lum )
          read ( txt, 1010 ) nd, x(nd,1), x(nd,2)
          write( luo, 2030 ) nd, x(nd,1), x(nd,2)
        end do
      else if ( ndim .eq. 3 ) then
C ..... Read 3-D coordinates
        write( luo, 2040 )
        do i = 1, nnp
          call gttxsg( txt, nl, lum )
          read ( txt, 1020 ) nd, x(nd,1), x(nd,2), x(nd,3)
          write( luo, 2040 ) nd, x(nd,1), x(nd,2), x(nd,3)
        end do
      endif

      return
 1000 format(i8,  e20.0 )
 1010 format(i8,2(e20.0))
 1020 format(i8,3(e20.0))
 2000 format(//,
     2 5x,'N o d a l  P o i n t  C o o r d i n a t e s '//,
     3 8x,'nd',9x,'x-ord'/)
 2010 format(1x,i9,e15.4)
 2020 format(//,
     2 5x,'N o d a l  P o i n t  C o o r d i n a t e s '//,
     3 8x,'nd',9x,'x-ord',11x,'y-ord',/)
 2030 format(1x,i9,2e15.4)
 2040 format(//,
     2 5x,'N o d a l  P o i n t  C o o r d i n a t e s '//,
     3 8x,'nd',9x,'x-ord',11x,'y-ord',11x,'z-ord',/)
 2050 format(1x,i9,3e15.4)
      end

      subroutine rwnset( ndset,nndset,   lum,   luo,   lus)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : rwnset - read/write nodeset lists
C
C Argument |    Meaning
C =========|============================================================
C ndset    | array containing nodeset id, # of nodes, node list pointer
c nndset   | # of nodesets
C lum      | mesh   LUN
C luo      | output LUN
C lus      | screen LUN
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
#include "cntlvec.h"
      dimension ndset(nndset,3)
      character *80 txt

C ... Part 1 of Nodeset data
      call gttxsg( txt, nl, lum )
      read( txt, 1000 ) ijunk
      if( ijunk .ne. nndset ) then
        write( lus, 1010 ) nndset, ijunk
        stop
      endif

C ... Part 2 of Nodeset data
      do i=1,nndset
        call gttxsg( txt, nl, lum )
        read( txt, 1020 ) ndset(i,1),ndset(i,2)
      end do

C ... Part 3 of Nodeset data
      do i=1,nndset
C ..... allocate from ia()
        ndset(i,3) = locind
        locind = locind + ndset(i,2)
        if (locind .gt. MAXMEM ) then
          write( lus, 2000 ) locind, MAXMEM
          stop
        endif
        call rdlist( ia(ndset(i,3)), ndset(i,2), lum )
      end do

C ... Echo the nodeset data
      write( luo, 3000 ) nndset
      do i=1,nndset
        write( luo, 3010 ) ndset(i,1), ndset(i,2)
c       call wrlist( ia(ndset(i,3)), ndset(i,1), ndset(i,2), luo )
      end do

      return
 1000 format(i10)
 1010 format(//,
     2  5x,'Error reading nodesets ...',/,
     3  5x,'Nnd_sets in header: ',i10,/,
     4  5x,'Nnd_sets in footer: ',i10,/)
 1020 format(2i10)
 2000 format(//,
     2  5x,'!!!!! Memory Exceeded !!!!!',/,
     3  5x,'Allocation of nodeset lists in rwnset',/
     4  5x,'Requested: ',i10,/,
     5  5x,'Available: ',i10,/,
     6  5x,'!!!!! Memory Exceeded !!!!!',//)
 3000 format(//,
     2 5x,'Nodeset Parameters:',/,
     3 5x,'===================',//,
     4 5x,'Number of Node Sets: ',i5,//,
     5 5x,'Node Set# Node Set Id',/,
     6 5x,'========= ===========')
 3010 format(5x,i9,3x,i9)
      end

      subroutine rdlist(ndlist,   nnd,   lum)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : rdlist - read a nodeset node list
C
C Argument |    Meaning
C =========|============================================================
C ndlist   | node list
C nnd      | # of nodes in the lsit
C lum      | mesh LUN
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
      dimension ndlist(nnd)
      character *80 txt
      do i=1,nnd
        call gttxsg( txt, nl, lum )
        read ( txt, 1000 ) nd, ndlist(i)
      end do
      return
 1000 format(2i10)
      end

      subroutine wrlist(ndlist,    id,   nnd,   luo)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : wrlist - write a nodeset node list
C
C Argument |    Meaning
C =========|============================================================
C ndlist   | node list
C id       | nodeset id
C nnd      | # of nodes in the lsit
C luo      | output LUN
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
      dimension ndlist(nnd)

      write( luo, 1000 ) id, nnd
      do i=1,nnd
        write( luo, 1010 ) ndlist(i)
      end do

      return
 1000 format(//,
     2 5x,'Nodeset id: ',i5,/,
     3 5x,'Nnp in set: ',i5,//,
     4 5x,'Nodes',/,
     5 5x,'=====')
 1010 format(5x,i5)
      end

      subroutine bandwd(    np,   nel,  nnpe,  nwix,   nbw)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : bandwd - compute the 1/2-bandwidth
C
C Argument |    Meaning
C =========|============================================================
C np       | connectivity array
C nel      | # of elements
C nnpe     | # of nodes per elements
C nwix     | size of connectivity
C nbw      | 1/2-bandwidth
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      dimension np(nel,nwix)

      nbw=0
      do i = 1, nnpe
        do j = i+1, nnpe
          do k=1,nel
            ibw = np(k,j) - np(k,i) + 1
            if( ibw .ge. nbw ) nbw = ibw
          end do
        end do
      end do

      return
      end

      subroutine formkf(     T, ltype,     x,    rk,     f,    np,
     2                     nnp,   nel,  nwix,   nbw)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : formkf - form & assemble the K and F matrices (UDU format)
C
C Argument |    Meaning
C =========|============================================================
C T        | wire tension
C ltype    | load type =0 constant, =1 linear, =2 split 
C x        | nodal coordinates
C rk       | stiffness
C f        | rhs (force)
C np       | connectivity
C nnp      | # of nodes
C nel      | # of elements
C nwix     | size of np
C nbw      | 1/2-bandwidth
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      dimension x(nnp), rk(nnp,nbw), f(nnp)
      dimension np(nel,nwix)

C ... zero first
      do j=1,nbw
        do i=1,nnp
          rk(i,j)=0.0
        end do
      end do
      do i=1,nnp
        f(i)=0.0
      end do

C ... Loop over elements to assemble
      do i=1, nel
        n1=np(i,1)
        n2=np(i,2)
        ibo=abs(n2-n1)+1

C ..... Setup stiffness first
        rle = abs( x(n2) - x(n1) )
        rkij = T/rle

C ..... Only assemble the diagonal & upper-diagonal terms
        rk(n1,1) = rk(n1,1) + rkij
        rk(n2,1) = rk(n2,1) + rkij
        if ( n1 > n2 ) then
          rk(n2,ibo) = rk(n2,ibo) - rkij
        else 
          rk(n1,ibo) = rk(n1,ibo) - rkij
        endif

C ..... Now compute the load
        if ( ltype .eq. 0 ) then
C ....... constant weight per unit length
          w0  = 10.0
          fnd = 0.5*w0*rle
          f(n1) = f(n1) + fnd
          f(n2) = f(n2) + fnd
        else if ( ltype .eq. 1 ) then
C ....... linear weight per unit length

C         BLEEPED CODE GOES HERE 

        else if ( ltype .eq. 2 ) then
C ....... split constant-linear weight per unit length

C         BLEEPED CODE GOES HERE 

        endif
      end do
      return
      end

      subroutine blstbc(    rk,     f,   nnp)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : blstbc - use penalty method to enforce bc's
C
C Argument |    Meaning
C =========|============================================================
C rk       | stiffness
C f        | force
C nnp      | # of nodes per element
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      dimension rk(nnp), f(nnp)

C ... All homogeneous EBC's for wire problem

      f(  1) = 0.0
      f(nnp) = 0.0
      rk(  1) = rk(  1)*1.0e+10
      rk(nnp) = rk(nnp)*1.0e+10

      return
      end

      subroutine writeu(     x,     u,   nnp,   luo,   lup)
C-----------------------------------------------------------------------
C                                                                       
C* Routine : writeu - write output and plot files
C
C Argument |    Meaning
C =========|============================================================
C x        | nodal coordinates
C u        | solution vector
C nnp      | # of ndoes
C luo      | output LUN
C lup      | plot   LUN
C
C* Author  : Mark A. Christon 
C
C-----------------------------------------------------------------------
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      dimension x(nnp), u(nnp)

      write( luo, 1000 )
      write( lup, 1010 ) nnp
      do i=1,nnp
        write(luo,1020)   i ,u(i)
        write(lup,1030) x(i),u(i)
      end do

      return
 1000 format(//,
     2  5x,'Nodal Displacement',/,
     3  5x,'==================',//,
     4  5x,'Node  Displacement',/,
     5  5x,'====  ============')
 1010 format(
     2  '#',/,
     3  '# Nnp=',i4,/,
     5  '#',/,
     6  '# Coord.  Displacement',/,
     7  '#')
 1020 format(5x,i5,2x,1pe10.4)
 1030 format(5x,1pe10.4,1x,1pe10.4)
      end