!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              FIND_VOID_NODES    Nov. 2006                                    |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine find_void_nodes (params,osolve,vo,istep,ov,threadinfo)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine finds all nodes that are cmpletely in the void
! by this I mean that they are not connected by any element to
! a node that is not in the void
! those nodes are given a vo%node=1 flag; all others have vo%node=0
! there are vo%nnode nodes that are not in the void

! it also calculates the nodes that are in the fluid vo%influid=.true.

! it also caculates the number of elements (vo%nleaves) and
! bad faces (vo%nface) that are not in the void
! and corresponding flag arrays vo%leaf and vo%flace

! finally it computes two "correspondence" arrays that tell us the connection
! between dofs in the global set and dofs in the subset that is not in the void

! nleaves is the number of leaves/fe
! nface is the number of bad faces
! nnode is the number of nodes
! icon is the element-node connectivity matrix
! iface is the bad face connectivity matrix
! mpe is the number of nodes per elements
! vo is the structure containing the information on where the void is
! it is constructed in this routine
! lsf contains the nodal values of the nlsf level set functions
! materialn contains the material number associated to each lsf

! the algorithms in this routine are quite straightforward
! and based on the value of the first (zeroth) lsf

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use threads
use definitions

implicit none

type (parameters) params
type (octreesolve) osolve
type (void) vo
integer istep
type (octreev) ov
type (thread) threadinfo

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i,k,kk,ilsf
integer iproc,nproc,ierr
character(len=2) cistep
character*72 shift

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

shift=' '

vo%node(1:osolve%nnode)=0
vo%influid(1:osolve%nnode)=.true.
!vo%leaf(1:nleaves)=0
vo%leaf(1:osolve%nnode)=0

vo%face(1:osolve%nface)=0

vo%nleaves = osolve%nleaves
vo%nnode   = osolve%nnode
vo%nface   = osolve%nface
        
do i=1,osolve%nnode
   vo%rtf(i)=i
   vo%ftr(i)=i
enddo

if (params%materialn(0).gt.0) then
   ov%whole_leaf_in_fluid=.true.
   return
end if

if (.not. params%excl_vol) then

   do i=1,osolve%nnode
      if (osolve%lsf(i,1).gt.0.d0) vo%influid(i)=.false.
   enddo

else

   vo%influid=.false.
   do i=1,osolve%nnode
      do ilsf=1,osolve%nlsf
         if(osolve%lsf(i,ilsf) .le. 0.d0) vo%influid(i)=.true.
      end do
   enddo

end if

write(threadinfo%Logunit,*) '[same] find_void_nodes:  count(vo%influid)',count(vo%influid) 

!if (iproc.eq.0) then
!   write(cistep,'(i2)') istep
!   if (istep<10) cistep(1:1)='0'
!   open(unit=9999,file='influid_'//cistep//'.vtk')
!   write(9999,'(a)')'# vtk DataFile Version 3.0'
!   write(9999,'(a)')'influid'
!   write(9999,'(a)')'ASCII'
!   write(9999,'(a)')'DATASET UNSTRUCTURED_GRID'
!   write(9999,'(a7,i10,a6)')'POINTS ',count(vo%influid),' float'
!   do i=1,nnode
!      if (vo%influid(i)) write(9999,'(3f8.4)') x(i),y(i),z(i)
!   end do
!   close(9999)
!end if

ov%whole_leaf_in_fluid=.false.
do i=1,osolve%nleaves
   do ilsf=1,osolve%nlsf
      if (osolve%lsf(osolve%icon(1,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(2,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(3,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(4,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(5,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(6,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(7,i),ilsf)<0.d0 .and. &
          osolve%lsf(osolve%icon(8,i),ilsf)<0.d0 ) then
         ov%whole_leaf_in_fluid(i)=.true.
      end if
   end do
end do

if(iproc.eq.0 .and. params%debug>=1) write(*,'(a,i8)') shift//'whole leaf in fluid ',count(ov%whole_leaf_in_fluid)

write(threadinfo%Logunit,*) '[same] find_void_nodes: count(ov%whole_leaf_in_fluid)',count(ov%whole_leaf_in_fluid) 

!=====[]=====

vo%node(1:osolve%nnode)=1

do ilsf=1,osolve%nlsf
do i=1,osolve%nleaves
   do k=1,params%mpe
      if (osolve%lsf(osolve%icon(k,i),ilsf).le.0.d0 ) then
         do kk=1,params%mpe
            vo%node(osolve%icon(kk,i))=0
         enddo
         goto 111
      endif
   enddo
111   continue
enddo
enddo


vo%nnode=0

do i=1,osolve%nnode
   vo%ftr(i)=0
   if (vo%node(i).eq.0) then
      vo%nnode=vo%nnode+1
      vo%rtf(vo%nnode)=i
      vo%ftr(i)=vo%nnode
   endif
enddo

vo%nleaves=0
vo%leaf(1:osolve%nleaves)=1

do i=1,osolve%nleaves
   if (sum(vo%node(osolve%icon(:,i))).eq.0) then
      vo%nleaves=vo%nleaves+1
      vo%leaf(i)=0
   endif
enddo

vo%nface=0
vo%face(1:osolve%nface)=1

do i=1,osolve%nface
   if (sum(vo%node(osolve%iface(:,i))).eq.0) then
      vo%nface=vo%nface+1
      vo%face(i)=0
   endif
enddo

if (iproc.eq.0) then
   write (8,*) 'total nnode , nodes out of void ',osolve%nnode,vo%nnode
   write (8,*) 'total nleaves , leaves out of void ',osolve%nleaves,vo%nleaves
   write (8,*) 'total nface, faces out of void ',osolve%nface,vo%nface
end if

write(threadinfo%Logunit,*) '[same] find_void_nodes: osolve%nnode  ',osolve%nnode
write(threadinfo%Logunit,*) '[same] find_void_nodes: osolve%nleaves',osolve%nleaves
write(threadinfo%Logunit,*) '[same] find_void_nodes: osolve%nface',osolve%nface
write(threadinfo%Logunit,*) '[same] find_void_nodes: vo%nnode      ',vo%nnode
write(threadinfo%Logunit,*) '[same] find_void_nodes: vo%nleaves    ',vo%nleaves
write(threadinfo%Logunit,*) '[same] find_void_nodes: vo%nface      ',vo%nface

call flush(threadinfo%Logunit) 

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
