c*************************************************************************
c*************************************************************************
cc                                                                      cc
cc                           T P L A T E                                cc
cc                                                                      cc
cc            Program for the Finite Difference approximation           cc
cc                of thin plate bending problems using the              cc
cc                       Kirchhoff-Love theory                          cc
cc                                                                      cc
cc                       Ramon Codina, ETSECCP                          cc
cc                                                                      cc
cc                              9/1/92                                  cc
cc                                                                      cc
c*************************************************************************
c*************************************************************************
      implicit real*8 (a-h,o-z)
      parameter (maxunk=2500, maxban=100)
      common/deriva/nunk,hxsi,hysi,nban,dest,pois
      real*8     amat(maxban+1,maxunk), disp(maxunk), forc(maxunk),
     .           xxm(maxunk),           yym(maxunk),  xym(maxunk)
c
c***  This routine reads the problem data
c
      call inputd
c
c***  Force terms are now created
c
      call loadte(forc)
c
c***  This routine constructs the matrix of the linear system 
c
      call comatr(amat)
c
c***  The linear system is solved
c
      call solves(nunk,nban,amat,amat,1,forc,disp)
c
c***  Output of results
c
      call result(disp,xxm,yym,xym)
c
      stop
      end
      subroutine inputd
c*******************************************************************
c
c  This routine opens the result file and asks for data
c
c*******************************************************************
      implicit real*8 (a-h,o-z)
      character*20 arx_escr
      character*10 bounco(2)
      common/contro/xlon,ylon,nxdi,nydi,ncx0,ncxa,ncy0,ncyb
      common/deriva/nunk,hxsi,hysi,nban,dest,pois
c
c***  Open result file
c
      write(6,'(a,$)') ' > Result file name :'
      read(5,'(a)') arx_escr
      open(9,file=arx_escr,status='unknown')
c
c***  Read data
c
      write(6,'(a,$)') ' > Length in the x-direction (a) :'
      read(5,*) xlon
      write(6,'(a,$)') ' > Length in the y-direction (b) :'
      read(5,*) ylon
      write(6,'(a,$)') ' > Number of x-divisions (min.=6, max.=51) :'
      read(5,*) nxdi
      write(6,'(a,$)') ' > Number of y-divisions (min.=6, max.=51) :'
      read(5,*) nydi
      write(6,5)
   5  format(' > The possible types of boundary conditions are :'
     . ,/,5x,
     .     ' > (1) Clamped edge  ',/,5x,
     .     ' > (2) Simply supported edge ')
      write(6,'(a,$)') ' > Boundary conditions for x=0 :' 
      read(5,*) ncx0
      write(6,'(a,$)') ' > Boundary conditions for x=a :' 
      read(5,*) ncxa
      write(6,'(a,$)') ' > Boundary conditions for y=0 :' 
      read(5,*) ncy0
      write(6,'(a,$)') ' > Boundary conditions for y=b :' 
      read(5,*) ncyb
      write(6,'(a,$)') ' > Young modulus :' 
      read(5,*) eyou
      write(6,'(a,$)') ' > Poisson coefficient :' 
      read(5,*) pois
      write(6,'(a,$)') ' > Plate thickness :' 
      read(5,*) thic
c
c***  Set some variables
c
      nunk = (nxdi-1)*(nydi-1)
      hxsi = xlon/float(nxdi)
      hysi = ylon/float(nydi)
      nban = 2*(nydi-1)
      dest = eyou*(thic**3)/(12.0*(1.0-pois*pois))
c
c***  Output data
c
      bounco(1)='clampled  '
      bounco(2)='supported '
      write(9,10) xlon,ylon,nxdi,nydi,dest,
     .            bounco(ncx0),bounco(ncxa),bounco(ncy0),bounco(ncyb)
   10 format(//,5x,'RESULT FILE OF PROGRAM TPLATE',/,
     .          5x,'_____________________________',///,
     .          5x,'PROBLEM DATA  :',//,
     .          5x,'Size of the rectangle : ',f5.2,'  x  ',f5.2,/,
     .          5x,'Number of cells : ',i2,' x ',i2,/,
     .          5x,'Plate stiffness : ',e15.5,/,
     .          5x,'Boundary prescriptions:',/,
     .               10x,'x=0 :       ',a10,/,
     .               10x,'x=a :       ',a10,/,
     .               10x,'y=0 :       ',a10,/,
     .               10x,'y=b :       ',a10,//)
c
      return
      end
      subroutine loadte(forc)
c*******************************************************************
c
c  This routine creates the load terms
c
c*******************************************************************
      implicit real*8 (a-h,o-z)
      common/contro/xlon,ylon,nxdi,nydi,ncx0,ncxa,ncy0,ncyb
      common/deriva/nunk,hxsi,hysi,nban,dest,pois
      real*8  forc(nunk)
c
c***  Read load data 
c
      write(6,'(a,$)') ' > Number of loaded sub-rectangles :' 
      read(5,*) nlod
      write(6,5)
   5  format(' > The position of each sub-rectangle must be given',/,
     .       '   in grid-spacing units, introducing x1,x2,y1 and y2,',/,
     .       '   in this order. The corner loads must be given in',/,
     .       '   anti-clockwise order')
      do ilod=1,nlod
        write(6,'((a,i5,a),$)') ' > Position of rectangle ',ilod,'  :' 
        read(5,*) nx1,nx2,ny1,ny2
        if(nx2.ge.nxdi) nx2 = nxdi - 1
        if(ny2.ge.nydi) ny2 = nydi - 1
        write(6,'((a,i5,a),$)') ' > Corner loads for rectangle ',
     .  ilod,'  :' 
        read(5,*) q1,q2,q3,q4
        call noload(forc,nx1,nx2,ny1,ny2,q1,q2,q3,q4)
      end do
c
      return
      end
c
c***  Create nodal loads
c
      subroutine noload(forc,nx1,nx2,ny1,ny2,q1,q2,q3,q4)
      implicit real*8 (a-h,o-z)
      common/contro/xlon,ylon,nxdi,nydi,ncx0,ncxa,ncy0,ncyb
      common/deriva/nunk,hxsi,hysi,nban,dest,pois
      real*8  forc(nunk)
c
c***  Initialization
c
      do i=1,nunk
        forc(i) = 0.0
      end do
c
c***  Function for the nodal numbering (column-based) 
c
      nn(i,j) = (i-1)*(nydi-1) + j 
c
c***  Coordinates of the sub-rectangle
c
      x1 = hxsi*float(nx1)
      x2 = hxsi*float(nx2)
      y1 = hysi*float(ny1)
      y2 = hysi*float(ny2)
c
c***  Punctual load
c
      if(x1.eq.x2. and .y1.eq.y2) then
        nod = nn(nx1,ny1)
        forc(nod) = q1/(dest*hxsi*hysi)
        return
      end if
c
c*** `Knife' load
c
      if(x1.eq.x2) then  
        coef = 1.0/(y2-y1)
        do j = ny1,ny2
          nod = nn(nx1,j)
          y = hysi*float(j)
          q = q1*(y2-y) + q3*(y-y1) 
          forc(nod) = q*coef/(dest*hxsi)
        end do
      else if(y1.eq.y2) then  
        coef = 1.0/(x2-x1)
        do i = nx1,nx2
          nod = nn(i,ny1)
          x = hxsi*float(i)
          q = q1*(x2-x) + q2*(x-x1) 
          forc(nod) = q*coef/(dest*hysi)
        end do
        return
      end if
c
c***  Surface load
c
      if(x1.ne.x2. and .y1.ne.y2) then
        coef = 1.0/((x2-x1)*(y2-y1))
        do i = nx1,nx2
          do j = ny1,ny2
            nod = nn(i,j)
            x = hxsi*float(i)
            y = hysi*float(j)
            q = q1*(x2-x)*(y2-y) + q2*(x-x1)*(y2-y)
     .        + q3*(x-x1)*(y-y1) + q4*(x2-x)*(y-y1)
            forc(nod) = q*coef/dest
          end do
        end do
        return
      end if
c
      return
      end
      subroutine comatr(amat)
c*******************************************************************
c
c  This routine constructs the matrix of the linear system incorpo-
c  rating boundary conditions.
c
c*******************************************************************
      implicit real*8 (a-h,o-z)
      common/contro/xlon,ylon,nxdi,nydi,ncx0,ncxa,ncy0,ncyb 
      common/deriva/nunk,hxsi,hysi,nban,dest,pois
      real*8 amat(nban+1,nunk), sigc(2)
c
c***  Function for the nodal numbering (column-based) 
c
      nn(i,j) = (i-1)*(nydi-1) + j 
c
c***  Factors of the finite difference stencil
c
      hxx4 =  hxsi**4
      hyy4 =  hysi**4
      hxy2 =  hxsi*hxsi*hysi*hysi
      co00 =  6.0/hxx4 + 8.0/hxy2 + 6.0/hyy4
      co10 = -4.0/hxx4 - 4.0/hxy2
      co01 = -4.0/hxy2 - 4.0/hyy4
      co11 =  2.0/hxy2
      co20 =  1.0/hxx4
      co02 =  1.0/hyy4
      sigc(1) =  1.0
      sigc(2) = -1.0
c
c***  Initialization of the system matrix
c
      do i = 1,nban+1
        do j = 1,nunk
          amat(i,j) = 0.0
        end do
      end do
c
c***  Difference equations for the inner rectangle 
c
      do i = 2,nxdi-2
        do j = 2,nydi-2
          nod = nn(i,j)
          amat(nn(i  ,j  )-nod+1,nod) = co00 
          amat(nn(i  ,j+1)-nod+1,nod) = co01 
          if(j.ne.nydi-2) 
     .    amat(nn(i  ,j+2)-nod+1,nod) = co02
          amat(nn(i+1,j-1)-nod+1,nod) = co11
          amat(nn(i+1,j  )-nod+1,nod) = co10
          amat(nn(i+1,j+1)-nod+1,nod) = co11
          if(i.ne.nxdi-2) 
     .    amat(nn(i+2,j  )-nod+1,nod) = co20 
        end do
      end do
c
c***  Equations for one-grid-spacing-distant nodes
c
      do i = 2,nxdi-2
c
c***  Bottom row
c
        j = 1
        nod = nn(i,j)
        amat(nn(i  ,j  )-nod+1,nod) = co00 + sigc(ncy0)*co02
        amat(nn(i  ,j+1)-nod+1,nod) = co01 
        amat(nn(i  ,j+2)-nod+1,nod) = co02
        amat(nn(i+1,j  )-nod+1,nod) = co10
        amat(nn(i+1,j+1)-nod+1,nod) = co11
        if(i.ne.nxdi-2) 
     .  amat(nn(i+2,j  )-nod+1,nod) = co20 
c
c***  Top row
c
        j = nydi-1
        nod = nn(i,j)
        amat(nn(i  ,j  )-nod+1,nod) = co00 + sigc(ncyb)*co02
        amat(nn(i+1,j-1)-nod+1,nod) = co11
        amat(nn(i+1,j  )-nod+1,nod) = co10
        if(i.ne.nxdi-2) 
     .  amat(nn(i+2,j  )-nod+1,nod) = co20 
      end do
      do j = 2,nydi-2
c
c***  Left column
c
        i = 1
        nod = nn(i,j)
        amat(nn(i  ,j  )-nod+1,nod) = co00 +sigc(ncx0)*co20
        amat(nn(i  ,j+1)-nod+1,nod) = co01
        if(j.ne.nydi-2) 
     .  amat(nn(i  ,j+2)-nod+1,nod) = co02 
        amat(nn(i+1,j-1)-nod+1,nod) = co11
        amat(nn(i+1,j  )-nod+1,nod) = co10 
        amat(nn(i+1,j+1)-nod+1,nod) = co11
        amat(nn(i+2,j  )-nod+1,nod) = co20
c
c***  Right column
c
        i = nxdi-1
        nod = nn(i,j)
        amat(nn(i  ,j  )-nod+1,nod) = co00 + sigc(ncxa)*co20
        amat(nn(i  ,j+1)-nod+1,nod) = co01
        if(j.ne.nydi-2) 
     .  amat(nn(i  ,j+2)-nod+1,nod) = co02 
      end do
c
c***  Corner nodes
c
      nc1 = nn(1,1)
      amat(nn(1,1)-nc1+1,nc1) = co00 + sigc(ncx0)*co20 
     .                               + sigc(ncy0)*co02 
      amat(nn(1,2)-nc1+1,nc1) = co01
      amat(nn(1,3)-nc1+1,nc1) = co02
      amat(nn(2,1)-nc1+1,nc1) = co10
      amat(nn(2,2)-nc1+1,nc1) = co11
      amat(nn(3,1)-nc1+1,nc1) = co20
c
      nc2 = nn(nxdi-1,1)
      amat(nn(nxdi-1,1)-nc2+1,nc2) = co00 + sigc(ncxa)*co20 
     .                                    + sigc(ncy0)*co02 
      amat(nn(nxdi-1,2)-nc2+1,nc2) = co01
      amat(nn(nxdi-1,3)-nc2+1,nc2) = co02
c
      nc3 = nn(nxdi-1,nydi-1)
      amat(nn(nxdi-1,nydi-1)-nc3+1,nc3) = co00 + sigc(ncxa)*co20 
     .                                         + sigc(ncyb)*co02 
c
      nc4 = nn(1     ,nydi-1)
      amat(nn(1,nydi-1)-nc4+1,nc4) = co00 + sigc(ncx0)*co20 
     .                                    + sigc(ncyb)*co02 
      amat(nn(2,nydi-1)-nc4+1,nc4) = co10
      amat(nn(3,nydi-1)-nc4+1,nc4) = co20
      amat(nn(2,nydi-2)-nc4+1,nc4) = co11
c
      return
      end
      subroutine solves(necua,nsban,amatr,cmatr,nsist,bvect,xvect)
c*******************************************************************
c
c  This routine solves a linear system A x = b using the Cholesky
c  decomposition and a forward and a backward substitution
c
c*******************************************************************
      implicit real*8(a-h,o-z)
      real*8  amatr(nsban+1,necua),bvect(necua,nsist),
     .        cmatr(nsban+1,necua),xvect(necua,nsist)

      call cholesky(necua,nsban,amatr,cmatr)
      call sub1(necua,nsban,cmatr,nsist,bvect,xvect)
      call sub2(necua,nsban,cmatr,nsist,xvect,xvect)
c
      return
      end
c
c***  Cholesky decomposition
c
      subroutine cholesky(n,l,a,c)
      implicit real*8(a-h,o-z)
      real*8  a(l+1,n),c(l+1,n)
c
      do i = 1,n
        c(1,i) = a(1,i)
        k1 = i-l
        if (k1.lt.1) k1 = 1
        if (k1.le.i-1) then 
          do k = k1,i-1
            c(1,i) = c(1,i) - c(i-k+1,k)*c(i-k+1,k)
          end do
        end if
        if (c(1,i).le.0.) call error(' non-positive element in row ',i)
        c(1,i) = dsqrt(c(1,i))
        j2 = i + l
        if (j2.gt.n) j2 = n
        if (i+1.le.j2) then 
          do j = i+1,j2
            c(j-i+1,i) = a(j-i+1,i)
            k1 = j - l
            if (k1.lt.1) k1 = 1
            if (k1.le.i-1) then 
              do k = k1,i-1
                c(j-i+1,i) = c(j-i+1,i) - c(i-k+1,k)*c(j-k+1,k)
              end do
            end if
            c(j-i+1,i) = c(j-i+1,i)/c(1,i)
          end do
        end if
      end do
c
      return
      end
c
c***  Forward substitution
c
      subroutine sub1(n,l,c,m,b,z)
      implicit real*8(a-h,o-z)
      real*8 c(l+1,n),b(n,m),z(n,m)
c
      do is = 1,m
        z(1,is) = b(1,is)/c(1,1)
        do i = 2,n
          z(i,is) = b(i,is)
          k1 = i-l
          if (k1.lt.1) k1 = 1
          if (k1.le.i-1) then 
            do k = k1,i-1
              z(i,is) = z(i,is)-c(i-k+1,k)*z(k,is)
            end do
          end if
          if (c(1,i).le.0.) call error(' non-positive pivot in row',i)
          z(i,is) = z(i,is)/c(1,i)
        end do
      end do
c 
      return
      end
c
c***  Backward substitution
c
      subroutine sub2(n,l,c,m,z,x)
      implicit real*8(a-h,o-z)
      real*8  c(l+1,n),z(n,m),x(n,m)

      do is = 1,m
        x(n,is) = z(n,is)/c(1,n)
        do i = n-1,1,-1
          x(i,is) = z(i,is)
          k2 = i+l
          if (k2.gt.n) k2 = n
          if (i+1.le.k2) then 
            do k = i+1,k2
              x(i,is) = x(i,is)-c(k-i+1,i)*x(k,is)
            end do
          end if
          if (c(1,i).le.0.) call error(' non-positive pivot in row',i)
          x(i,is) = x(i,is)/c(1,i)
        end do
      end do
c 
      return
      end
c
c***  Output of possible errors
c
      subroutine error(alfa,i)
      implicit real*8(a-h,o-z)
      character alfa*(*)
      write(6,'(a,i5)') alfa,i
      stop
      end
      subroutine result(disp,xxm,yym,xym)
c*******************************************************************
c
c  This routine outputs the computed displacements
c
c*******************************************************************
      implicit real*8 (a-h,o-z)
      common/contro/xlon,ylon,nxdi,nydi,ncx0,ncxa,ncy0,ncyb
      common/deriva/nunk,hxsi,hysi,nban,dest,pois
      character*1  con
      character*20 arx_post, xtitle, ytitle
      character*35 gtitle
      character*5  stitle
      real*8  disp(nunk),xxm(nunk),yym(nunk),xym(nunk)
c
c***  Function for the nodal numbering (column-based) 
c
      nn(i,j) = (i-1)*(nydi-1) + j 
c
c***  Bending moments in interior points
c
      do i = 2,nxdi-2
        do j = 2,nydi-2
          nod = nn(i,j)
          xxm(nod) = -dest*((disp(nn(i+1,j))-2.0*disp(nod)
     .                +disp(nn(i-1,j)))/(hxsi*hxsi)
     .                +(disp(nn(i,j+1))-2.0*disp(nod)
     .                +disp(nn(i,j-1)))*pois/(hysi*hysi))
          yym(nod) = -dest*((disp(nn(i+1,j))-2.0*disp(nod)
     .                +disp(nn(i-1,j)))*pois/(hxsi*hxsi)
     .                +(disp(nn(i,j+1))-2.0*disp(nod)
     .                +disp(nn(i,j-1)))/(hysi*hysi))
          xym(nod) = -dest*(1.0-pois)*(disp(nn(i+1,j+1))
     .                -disp(nn(i+1,j-1))-disp(nn(i-1,j+1))
     .                +disp(nn(i-1,j-1)))/(4.0*hxsi*hysi)
        end do
      end do
c
c***  Writes to result file
c
      write(9,10)
      do i = 1,nxdi-1
        do j = 1,nydi-1
          nod = nn(i,j)
          write(9,20) i,j,disp(nod)
        end do
      end do
   10 format(/,5x,'NODAL DISPLACEMENTS',//,
     .         5x,'Row number   Column number    Displacement',/)
   20 format(5x,i5,10x,i5,5x,f17.10)
c
      write(9,11)
      do i = 2,nxdi-2
        do j = 2,nydi-2
          nod = nn(i,j)
          write(9,21) i,j,xxm(nod),yym(nod),xym(nod)
        end do
      end do
   11 format(/,5x,'NODAL MOMENTS',//,
     .         5x,'Row No.   Col. No.   x-Moment    y-Moment',
     .                                         '   xy-Moment',/)
   21 format(5x,i7,6x,i5,1x,f10.5,2x,f10.5,3x,f10.5)
c
c***  Writes for post-processor
c
      write(6,'(a,$)') ' > Do you want to plot some results (y,n) ? :'
      read(5,'(a)') con
      if(con.eq.'n') return
      write(6,'(a,$)') ' > Plotting file name :'
      read(5,'(a)') arx_post
      open(10,file=arx_post,status='unknown')
      write(6,'(a,$)') ' > Global title for the figure :'
      read(5,'(a)') gtitle
      write(6,'(a,$)') ' > Title for the x-axis (a20) :'
      read(5,'(a)') xtitle
      write(6,'(a,$)') ' > Title for the y-axis (a20) :'
      read(5,'(a)') ytitle
      write(6,'(a,$)') ' > Number of x=const. lines :'
      read(5,*) nxli
      write(6,'(a,$)') ' > Number of y=const. lines :'
      read(5,*) nyli
      ntli = nxli + nyli
      write(10,30) ntli,'(f20.7,f20.7)',gtitle
c
      npox = nxdi + 1
      npoy = nydi + 1
      ntyc = -80
      icur = 0
      if(nxli.gt.0) then
        stitle = ' x = '
        do ixli = 1,nxli
          icur = icur + 1
          ntyc = ntyc - 1
          write(6,'((a,i5,a),$)') ' > Position of line x=c. No.',
     .                                ixli,'  :' 
          read(5,*) kx
          xco = hxsi*float(kx)
          write(10,40) npoy, ntyc, 0, 0, 1, 0.0, 0.0, 0.25, stitle, xco
          if(icur.eq.1) then
            write(10,50) xtitle,ytitle
          else
            write(10,50) '          ', '          '
          end if
          write(10,60) 0.0, 0.0
          do j = 1,nydi - 1
            y = hysi*float(j)
            nod = nn(kx,j)
            ade = -disp(nod)
            write(10,60) y, ade
          end do
          write(10,60) ylon, 0.0
        end do
      end if
c
      if(nyli.gt.0) then
        stitle = ' y = '
        do iyli = 1,nyli
          icur = icur + 1
          ntyc = ntyc - 1
          write(6,'((a,i5,a),$)') ' > Position of line y=c. No.',
     .                                iyli,'  :' 
          read(5,*) ky
          yco = hysi*float(ky)
          write(10,40) npox, ntyc, 0, 0, 1, 0.0, 0.0, 0.25, stitle, yco
          if(icur.eq.1) then
            write(10,50) xtitle,ytitle
          else
            write(10,50) '          ', '          '
          end if
          write(10,60) 0.0, 0.0
          do i = 1,nxdi - 1
            x = hxsi*float(i)
            nod = nn(i,ky)
            ade = -disp(nod)
            write(10,60) x, ade
          end do
          write(10,60) xlon, 0.0
        end do
      end if
c
   30 format(i5,a40,a35)
   40 format(5i5,3f10.2,/,a5,f5.2)
   50 format(a20,/,a20)
   60 format(f20.7,f20.7)
c
      return
      end
