! abc_la.f90  User problem definitions for a specific PDE to be solved
!          this includes the user choice of degree and discretization
!
!          The problem must be well posed and the solution must be
!          continuous and continuously differentiable
!
! The general PDE to be solved is:
!  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)
!

module abc_la

! Define the region for the solution:
! xmin <= x <= xmax, ymin <= y <= ymax  using  nx by ny points

  double precision, parameter :: xmin =  -1.0
  double precision, parameter :: xmax =   1.0
  double precision, parameter :: ymin =  -1.0
  double precision, parameter :: ymax =   1.0
  integer, parameter :: nx = 6
  integer, parameter :: ny = 5
  integer, parameter :: ifcheck = 1  !set ifcheck = 1 or more, else ifcheck = 0

  public :: a1
  public :: b1
  public :: c1
  public :: d1
  public :: e1
  public :: f1
  public :: u
  public :: c

contains

! the user must provide all functions, even if only  return 0.0
! A specific set of test functions is covered in abc.txt

  function a1(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = exp(x/2.0)*exp(y)/2.0
  end function a1

  function b1(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = 0.7/(x*x*y*y+0.5)
  end function b1

  function c1(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = (4.0 - exp(x) - exp(y/2.0))*2.0
  end function c1

  function d1(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = x*x+y
  end function d1

  function e1(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = x*y*y
  end function e1

  function f1(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = 3.0*x + 2.0*y
  end function f1

  function u(x, y) result(z)
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z
                     ! must compute boundary, if solution, then set ifcheck = 1
    z = x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*x*y*y + 5.0*x*y + &
        6.0*x + 7.0*y + 8.0
  end function u

  function c(x, y) result(z)  ! this becomes defined by u and a1...f1
    implicit none
    double precision, intent(in) :: x, y
    double precision :: z

    z = 0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)                  +  &
        0.7*(6.0*x + 8.0*y + 5.0)/(x*x*y*y+0.5)              +  &
        (8.0 - 2.0*exp(x) - 2.0*exp(y/2.0))*(12.0*y + 8.0*x) +  &
        (x*x+y)*(3.0*x*x + 6.0*x*y + 4.0*y*y + 5.0*y + 6.0)  +  &
        x*y*y*(6.0*y*y + 3.0*x*x + 8.0*x*y + 5.0*x +7.0)     +  &
        (3.0*x + 2.0*Y)*(x*x*x + 2.0*y*y*y + 3.0*x*x*y +        &
         4.0*x*y*y + 5.0*x*y + 6.0*x + 7.0*y + 8.0)
  end function c

end module abc_la