Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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, &
!use mpi
implicit none
include 'mpif.h'
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
!------------------------------------------------------------------------------|
!(((((((((((((((( 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))
Dave Whipp
committed
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
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|