! fem_check44_la.f90  Galerkin FEM
! solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
!     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
!     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
!     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
!     F(x,y,z,t)
!     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
!                  180 z^2*t + 200 y^2z*t + 44 t^3 +
!                  110 y^2z^2t + 66 xy + 12t^4 +
!                  24 z^4 + 36 y^4 + 48 x^4 +
!                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t 
!
! boundary conditions computed using u(x,y,z,t)
! analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
!                                   5 y^2 z^2 t^2 + 6 x y t + 
!                                   7 x z + 8 t + 9
! gauss-legendre integration used
! solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
! kg * ug = fg, solve simultaneous equations for U

module global
  implicit none
  integer, parameter :: nx=5
  integer, parameter :: ny=5
  integer, parameter :: nz=5
  integer, parameter :: nt=5
  double precision, dimension(nx*ny*nz*nt,nx*ny*nz*nt) :: kg 
     ! (i*ny*nz*nt+ii*nz*nt+iii*nz+iiii),(j*ny*nz*nt+jj*nz*nt+jjj*nz+jjjj)
  double precision, dimension(nx) :: xg
  double precision, dimension(ny) :: yg
  double precision, dimension(nz) :: zg
  double precision, dimension(nt) :: tg
  double precision, dimension(nx*ny*nz*nt) :: fg
  double precision, dimension(nx*ny*nz*nt) :: ug
  double precision, dimension(nx*ny*nz*nt) :: 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 :: tmin = 0.0
  double precision :: tmax = 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, dimension(100) :: tt
  double precision, dimension(100) :: wt
  double precision :: val, err, avgerr, maxerr
  integer :: px, py, pz, pt, npx, npy, npz, npt

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

  end interface

end module global

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

    v= 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t + &
       200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y + &
       12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x + &
       60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t
  end function FF

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

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

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

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

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

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

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


program fem_check44_la
  use global
  use laphi
  implicit none
  integer :: i, j, ii, jj, iii, jjj, iiii, jjjj
  double precision :: x, y, z, t, hx, hy, hz, ht

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

  print *, 'fem_check44_la.f90 running'
  print *,' solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +'
  print *,'     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +'
  print *,'     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +'
  print *,'     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) ='
  print *,'     F(x,y,z,t)'
  print *,'     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +'
  print *,'                  180 z^2*t + 200 y^2z*t + 44 t^3 +'
  print *,'                  110 y^2z^2t + 66 xy + 12t^4 +'
  print *,'                  24 z^4 + 36 y^4 + 48 x^4 +'
  print *,'                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t' 
  print *,' '
  print *,' boundary conditions computed using u(x,y,z,t)'
  print *,' analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +'
  print *,'                                   5 y^2 z^2 t^2 + 6 x y t + '
  print *,'                                   7 x z + 8 t + 9'
  print *, 'Gauss-Legendre integration used'
  print *, ' '
  print *, 'xmin=',xmin,', ymin =',ymin,', zmin=',zmin,', tmin=',tmin
  print *, 'xmax=',xmax,', ymax =',ymax,', zmax=',zmax,', tmax=',tmax
  print *, 'nx=',nx,', ny =',ny,', nz =',nz,', nt=',nt 

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

  ht = (tmax-tmin)/(nt-1) ! does not have to be uniform spacing
  print *, 't grid and analytic solution, at xmin,ymin,zmin'
  do iiii=1,nt 
    tg(iiii) = tmin+(iiii-1)*ht
    Ua(iiii) = uana(xmin,ymin,zmin,tg(iiii))
    print *, 'iiii=',iiii,', Ua(',tg(iiii),')=',Ua(iiii)
  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
        do iiii=1,nt
          t = tmin+(iiii-1)*ht
          Ua(s(i,ii,iii,iiii)) = uana(x,y,z,t)
        end do ! iiii
      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 
              do iiii=1,nt 
                do jjjj=1,nt 
                  kg(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj)) = 0.0
                end do ! jjjj
              end do ! iiii
            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 
        do iiii=1,nt 
          fg(s(i,ii,iii,iiii)) = 0.0
        end do ! iiii
      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)
        do iiii=1,nt
          t = tg(iiii)
          kg(s(i,ii,iii,iiii),s(i,ii,iii,iiii)) = 1.0
          val = uana(x,y,z,t)
          fg(s(i,ii,iii,iiii)) = val
        end do ! iiii
      end do ! iii
    end do ! ii
    print "('boundary i=',I1,', x=',f6.3,', ii=',I1,', y=',f6.3,', iii=',I1,', z=',f6.3,', iiii=',I1,', t=',f6.3,' is ',f10.5)", &
            i, x, ii, y, iii, z, iiii, t, val
  end do ! i


  npx = 12
  npy = 12
  npz = 12
  npt = 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 *, 'calling gaulegf npt=',npt
  call gaulegf(tmin, tmax, tt, wt, npt)
  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)
  print *, 'tt(1)=',tt(1),', tt(2)=',tt(2),', wt(1)=',wt(1),', wt(2)=',wt(2)
  val = galk(xx(3),yy(3),zz(3),tt(3),2,2,2,2,2,2,2,2) ! one more than C
  print *, 'galk(xx(3),yy(3),zz(3),tt(3),2,2,2,2,2,2,2,2)=',val
  val = galf(xx(3),yy(3),zz(3),tt(3),2,2,2,2)
  print *, 'galf(xx(3),yy(3),zz(3),tt(3),2,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
              do iiii=2,nt-1 
                do jjjj=1,nt
                  ! compute integral for k(i,j)
                  val = 0.0
                  do px=1,npx 
                    do py=1,npy 
                      do pz=1,npz 
                        do pt=1,npt 
                          val = val + wx(px)*wy(py)*wz(pz)*wt(pt)* &
                                      galk(xx(px),yy(py),zz(pz),tt(pt), &
                                      i,j,ii,jj,iii,jjj,iiii,jjjj)
                        end do ! pt
                      end do ! pz
                    end do ! py
                  end do ! px
                  kg(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj)) = val
                end do ! jjjj
              end do ! iiii
            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,',iiii=',I1,', jjjj=',I1)"&
            , val, i, j, ii, jj, iii, jjj, iiii, jjjj
    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 
        do iiii=2,nt-1 
          val = 0.0
          do px=1,npx 
            do py=1,npy 
              do pz=1,npz 
                do pt=1,npt 
                  val = val + wx(px)*wy(py)*wz(pz)*wt(pt)* &
                              galf(xx(px),yy(py),zz(pz),tt(pt),i,ii,iii,iiii)
                end do ! pt
              end do ! pz
            end do ! py
          end do ! px
          fg(s(i,ii,iii,iiii)) = val
        end do ! iiii
      end do ! iii
    end do ! ii
    print "('Legendre integration=',f12.5,' F at i=',I1,', ii=',I1,', iii=',I1,', iiii=',I1)", &
            val, i, ii, iii, iiii
  end do ! i

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

  call simeq(nx*ny*nz*nt, nx*ny*nz*nt, 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 
        do iiii=1,nt 
          err = abs(ug(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii)))
          if(err>maxerr) then
            maxerr = err
          end if
          avgerr = avgerr + err
          print "('ug(',I1,',',I1,','I1,',',I1')=',f7.4,', Ua=',f7.4,', err=',e15.7)", &
              i, ii, iii, iiii, ug(s(i,ii,iii,iiii)), Ua(s(i,ii,iii,iiii)), &
              (ug(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii)))
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i
  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nx*ny*nz*nt)
  print *, 'end fem_check44_la.f90 '

end program fem_check44_la