!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              COMPUTE_PRESSURE    Oct. 2008                                   |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine compute_pressure (params,osolve,ov,vo,mat,cl)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine is used to calculate the pressure
! It is an equivalent to build_system but in a stripped down version

! 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
! mpe is the number of nodes per elements
! x,y,z are the coordinates of the nodes of the fe grid
! lsf contains the nodal values of the nlsf level set functions
! mat is the material matrix for the nmat materials
! materialn contains the material number associated to each lsf
! u,v,w are the nodal velocities from the previous time step/iteration
! temp and pressure are the nodal temperature and pressure interpolated from the
!   previous time step/iteration
! vo is the structure containing the information on where the void is
! levelcut is the maximum number of levels used for dividing cut cells
! levelapprox is used to improve the estimate of the positive volume
!   in cut cells by another subdivision

! Like in build_system, we go through a routine called pressure_cut
! that deals with the integration of cut cells with precision
! It is necessary to go through this lengthy process of checking for cut cells
! in case the penalty parameter is not uniform in the fe grid

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

use definitions
!use mpi

implicit none

include 'mpif.h'

type (parameters) params
type (octreesolve) osolve
type (octreev) ov
type (void) vo
type (material) mat(0:params%nmat)
type (cloud) cl

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

integer,dimension(:,:),allocatable :: cloud_elem_mat_bins
integer,dimension(:),allocatable :: leaf_mat_bin
integer err,iproc,nproc,ierr,ileaves,level,i,k,levelmax,icut
double precision,dimension(:,:),allocatable::lsf_el
double precision r0,s0,t0,rst
double precision,dimension(:),allocatable::press
!integer,dimension(:),allocatable::levs
character*72  :: shift

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

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

shift=' '

allocate (press(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc press in compute_pressure$')

press=0.d0

allocate (lsf_el(params%mpe,osolve%nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsf_el in compute_pressure$')
! Allocate array to store cloud particle material types for each element
allocate (cloud_elem_mat_bins(osolve%nleaves,0:params%nmat),stat=err)
if (err.ne.0) call stop_run ('Error alloc cloud_elem_mat_bins in build_system$')



!write(*,*) 'After alloc cloud_elem_mat_bins'



! Calculate number of cloud particles of each material in all elements
call find_mat_numbers_from_cloud(params,cl,osolve,cloud_elem_mat_bins)


!write(*,*) 'After find_mat_numbers'


! Allocate small array to store cloud particle material types for individual
! elements
allocate (leaf_mat_bin(0:params%nmat),stat=err)
if (err.ne.0) call stop_run ('Error alloc leaf_mat_bin in build_system$')


!write(*,*) 'After alloc leaf_mat_bin'


do ileaves=1+iproc,osolve%nleaves,nproc
   ! first test if in the void
   if (vo%leaf(ileaves).eq.0) then
      level=0
      levelmax=params%levelcut
      do k=1,params%mpe
        do i=1,osolve%nlsf
          lsf_el(k,i)=osolve%lsf(osolve%icon(k,ileaves),i)
        enddo
      enddo
      ! Grab material bins for current element
      leaf_mat_bin=cloud_elem_mat_bins(ileaves,:)
      
      
      !write(*,*) 'After fill leaf_mat_bin'
      !write(*,*) 'osolve%matnum for element ',ileaves,': ',osolve%matnum(ileaves)
      
      r0=-1.d0
      s0=-1.d0
      t0=-1.d0
      rst=2.d0
      call pressure_cut (params,level,levelmax,osolve%icon(1,ileaves),         &
                         osolve%x,osolve%y,osolve%z,mat,params%materialn,      &
                         ov%unode,ov%vnode,ov%wnode,osolve%temp,press(ileaves),&
                         osolve%nnode,lsf_el,osolve%nlsf,r0,s0,t0,rst,icut,    &
                         ileaves,osolve%eviscosity(ileaves),                   &
                         osolve%matnum(ileaves),leaf_mat_bin,                  &
                         osolve%matnump(ileaves))


!write(*,*) 'After pressure_cut'



   endif
enddo

deallocate (lsf_el)
! Deallocate material bin arrays
deallocate(cloud_elem_mat_bins,leaf_mat_bin)

osolve%pressure=0.d0
call mpi_allreduce (press,osolve%pressure,osolve%nleaves,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

deallocate (press)

if (iproc.eq.0 .and. params%debug>=1) then
   write(*,'(a,2e13.5)') shift//'raw pressures :',minval(osolve%pressure),maxval(osolve%pressure)
end if


!12/02/2008
! modification by Jean (no need to rescale pressures by element size
! as it is not done at element level (in make_pressure) any longer
!allocate (levs(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc levs in compute_pressure$')
!call octree_find_element_level (octree,noctree,levs,nleaves)
!pressure=pressure*(8.d0**levs)
!deallocate (levs)

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