! fem_checke_la.f90  Galerkin FEM
! solve x^3 u'''' + x^2 u''' + u'' + u' + u = x^4 + 64 x^3 - x^2 - 4x + 1
! analytic solution u(x) = x^4 - x^2 + 1
! Gauss-Legendre integration
! gauss-Legendre
! 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) :: Ua
  double precision :: xmin = 0.0
  double precision :: xmax = 1.0
  double precision :: h
  integer, parameter :: np = 24
  double precision, dimension(np) :: xx ! for Gauss-Legendre
  double precision, dimension(np) :: ww
  integer :: p
  double precision :: val, err, avgerr, maxerr

  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
  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 = x**4 + 64.0*x**3 - x**2 - 4.0*x + 1.0
  end function FF

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

    y = x**4 - x**2 + 1.0
  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 = (x**3 *phipppp(x,j,n,xg) + &
         x**2 * phippp(x,j,n,xg) + &
         x*      phipp(x,j,n,xg) + &
                  phip(x,j,n,xg) + &
                   phi(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 x^3 u'''' + x^2 u''' + u'' + u' + u = x^4 + 64 x^3 - x^2 - 4x + 1"
  print *, '0 <= x <= 1  Boundary u(0) = 1, u(1) = 1'
  print *, 'Analytic solution u(x) = x^4 - x^2 + 1'

  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
    end do
    f(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) = 1.0

  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
    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
  end do ! i

  print *, 'k computed stiffness matrix, see above'
  print *, 'f computed forcing function, see above'
  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(',I2,')=',f7.4,', Ua(i)=',f7.4,', err=',e15.7)", & 
          i, u(i), Ua(i), (u(i)-Ua(i))
  end do
  print *, '                   maxerr=',maxerr,', avgerr=',avgerr/n
  print *, ' '

end program fem_checke_la