!------------------------------------------------------------------
!    FVCA5 Benchmark on anisotropic diffusion problems

! Fortran subroutines for the diffusion tensors, right hand sides and 
! exact solutions (when available). 
!------------------------------------------------------------------

!------------------------------------------------------------------
!  TEST 1.1 regular solution1_1 
!------------------------------------------------------------------
subroutine  solution1_1(x,y,kxx,kxy,kyy,u,ux,uy,f)

  double precision  x,y,kxx,kxy,kyy,u,ux,uy,f,uxx,uxy,uyy

  !     input:   x,y
  !     output : kxx,kxy,kyy,u,ux,uy,f 

  u  = x*(1-x)*y*(1-y)*16. 

  ux = (-2.*x+1)*y*(1-y)*16. 
  uy = (-2.*y+1)*x*(1-x)*16. 

  uxx = -2.*y*(1-y)*16.
  uxy = (-2.*x+1)*(-2.*y+1)*16. 
  uyy = -2.*x*(1-x)*16. 

  kxx = 1.5 
  kxy = .5 
  kyy = 1.5 

  f = -(kxx*uxx + 2.*kxy*uxy + kyy*uyy)  

end  subroutine  solution1_1

!------------------------------------------------------------------
!  TEST 1.2  increasing solution
!------------------------------------------------------------------
subroutine  solution1_2(delta,x,y,kxx,kxy,kyy,u,ux,uy,f)

  double precision   delta,x,y,kxx,kxy,kyy,u,ux,uy,uxy,f,x1,y1 

 !      input :  delta,x,y 
 !      output   :  kxx,kxy,kyy,u,ux,uy,f 

  x1 = 1.-x
  y1 = 1.-y
  
  u  = sin(x1*y1) + x1**3. * y1**2.

  ux = -y1*cos(x1*y1) - 3.* (x1*y1)**2.
  uy = -x1*cos(x1*y1) - 2.*y1* x1**3.

  uxx = -y1*y1*sin(x1*y1) + 6.*x1* y1**2.
  uxy = -x1*y1*sin(x1*y1) + cos(x1*y1) + 6.*y1* x1**2.
  uyy = -x1*x1*sin(x1*y1) + 2.* x1**3.

  kxx = 1.5
  kxy = .5
  kyy = 1.5

  f = -(kxx*uxx + 2.*kxy*uxy+ kyy*uyy)

end subroutine solution1_2

!------------------------------------------------------------------
!  TEST 2  numerical locking
!------------------------------------------------------------------

subroutine solution2(delta,x,y,kxx,kxy,kyy,u,ux,uy,f) 
  double  precision   delta,x,y,kxx,kxy,kyy,u,ux,uy,f,pi,x1,x2 

  !      input :  delta,x,y 
  !      output:  kxx,kxy,kyy,u,ux,uy,f 

  pi = 4.*atan(1.) 
  x1 = 2.*pi  
  x2 = x1/sqrt(delta) 
  u  = sin(x1*x)*exp(-x2*y) 
  ux = x1*cos(x1*x)*exp(-x2*y) 
  uy = -x2*u 
  kxx = 1. 
  kxy = 0. 
  kyy = delta 

  f = 0.  

end subroutine  solution2


!------------------------------------------------------------------
!  TEST 3 oblique flow
!------------------------------------------------------------------
subroutine solution3(delta,kxx,kxy,kyy,f) 
  double  precision   delta,kxx,kxy,kyy,f,pi,cost,sint

  !      input :  delta 
  !      output:  kxx,kxy,kyy,f 

  pi = 4.*atan(1.) 
  cost = cos(40.*pi/180.) 
  sint = sqrt(1. - cost*cost) 

  kxx = cost*cost+delta*sint*sint  
  kxy = cost*sint*(1.-delta) 
  kyy = sint*sint+delta*cost*cost 

  f = 0.  

end subroutine solution3

!------------------------------------------------------------------
!  TEST 4 Vertical fault
!------------------------------------------------------------------


subroutine solution4(x,y,kxx,kxy,kyy,f)  

  double precision x,y,kxx,kxy,kyy,f
  integer intx, inty

  !      input  :  x,y 
  !      output :  kxx,kxy,kyy,f 

  if (x<.5) then
     inty = int(10. * (y+.15))
  !  inty = int(10. * (y-.05))
     if (inty - 2*(inty/2) ==0) then
        kxx = 100.
        kxy = 0.
        kyy = 10.
     else
        kxx = .01
        kxy = 0.
        kyy = .001
     endif
  else
     inty = int(10. * y)
     if (inty - 2*(inty/2) ==0) then
        kxx = 100.
        kxy = 0.
        kyy = 10.
     else
        kxx = .01
        kxy = 0.
        kyy = .001
     endif
  endif

  f = 0.

end subroutine solution4

!------------------------------------------------------------------
!  TEST 5 Heterogeneous rotating anisotropy
!------------------------------------------------------------------

subroutine solution5(delta,x,y,kxx,kxy,kyy,u,ux,uy,f) 
  double  precision   delta,x,y,kxx,kxy,kyy,u,ux,uy,f,pi,rt,f0

  !     input   :  delta,x,y 
  !     output  :  kxx,kxy,kyy,u,ux,uy,f 

  pi =  4.*atan(1.) 

  u  = sin(pi*x)*sin(pi*y) 

  ux =  pi * cos(pi*x)*sin(pi*y)  
  uy =  pi * cos(pi*y)*sin(pi*x) 

  rt =  x*x+y*y 

   kxx =  (x*x+delta*y*y)/rt 
   kxy =  -(1-delta)*x*y/rt  
   kyy =  (x*x+delta*y*y)/rt 

  f0 =   sin(pi*x)*sin(pi*y)*pi*pi*(1+delta)*(x*x+y*y)  &
       + cos(pi*x)*sin(pi*y)*pi*(1.-3.*delta)*x   &
       + cos(pi*y)*sin(pi*x)*pi*(1.-3.*delta)*y  & 
       + cos(pi*y)*cos(pi*x)*2.*pi*pi*(1.-delta)*x*y    

  f =  (f0+2.*(x*(kxx*ux+kxy*uy)+y*(kxy*ux+kyy*uy)))/rt  

end  subroutine  solution5

!------------------------------------------------------------------
!  TEST 6 Oblique drain
!------------------------------------------------------------------

subroutine  solution6(delta,x,y,kxx,kxy,kyy,u,ux,uy,f) 
  double  precision  delta,x,y,kxx,kxy,kyy,u,ux,uy,f,phi1,phi2,alpha,beta  

  !      input  :  delta,x,y 
  !      output :  kxx,kxy,kyy,u,ux,uy,f 
  u  = - x - y * delta 
  ux = - 1. 
  uy = - delta 
  phi1 = y - delta * (x - .5) - .475 
  phi2 = phi1 - .05 
  if (phi1<0  .or. phi2>0) then 
     alpha = 1. 
     beta = .1
  else 
     alpha = 100.
     beta = 10.
  endif
  cost = 1./sqrt(1.+delta*delta)
  sint = delta*cost
  kxx = alpha*cost*cost+beta*sint*sint
  kxy = cost*sint*(alpha-beta)
  kyy = alpha*sint*sint+beta*cost*cost

  f = 0.  

end subroutine solution6

!------------------------------------------------------------------
!  TEST 7 Oblique barrier
!------------------------------------------------------------------

subroutine solution7(delta,x,y,kxx,kxy,kyy,u,ux,uy,f)  
  double precision delta,x,y,kxx,kxy,kyy,u,ux,uy,f,phi1,phi2   

  !      input  :  delta,x,y 
  !      output :  kxx,kxy,kyy,u,ux,uy,f 


  phi1 = y - delta * (x - .5) - .475 
  phi2 = phi1 - .05  
  if (phi1<0) then  
     u  = -phi1  
     ux = delta 
     uy = -1  
     kxx = 1. 
     kxy = 0. 
     kyy = 1.  
  else if (phi2<0) then   
     u  = - phi1/.01   
     ux = delta/.01 
     uy = - 1./.01  
     kxx = .01 
     kxy = 0. 
     kyy = .01  
  else   
     u  = - phi2 - 5.  
     ux = delta 
     uy = - 1.  
     kxx = 1. 
     kxy = 0. 
     kyy = 1. 
  endif

  f = 0.    

end subroutine solution7