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, &
plasticity_type,plasticity_parameters,eviscosity)
Dave Whipp
committed
use mpi
25
26
27
28
29
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
111
!------------------------------------------------------------------------------|
!(((((((((((((((( 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 ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
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
character (len=8) plasticity_type
double precision plasticity_parameters(9)
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 :: rr,ss,tt,ww,h
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|