! fem_bihar2dps_la.f90  Galerkin FEM biharmonic 4th order
!                                                                   
! solve  uxxxx(x,y) + uyyyy(x,y) + 2*uxxyy(x,y) = f(x,y)
!
!  f(x,y) := 4.0*sin(x)*sin(y)
!
!  in the rectangle  xmin< x <xmax and  ymin< y <ymax  ,
!
!     subject to boundary conditions of the first
!     kind, Dirichlet
!
!  Dirichlet boundary conditions are computed using:
!  ub(x,y) := sin(x)*sin(y)
!
!  uses high order shape functions
!  uses high order Gauss-Legendre quadrature
!  uses few degrees of freedom
!  uses laphi.f90  simeq.f90  deriv.f90  gaulegf.f90 

module global
  implicit none
  integer, parameter :: nx=9
  integer, parameter :: ny=9
  integer, parameter :: nxy = nx*ny
  integer, parameter :: npx = 9
  integer, parameter :: npy = 9
  double precision, dimension(nx*ny,nx*ny) :: kg ! ((i-1)*ny+ii,(j-1)*ny+jj)
  double precision, dimension(nx) :: xg
  double precision, dimension(ny) :: yg
  double precision, dimension(nx*ny) :: fg
  double precision, dimension(nx*ny) :: ug ! solution
  double precision, dimension(nx*ny) :: Ua
  double precision :: xmin = 0.0
  double precision :: xmax = 3.14159
  double precision :: ymin = 0.0
  double precision :: ymax = 3.14159
  double precision :: hx, hy 
  double precision, dimension(npx) :: xx ! for Gauss-Legendre
  double precision, dimension(npx) :: wx
  double precision, dimension(npy) :: yy
  double precision, dimension(npy) :: wy
  integer :: px, py

  interface
    function FF(x,y) result(z) ! forcing function
      double precision, intent(in) :: x,y
      double precision :: z
    end function FF
    function ub(x,y) result(z) ! analytic solution and boundary
      double precision, intent(in) :: x,y
      double precision :: z
    end function ub
    subroutine simeq(n, sz, A, Y, X)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz), intent(in) :: Y
      double precision, dimension(sz), intent(out) :: X
    end subroutine simeq
    subroutine gaulegf(x1, x2, x, w, n)
      integer, intent(in) :: n
      double precision, intent(in) :: x1, x2
      double precision, dimension(n), intent(out) :: x, w
    end subroutine gaulegf
    subroutine deriv(order, npoints, point, a, bb)
      integer, intent(in) :: order, npoints, point
      integer, dimension(0:*), intent(out) :: a
      integer, intent(out) :: bb
    end subroutine deriv 
    subroutine rderiv(order, npoints, point, hx, cx)
      integer, intent(in) :: order, npoints, point
      double precision, intent(in) :: hx
      double precision, dimension(0,*), intent(out) :: cx
    end subroutine rderiv 
  end interface

end module global

! functions called by program
! PDE problem definition functions  f and ub 
  function FF(x,y) result(z) ! forcing function
    implicit none
    double precision, intent(in) :: x,y
    double precision :: z

    z = 4.0*sin(x)*sin(y)
  end function FF

  function ub(x,y) result(z) ! analytic solution and boundary
    implicit none
    double precision, intent(in) :: x,y
    double precision :: z

    z = sin(x)*sin(y)
  end function ub

  function galk(x,y,  i,j, ii,jj) result(z)
                     ! Galerkin k stiffness function for this specific problem
    use global
    use laphi
    implicit none
    double precision, intent(in) :: x,y
    integer, intent(in) :: i, j, ii, jj
    double precision :: z

    z = (    phipppp(x,j,nx,xg)*     phi(y,jj,ny,yg) +  &
         2.0 * phipp(x,j,nx,xg)*   phipp(y,jj,ny,yg) +  &
                 phi(x,j,nx,xg)* phipppp(y,jj,ny,yg) )* &
                 phi(x,i,nx,xg)*     phi(y,ii,ny,yg)

  end function galk

  function galf(x, y, i, ii) result(z) ! Galerkin f function
    use global
    use laphi
    implicit none
    double precision, intent(in) :: x,y
    integer, intent(in) :: i, ii
    double precision :: z

    z = FF(x,y)*phi(x,i,nx,xg)*phi(y,ii,ny,yg)
  end function galf

  subroutine write_soln()
    use global
    implicit none
    integer :: i, ii, iostat

    open(unit=11, file="fem_bihar2dps_la_f90.dat", status="replace",  &
         form="formatted", action="readwrite", iostat=iostat)
    print *, "iostat=", iostat
    if(iostat/=0) then
      print *,'file fem_bihar2dps_la_f90.dat can not be opened for writing'
      return
    end if
 
    print *,'writing fem_bihar2dps_la_f90.dat'
    do i=1,nx
      do ii=1,ny
        write(unit=11, fmt='(3f11.6)', iostat=iostat) &
              xg(i),yg(ii),ug((i-1)*ny+ii)
      end do 
        write(unit=11, fmt='(1h )', iostat=iostat)
    end do 
    close(unit=11, status="keep")
    print *,'finished writing fem_bihar2dps_la_f90.dat'
  end subroutine write_soln 

  subroutine check_soln()
               ! for all interior points 
    use global
    implicit none
    double precision :: x, y, err, maxerr, avgerr, rmserr
    double precision :: sumerr, sumsqerr, eval_deriv
    double precision, dimension(nx,nx) :: cxx
    double precision, dimension(nx,nx) :: cxxxx
    double precision, dimension(ny,ny) :: cyy
    double precision, dimension(ny,ny) :: cyyyy
    double precision, dimension(ny) :: tcxxyy
    double precision, dimension(nx) :: tempx
    double precision, dimension(ny) :: tempy
    integer :: cnt = 0
    integer :: i, j, ii, jj, k, m, p

    do i=2,nx-1
      call rderiv(2,nx,i,hx,tempx)
      do p=1,nx
        cxx(i,p) = tempx(p)
      end do
      call rderiv(4,nx,i,hx,tempx)
      do p=1,nx
        cxxxx(i,p) = tempx(p)
      end do
    end do
    do ii=2,ny-1
      call rderiv(2,ny,ii,hy,tempy)
      do p=1,ny
        cyy(i,p) = tempy(p)
      end do
      call rderiv(4,ny,ii,hy,tempy)
      do p=1,ny
        cyyyy(i,p) = tempy(p)
      end do
    end do

    sumerr = 0.0
    sumsqerr = 0.0
    maxerr = 0.0
    do i=2,nx-1 ! full equation 
      x = xmin + i*hx
      do j=2,ny-1
        y = ymin + j*hy
        eval_deriv = 0.0 ! no PDE u term 
        do k=1,nx
	  eval_deriv = eval_deriv+cxxxx(i,k)*ug((k-1)*ny+j) ! PDE  uxxxx term 
        end do
        do k=1,ny
	  eval_deriv = eval_deriv+cyyyy(j,k)*ug((i-1)*ny+k)  ! PDE uyyyy term 
        end do
        !  + 2*uxxyy 
        do jj=1,ny ! walk jj in y, computing uxx
 	  tcxxyy(jj) = 0.0 ! numeric  uxx part of uxxyy 
          do k=1,nx
	    tcxxyy(jj) = tcxxyy(jj)+cxx(i,k)*ug((k-1)*ny+jj)
          end do
        end do
        do k=1,ny  ! uyy part of uxxyy
          eval_deriv = eval_deriv+2.0*cyy(j,k)*tcxxyy(k)  ! PDE uxxyy term 
        end do

        err = eval_deriv-FF(x,y)
        err = abs(err)
        sumerr = sumerr+err
        sumsqerr = sumsqerr+err*err
        maxerr = max(err,maxerr)
        cnt = cnt+1
      end do ! end y 
    end do ! end x 
    rmserr = sqrt(sumsqerr/cnt)
    avgerr = sumerr/cnt
    print *,'check_soln values '
    print *,'maxerr=',maxerr,', rmserr=',rmserr,', avgerr=',avgerr
  end subroutine check_soln  

program fem_bihar2dps_la
  use global
  use laphi
  implicit none
  integer :: i, j, ii, jj
  double precision :: x, y
  double precision :: val, err, avgerr, maxerr
  integer :: start, now, prev ! timing 
  integer :: count, count_rate, count_max  ! for system_clock 

  interface
    function galk(x,y, i,j, ii,jj) result(z)
      double precision, intent(in) :: x,y
      integer, intent(in) :: i, j, ii, jj
      double precision :: z
    end function galk
    function galf(x,y, i, ii) result(z)
      double precision, intent(in) :: x,y
      integer, intent(in) :: i, ii
      double precision :: z
    end function galf
  end interface

  print *, 'fem_bihar2dps_la.f90 running'
  print *, 'solve uxxxx(x,y) + 2 uxxyy(x,y) + uyyyy(x,y) = F(x,y)'
  print *, 'F(x,y) = 4*sin(x)*sin(y)'
  print *, 'boundary conditions computed using ub(x,y)'
  print *, 'ub(x,y) = sin(x)*sin(y)'
  print *, 'high order shape functions used'
  print *, 'high order Gauss-Legendre quadrature used'
  print *, ((nx-2)*(ny-2)),'  Degrees of Freedom'
  print *, ' '
  print *, 'xmin=',xmin,', ymin =',ymin
  print *, 'xmax=',xmax,', ymax =',ymax
  print *, 'nx=',nx,', ny =',ny

  call system_clock(count, count_rate, count_max)
  start = count
  prev = count

  hx = (xmax-xmin)/(nx-1) ! does not have to be uniform spacing
  print *, 'x grid and analytic solution, at ymax'
  do i=1,nx 
    xg(i) = xmin+(i-1)*hx
    print *, 'i=',i,', ub(',xg(i),',',ymax,')=',ub(xg(i),ymax)
  end do
  print *, ' '

  hy = (ymax-ymin)/(ny-1) ! does not have to be uniform spacing
  print *, 'y grid and analytic solution, at xmax'
  do ii=1,ny 
    yg(ii) = ymin+(ii-1)*hy
    print *, 'ii=',ii,', ub(',xmax,',',yg(ii),')=',ub(xmax,yg(ii))
  end do
  print *, ' '

  ! put solution in Ua vector
  do i=1,nx
    x = xmin+(i-1)*hx
    do ii=1,ny
      y = ymin+(ii-1)*hy
      Ua((i-1)*ny+ii) = ub(x,y)
      ! print "('solution at i=',I1,', x=',f6.3,', ii=',I1, &
      !         ', y=',f6.3,' is ',f7.3)", &
      !         i, x, ii, y, Ua((i-1)*ny+ii)
    end do
  end do

  ! initialize k and f
  do i=1,nx
    do j=1,nx 
      do ii=1,ny 
        do jj=1,ny 
          kg((i-1)*ny+ii,(j-1)*ny+jj) = 0.0
        end do
      end do
    end do
  end do
  do i=1,nx
    do ii=1,ny 
      fg((i-1)*ny+ii) = 0.0
    end do
  end do

  ! set boundary conditions, four sides
  do i=1,nx
    x = xg(i)
    ii = 1
    y = yg(ii)
    kg((i-1)*ny+ii,(i-1)*ny+ii) = 1.0
    fg((i-1)*ny+ii) = ub(x,y)
    ! print "('boundary i=',I1,', x=',f6.3,', ii=',I1,', &
    !          y=',f6.3,' is ',f10.5)", &
    !          i, x, ii, y, fg((i-1)*ny+ii)

    ii = ny
    y = yg(ii)
    kg((i-1)*ny+ii,(i-1)*ny+ii) = 1.0
    fg((i-1)*ny+ii) = ub(x,y)
    ! print "('boundary i=',I1,', x=',f6.3,', ii=',I1,', &
    !         y=',f6.3,' is ',f10.5)", &
    !         i, x, ii, y, fg((i-1)*ny+ii)
  end do ! nx
  do ii=1,ny
    y = yg(ii)
    i = 1
    x = xg(i)
    kg((i-1)*ny+ii,(i-1)*ny+ii) = 1.0
    fg((i-1)*ny+ii) = ub(x,y)
    ! print "('boundary i=',I1,', x=',f6.3,', ii=',I1,', &
    !         y=',f6.3,' is ',f10.5)", &
    !         i, x, ii, y, fg((i-1)*ny+ii)

    i = nx
    x = xg(i)
    kg((i-1)*ny+ii,(i-1)*ny+ii) = 1.0
    fg((i-1)*ny+ii) = ub(x,y)
    ! print "('boundary i=',I1,', x=',f6.3,', ii=',I1,', &
    !         y=',f6.3,' is ',f10.5)", &
    !         i, x, ii, y, fg((i-1)*ny+ii)
  end do ! ny

  print *, 'calling gaulegf npx=',npx
  call gaulegf(xmin, xmax, xx, wx, npx) ! xx and ww set up for integrations
  print *, 'calling gaulegf npy=',npy
  call gaulegf(ymin, ymax, yy, wy, npy)
  print *, 'xx(1)=',xx(1),', xx(2)=',xx(2),', wx(1)=',wx(1),', wx(2)=',wx(2)
  print *, 'yy(1)=',yy(1),', yy(2)=',yy(2),', wy(1)=',wy(1),', wy(2)=',wy(2)
  val = galk(xx(2),yy(2),2,2,2,2)
  print *, 'galk(xx(2),yy(2),2,2,2,2)=',val
  val = galf(xx(2),yy(2),2,2)
  print *, 'galf(xx(2),yy(2),2,2)=',val

  call system_clock(count, count_rate, count_max)
  now = count
  print *,'time to initialize = ',real(now-prev)/1000.0,' seconds'
  prev = now

  ! compute entries=stiffness matrix
  print *, 'compute stiffness matrix'
  do i=2,nx-1 
    do j=1,nx 
      do ii=2,ny-1 
        do jj=1,ny
          ! compute integral for k(i,j)
          val = 0.0
          do px=1,npx 
            do py=1,npy 
              val = val + wx(px)*wy(py)*galk(xx(px),yy(py),i,j,ii,jj)
            end do
          end do
          ! print "('Legendre integration=',f12.5,' at i=',I1,', j=', &
          !         I1,', ii=',I1,', jj=',I1)", &
          !         val, i, j, ii, jj
          kg((i-1)*ny+ii,(j-1)*ny+jj) = val
        end do
      end do
    end do
  end do

  ! compute integral for f(i)
  do i=2,nx-1 
    do ii=2,ny-1 
      val = 0.0
      do px=1,npx 
        do py=1,npy 
          val = val + wx(px)*wy(py)*galf(xx(px),yy(py),i,ii)
        end do
      end do
      ! print "('Legendre integration=',f12.5,' at i=',I1,', ii=',I1)", &
      !          val, i, ii
      fg((i-1)*ny+ii) = val
    end do
  end do

  print *, 'k computed stiffness matrix, see above'
  print *, 'f computed forcing function, see above'
  print *, ' '


  call system_clock(count, count_rate, count_max)
  now = count
  print *,'time for stiffness and forcing = ',real(now-prev)/1000.0,' seconds'
  prev = now

  call simeq(nx*ny, nx*ny, kg, fg, ug)

  call system_clock(count, count_rate, count_max)
  now = count
  print *,'time to solve equations = ',real(now-prev)/1000.0,' seconds'
  prev = now

  call write_soln() ! to file x y soln 

  call check_soln() ! independent of nowing analytic solution 


  call system_clock(count, count_rate, count_max)
  now = count
  print *,'time write and check solution = ',real(now-prev)/1000.0,' seconds'
  prev = now

  avgerr = 0.0
  maxerr = 0.0
  print *, 'ug computed Galerkin, Ua analytic, error'
  do i=1,nx 
    do ii=1,ny 
      err = abs(ug((i-1)*ny+ii)-Ua((i-1)*ny+ii))
      if(err>maxerr) then
        maxerr = err
      end if
      avgerr = avgerr + err
      print "('ug(',I1,',',I1,')=',f7.4,', Ua=',f7.4,', err=',e15.7)", &
            i, ii, ug((i-1)*ny+ii), Ua((i-1)*ny+ii), &
            (ug((i-1)*ny+ii)-Ua((i-1)*ny+ii))
    end do
  end do
  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nx*ny)
  print *,'xmax=',xmax,', ymax=',ymax,', nx=',nx,', ny=',ny,', npx=',npx, &
          ', npy=',npy

  call system_clock(count, count_rate, count_max)
  now = count
  print *,'total wall time = ',real(now-start)/1000.0,' seconds'
  print *, 'end fem_bihar2dps_la.f90 '

end program fem_bihar2dps_la