!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              MAKE_PRESSURE    Apr. 2007                                      |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine make_pressure (params,icon,xg,yg,zg,viscosity0,penalty,expon,unode, &
                         vnode,wnode,temp,pressure,strain,nnode,r0,s0,t0,rst,  &
                         eviscosity)

use definitions
use invariants
!use mpi

implicit none

include 'mpif.h'

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine is a stripped down version of make_matrix but is used in the
! calculation of the pressure

! icon is connectivity array for the current element
! xg,yg,zg are the global coordinate arrays of length nnode
! viscosity0 is viscosity
! penalty is incompressibility penalty factor
! expon is nonlinear viscosity exponent
! unode,vnode,wnode is the velocity array (obtained from previous time step
! or at least containing the proper velocity at the fixed dofs)
! temp and pressure anre temperature and pressure
! strain is total accumulated strain as interpolated from the 3D cloud
! nnode is number of nodes
! r0,s0,t0 are the bottom left back corner of the part of the element
!   we are now integrating (in local r,s,t coordinates)
! rst is the size of the part of the element we are integrating

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

type(parameters) params
integer icon(params%mpe)
double precision xg(nnode),yg(nnode),zg(nnode)
double precision viscosity0
double precision penalty
double precision expon
double precision unode(nnode),vnode(nnode),wnode(nnode)
double precision temp(nnode)
double precision pressure
double precision strain(nnode)
integer nnode
double precision r0,s0,t0,rst
double precision eviscosity

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
!integer nint
integer iproc,nproc,ierr,ic,k,iint,err
double precision viscosity,straintot
double precision exx,eyy,ezz,exy,eyz,ezx,E2d,E3d
double precision f1(3,3),defgrad(3,3)
double precision detf1,detf,e2dref,viscomean
double precision velox,veloy,veloz,tau,uvwnorm,xmin,xmax,ymin,ymax
double precision zmin,zmax,dx,dy,dz,alpha,div
double precision eps,r,s,t,w,volume
double precision jcb(3,3),jcbi(3,3),jcbp(3,3),jcbip(3,3)
double precision,dimension(:),allocatable :: x,y,z
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
logical is_plastic_dummy,flag_vrm_pb

double precision :: rr,ss,tt,ww

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

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

eps=tiny(eps)

allocate (x(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc x in make_pressure$')
allocate (y(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc y in make_pressure$')
allocate (z(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc z in make_pressure$')
allocate (dhdx(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdx in make_pressure$')
allocate (dhdy(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdy in make_pressure$')
allocate (dhdz(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdz in make_pressure$')

rr=0.d0
ss=0.d0
tt=0.d0
ww=8.d0

do k=1,params%mpe
   x(k)=xg(icon(k))
   y(k)=yg(icon(k))
   z(k)=zg(icon(k))*params%vex
enddo

! integrate over nint integration points
! here nint has been chosen to 1 because incompresibility
! is imposed in the FE matrix with a single integration point

   r=r0+rst*(rr+1.d0)/2.d0
   s=s0+rst*(ss+1.d0)/2.d0
   t=t0+rst*(tt+1.d0)/2.d0
   w=ww

   ! added by Jean 13/2/08
   r=rr
   s=ss
   t=tt
   w=ww
   !end addition

   call compute_dhdx_dhdy_dhdz (params%mpe,r,s,t,x,y,z,dhdx,dhdy,dhdz,volume)

   exx=0.d0
   eyy=0.d0
   ezz=0.d0
   exy=0.d0
   eyz=0.d0
   ezx=0.d0
   do k=1,params%mpe
      ic=icon(k)
      exx=exx+dhdx(k)*unode(ic)
      eyy=eyy+dhdy(k)*vnode(ic)
      ezz=ezz+dhdz(k)*wnode(ic)
   enddo

! Change by Jean 13/2/08
! divergence should not be scaled by volume or weight
!   div=(exx+eyy+ezz)*volume*w
   div=exx+eyy+ezz
! end of change

if (params%bulkvisc) then
  pressure=div*penalty
else
  pressure=div*penalty*eviscosity
endif

deallocate (x)
deallocate (y)
deallocate (z)
deallocate (dhdx)
deallocate (dhdy)
deallocate (dhdz)

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