! fem_check22_la.f90  Galerkin FEM
! solve uxx(x,y) + 2 uxy(x,y) + 3 uyy(x,y) +
!       4 ux(x,y)  + 5 uy(x,y)  + 6 u(x,y) = c(x,y)
! boundary conditions computed using u(x)
! analytic solution may be given by u(x) = x^3 + y^3 + xy + 1
! Gauss-Legendre integration used
! solve for Uij numerical approximation U(x_i,y_j) of u(x_i,y_j)
! k * U = f, solve simultaneous equations for U

module global
  implicit none
  integer, parameter :: nx=4
  integer, parameter :: ny=4
  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
  double precision, dimension(nx*ny) :: Ua
  double precision :: xmin = 0.0
  double precision :: xmax = 1.0
  double precision :: ymin = 0.0
  double precision :: ymax = 1.0
  double precision, dimension(100) :: xx ! for Gauss-Legendre
  double precision, dimension(100) :: wx
  double precision, dimension(100) :: yy
  double precision, dimension(100) :: wy
  double precision :: val, err, avgerr, maxerr
  double precision :: eps = 1.0e-6
  integer :: px, py, npx, npy

  interface
    function FF(x,y) result(z) ! forcing function
      double precision, intent(in) :: x,y
      double precision :: z
    end function FF
    function uana(x,y) result(z) ! analytic solution and boundary
      double precision, intent(in) :: x,y
      double precision :: z
    end function uana
    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
  end interface

end module global

! functions called by program
  function FF(x,y) result(z) ! forcing function
    implicit none
    double precision, intent(in) :: x,y
    double precision :: z

    z = 6.0*x*x*x+6.0*y*y*y+12.0*x*x+15.0*y*y+6.0*x*y+11.0*x+22.0*y+8.0
  end function FF

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

    z = x*x*x+y*y*y+x*y+1.0
  end function uana

  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 =    ( phipp(x,j,nx,xg)*  phi(y,jj,ny,yg) +  &
        2.0 * phip(x,j,nx,xg)* phip(y,jj,ny,yg) +  &
        3.0 *  phi(x,j,nx,xg)*phipp(y,jj,ny,yg) +  &
        4.0 * phip(x,j,nx,xg)*  phi(y,jj,ny,yg) +  &
        5.0 *  phi(x,j,nx,xg)* phip(y,jj,ny,yg) +  &
        6.0 *  phi(x,j,nx,xg)*  phi(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

program fem_check22_la
  use global
  use laphi
  implicit none
  integer :: i, j, ii, jj
  double precision :: x, y, hx, hy

  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_check22_la.f90 running'
  print *, 'solve uxx(x,y) + 2 uxy(x,y) + 3 uyy(x,y) +'
  print *, '      4 ux(x,y)  + 5 uy(x,y)  + 6 u(x,y) = c(x,y)'
  print *, 'boundary conditions computed using u(x)'
  print *, 'analytic solution may be given by u(x) = x^3 + y^3 + xy + 1'
  print *, 'Gauss-Legendre integration used'
  print *, ' '
  print *, 'xmin=',xmin,', ymin =',ymin
  print *, 'xmax=',xmax,', ymax =',ymax
  print *, 'nx=',nx,', ny =',ny

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

  hy = (ymax-ymin)/(ny-1) ! does not have to be uniform spacing
  print *, 'y grid and analytic solution, at xmin'
  do ii=1,ny 
    yg(ii) = ymin+(ii-1)*hy
    Ua(ii) = uana(xmin,yg(ii))
    print *, 'ii=',ii,', Ua(',yg(ii),')=',Ua(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) = uana(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) = uana(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) = uana(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) = uana(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) = uana(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

  npx = 48
  npy = 48
  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


  ! 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 simeq(nx*ny, nx*ny, kg, fg, ug)

  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 *, 'end fem_check22_la.f90 '

end program fem_check22_la