!> emiss boundary condition for PM
module emiss_bc
  use main_data
  !  use problem_oper
  use eval_sol
  use set_solution
  use stdgm_mod
  use mesh_mod


  implicit none
  
  public:: SetEmissionBC
  public:: Elem_Flux_edge

contains
  
  !>  seeking elements for emission BC
  subroutine SetEmissionBC( FinalStep, file_name )
    logical, intent(in) :: FinalStep ! if true, solution from t_m^-
    character(len=50), intent(in) :: file_name
    class(element), pointer :: elem
    integer :: i, k, l, ie, Qdof
    real, dimension(:,:), allocatable :: wi   ! w recomputed  in integ nodes
    real, dimension(:,:), allocatable :: PMflows         ! porous media
    real, dimension(:,:), allocatable :: Re_1
    real :: flux, head, PM_min, sig_h, sig_fl, tol_h, tol_fl, aver, val
    integer :: ifile, is, inum, idat

    
    if(FinalStep) then
       ifile = 11
       open(ifile, file = file_name, status='UNKNOWN')  !, position='APPEND')
    endif
    
    
    PM_min = 1E-5

    do i=1,grid%nelem
       elem => grid%elem(i)

       do ie = 1, elem%flen
          k = elem%face(neigh,ie)

          !  state%BC(elem%iBC(ie))%inout   -- original value of type BC
          !        elem%tBC(ie)             -- new setting

          if( k <= 0  .and. elem%iBC(ie)> 0 ) then

             if( state%BC(elem%iBC(ie))%inout == 3 ) then ! emission boundary condition

                !! seting of degree of the Gauss quadrature
                Qdof = elem%face(fGdof,ie)

                ! solution at the last time level
                !if(FinalStep) & 
                !     elem%w(0,1:elem%dof) =  elem%wSTfin(1,1:elem%dof)

                ! evaluation of w_ie in integ nodes
                allocate(wi(1:Qdof,1:ndim) )
                call Eval_w_Edge(elem, ie, wi, .false.)

                allocate(PMflows(1:3, 1:Qdof), Re_1(1:iRe, 1:Qdof) )
                Re_1(2:iRe, 1:Qdof) = transpose( elem%xi(ie, 1:Qdof,  2+1:2+iRe-1) )
                
                Re_1(1, 1:Qdof) = 0.
                do l=1,3
                   Re_1(1, 1:Qdof) = Re_1(1, 1:Qdof) +  Re_1(1+l, 1:Qdof) * soilpar(l)%alpha
                enddo
                
                ! preassure head
                PMflows(1, 1:Qdof) =  wi(1:Qdof,1 ) -  grid%b_edge(-elem%face(neigh,ie))%x_div(1:Qdof,2)
                call IntegrateFunctionEdge_MeanValue(elem, ie,PMflows(1, 1:Qdof), head)
                !print*,'###S#D#:',Qdof, head, PMflows(1, 1:4)
                
                ! flux through the boundary element edge
                call Elem_Flux_edge(elem, ie, flux, 0, PMflows(2:3, 1:Qdof) )

             !!!deallocate(wi, PMflows, Re_1)
             !!!if( state%BC(elem%iBC(ie))%inout == 3 ) then ! emission boundary condition


                !write(*,'(a8, 3i5, 30es12.4)') &
                !     'Emiss:',state%time%iter, state%nlSolver%iter, elem%i, &
                !     elem%xi(ie, 1:elem%face(fGdof, ie), 2+iRe)


                !tol_h =  1E-1
                !tol_h =  2E-1
                !tol_h =  5E-1
                !tol_h = elem%dn(ie) /5.   ! WORKS
                tol_h = elem%dn(ie) /10.
                
                !tol_fl =  2E-1
                !tol_fl = 1E-3
                tol_fl = elem%dn(ie) * elem%dn(ie)

                
                ! setting of the local "sigma"parameter for smooth changing of the
                ! emission BC
                do l =1,Qdof
                   !tol_h = 0.5 / Re_1(1,l)  *1E-2
                   !tol_fl = PMflows(3, l) * 1E-2

                   ! setting of local sigma for pressure head
                   !!sig_h = smooth_Heaviside( head, -tol_h, 0., 0., 1.)
                   !sig_h = smooth_Heaviside( PMflows(1, l), -tol_h, 0., 0., 1.)
                   sig_h = smooth_Heaviside( PMflows(1, l), 0., tol_h,  0., 1.)
                   
                   ! setting of local sigma for pressure flux
                   !!sig_fl = smooth_Heaviside( flux, 0., tol_fl,  1., 0.)
                   !sig_fl = smooth_Heaviside( PMflows(2, l), 0., tol_fl,  1., 0.)
                   !!!sig_fl = smooth_Heaviside( PMflows(2, l), -tol_fl,  1.,  0., 1.)
                   sig_fl = smooth_Heaviside( PMflows(2, l), 0., tol_fl,   1., 0.)


                   ! version (1)
                   elem%xi(ie, l, 2+iRe) = sig_fl *  sig_h
                   !elem%xi(ie, l, 2+iRe) = sig_h  ! WORKS

                   !if(sig_fl < 0.9) &
                   !     write(*,'(3i5, 20es12.4)') &
                   !     state%time%iter, state%Nlsolver%iter, elem%i, &
                   !     PMflows(1:2, l) , sig_h, sig_fl, elem%xi(ie, l, 1:2), elem%xi(ie, l, 2+iRe)
                   
                   ! version (2),
                   !elem%xi(ie, l, 2+iRe) = max(0., PMflows(1, l)) * 20. / elem%diam
                   
                   !idat = 100*(state%time%iter+1) +  state%Nlsolver%iter
                   !if(sig_fl * sig_h > 0.) &
                   !write(idat,*) state%time%iter, state%Nlsolver%iter, elem%i, &
                   !     PMflows(1:2, l) , sig_h, sig_fl, elem%xi(ie, l, 1:2), elem%xi(ie, l, 2+iRe)
                   if(FinalStep) then
                      write(ifile,'(30es12.4)') state%time%ttime,  &
                           grid%b_edge(-elem%face(neigh,ie))%x_div(l,1:2) , &  !2:3 x_i 
                           PMflows(1:2, l),  PMflows(1, l) * PMflows(2, l),& !4:6  h, flux
                           sig_h, sig_fl, sig_h*sig_fl, elem%xi(ie, l, 2+iRe), &  ! 7:10
                           -PMflows(2, l)/elem%dn(ie), elem%dn(ie)
                      if(l == Qdof) write(ifile,'(x)') 
                   endif

                enddo

                ! averaging
                !val = 0.5
                !aver = sum(elem%xi(ie, 1:Qdof, 2+iRe) ) /Qdof
                !elem%xi(ie, 1:Qdof, 2+iRe) = val*aver + (1.-val)*elem%xi(ie, 1:Qdof, 2+iRe)

                
!                if(FinalStep) then
!                   do l=1,Qdof
!                      write(ifile,'(30es12.4)') state%time%ttime,  &
!                           grid%b_edge(-elem%face(neigh,ie))%x_div(l,1:2) , &  !2:3 x_i 
!                           PMflows(1:2, l),  PMflows(1, l) * PMflows(2, l),& !4:6  h, flux
!                           sig_h, sig_fl, sig_h*sig_fl, elem%xi(ie, l, 2+iRe)  ! 7:10
!                   enddo
!                   write(ifile,'(x)') 
!                endif

                

                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                !elem%xi(ie, 1:Qdof, 2+iRe) = 1.
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                
                !if(minval(elem%xi(ie, 1:Qdof, 2+iRe))  == 0.) then
                if(maxval(elem%xi(ie, 1:Qdof, 2+iRe))  == 0.) then
                   elem%tBC(ie) = -3  ! acts as the Neumann BC
                else
                   elem%tBC(ie) = 3  ! acts as the Dirichlet BC (with possibly no influence)
                   ! write(*,'(a10, 4i6, 300es12.4)') &
                   !      'emisc deg5:',state%time%iter, state%Nlsolver%iter, elem%i, &
                   !      elem%tBC(ie), minval(elem%xi(ie, 1:Qdof, 2+iRe)) , &
                   !      maxval(elem%xi(ie, 1:Qdof, 2+iRe)),elem%xi(ie, 1:Qdof, 2+iRe) 
                   ! write(*,'(a10, 4i6, 300es12.4)') &
                   !      'emisc deg5:',state%time%iter, state%Nlsolver%iter, elem%i, &
                   !      elem%tBC(ie), minval(PMflows(1, 1:Qdof)) , &
                   !      maxval(PMflows(1, 1:Qdof)),PMflows(1, 1:Qdof)
                   ! print*
                endif
                
                ! output for verification of emmision BC
                ! if(FinalStep) then
                !    !if( state%BC(elem%iBC(ie))%inout == 3 ) then
                !       do l=1,Qdof
                !          write(ifile,'(30es12.4)') state%time%ttime,  &
                !               grid%b_edge(-elem%face(neigh,ie))%x_div(l,1:2) , &  !2:3 x_i 
                !               PMflows(1:2, l),  PMflows(1, l) * PMflows(2, l),& !4:6  h, flux
                !               sig_h, sig_fl, sig_h*sig_fl  ! 7:9
                !       enddo
                !       write(ifile,'(x)') 
                !    !end if
                ! endif


                
                !write(*,'(a8, 4i5, 30es12.4)') 'emiss:',state%time%iter, state%nlSolver%iter, &
                !     elem%i,  elem%tBC(ie) , head, flux,  elem%xi(ie, 1:Qdof, 2+iRe)

                ! ! new implementation of the OLD variant
                ! if(head  > PM_min .and. flux < 0.) then
                !    ! Dirichlet BC, we go to call ElementViscousBoundEdge(elem, ie , ...)

                !    elem%tBC(ie) = 3  !1  !
                !    elem%xi(ie, :, 2+iRe) = 1.
                   
                   
                !    write(*,'(a8, 3i5, 30es12.4)') 'emiss:',state%time%iter, state%nlSolver%iter, &
                !         elem%i, head, flux

                !    open(199, file ='emiss_BC', status='unknown', position='append')
                !    write(199, *) state%time%iter,  elem%i, &
                !         grid%b_edge(-elem%face(neigh,ie))%x_div(Qdof/2, 1:nbDim),&
                !         head, flux
                !    close(199)
                   
                !    !do l=1,Qdof
                !    !   write(200+ state%time%iter, *) &
                !    !        grid%b_edge(-elem%face(neigh,ie))%x_div(l, 1:nbDim),PMflows(1:2, l)
                !    !enddo

                ! else ! homogeneous Neuman BC

                !    elem%xi(ie, :, 2+iRe) = 0.
                !    elem%tBC(ie) = -3 !!!    0  !3 

                !    !do l=1,Qdof
                !    !   write(100+ state%time%iter, *) &
                !    !        grid%b_edge(-elem%face(neigh,ie))%x_div(l, 1:nbDim),PMflows(1:2, l)
                !    !enddo



                ! endif

                
                ! ! OLD variant
                ! if(head  > PM_min .and. flux < 0.) then
                !    ! Dirichlet BC, we go to call ElementViscousBoundEdge(elem, ie , ...)

                !    elem%tBC(ie) = 1  !

                !    write(*,'(a8, 3i5, 30es12.4)') 'emiss:',state%time%iter, state%nlSolver%iter, &
                !         elem%i, head, flux

                !    open(199, file ='emiss_BC', status='unknown', position='append')
                !    write(199, *) state%time%iter,  elem%i, &
                !         grid%b_edge(-elem%face(neigh,ie))%x_div(Qdof/2, 1:nbDim),&
                !         head, flux
                !    close(199)
                   
                !    !do l=1,Qdof
                !    !   write(200+ state%time%iter, *) &
                !    !        grid%b_edge(-elem%face(neigh,ie))%x_div(l, 1:nbDim),PMflows(1:2, l)
                !    !enddo

                ! else ! homogeneous Neuman BC


                !    elem%tBC(ie) = 3

                !    !do l=1,Qdof
                !    !   write(100+ state%time%iter, *) &
                !    !        grid%b_edge(-elem%face(neigh,ie))%x_div(l, 1:nbDim),PMflows(1:2, l)
                !    !enddo



                ! endif

                deallocate(wi, PMflows, Re_1)

             endif  ! BC == 3

          endif  ! if (k <= 0)

       enddo  ! ie =1,elem%flen

    end do  ! i =1,grid%nelem

    if(FinalStep) close(ifile)

  end subroutine SetEmissionBC


  !> compute the flux throught the edge of an element
  subroutine Elem_Flux_edge(elem, ie, val, alpha, flux_nodes)
    class(element), intent(inout) :: elem
    integer, intent(in) :: ie  ! index of the edge
    integer, intent(in) :: alpha ! time integ nodes
    real, intent(inout) :: val  ! integral value
    real, dimension(1:2, 1: elem%face(fGdof,ie)), optional :: flux_nodes ! flux in ineg nodes
    real, dimension(:,:,:), allocatable :: Dwi
    real, dimension(:,:,:), allocatable :: K_sk
    real, dimension(:), allocatable :: Re_1, weights
    integer :: Qdof, Qnum, l

    Qdof = elem%face(fGdof,ie)
    Qnum = elem%face(fGnum,ie)

    allocate(Dwi(1:Qdof, 1:ndim, 0:nbDim), K_sk(1:Qdof, 1:nbDim, 1:nbDim) )
    allocate(Re_1(1:iRe), weights(1:Qdof))

    ! w in integ nodes
    call Eval_w_Edge(elem, ie, Dwi(1:Qdof, 1:ndim, 0), .false.) !!elem%n(ie,:) !

    ! Dw ininteg nodes
    call Eval_Dw_Edge(elem, ie, Dwi(1:Qdof, 1:ndim, 1:nbDim), .false.) !!elem%n(ie,:) !


    weights(1:Qdof) = state%space%G_rule(Qnum)%weights(1:Qdof) !* elem%dn(ie)

    val = 0.
    ! evaluation of the water content in integ nodes
    do l=1,Qdof
       Re_1(1) = 1.
       Re_1(2:iRe) = elem%xi(ie, l,  2+1:2+iRe-1)

       call Eval_Diff_Porous_Coeffs(Dwi(l,1,0), Dwi(l,1,1:2), K_sk(l, 1:nbDim, 1:nbDim), &
            Re_1(1:iRe), 0, elem%xi(ie, l, 1:2) )

       !val = val + K_sk(l,1,1) * weights(l) *  dot_product( Dwi(l,1,1:2),  elem%n(ie,1:2) )
       val = val +  weights(l) &
            * dot_product( matmul( K_sk(l,1:2,1:2),  Dwi(l,1,1:2)),  elem%n(ie,1:2) )

       if( present(flux_nodes) ) then
          !flux_nodes(1, l) =  K_sk(l,1,1) * dot_product( Dwi(l,1,1:2),  elem%n(ie,1:2) )
          flux_nodes(1, l) =  dot_product( matmul( K_sk(l,1:2,1:2), Dwi(l,1,1:2)), elem%n(ie,1:2) )
          flux_nodes(2, l) =  max( K_sk(l,1,1) , K_sk(l,2,2) )
       endif

       !if(elem%i  == 391 ) then
       !!if( l == 1 .and. alpha == 1) print*,'__________________________'
       !!if(l == 1) print*
       !   write(*,'(a8, es12.4, 2i5 , 30es12.4)') &
       !        'flux =' ,val, ie, l,  Dwi(l,1,1:2),  elem%n(ie,1:2), &
       !        dot_product( matmul( K_sk(l,1:2,1:2),  Dwi(l,1,1:2)),  elem%n(ie,1:2) ), &
       !        elem%xi(ie, l, 1:2)
       !
       !endif
          

    enddo

    !if(elem%i == 88) &
    !if(state%time%iter >= 104) &
    !     write(*,'(a8, 3i5, 30es12.4)') 'flux:', &
    !     state%time%iter, state%nlSolver%iter,elem%i, &
    !     val, sum(K_sk(1:Qdof,1,1) )/Qdof, val /( sum(K_sk(1:Qdof,1,1) )/Qdof ), &
    !     sum(Dwi(1:Qdof, 1, 1)) / Qdof, sum(Dwi(1:Qdof, 1, 2)) / Qdof

    !if(elem%i  <= 4) then
    !   write(*,'(a8, es12.4, i5, 30es12.4)')'val =' ,val, elem%i, 1.* elem%i
    !   print*
    !endif

    !if(val < 0.) then
    !   write(19,'(a8, es12.4, i5, 30es12.4)')' FLUX:', val, elem%i, elem%xc(:), elem%n(ie,1:2),&
    !   elem%xi(ie, 1, 1:2)
    !endif


    !
    !call  IntegrateFunctionNormalEdge(elem, ie, transpose(Dwi(1:Qdof, 1, 1:nbdim)), K_sk(1:Qdof, 1, 1),&
    !     val)
    !print*,'val =' ,val
    !print*

    deallocate(Dwi, K_sk, Re_1, weights)

  end subroutine Elem_Flux_edge

  !> a smoothing of the Heaviside function within the interval [xL, xR]
  function smooth_Heaviside(x, xL, xR, fL, fR) result(f)
    real :: f
    real, intent(in) :: x, xL, xR, fL, fR
    real :: r, d

    if(x <= xL) then
       f = fL
    elseif (x >= xR) then
       f = fR
    else
       r = (x-xL)/(xR-xL)
       d = fR - fL
       f = fL + r*r*d*(3 - 2*r)
    endif

  end function smooth_Heaviside

end module emiss_bc
