! fem_check_abc_la.f90  Galerkin FEM
! solve a1(x,y) uxx(x,y) + b1(x,y) uxy(x,y) + c1(x,y) uyy(x,y) +
!       d1(x,y) ux(x,y)  + e1(x,y) uy(x,y)  + f1(x,y) u(x,y) = c(x,y)
! boundary conditions computed using u(x)
! analytic solution may be given by u(x)
! Gauss-Legendre integration used
! solve for Uij numerical approximation U(x_i,y_j) of u(x_i,y_j)
! kg * ug = fg, solve simultaneous equations for U

module global
  use abc_la
  implicit none
  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, 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
    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 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 = (a1(x,y)* phipp(x,j,nx,xg)*  phi(y,jj,ny,yg) +  &
         b1(x,y)*  phip(x,j,nx,xg)* phip(y,jj,ny,yg) +  &
         c1(x,y)*   phi(x,j,nx,xg)*phipp(y,jj,ny,yg) +  &
         d1(x,y)*  phip(x,j,nx,xg)*  phi(y,jj,ny,yg) +  &
         e1(x,y)*   phi(x,j,nx,xg)* phip(y,jj,ny,yg) +  &
         f1(x,y)*   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 = c(x,y)*phi(x,i,nx,xg)*phi(y,ii,ny,yg)
  end function galf

program fem_check_abc_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_check_abc_la.f90 running'
  print *, 'solve a1(x,y) uxx(x,y) + b1(x,y) uxy(x,y) + c1(x,y) uyy(x,y) +'
  print *, '     d1(x,y) ux(x,y)  + e1(x,y) uy(x,y)  + f1(x,y) u(x,y) = c(x,y)'
  print *, 'boundary conditions computed using u(x)'
  print *, 'analytic solution may be given by u(x) ifcheck>0'
  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) = u(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) = u(xmin,yg(ii))
    print *, 'ii=',ii,', Ua(',yg(ii),')=',Ua(ii)
  end do
  print *, ' '

  ! put solution in Ua vector
  if(ifcheck>0) then
    do i=1,nx
      x = xmin+(i-1)*hx
      do ii=1,ny
        y = ymin+(ii-1)*hy
        Ua((i-1)*ny+ii) = u(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
  end if ! ifcheck

  ! 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) = u(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) = u(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) = u(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) = u(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)

  ! 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
          if(ifcheck>1) then
            print "('Legendre integration=',f12.5,' at i=',I1,', j=',I1,', ii=',I1,', jj=',I1)", &
                    val, i, j, ii, jj
          end if
          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
      if(ifcheck>1) then
        print "('Legendre integration=',f12.5,' at i=',I1,', ii=',I1)", &
              val, i, ii
      end if
      fg((i-1)*ny+ii) = val
    end do
  end do

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

  if(ifcheck>0) then
    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)
  else ! not ifcheck
    print *, 'solution ug computed by Galerkin method'
    do i=1,nx 
      do ii=1,ny 
        print "('ug(',I1,',',I1,')=',e15.7)", i, ii, ug((i-1)*ny+ii)
      end do
    end do
  end if ! ifcheck
  print *, 'end fem_check_abc_la.f90 '

end program fem_check_abc_la