module weight_fun_mod
   use paramets

implicit none

public :: weightFun ! use in target functional

contains

!> weight function used in target functional
!> i.e. j_\om, j_D or j_N depending on the type of the target functional
!> used only in boundary now, in SCALAR case only
function weightFun( xi, i, Re1) result (w)
   real, dimension(1:nbDim), intent(in) :: xi
   integer, intent(in) :: i
   real, intent(in) :: Re1
   real :: w
   real :: pi
   real :: eps
   real :: diff, conv

   ! tolerance around zero
   eps = 1.E-12


   select case(i)
   case(1) ! constant 1
      w = 1.0
      ! TEST ONLY
      !w = ( (xi(1) - 2.5) *(3.5 - xi(1) )* (xi(2) - 2.5) *(3.5 - xi(2) )) ** 0.5
      ! Ajay
      !w = exp(-1000 * sqrt( (xi(1) - 1.3)**2 + (xi(2) - 0.3)**2 ) )
      !w = exp(-1000 * sqrt( max((xi(1) - 1.3)**2, 0.01) - 0.01  &
      !     + max((xi(2) - 0.3)**2, 0.01 )-0.01 ) )
      !if(w < 1E-20) w = 0.
      
      !if(w > 1E-10) print*,'#####', xi(1:2), w
   case(2) ! y
      w = xi(2)
   case(3) ! y * sin( (y-0.5)*4pi)
      pi = asin(1.0)*2.0
      w = xi(2) * sin( (xi(2) - 0.5)*4.*pi )
   case(4)
   ! RHS for Laplace problem with BC is zero Dirichlet and subdomain is the whole Omega
   !such that the solution z = 16x(1-x)y(1-y),
      w = 32.0* ( xi(1)*(1.0-xi(1)) + xi(2)*(1.0-xi(2)) )
   case(5) ! constant -1
      w = -1.0
   case(6) ! specially for icase=79 Carpio 2013, Lshaped domain
      ! j = (b*n) = (y,-x)*n
      ! n = (1,0)    on {4}*(0,4) right edge
      ! n = (0,-1)   on (2,4)*{0} bottom edge

      if ( xi(1) >= 4.0 - eps ) then ! right boundary
        w = xi(2)
      elseif ( xi(2) <= 0.0 + eps ) then ! bottom
        w = xi(1)
        !print*, 'weightFun on the bottom'
        !w = 0.
        !write(*, *) '^%E#^FE^#Y', xi(1:2), w
      else
        print*, 'Problem in weightFun for Boundary functional, xi=' , xi
        w = 0.0
     endif
     !write(38, *) xi(1:2), w

   case(7)
      ! RHS for problem with BC is zero Dirichlet and subdomain is the whole Omega
      !such that the solution z = 16x(1-x)y(1-y),
      diff = Re1 * 32.0* ( xi(1)*(1.0-xi(1)) + xi(2)*(1.0-xi(2)) )

      !if(xi(1) > 0) then  ! alternative (singular) test case, iexact = 78 
      !   diff = Re1 * (15*xi(1) -3 )/4/ sqrt(xi(1)) * xi(2) * ( 1- xi(2) ) &
      !        + 2.* sqrt(xi(1)**3)* (1-xi(1) )
      !endif
      
      ! b*grad(u), b = (1,1) so in DP we should have -b * grad(u) =
      ! iconv == 1 or 0
      conv = 0.0 !-16 * ( (1.-2.0*xi(1))*xi(2)*(1.0-xi(2)) + (1.-2.0*xi(2))*xi(1)*(1.0-xi(1)) )
      w = diff + conv

   case(8) ! specially for icase=76, Harriman
      if ( xi(2) > 1-eps .and. xi(1) >= 0.25 .and. xi(1) <= 0.625 ) then ! bottom
        w = 1.
        !print*, 'weightFun on the bottom'
        !w = 0.
        !write(*, *) '^%E#^FE^#Y', xi(1:2), w
      else
         print*, '?Problem in weightFun for Boundary functional, xi=' , xi
         w = 0.0
      endif
     !write(38, *) xi(1:2), w

   case(9) ! specially for icase=63, battery
      if ( xi(1) < eps .and. xi(2) >= -eps .and. xi(2) <= 24+eps ) then ! bottom
        w = 1.
        !print*, 'weightFun on the bottom'
        !w = 0.
        !write(*, *) '^%E#^FE^#Y', xi(1:2), w
      else
         print*, '?Problem in weightFun for Boundary functional, xi=' , xi
         w = 0.0
      endif
     !write(38, *) xi(1:2), w
   case(10) ! specially for icase=67, simplified battery
      if ( xi(1) < eps .and. xi(2) >= -eps .and. xi(2) <= 1+eps ) then ! bottom
        w = 1.
        !print*, 'weightFun on the bottom'
        !w = 0.
        !write(*, *) '^%E#^FE^#Y', xi(1:2), w
      else
         print*, '?Problem in weightFun for Boundary functional, xi=' , xi
         w = 0.0
      endif
     !write(38, *) xi(1:2), w


   case default
      print*, 'Should stop! unknown parameter i in weightFun!'
   end select

end function weightFun

end module weight_fun_mod
