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)
!------------------------------------------------------------------------------|
!(((((((((((((((( 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
Dave Whipp
committed
type(parameters),intent(in) :: params
integer,intent(in) :: level
integer,intent(in) :: levelmax
integer,intent(in) :: icon(params%mpe)
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(in) :: materialn(0:nlsf)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(inout) :: pressure
double precision,intent(in) :: strain(nnode)
integer,intent(in) :: nnode
double precision,intent(in) :: lsf(params%mpe,nlsf)
integer,intent(in) :: nlsf
double precision,intent(in) :: r0,s0,t0,rst
integer,intent(out) :: icut
integer,intent(in) :: ileaves
double precision,intent(in) :: eviscosity
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
Dave Whipp
committed
integer :: matel,matrule,i,j,k,err
double precision :: prod,vol_lsf0,eps,pressurep
double precision,allocatable :: vol_lsf(:)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
Dave Whipp
committed
eps=1.d-10
Dave Whipp
committed
icut=0
Dave Whipp
committed
Dave Whipp
committed
pressure=0.d0
Dave Whipp
committed
matrule=params%matrule
Dave Whipp
committed
if (prod*lsf(k,i).le.0.d0 .and. icut.eq.0) then
icut=1
allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in make_cut$')
call get_mat_volume_pc(params,nlsf,lsf,vol_lsf,vol_lsf0)
Dave Whipp
committed
if (vol_lsf0.gt.eps .or. level.gt.0) matrule=0
if (matrule.eq.1 .or. matrule.eq.2) then ! Use majority or minority rule
call set_elem_props_pc(params,mat,matrule,nlsf,nnode,icon,level, &
materialn,vol_lsf,vol_lsf0,x,y,z,u,v,w,temp, &
pressure,strain,r0,s0,t0,rst,eviscosity)
Dave Whipp
committed
else ! Use divFEM
if (level.eq.levelmax) then
call set_elem_props_pc(params,mat,matrule,nlsf,nnode,icon,level, &
materialn,vol_lsf,vol_lsf0,x,y,z,u,v,w,temp, &
pressure,strain,r0,s0,t0,rst,eviscosity)
Dave Whipp
committed
else
call divide_element_pc(params,mat,level,levelmax,nlsf,nnode,icon, &
ileaves,materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,&
w,temp,pressure,strain,eviscosity)
Dave Whipp
committed
endif
endif
deallocate (vol_lsf)
endif
Dave Whipp
committed
if (prod.lt.0.d0 .and. icut.eq.0) then
matel=materialn(i)
Dave Whipp
committed
if (icut.eq.0) then
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, &
Dave Whipp
committed
pressure=pressure+pressurep/(8.d0**level)
endif
Dave Whipp
committed
end subroutine pressure_cut
Dave Whipp
committed
!-------------------------------------------------------------------------------
subroutine get_mat_volume_pc (params,nlsf,lsf,vol_lsf,vol_lsf0)
Dave Whipp
committed
use definitions
Dave Whipp
committed
implicit none
Dave Whipp
committed
type(parameters),intent(in) :: params
integer,intent(in) :: nlsf
double precision,intent(in) :: lsf(params%mpe,nlsf)
double precision,intent(out) :: vol_lsf(nlsf)
double precision,intent(out) :: vol_lsf0
Dave Whipp
committed
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: i
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
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=nlsf-1,1,-1
vol_lsf(i)=max(vol_lsf(i),vol_lsf(i+1))
enddo
endif
Dave Whipp
committed
if (.not.params%excl_vol) then
do i=1,nlsf-1
vol_lsf(i)=vol_lsf(i)-vol_lsf(i+1)
enddo
Dave Whipp
committed
! this is a little fix that is necessary due to the finite precision with which
! we compute the positive volumes and could lead to a contributing volume being
! either marginally negative or greater than 1
do i=1,nlsf
vol_lsf(i)=max(vol_lsf(i),0.d0)
vol_lsf(i)=min(vol_lsf(i),1.d0)
enddo
vol_lsf0=1.d0-sum(vol_lsf)
end subroutine get_mat_volume_pc
Dave Whipp
committed
!-------------------------------------------------------------------------------
subroutine set_elem_props_pc (params,mat,matrule,nlsf,nnode,icon,level, &
materialn,vol_lsf,vol_lsf0,x,y,z,u,v,w,temp, &
pressure,strain,r0,s0,t0,rst,eviscosity)
Dave Whipp
committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
use definitions
implicit none
type(parameters),intent(in) :: params
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(in) :: matrule
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: level
integer,intent(in) :: materialn(0:nlsf)
double precision,intent(in) :: vol_lsf(nlsf)
double precision,intent(in) :: vol_lsf0
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: r0,s0,t0,rst
double precision,intent(in) :: eviscosity
double precision,intent(out) :: pressure
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: matel,i
Dave Whipp
committed
double precision :: eps,viscosity,penal,expon,pressurep,xmean,ymean,zmean
Dave Whipp
committed
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
eps=1.d-10
Dave Whipp
committed
xmean=0.d0
ymean=0.d0
Dave Whipp
committed
zmean=0.d0
do i=1,params%mpe
Dave Whipp
committed
xmean=xmean+x(icon(i))
ymean=ymean+y(icon(i))
Dave Whipp
committed
zmean=zmean+z(icon(i))*params%vex
enddo
Dave Whipp
committed
select case (matrule)
case (1) ! Assign material properties of volumetric majority material
matel=materialn(sum(maxloc(vol_lsf)))
Dave Whipp
committed
if (xmean.lt.mat(matel)%mattrans(1)) matel=mat(matel)%transnum(1)
if (xmean.gt.mat(matel)%mattrans(2)) matel=mat(matel)%transnum(2)
if (ymean.lt.mat(matel)%mattrans(3)) matel=mat(matel)%transnum(3)
if (ymean.gt.mat(matel)%mattrans(4)) matel=mat(matel)%transnum(4)
if (zmean.lt.mat(matel)%mattrans(5)) matel=mat(matel)%transnum(5)
if (zmean.gt.mat(matel)%mattrans(6)) matel=mat(matel)%transnum(6)
Dave Whipp
committed
viscosity=mat(matel)%viscosity
penal=mat(matel)%penalty
expon=mat(matel)%expon
case(2) ! Assign material properties of volumetric minority material
matel=materialn(sum(minloc(vol_lsf)))
Dave Whipp
committed
if (xmean.lt.mat(matel)%mattrans(1)) matel=mat(matel)%transnum(1)
if (xmean.gt.mat(matel)%mattrans(2)) matel=mat(matel)%transnum(2)
if (ymean.lt.mat(matel)%mattrans(3)) matel=mat(matel)%transnum(3)
if (ymean.gt.mat(matel)%mattrans(4)) matel=mat(matel)%transnum(4)
if (zmean.lt.mat(matel)%mattrans(5)) matel=mat(matel)%transnum(5)
if (zmean.gt.mat(matel)%mattrans(6)) matel=mat(matel)%transnum(6)
Dave Whipp
committed
viscosity=mat(matel)%viscosity
penal=mat(matel)%penalty
expon=mat(matel)%expon
case default
if (matrule.ne.0) write(*,*) 'Invalid matrule value, using divFEM'
viscosity=0.d0
penal=0.d0
expon=0.d0
do i=1,nlsf
matel=materialn(i)
Dave Whipp
committed
if (xmean.lt.mat(matel)%mattrans(1)) matel=mat(matel)%transnum(1)
if (xmean.gt.mat(matel)%mattrans(2)) matel=mat(matel)%transnum(2)
if (ymean.lt.mat(matel)%mattrans(3)) matel=mat(matel)%transnum(3)
if (ymean.gt.mat(matel)%mattrans(4)) matel=mat(matel)%transnum(4)
if (zmean.lt.mat(matel)%mattrans(5)) matel=mat(matel)%transnum(5)
if (zmean.gt.mat(matel)%mattrans(6)) matel=mat(matel)%transnum(6)
Dave Whipp
committed
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
matel=materialn(sum(maxloc(vol_lsf)))
end select
call make_pressure (params,icon,x,y,z,viscosity,penal,expon,u,v,w,temp, &
pressurep,strain,nnode,r0,s0,t0,rst,eviscosity)
Dave Whipp
committed
pressure=pressure+pressurep/(8.d0**level)
end subroutine set_elem_props_pc
Dave Whipp
committed
!-------------------------------------------------------------------------------
subroutine divide_element_pc(params,mat,level,levelmax,nlsf,nnode,icon,ileaves,&
materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp, &
pressure,strain,eviscosity)
Dave Whipp
committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
use definitions
implicit none
type(parameters),intent(in) :: params
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(in) :: level
integer,intent(in) :: levelmax
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: ileaves
integer,intent(in) :: materialn(0:nlsf)
double precision,intent(in) :: lsf(params%mpe,nlsf)
double precision,intent(in) :: r0,s0,t0,rst
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: eviscosity
double precision,intent(inout) :: pressure
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: err,i,j,k,ii,jj,kk,levelp,jcut
double precision,allocatable :: lsfp(:,:)
double precision :: r(params%mpe),s(params%mpe),t(params%mpe),h(params%mpe)
double precision :: r0p,s0p,t0p,rstp
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in make_cut$')
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
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
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)
end subroutine divide_element_pc
Dave Whipp
committed
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------