!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DO_LEAF_MEASUREMENTS Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine do_leaf_measurements (osolve,ov,params,istep,iter,iter_nl)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine calculates the second invariant of the strain rate,
! and a value for the criterion used by the improve_octree subroutine,
! inside each element by using the velocity information at the nodes.
! mpe is number of nodes per element (8)
! icon is connectivity array
! xnode,ynode,znode are the global coordinate arrays of length nnode
! unode,vnode,wnode are the velocity arrays of length nnode (obtained
! from previous time step or at least containing the proper velocity
! at the fixed dofs)
! nnode is number of nodes
! nleaves is the number of leaves/elements
! e2d is the second invariant of the strain rate
! crit is the criterion
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
use gauss
use invariants
use constants
implicit none
type (octreev) ov
type (octreesolve) osolve
type (parameters) params
integer istep,iter,iter_nl
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision,dimension(:),allocatable :: x,y,z,vx,vy,vz
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
double precision :: r,s,t,volume,w
double precision :: exx,eyy,ezz,exy,eyz,ezx,J2d,J3d
integer err,ierr,nproc,iproc,k,i,iint,nleaves,mpe
double precision, dimension(:), allocatable :: e2dtemp,e3dtemp,crittemp,lodetemp
double precision, dimension(:), allocatable :: qtemp,dilatrtemp
character*72 :: shift
logical cutcell
double precision :: prod
!=====[2 oct 2007]====
!double precision pw,pw2
!double precision x_P1,y_P1,z_P1
!double precision x_P2,y_P2,z_P2
!double precision x_P3,y_P3,z_P3
!double precision x_P4,y_P4,z_P4
!double precision x_P5,y_P5,z_P5
!double precision x0,y0,z0,dxyz
!integer np,iy
!double precision x_A,x_B,xxx,yyy,zzz,Vme2d,maxe2d
!logical, dimension(:), allocatable :: filter
!integer ny
!integer leaf,level,locc
!==========
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
mpe=params%mpe
nleaves=osolve%nleaves
allocate (x(mpe),y(mpe),z(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc x,y,z in compute_e2d_crit$')
allocate (vx(mpe),vy(mpe),vz(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vx,vy,vz in compute_e2d_crit$')
allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdx,dhdy,dhdz in compute_e2d_crit$')
allocate (e2dtemp(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc e2dtemp in compute_e2d_crit$')
allocate (e3dtemp(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc e2dtemp in compute_e3d_crit$')
allocate (crittemp(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc crittemp in compute_e2d_crit$')
allocate (lodetemp(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lodetemp in compute_e2d_crit$')
allocate (qtemp(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc qtemp in compute_e2d_crit$')
allocate (dilatrtemp(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dilatrtemp in compute_e2d_crit$')
e2dtemp=0.d0
e3dtemp=0.d0
crittemp=0.d0
lodetemp=0.d0
qtemp=0.d0
dilatrtemp=0.d0
osolve%e2d=0.d0
osolve%e3d=0.d0
osolve%crit=0.d0
osolve%lode=0.d0
osolve%q=0.d0
osolve%dilatr=0.d0
do i=1+iproc,nleaves,nproc
do k=1,mpe
x(k) =osolve%x(osolve%icon(k,i))
y(k) =osolve%y(osolve%icon(k,i))
z(k) =osolve%z(osolve%icon(k,i))
vx(k)=ov%unode(osolve%icon(k,i))
vy(k)=ov%vnode(osolve%icon(k,i))
vz(k)=ov%wnode(osolve%icon(k,i))
enddo
r=0.d0
s=0.d0
t=0.d0
call compute_dhdx_dhdy_dhdz (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,mpe
exx=exx+dhdx(k)*vx(k)
eyy=eyy+dhdy(k)*vy(k)
ezz=ezz+dhdz(k)*vz(k)
exy=exy+(dhdx(k)*vy(k)+dhdy(k)*vx(k))/2.d0
eyz=eyz+(dhdy(k)*vz(k)+dhdz(k)*vy(k))/2.d0
ezx=ezx+(dhdz(k)*vx(k)+dhdx(k)*vz(k))/2.d0
enddo
! Dilatation rate
dilatrtemp(i)=(exx+eyy+ezz)
call deviatoric (params%invariants_2d,exx,eyy,ezz,exy,eyz,ezx)
! if (.not.osolve%is_plastic(i)) then
! qtemp(i)=0.d0
! else
qtemp(i)=2.d0*osolve%eviscosity(i)*sqrt(second_invariant(exx,eyy,ezz,exy,eyz,ezx))
! end if
J2d=second_invariant(exx,eyy,ezz,exy,eyz,ezx)
J3d=third_invariant (exx,eyy,ezz,exy,eyz,ezx)
e2dtemp(i)=sqrt(J2d)
if (J2d/=0.d0) lodetemp(i)=lode_angle (J2d,J3d)/pi
if (J3d<0.d0) then
e3dtemp(i)=-abs(J3d)**(1./3.d0)
else
e3dtemp(i)=J3d**(1./3.d0)
end if
select case (params%refine_criterion)
case(1)
crittemp(i)=e2dtemp(i)
case(2)
crittemp(i)=sqrt(exx**2+eyy**2+ezz**2)
case(3)
crittemp(i)=sqrt(J2d)*abs(x(1)-x(2))
case default
crittemp(i)=0.d0
end select
end do
call mpi_reduce (e2dtemp, osolve%e2d, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (e3dtemp, osolve%e3d, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (crittemp,osolve%crit,nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (lodetemp,osolve%lode,nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (qtemp, osolve%q, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
call mpi_reduce (dilatrtemp, osolve%dilatr, nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
deallocate (x,y,z)
deallocate (vx,vy,vz)
deallocate (dhdx,dhdy,dhdz)
deallocate (crittemp)
deallocate (e2dtemp)
deallocate (e3dtemp)
deallocate (lodetemp)
deallocate (qtemp)
deallocate (dilatrtemp)
if (iproc.eq.0 .and. params%debug>=1 .and. params%compute_qpgram) then
write(*,'(a,f11.5)') shift//'min(q) =', minval(osolve%q)
write(*,'(a,f11.5)') shift//'max(q) =', maxval(osolve%q)
end if
call DoRuRe_vel_stats (params%doDoRuRe,istep,iter,iter_nl,ov%nnode,ov%unode,ov%vnode,ov%wnode)
call DoRuRe_leaf_stats (params%doDoRuRe,osolve%nleaves,osolve%e2d,osolve%e3d,osolve%pressure,osolve%q,osolve%dilatr,istep,iter,iter_nl)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|