Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! PRESSURE_CUT Nov. 2006 I |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
Dave Whipp
committed
recursive subroutine pressure_cut (params,level,levelmax,icon,x,y,z,mat, &
materialn,u,v,w,temp,pressure,strain,nnode, &
lsf,nlsf,r0,s0,t0,rst,icut,ileaves, &
eviscosity)
23
24
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
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this is an intermediary routine between find_pressure and make_pressure
! in the same way as make_cut is between build_system and make_matrix
! level is current level in cut cell algortihm
! it varies between 0 and levelmax
! icon is connectivity array for the current element
! x,y,z are the global coordinate arrays of length nnode
! mat is the material matrix for the nmat materials
! materialn contains the material number associated to each lsf
! u,v,w is the velocity array (obtained from previous time step
! or at least containing the proper velocity at the fixed dofs)
! temp, pressure and strain are temperature, pressure and strain
! nnode is number of nodes
! lsf is global array of level set functions defining the surfaces
! nlsf is number of lsfs
! 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
! icut is returned (0 if homogeneous element - 1 if cut element)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
implicit none
type (parameters) params
integer level
integer levelmax
integer icon(params%mpe)
double precision x(nnode),y(nnode),z(nnode)
type (material) mat(0:params%nmat)
integer materialn(0:nlsf)
double precision u(nnode),v(nnode),w(nnode)
double precision temp(nnode)
double precision pressure
double precision strain(nnode)
integer nnode
double precision lsf(params%mpe,nlsf)
integer nlsf
double precision r0,s0,t0,rst
integer icut
integer ileaves
double precision eviscosity
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer err,i,j,k,ii,jj,kk,levelp,matel,jcut
double precision r(params%mpe),s(params%mpe),t(params%mpe)
double precision h(params%mpe),vol_lsf0
Dave Whipp
committed
double precision r0p,s0p,t0p,rstp,dt,eps
double precision viscosity,density,penal,expon,diffusivity,heat
double precision pressurep,prod,volmax
double precision,dimension(:,:),allocatable::lsfp
double precision,dimension(:),allocatable::vol_lsf
character (len=8) plasticity_type
double precision plasticity_parameters(9)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
Dave Whipp
committed
eps=1.d-10
if (level.eq.0) then
pressure=0.d0
endif
matel=materialn(0)
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
if (prod*lsf(k,i).le.0.d0) goto 222
enddo
if (prod.lt.0.d0) then
matel=materialn(i)
endif
enddo
Dave Whipp
committed
call make_pressure (params,icon,x,y,z,mat(matel)%viscosity,mat(matel)%penalty, &
mat(matel)%expon,u,v,w,temp,pressurep,strain,nnode,r0,s0,t0,&
rst,mat(matel)%plasticity_type,&
mat(matel)%plasticity_parameters,eviscosity)
pressure=pressure+pressurep/(8.d0**level)
icut=0
return
222 continue
icut=1
! when we get to the bottom of the division
! we use an approximate algorithm
if (level.eq.levelmax) then
Dave Whipp
committed
allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in pressure_cut$')
do i=1,nlsf
call compute_positive_volume (lsf(1,i),vol_lsf(i),params%levelapprox)
enddo
vol_lsf=1.d0-vol_lsf
if (.not.params%excl_vol) then
do i=1,nlsf-1
vol_lsf(i)=vol_lsf(i)-vol_lsf(i+1)
enddo
end if
do i=1,nlsf
vol_lsf(i)=max(vol_lsf(i),0.d0)
vol_lsf(i)=min(vol_lsf(i),1.d0)
enddo
Dave Whipp
committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
vol_lsf0=1.d0-sum(vol_lsf)
if (vol_lsf0.gt.eps) then ! Element contains void, always use divFEM
viscosity=0.d0
penal=0.d0
expon=0.d0
do i=1,nlsf
matel=materialn(i)
viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
penal=penal+vol_lsf(i)*mat(matel)%penalty
expon=expon+vol_lsf(i)*mat(matel)%expon
enddo
matel=materialn(0)
viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
penal=penal+vol_lsf0*mat(matel)%penalty
expon=expon+vol_lsf0*mat(matel)%expon
elseif (params%matrule.eq.1) then ! Assign material properties of volumetric majority material
matel=params%materialn(sum(maxloc(vol_lsf)))
viscosity=mat(matel)%viscosity
penal=mat(matel)%penalty
expon=mat(matel)%expon
elseif (params%matrule.eq.2) then ! Assign material properties of volumetric minority material
matel=params%materialn(sum(minloc(vol_lsf)))
viscosity=mat(matel)%viscosity
penal=mat(matel)%penalty
expon=mat(matel)%expon
else
if (params%matrule.ne.0) print *,'Invalid matrule value, using divFEM'
viscosity=0.d0
penal=0.d0
expon=0.d0
do i=1,nlsf
matel=materialn(i)
viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
penal=penal+vol_lsf(i)*mat(matel)%penalty
expon=expon+vol_lsf(i)*mat(matel)%expon
enddo
matel=materialn(0)
viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
penal=penal+vol_lsf0*mat(matel)%penalty
expon=expon+vol_lsf0*mat(matel)%expon
endif
if (maxval(vol_lsf).gt.0.d0) then
if (params%matrule.eq.2) then
matel=params%materialn(sum(minloc(vol_lsf)))
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
else
matel=params%materialn(sum(maxloc(vol_lsf)))
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
Dave Whipp
committed
else
matel=params%materialn(0)
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
endif
call make_pressure (params,icon,x,y,z,viscosity,penal,expon,u,v,w,temp, &
pressurep,strain,nnode,r0,s0,t0,rst,plasticity_type, &
plasticity_parameters,eviscosity)
pressure=pressure+pressurep/(8.d0**level)
deallocate (vol_lsf)
return
endif
! If we are not at the bottom level, we keep dividing
Dave Whipp
committed
allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in pressure_cut$')
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
do kk=1,2
t(1:4)=-1.d0+float(kk-1)
t(5:8)=float(kk-1)
do jj=1,2
s(1:2)=-1.d0+float(jj-1)
s(3:4)=float(jj-1)
s(5:6)=-1.d0+float(jj-1)
s(7:8)=float(jj-1)
do ii=1,2
r(1)=-1.d0+float(ii-1)
r(2)=float(ii-1)
r(3)=-1.d0+float(ii-1)
r(4)=float(ii-1)
r(5)=-1.d0+float(ii-1)
r(6)=float(ii-1)
r(7)=-1.d0+float(ii-1)
r(8)=float(ii-1)
do k=1,8
h(1)=(1.d0-r(k))*(1.d0-s(k))*(1.d0-t(k))/8.d0
h(2)=(1.d0+r(k))*(1.d0-s(k))*(1.d0-t(k))/8.d0
h(3)=(1.d0-r(k))*(1.d0+s(k))*(1.d0-t(k))/8.d0
h(4)=(1.d0+r(k))*(1.d0+s(k))*(1.d0-t(k))/8.d0
h(5)=(1.d0-r(k))*(1.d0-s(k))*(1.d0+t(k))/8.d0
h(6)=(1.d0+r(k))*(1.d0-s(k))*(1.d0+t(k))/8.d0
h(7)=(1.d0-r(k))*(1.d0+s(k))*(1.d0+t(k))/8.d0
h(8)=(1.d0+r(k))*(1.d0+s(k))*(1.d0+t(k))/8.d0
lsfp(k,:)=0.d0
do i=1,nlsf
do j=1,8
lsfp(k,i)=lsfp(k,i)+h(j)*lsf(j,i)
enddo
enddo
enddo
r0p=r0+rst*(r(1)+1.d0)/2.d0
s0p=s0+rst*(s(1)+1.d0)/2.d0
t0p=t0+rst*(t(1)+1.d0)/2.d0
rstp=rst/2.d0
levelp=level+1
Dave Whipp
committed
call pressure_cut (params,levelp,levelmax,icon,x,y,z,mat,materialn,u, &
v,w,temp,pressure,strain,nnode,lsfp,nlsf,r0p,s0p, &
t0p,rstp,jcut,ileaves,eviscosity)
enddo
enddo
enddo
deallocate (lsfp)
return
end
!------------------------------------------------------------------------------|
Dave Whipp
committed
!------------------------------------------------------------------------------|