! fem_check33_la.f90  Galerkin FEM
! solve uxxx(x,y,z) + 2 uyyy(x,y,z) + 3 uzzz(x,y,z) +
!     4 uxxy(x,y,z) + 5 uxxz(x,y,z) + 6 uyyx(x,y,z) +
!     7 uyyz(x,y,z) + 8 uzzx(x,y,z) + 9 uzzy(x,y,z) +
!     10 uxx(x,y,z) + 11 uyy(x,y,z) + 12 uzz(x,y,z) +
!     13 uxy(x,y,z) + 14 uxz(x,y,z) + 15 uyz(x,y,z) +
!     16  ux(x,y,z) + 17  uy(x,y,z) + 18  uz(x,y,z) +
!     19   u(x,y,z) = F(x,y,z)
! F(x,y,z) = 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z +
!            136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y +
!            80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y +
!            224.0*x*y*y*z + 76.0*x*x*y*y*z*z + 144.0*x*x*y +
!            80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y +
!            112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z +
!            19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x +
!            136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z
!
! boundary conditions computed using u(x,y,z)
! analytic solution may be given by u(x,y,z) = x**4 + 2 y**4 + 3 z**4 +
!                   4 x**2 y**2 z**2 + 5 x y + 6 x z + 7 y z + 8
! gauss-legendre integration used
! solve for Uijk numerical approximation U(x_i,y_j,z_k) of u(x_i,y_j,z_k)
! kg * ug = fg, solve simultaneous equations for U

module global
  implicit none
  integer, parameter :: nx=5
  integer, parameter :: ny=5
  integer, parameter :: nz=5
  double precision, dimension(nx*ny*nz,nx*ny*nz) :: kg 
                 ! (i*ny*nz+ii*nz+iii)*ny*ny*nz + (j*ny*nz+jj*nz+jjj)
  double precision, dimension(nx) :: xg
  double precision, dimension(ny) :: yg
  double precision, dimension(nz) :: zg
  double precision, dimension(nx*ny*nz) :: fg
  double precision, dimension(nx*ny*nz) :: ug
  double precision, dimension(nx*ny*nz) :: Ua
  double precision :: xmin = 0.0
  double precision :: xmax = 1.0
  double precision :: ymin = 0.0
  double precision :: ymax = 1.0
  double precision :: zmin = 0.0
  double precision :: zmax = 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, dimension(100) :: zz
  double precision, dimension(100) :: wz
  double precision :: val, err, avgerr, maxerr
  double precision :: eps = 1.0e-6
  integer :: px, py, pz, npx, npy, npz

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

  end interface

end module global

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

  v = 406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z +      &
      136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y +     &
      80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y +             &
      224.0*x*y*y*z + 76.0*x*x*y*y*z*z + 144.0*x*x*y +            &
      80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y +         &
      112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z + &
      19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x +   &
      136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z
  end function FF

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

    v = x*x*x*x + 2.0*y*y*y*y + 3.0*z*z*z*z + 4.0*x*x*y*y*z*z + &
        5.0*x*y + 6.0*x*z + 7.0*y*z + 8.0
  end function uana

  function galk(x,y,z,  i,j, ii,jj, iii,jjj) result(v)
                     ! Galerkin k stiffness function for this specific problem
    use global
    use laphi
    implicit none
    double precision, intent(in) :: x,y,z
    integer, intent(in) :: i, j, ii, jj, iii, jjj
    double precision :: v
    v =                                                              &
    (1.0*phippp(x,j,nx,xg)*   phi(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
     2.0*   phi(x,j,nx,xg)*phippp(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
     3.0*   phi(x,j,nx,xg)*   phi(y,jj,ny,yg)* phippp(z,jjj,nz,zg) + &
     4.0* phipp(x,j,nx,xg)*  phip(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
     5.0* phipp(x,j,nx,xg)*   phi(y,jj,ny,yg)*   phip(z,jjj,nz,zg) + &
     6.0*  phip(x,j,nx,xg)* phipp(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
     7.0*   phi(x,j,nx,xg)* phipp(y,jj,ny,yg)*   phip(z,jjj,nz,zg) + &
     8.0*  phip(x,j,nx,xg)*   phi(y,jj,ny,yg)*  phipp(z,jjj,nz,zg) + &
     9.0*   phi(x,j,nx,xg)*  phip(y,jj,ny,yg)*  phipp(z,jjj,nz,zg) + &
    10.0* phipp(x,j,nx,xg)*   phi(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
    11.0*   phi(x,j,nx,xg)* phipp(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
    12.0*   phi(x,j,nx,xg)*   phi(y,jj,ny,yg)*  phipp(z,jjj,nz,zg) + &
    13.0*  phip(x,j,nx,xg)*  phip(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
    14.0*  phip(x,j,nx,xg)*   phi(y,jj,ny,yg)*   phip(z,jjj,nz,zg) + &
    15.0*   phi(x,j,nx,xg)*  phip(y,jj,ny,yg)*   phip(z,jjj,nz,zg) + &
    16.0*  phip(x,j,nx,xg)*   phi(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
    17.0*   phi(x,j,nx,xg)*  phip(y,jj,ny,yg)*    phi(z,jjj,nz,zg) + &
    18.0*   phi(x,j,nx,xg)*   phi(y,jj,ny,yg)*   phip(z,jjj,nz,zg) + &
    19.0*   phi(x,j,nx,xg)*   phi(y,jj,ny,yg)*    phi(z,jjj,nz,zg))* &
           (phi(x,i,nx,xg)*   phi(y,ii,ny,yg)*    phi(z,iii,nz,zg))
  end function galk

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

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

  function s(i, ii, iii) result(v)
    use global
    implicit none
    integer, intent(in) :: i, ii, iii
    integer :: v

    v = (i-1)*ny*nz+(ii-1)*nz+iii
  end function s


program fem_check33_la
  use global
  use laphi
  implicit none
  integer :: i, j, ii, jj, iii, jjj
  double precision :: x, y, z, hx, hy, hz

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

  print *, 'fem_check33_la.f90 running'
  print *, 'solve uxxx(x,y,z) + 2 uyyy(x,y,z) + 3 uzzz(x,y,z) +'
  print *,'     4 uxxy(x,y,z) + 5 uxxz(x,y,z) + 6 uyyx(x,y,z) +'
  print *,'     7 uyyz(x,y,z) + 8 uzzx(x,y,z) + 9 uzzy(x,y,z) +'
  print *,'     10 uxx(x,y,z) + 11 uyy(x,y,z) + 12 uzz(x,y,z) +'
  print *,'     13 uxy(x,y,z) + 14 uxz(x,y,z) + 15 uyz(x,y,z) +'
  print *,'     16  ux(x,y,z) + 17  uy(x,y,z) + 18  uz(x,y,z) +'
  print *,'     19   u(x,y,z) = F(x,y,z)'
  print *,' F(x,y,z) = '
  print *,'    406.0 + 208.0*x*y*z*z + 240.0*x*x*y*z + 88.0*x*x*z*z +'
  print *,'    136.0*x*x*y*z*z + 128.0*x*y*y*z*z + 217.0*x + 302.0*y +'
  print *,'    80.0*y*y*z*z + 144.0*x*x*y*y*z + 96.0*x*x*y*y +'
  print *,'    224.0*x*y*y*z + 76.0*x*x*y*y*z*z + 144.0*x*x*y +'
  print *,'    80.0*y*y*z + 64.0*y*z*z + 133.0*y*z + 128.0*x*y*y +'
  print *,'    112.0*x*x*z + 96.0*x*z*z + 114.0*x*z + 95.0*x*y + 431.0*z +'
  print *,'    19.0*x*x*x*x + 38.0*y*y*y*y + 57.0*z*z*z*z + 64.0*x*x*x +'
  print *,'    136.0*y*y*y + 216.0*z*z*z + 120.0*x*x + 264.0*y*y + 432.0*z*z'
  print *, ' '
  print *, 'boundary conditions computed using u(x,y,z)'
  print *, 'analytic solution may be given by u(x,y,z) = '
  print *, '  x**4 + 2 y**4 + 3 z**4 + 4 x**2 y**2 z**2 +' 
  print *, '  5 x y + 6 x z + 7 y z + 8'
  print *, 'Gauss-Legendre integration used'
  print *, ' '
  print *, 'xmin=',xmin,', ymin =',ymin,', zmin=',zmin
  print *, 'xmax=',xmax,', ymax =',ymax,', zmax=',zmax
  print *, 'nx=',nx,', ny =',ny,', nz =',nz

  hx = (xmax-xmin)/(nx-1) ! does not have to be uniform spacing
  print *, 'x grid and analytic solution, at ymin,zmin'
  do i=1,nx 
    xg(i) = xmin+(i-1)*hx
    Ua(i) = uana(xg(i),ymin,zmin) ! not saved
    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,zmin'
  do ii=1,ny 
    yg(ii) = ymin+(ii-1)*hy
    Ua(ii) = uana(xmin,yg(ii),zmin)
    print *, 'ii=',ii,', Ua(',yg(ii),')=',Ua(ii)
  end do
  print *, ' '

  hz = (zmax-zmin)/(nz-1) ! does not have to be uniform spacing
  print *, 'z grid and analytic solution, at xmin,ymin'
  do iii=1,nz 
    zg(iii) = zmin+(iii-1)*hz
    Ua(iii) = uana(xmin,ymin,zg(iii))
    print *, 'iii=',iii,', Ua(',zg(iii),')=',Ua(iii)
  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
      do iii=1,nz
        z = zmin+(iii-1)*hz
        Ua(s(i,ii,iii)) = uana(x,y,z)
        ! print "('solution at i=',I1,', x=',f6.3,', ii=',I1,', y=',f6.3,', iii=',I1,' z=',f7.3,' is ',f7.3)",&
        !      i, x, ii, y, iii, z, Ua(s(i,ii,iii))
      end do ! iii
    end do ! ii
  end do ! i

  ! initialize k and f
  do i=1,nx
    do j=1,nx 
      do ii=1,ny 
        do jj=1,ny 
          do iii=1,nz 
            do jjj=1,nz 
              kg(s(i,ii,iii),s(j,jj,jjj)) = 0.0
            end do ! jjj
          end do ! iii
        end do ! jj
      end do ! ii
    end do ! j
  end do ! i

  do i=1,nx
    do ii=1,ny 
      do iii=1,nz 
        fg(s(i,ii,iii)) = 0.0
      end do ! iii
    end do ! ii
  end do ! i

  ! set boundary conditions, overkill
  do i=1,nx
    x = xg(i)
    do ii=1,ny
      y = yg(ii)
      do iii=1,nz
        z = zg(iii)
        kg(s(i,ii,iii),s(i,ii,iii)) = 1.0
        fg(s(i,ii,iii)) = uana(x,y,z)
      end do ! iii
    end do ! ii
  end do ! i


  npx = 12
  npy = 12
  npz = 12
  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 *, 'calling gaulegf npz=',npz
  call gaulegf(zmin, zmax, zz, wz, npz)
  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)
  print *, 'zz(1)=',zz(1),', zz(2)=',zz(2),', wz(1)=',wz(1),', wz(2)=',wz(2)
  val = galk(xx(3),yy(3),zz(3),2,2,2,2,2,2)
  print *, 'galk(xx(3),yy(3),zz(3),2,2,2,2,2,2)=',val
  val = galf(xx(3),yy(3),zz(3),2,2,2)
  print *, 'galf(xx(3),yy(3),zz(3),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
          do iii=2,nz-1 
            do jjj=1,nz
              ! compute integral for k(i,j)
              val = 0.0
              do px=1,npx 
                do py=1,npy 
                  do pz=1,npz 
                    val = val + wx(px)*wy(py)*wz(pz)* &
                          galk(xx(px),yy(py),zz(pz), &
                          i,j,ii,jj,iii,jjj)
                  end do !pz
                end do ! py
              end do ! px
              kg(s(i,ii,iii),s(j,jj,jjj)) = val
            end do ! jjj
          end do ! iii
        end do ! jj
      end do ! ii
      print "('Legendre integration=',f12.5,' at i=',I1,', j=',I1,', ii=',I1,', jj=',I1,', iii=',I1,',jjj=',I1)", &
             val, i, j, ii, jj, iii, jjj
    end do ! j
  end do ! i

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

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

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

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

end program fem_check33_la