!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! MAKE_CUT Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
recursive subroutine make_cut (level,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,&
u,v,w,temp,pressure,strain,is_plastic,nnode,f, &
lsf,nlsf,r0,s0,t0,rst,icut,ileaves,eviscosity, &
vbounded,params,threadinfo,weightel,matnum)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this subroutine is a intermediary routine between build_system and make_matrix
! to take into account the complex geometry of cut cells
! if we are in a non cut cell, make-matrix is called
! if we are in a cut cell but at a level that it smaller than levelmax, the
! cell is further cut and make_cut is recursively called
! if we are in a cut cell and level is equal to levelmax, we call make_matrix
! with material properties that have been interpolated from the
! various material properties contribnuting to the cut cell
! level : current level in cut cell algorithm. It varies between 0 and levelmax.
! levelapprox : used to improve the postitive volutme calculation by further division
! mpe : number of nodes per element (8)
! ndof : number of degrees of freedom per node (3)
! ael : computed finite element matrix
! bel : computed rhs vector
! icon : connectivity array for the current element
! xg,yg,zg : global coordinate arrays of length nnode
! penalty : penalty factor used to impose the boundary conditions
! tempscale : temperature scaling parameter
! kfix : bc array of length ndof*nnode (kfix=1 means the dof is fixed to the value stored in velo)
! mat : material matrix for the nmat materials
! materialn : contains the material number associated to each lsf
! dt : time step length (only needed for the temperature calculations)
! u,v,w : velocity array (obtained from previous time step or at least containing the proper velocity at the fixed dofs)
! nnode : number of nodes
! f : global rhs vector
! lsf : global array of level set functions defining the surfaces
! nlsf : number of lsfs
! r0,s0,t0 : bottom left back corner of the part of the element. (we are now integrating in local r,s,t coordinates)
! rst : size of the part of the element we are integrating
! icut : returned, 0 if homogeneous element - 1 if cut element
! ileaves : current leaf number (useless except for debugging)
! matnum : material number, stored for output
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use threads
use definitions
implicit none
type(parameters),intent(in) :: params
integer,intent(in) :: level
integer,intent(in) :: levelmax
integer,intent(in) :: ndof
double precision,intent(inout) :: ael(params%mpe*ndof,params%mpe*ndof)
double precision,intent(inout) :: bel(params%mpe*ndof)
integer,intent(in) :: icon(params%mpe)
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
integer,intent(in) :: kfix(nnode*ndof)
type(material),intent(in) :: mat(0:params%nmat)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(in) :: pressure
double precision,intent(in) :: strain(nnode)
logical,intent(inout) :: is_plastic
integer,intent(in) :: nnode
double precision,intent(inout) :: f(nnode*ndof)
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(inout) :: eviscosity
logical,intent(inout) :: vbounded
type(thread),intent(inout) :: threadinfo
double precision,intent(inout) :: weightel
integer,intent(inout) :: matnum
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer :: matel,matrule,i,j,k,err
double precision,allocatable :: vol_lsf(:)
double precision :: aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision :: prod,vol_lsf0,eviscosityp,eps,weight,zmean
logical :: is_plastic_temp,vbounded_temp
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
eps=1.d-10
icut=0
if (level.eq.0) then
ael=0.d0
bel=0.d0
eviscosity=0.d0
if (ndof.eq.3) weightel=0.d0
endif
matel=params%materialn(0)
matrule=params%matrule
zmean=0.d0
do i=1,params%mpe
zmean=zmean+z(icon(i))*params%vex
enddo
zmean=zmean/params%mpe
if (.not.params%excl_vol) then
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
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(params,nlsf,lsf,vol_lsf,vol_lsf0)
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(params,mat,matrule,ndof,nlsf,nnode,icon,kfix, &
ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,&
w,temp,pressure,strain,r0,s0,t0,rst,zmean,f,ael,&
bel,eviscosity,weightel,is_plastic,vbounded, &
threadinfo)
else ! Use divFEM
if (level.eq.levelmax) then
call set_elem_props(params,mat,matrule,ndof,nlsf,nnode,icon,kfix,&
ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,&
v,w,temp,pressure,strain,r0,s0,t0,rst,zmean,f,&
ael,bel,eviscosity,weightel,is_plastic, &
vbounded,threadinfo)
else
call divide_element(params,mat,level,levelmax,ndof,nlsf,nnode, &
icon,kfix,ileaves,matnum,lsf,r0,s0,t0,rst,x,y,&
z,u,v,w,temp,pressure,strain,eviscosity, &
weightel,ael,bel,f,is_plastic,vbounded, &
threadinfo)
endif
endif
deallocate (vol_lsf)
endif
enddo
if (prod.lt.0.d0 .and. icut.eq.0) then
matel=params%materialn(i)
if (zmean.lt.mat(matel)%ztrans) matel=mat(matel)%transnum
matnum=matel
endif
enddo
else
!check whether this is a cut cell
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
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(params,nlsf,lsf,vol_lsf,vol_lsf0)
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(params,mat,matrule,ndof,nlsf,nnode,icon,kfix, &
ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,&
w,temp,pressure,strain,r0,s0,t0,rst,zmean,f,ael,&
bel,eviscosity,weightel,is_plastic,vbounded, &
threadinfo)
else ! Use divFEM
if (level.eq.levelmax) then
call set_elem_props(params,mat,matrule,ndof,nlsf,nnode,icon,kfix,&
ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,&
v,w,temp,pressure,strain,r0,s0,t0,rst,zmean,f,&
ael,bel,eviscosity,weightel,is_plastic, &
vbounded,threadinfo)
else
call divide_element(params,mat,level,levelmax,ndof,nlsf,nnode, &
icon,kfix,ileaves,matnum,lsf,r0,s0,t0,rst,x,y,&
z,u,v,w,temp,pressure,strain,eviscosity, &
weightel,ael,bel,f,is_plastic,vbounded, &
threadinfo)
endif
endif
deallocate (vol_lsf)
endif
enddo
enddo
if (icut.eq.0) then
!assign material to plain cell, since at that point we know that this is a plain cell
do i=1,nlsf
if (lsf(1,i).lt.0.d0) then
matel=params%materialn(i)
if (zmean.lt.mat(matel)%ztrans) matel=mat(matel)%transnum
matnum=matel
endif
enddo
endif
endif
if (icut.eq.0) then
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,mat(matel)%viscosity,&
mat(matel)%density,mat(matel)%penalty,mat(matel)%expon, &
mat(matel)%activationenergy,mat(matel)%expansion, &
mat(matel)%diffusivity,mat(matel)%heat, &
mat(matel)%plasticity_type,mat(matel)%plasticity_parameters,&
u,v,w,temp,pressure,strain,is_plastic_temp,nnode,f,r0,s0,t0,&
rst,ileaves,eviscosityp,vbounded_temp,threadinfo,weight)
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
eviscosity=eviscosity+eviscosityp/(8.d0**level)
is_plastic=(is_plastic.or.is_plastic_temp)
vbounded=(vbounded.or.vbounded_temp)
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
endif
end subroutine make_cut
!-------------------------------------------------------------------------------
subroutine get_mat_volume (params,nlsf,lsf,vol_lsf,vol_lsf0)
use definitions
implicit none
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
!------------------------------------------------------------------------------|
!(((((((((((((((( 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=1,nlsf-1
vol_lsf(i)=vol_lsf(i)-vol_lsf(i+1)
enddo
endif
! 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
!-------------------------------------------------------------------------------
subroutine set_elem_props (params,mat,matrule,ndof,nlsf,nnode,icon,kfix, &
ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,w, &
temp,pressure,strain,r0,s0,t0,rst,zmean,f,ael,bel, &
eviscosity,weightel,is_plastic,vbounded,threadinfo)
use threads
use definitions
implicit none
type(parameters),intent(in) :: params
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(in) :: matrule
integer,intent(in) :: ndof
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: kfix(nnode*ndof)
integer,intent(in) :: ileaves
integer,intent(in) :: level
integer,intent(out) :: matnum
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) :: pressure
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: r0,s0,t0,rst
double precision,intent(in) :: zmean
double precision,intent(inout) :: f(nnode*ndof)
double precision,intent(out) :: ael(params%mpe*ndof,params%mpe*ndof)
double precision,intent(out) :: bel(params%mpe*ndof)
double precision,intent(out) :: eviscosity
double precision,intent(out) :: weightel
logical,intent(out) :: is_plastic
logical,intent(out) :: vbounded
type(thread),intent(inout) :: threadinfo
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision :: eps,viscosity,density,penal,expon,activ,expan,diffusivity
double precision :: heat,eviscosityp,weight
double precision :: plasticity_parameters(9)
double precision :: aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
integer :: matel,i
character(len=8) :: plasticity_type
logical :: is_plastic_temp,vbounded_temp
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
eps=1.d-10
matnum=0
select case (matrule)
case (1) ! Assign material properties of volumetric majority material
matel=params%materialn(sum(maxloc(vol_lsf)))
if (zmean.lt.mat(matel)%ztrans) matel=mat(matel)%transnum
matnum=matel
viscosity=mat(matel)%viscosity
density=mat(matel)%density
penal=mat(matel)%penalty
expon=mat(matel)%expon
activ=mat(matel)%activationenergy
expan=mat(matel)%expansion
diffusivity=mat(matel)%diffusivity
heat=mat(matel)%heat
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
case(2) ! Assign material properties of volumetric minority material
matel=params%materialn(sum(minloc(vol_lsf)))
if (zmean.lt.mat(matel)%ztrans) matel=mat(matel)%transnum
matnum=matel
viscosity=mat(matel)%viscosity
density=mat(matel)%density
penal=mat(matel)%penalty
expon=mat(matel)%expon
activ=mat(matel)%activationenergy
expan=mat(matel)%expansion
diffusivity=mat(matel)%diffusivity
heat=mat(matel)%heat
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
case default
if (matrule.ne.0) write(*,*) 'Invalid matrule value, using divFEM'
viscosity=0.d0
density=0.d0
penal=0.d0
expon=0.d0
activ=0.d0
expan=0.d0
diffusivity=0.d0
heat=0.d0
do i=1,nlsf
matel=params%materialn(i)
if (zmean.lt.mat(matel)%ztrans) matel=mat(matel)%transnum
if (vol_lsf(i).ge.0.015625-eps) matnum=matnum+matel
viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
density=density+vol_lsf(i)*mat(matel)%density
penal=penal+vol_lsf(i)*mat(matel)%penalty
expon=expon+vol_lsf(i)*mat(matel)%expon
activ=activ+vol_lsf(i)*mat(matel)%activationenergy
expan=expan+vol_lsf(i)*mat(matel)%expansion
diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
heat=heat+vol_lsf(i)*mat(matel)%heat
enddo
matel=params%materialn(0)
if (vol_lsf0.ge.0.015625-eps) matnum=matnum+100
viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
density=density+vol_lsf0*mat(matel)%density
penal=penal+vol_lsf0*mat(matel)%penalty
expon=expon+vol_lsf0*mat(matel)%expon
activ=activ+vol_lsf0*mat(matel)%activationenergy
expan=expan+vol_lsf0*mat(matel)%expansion
heat=heat+vol_lsf0*mat(matel)%heat
diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
diffusivity=1.d0/diffusivity ! Note that some properties add geometrically not algebraically
matel=params%materialn(sum(maxloc(vol_lsf)))
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
end select
if (maxval(vol_lsf).lt.eps) then
matel=params%materialn(0)
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
endif
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,viscosity,density, &
penal,expon,activ,expan,diffusivity,heat,plasticity_type, &
plasticity_parameters,u,v,w,temp,pressure,strain, &
is_plastic_temp,nnode,f,r0,s0,t0,rst,ileaves,eviscosityp, &
vbounded_temp,threadinfo,weight)
is_plastic=(is_plastic.or.is_plastic_temp)
vbounded=(vbounded.or.vbounded_temp)
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
eviscosity=eviscosity+eviscosityp/(8.d0**level)
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
end subroutine set_elem_props
!-------------------------------------------------------------------------------
subroutine divide_element(params,mat,level,levelmax,ndof,nlsf,nnode,icon,kfix, &
ileaves,matnum,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp, &
pressure,strain,eviscosity,weightel,ael,bel,f, &
is_plastic,vbounded,threadinfo)
use threads
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) :: ndof
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: kfix(nnode*ndof)
integer,intent(in) :: ileaves
integer,intent(inout) :: matnum
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) :: pressure
double precision,intent(in) :: strain(nnode)
double precision,intent(inout) :: eviscosity
double precision,intent(inout) :: weightel
double precision,intent(inout) :: ael(params%mpe*ndof,params%mpe*ndof)
double precision,intent(inout) :: bel(params%mpe*ndof)
double precision,intent(inout) :: f(nnode*ndof)
logical,intent(inout) :: is_plastic
logical,intent(inout) :: vbounded
type(thread),intent(inout) :: threadinfo
!------------------------------------------------------------------------------|
!(((((((((((((((( 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$')
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 make_cut (levelp,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,u,v,w,&
temp,pressure,strain,is_plastic,nnode,f,lsfp,nlsf,r0p, &
s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded,params, &
threadinfo,weightel,matnum)
enddo
enddo
enddo
deallocate (lsfp)
end subroutine divide_element
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------