! fem_checke_la.f90  Galerkin FEM
!  solve  u'(x)=e^x  u(0)=1, u(1)=e  0 < x < 1
! initial try 10 points
! investigate gauss-legendre, adaptive, and trapezoidal integration
! solve for Uj
! sum j=0,9 int x=0,1 phip(x,j)*phi(x,i) dx = int x=0,1 f(x)*phi(x,i)dx
! k(i,j) = int x=0,1 phip(x,j)*phi(x,i) dx  (not i==0 or i==9, boundary)
! f(i)   = int x=0,1 exp(x)*phi(x,i) dx
! k * U = f, solve simultaneous equations for U

module global
  implicit none
  integer, parameter :: n=10
  double precision, dimension(n,n) :: k ! (i,j)  Gauss-Legendre integration
  double precision, dimension(n,n) :: ks ! simple trapazoidal integration
  double precision, dimension(n,n) :: kd! adaptive integration
  double precision, dimension(n) :: xg
  double precision, dimension(n) :: f
  double precision, dimension(n) :: u
  double precision, dimension(n) :: fs
  double precision, dimension(n) :: us
  double precision, dimension(n) :: fd
  double precision, dimension(n) :: ud
  double precision, dimension(n) :: Ua
  double precision :: xmin = 0.0
  double precision :: xmax = 1.0
  double precision :: h, a, b
  double precision, dimension(100) :: xx ! for Gauss-Legendre
  double precision, dimension(100) :: ww
  double precision :: val, err, avgerr, maxerr
  double precision :: eps = 1.0e-6
  integer :: p, np

  interface
    function FF(x) result(y) ! forcing function
      double precision, intent(in) :: x
      double precision :: y
    end function FF
    function uana(x) result(y) ! analytic solution and boundary
      double precision, intent(in) :: x
      double precision :: y
    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
    function aquad3e(f, xmin, xmax, eps, i, j) result(integral)
      double precision, intent(in) :: xmin, xmax, eps
      integer, intent(in) :: i,j
      double precision :: integral
      interface
        function f(x,i,j) result(y)
          double precision, intent(in) :: x
          integer, intent(in) :: i,j
          double precision :: y
        end function f
      end interface
    end function aquad3e
  end interface

end module global

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

    y = exp(x)
  end function FF

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

    y = exp(x)
  end function uana

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

    y =  phip(x,j,n,xg)*phi(x,i,n,xg)
  end function galk

  function galf(x, i, j) result(y) ! Galerkin f function
    use global
    use laphi
    implicit none
    double precision, intent(in) :: x
    integer, intent(in) :: i, j ! j unused yet needed due to galk
    double precision :: y

    y = FF(x)*phi(x,i,n,xg)
  end function galf

program fem_checke_la
  use global
  use laphi
  implicit none
  integer :: i, j

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

  print *, 'fem_checke_la.f90 running'
  print *, 'Given du/dx = exp(x)  0<x<1  Boundary u(0) = 1, u(1) = e'
  print *, 'Analytic solution u(x) = exp(x)'

  print *, 'Try 10 points'
  h = (xmax-xmin)/(n-1)
  print *, 'x grid and analytic solution'
  do i=1,n 
    xg(i) = xmin+(i-1)*h
    Ua(i) = uana(xg(i))
    print *, 'i=',i,', Ua(',xg(i),')=',Ua(i)
  end do
  print *, ' '

  ! initialize k and f
  do i=1,n 
    do j=1,n 
      k(i,j) = 0.0
      ks(i,j) = 0.0
      kd(i,j) = 0.0
    end do
    f(i) = 0.0
    fs(i) = 0.0
    fd(i) = 0.0
  end do
  ! set boundary conditions
  k(1,1) = 1.0
  f(1) = uana(xmin)        ! u(x=0.0) = 1.0
  k(n,n) = 1.0
  f(n) = uana(xmax)        ! u(x=1.0) = exp(1)
  ks(1,1) = 1.0
  fs(1) = uana(xmin)       ! u(x=0.0) = 1.0
  ks(n,n) = 1.0
  fs(n) = uana(xmax)       ! u(x=1.0) = exp(1)
  kd(1,1) = 1.0
  fd(1) = uana(xmin)       ! u(x=0.0) = 1.0
  kd(n,n) = 1.0
  fd(n) = uana(xmax)       ! u(x=1.0) = exp(1)

  np = 98
  print *, 'calling gaulegf np=',np
  call gaulegf(xmin, xmax, xx, ww, np) ! xx and ww set up for integrations
  h = (xmax-xmin)/np ! for trapezoidal integration

  ! compute entries=stiffness matrix
  print *, 'compute stiffness matrix'
  do i=2,n-1 
    do j=1,n 
      ! compute integral for k(i,j)
      val = 0.0
      do p=1,np 
        val = val + ww(p)*galk(xx(p),i,j)
      end do
      print *, 'Legendre integration=',val,', at i=',i,', j=',j
      k(i,j) = val

      val = aquad3e(galk, xmin, xmax, eps, i, j)
      print *, 'adaptive integration=',val,', at i=',i,', j=',j
      kd(i,j) = val

      val = 0.0
      do p=1,np-1 
        val = val + galk(xmin+(p)*h,i,j)*h
      end do
      val = val + (galk(xmin,i,j)+ galk(xmax,i,j))*h/2.0
      print *, 'simple   integration=',val,', at i=',i,', j=',j
      ks(i,j) = val
    end do ! j

    ! compute integral for f(i)
    val = 0.0
    do p=1,np 
      val = val + ww(p)*galf(xx(p),i,0)
    end do
    print *, 'Legendre integration=',val,', i=',i
    f(i) = val

    val = aquad3e(galf, xmin, xmax, eps, i, 0)
    print *, 'adaptive integration=',val,', i=',i
    fd(i) = val

    val = 0.0
    do p=1,np-1 
      val = val + galf(xmin+(p)*h,i,0)*h
    end do
    val = val + (galf(xmin,i,0) + galf(xmax,i,0))*h/2.0
    print *, 'simple   integration=',val,', i=',i
    fs(i) = val
  end do ! i

  print *, 'k computed stiffness matrix'
  do i=1,n 
    do j=1,n 
      print *, 'i=',i,', j=',j,', k(i,j)=',k(i,j)
    end do
  end do
  print *, ' '
  print *, 'f computed forcing function'
  do i=1,n 
    print *, 'f(',i,')=',f(i)
  end do
  print *, ' '
  call simeq(n, n, k, f, u)

  avgerr = 0.0
  maxerr = 0.0
  print *, 'u computed Galerkin, Ua analytic, error'
  do i=1,n 
    err = abs(u(i)-Ua(i))
    if(err>maxerr) then
      maxerr = err
    end if
    avgerr = avgerr + err
    print *, 'u(',i,')=',u(i),', Ua=',Ua(i),', err=',(u(i)-Ua(i))
  end do
  print *, '                   maxerr=',maxerr,', avgerr=',avgerr/n
  print *, ' '



  print *, 'kd computed adaptive stiffness matrix'
  do i=1,n 
    do j=1,n 
      print *, 'i=',i,', j=',j,', kd(i,j)=',kd(i,j)
    end do
  end do
  print *, ' '
  print *, 'fd computed forcing function'
  do i=1,n 
    print *, 'fd(',i,')=',fd(i)
  end do
  print *, ' '
  call simeq(n, n, kd, fd, ud)

  avgerr = 0.0
  maxerr = 0.0
  print *, 'ud computed Galerkin, Ua analytic, error'
  do i=1,n 
    err = abs(ud(i)-Ua(i))
    if(err>maxerr) then
      maxerr = err
    end if
    avgerr = avgerr + err
    print *, 'ud(',i,')=',ud(i),', Ua=',Ua(i),', err=',(ud(i)-Ua(i))
  end do
  print *, '                   maxerr=',maxerr,', avgerr=',avgerr/n
  print *, ' '

  print *, 'ks computed simple stiffness matrix'
  do i=1,n 
    do j=1,n 
      print *, 'i=',i,', j=',j,', ks(i,j)=',ks(i,j)
    end do
  end do
  print *, ' '
  print *, 'fs computed forcing function'
  do i=1,n 
    print *, 'fs(',i,')=',fs(i)
  end do
  print *, ' '
  call simeq(n, n, ks, fs, us)

  avgerr = 0.0
  maxerr = 0.0
  print *, 'us computed Galerkin, Ua analytic, error'
  do i=1,n 
    err = abs(us(i)-Ua(i))
    if(err>maxerr) then
      maxerr = err
    end if
    avgerr = avgerr + err
    print *, 'us(',i,')=',us(i),', Ua=',Ua(i),', err=',(us(i)-Ua(i))
  end do
  print *, '                   maxerr=',maxerr,', avgerr=',avgerr/n
  print *, ' '

end program fem_checke_la