!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! PRESSURE_CUT Nov. 2006 I |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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
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
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)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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
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
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
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
endif
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
allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in pressure_cut$')
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
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|