Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! MAKE_CUT Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
recursive subroutine make_cut (level,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,&
Dave Whipp
committed
u,v,w,temp,pressure,strain,e2dp,is_plastic, &
nnode,lsf,nlsf,r0,s0,t0,rst,icut,ileaves, &
Dave Whipp
committed
eviscosity,vbounded,yield_ratio,frict_angle, &
compaction_density,params,threadinfo,weightel, &
matnum,leaf_mat_bin,ematnump)
!------------------------------------------------------------------------------|
!(((((((((((((((( 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
Dave Whipp
committed
! 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
Dave Whipp
committed
! 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
! 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)
! compaction_density: the compacting sediment density of calculated for element ileaves
! matnum : material number, stored for output
! leaf_mat_bin: material bins for current leaf
! ematnump : material number of leaf ileaves from previous time step (derived from cloud)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use threads
use definitions
implicit none
Dave Whipp
committed
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)
Dave Whipp
committed
double precision,intent(in) :: e2dp(nnode)
Dave Whipp
committed
logical,intent(inout) :: is_plastic
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(inout) :: eviscosity
logical,intent(inout) :: vbounded
double precision,intent(inout) :: yield_ratio
double precision,intent(inout) :: frict_angle
double precision,intent(in) :: compaction_density
Dave Whipp
committed
type(thread),intent(inout) :: threadinfo
double precision,intent(inout) :: weightel
integer,intent(in) :: leaf_mat_bin(0:params%nmat)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
Dave Whipp
committed
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)
Dave Whipp
committed
double precision :: prod,vol_lsf0,eviscosityp,eps,weight,yield_ratio_temp
double precision :: xmean,ymean,zmean,frict_angle_temp,density
double precision :: viscosity,expon,activ,fviscosity,fvisc_weak_onset
character(len=8) :: plasticity_type
logical :: is_plastic_temp,vbounded_temp,cloud_for_mat_props
Dave Whipp
committed
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
Dave Whipp
committed
eps=1.d-10
Dave Whipp
committed
icut=0
cloud_for_mat_props=.false.
Dave Whipp
committed
if (level.eq.0) then
ael=0.d0
bel=0.d0
eviscosity=0.d0
frict_angle=0.d0
if (ndof.eq.3) weightel=0.d0
Dave Whipp
committed
matrule=params%matrule
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
xmean=xmean/params%mpe
ymean=ymean/params%mpe
Dave Whipp
committed
zmean=zmean/params%mpe
if (.not.params%excl_vol) then
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
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(params,nlsf,lsf,vol_lsf,vol_lsf0)
! Added option to use cloud particles for material property assignment
! dwhipp/jallen 11/2012
if (params%materials_on_cloud) then
if (vol_lsf0.gt.eps .or. level.gt.0) then
matrule=0
else
if (params%surface_for_mat_props(i)) then
if (params%dommat.gt.0) then
if (vol_lsf(params%dommat).gt.params%domvol) then ! If the volume of dominant material within the element exceeds domvol,
matrule=1 ! use majority rule with the dominant material
vol_lsf(params%dommat) = 1.d0
endif
endif
else
! Assign material properties for given element from cloud
Dave Whipp
committed
! if (ileaves == 1673) then
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) 'calling get_elem_mat_number_from_cloud at location 1'
! write (*,*) 'cloud_for_mat_props: ',cloud_for_mat_props
! write (*,*) 'icut: ',icut
! write (*,*) 'ileaves: ',ileaves
! write (*,*) 'i: ',i
! write (*,*) 'k: ',k
! write (*,*) 'vol_lsf: ',vol_lsf
! write (*,*) 'vol_lsf0: ',vol_lsf0
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) ''
! endif
call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
cloud_for_mat_props=.true.
endif
endif
else
if (vol_lsf0.gt.eps .or. level.gt.0) then
matrule=0
elseif (params%dommat.gt.0) then
if (vol_lsf(params%dommat).gt.params%domvol) then ! If the volume of dominant material within the element exceeds domvol,
matrule=1 ! use majority rule with the dominant material
vol_lsf(params%dommat) = 1.d0
endif
Dave Whipp
committed
endif
endif
! Skip the normal divFEM stuff if cloud was used for material props
if (.not.cloud_for_mat_props) then
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,e2dp,r0,s0,t0,rst, &
xmean,ymean,zmean,compaction_density,ael,bel, &
eviscosity,weightel,yield_ratio,frict_angle, &
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,e2dp,r0,s0,&
t0,rst,xmean,ymean,zmean,compaction_density,&
ael,bel,eviscosity,weightel,yield_ratio, &
frict_angle,is_plastic,vbounded,threadinfo)
call divide_element(params,mat,level,levelmax,ndof,nlsf,nnode, &
icon,kfix,ileaves,leaf_mat_bin,ematnump, &
matnum,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp, &
pressure,strain,e2dp,compaction_density, &
eviscosity,weightel,ael,bel,yield_ratio, &
frict_angle,is_plastic,vbounded,threadinfo)
Dave Whipp
committed
endif
deallocate (vol_lsf)
Dave Whipp
committed
endif
endif
Dave Whipp
committed
if (params%materials_on_cloud .and. level.eq.0) then ! WILL ONLY OCCUR FOR UNCUT ELEMENTS
if (.not.cloud_for_mat_props .and. icut.eq.0 .and. i.eq.nlsf) then
! Assign material properties for given element from cloud
Dave Whipp
committed
! if (ileaves == 1673) then
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) 'calling get_elem_mat_number_from_cloud at location 2'
! write (*,*) 'cloud_for_mat_props: ',cloud_for_mat_props
! write (*,*) 'icut: ',icut
! write (*,*) 'ileaves: ',ileaves
! write (*,*) 'i: ',i
! write (*,*) 'k: ',k
! write (*,*) 'vol_lsf: ',vol_lsf
! write (*,*) 'vol_lsf0: ',vol_lsf0
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) ''
! endif
call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
cloud_for_mat_props=.true.
endif
else
if (prod.lt.0.d0 .and. icut.eq.0) matel=params%materialn(i)
endif
enddo
else
!check whether this is a cut cell
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
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(params,nlsf,lsf,vol_lsf,vol_lsf0)
! Added option to use cloud particles for material property assignment
! dwhipp/jallen 11/2012
if (params%materials_on_cloud) then
if (vol_lsf0.gt.eps .or. level.gt.0) then
matrule=0
else
if (params%surface_for_mat_props(i)) then
if (params%dommat.gt.0) then
if (vol_lsf(params%dommat).gt.params%domvol) then ! If the volume of dominant material within the element exceeds domvol,
matrule=1 ! use majority rule with the dominant material
vol_lsf(params%dommat) = 1.d0
endif
endif
else
! Assign material properties for given element from cloud
Dave Whipp
committed
! if (ileaves == 1673) then
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) 'calling get_elem_mat_number_from_cloud at location 3'
! write (*,*) 'cloud_for_mat_props: ',cloud_for_mat_props
! write (*,*) 'icut: ',icut
! write (*,*) 'ileaves: ',ileaves
! write (*,*) 'i: ',i
! write (*,*) 'k: ',k
! write (*,*) 'vol_lsf: ',vol_lsf
! write (*,*) 'vol_lsf0: ',vol_lsf0
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) ''
! endif
call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
cloud_for_mat_props=.true.
endif
endif
else
if (vol_lsf0.gt.eps .or. level.gt.0) then
matrule=0
elseif (params%dommat.gt.0) then
if (vol_lsf(params%dommat).gt.params%domvol) then ! If the volume of dominant material within the element exceeds domvol,
matrule=1 ! use majority rule with the dominant material
vol_lsf(params%dommat) = 1.d0
endif
! Skip the normal divFEM stuff if cloud was used for material props
if (.not.cloud_for_mat_props) then
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,e2dp,r0,s0,t0,rst, &
xmean,ymean,zmean,compaction_density,ael,bel, &
eviscosity,weightel,yield_ratio,frict_angle, &
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,e2dp,r0,s0,&
t0,rst,xmean,ymean,zmean,compaction_density,&
ael,bel,eviscosity,weightel,yield_ratio, &
frict_angle,is_plastic,vbounded,threadinfo)
call divide_element(params,mat,level,levelmax,ndof,nlsf,nnode, &
icon,kfix,ileaves,leaf_mat_bin,ematnump, &
matnum,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp, &
pressure,strain,e2dp,compaction_density, &
eviscosity,weightel,ael,bel,yield_ratio, &
frict_angle,is_plastic,vbounded,threadinfo)
Dave Whipp
committed
endif
deallocate (vol_lsf)
Dave Whipp
committed
endif
endif
Dave Whipp
committed
enddo
Dave Whipp
committed
if (params%materials_on_cloud .and. level.eq.0) then ! WILL ONLY OCCUR FOR UNCUT ELEMENTS
if (.not.cloud_for_mat_props .and. icut.eq.0) then
! Assign material properties for given element from cloud
Dave Whipp
committed
! if (ileaves == 1673) then
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) 'calling get_elem_mat_number_from_cloud at location 4'
! write (*,*) 'cloud_for_mat_props: ',cloud_for_mat_props
! write (*,*) 'icut: ',icut
! write (*,*) 'ileaves: ',ileaves
! write (*,*) 'i: ',i
! write (*,*) 'k: ',k
! write (*,*) 'vol_lsf: ',vol_lsf
! write (*,*) 'vol_lsf0: ',vol_lsf0
! write (*,*) '-----------------------------------------------------------------------'
! write (*,*) ''
! endif
call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
cloud_for_mat_props=.true.
endif
else
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) matel=params%materialn(i)
enddo
endif
Dave Whipp
committed
endif
endif
! Check for overriding material transitions
!write (*,*) 'Calling mat_trans for an uncut element'
call mat_trans(params,mat,matel,xmean,ymean,zmean)
!write (*,*) 'Finished calling mat_trans for an uncut element'
matnum=matel
! If the element is uncut and using divFEM, or the cloud was used for material
! property assignment, call make_matrix here
if (icut.eq.0 .or. cloud_for_mat_props) then
if (params%sstemp_spinup .and. params%first_step .and. params%sstemp_viscosity_spinup>=eps) then
! Use a linear viscous material for the first time step if using a
! steady-state thermal calculation
viscosity = params%sstemp_viscosity_spinup
plasticity_type = 'No'
fviscosity = 1.d0
fvisc_weak_onset = -1.d0
activ = 0.d0
expon = 1.d0
else
viscosity = mat(matel)%viscosity
plasticity_type = mat(matel)%plasticity_type
fviscosity = mat(matel)%fviscosity
fvisc_weak_onset = mat(matel)%fvisc_weak_onset
activ = mat(matel)%activationenergy
expon = mat(matel)%expon
endif
if (params%compaction) then
density = compaction_density
else
density = mat(matel)%density
endif
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,viscosity, &
density,mat(matel)%penalty,expon,activ, &
mat(matel)%expansion,mat(matel)%diffusivity,mat(matel)%heat,&
plasticity_type,mat(matel)%plasticity_ss_type_coh, &
mat(matel)%plasticity_parameters,fviscosity, &
mat(matel)%fvisc_weak_type,fvisc_weak_onset, &
mat(matel)%fvisc_weak_end,mat(matel)%fvisc_weak_final,u,v, &
w,temp,pressure,strain,e2dp,is_plastic_temp,nnode,r0,s0,t0, &
rst,ileaves,eviscosityp,vbounded_temp,yield_ratio_temp, &
frict_angle_temp,threadinfo,weight)
Dave Whipp
committed
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)
yield_ratio=yield_ratio+yield_ratio_temp/(8.d0**level)
frict_angle=frict_angle+frict_angle_temp/(8.d0**level)
Dave Whipp
committed
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
endif
Dave Whipp
committed
end subroutine make_cut
Dave Whipp
committed
!-------------------------------------------------------------------------------
Dave Whipp
committed
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
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=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
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, &
Dave Whipp
committed
ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,w, &
Dave Whipp
committed
temp,pressure,strain,e2dp,r0,s0,t0,rst,xmean,ymean, &
zmean,compaction_density,ael,bel,eviscosity,weightel,&
yield_ratio,frict_angle,is_plastic,vbounded, &
threadinfo)
Dave Whipp
committed
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
Dave Whipp
committed
integer,intent(out) :: matnum
Dave Whipp
committed
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)
Dave Whipp
committed
double precision,intent(in) :: e2dp(nnode)
Dave Whipp
committed
double precision,intent(in) :: r0,s0,t0,rst
Dave Whipp
committed
double precision,intent(in) :: xmean,ymean,zmean
double precision,intent(in) :: compaction_density
Dave Whipp
committed
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
double precision,intent(out) :: yield_ratio
double precision,intent(out) :: frict_angle
Dave Whipp
committed
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,yield_ratio_temp,frict_angle_temp
Dave Whipp
committed
double precision :: plasticity_parameters(9)
double precision :: aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision :: fviscosity,fvisc_weak_onset,fvisc_weak_end,fvisc_weak_final
Dave Whipp
committed
integer :: matel,i
character(len=8) :: plasticity_type
Dave Whipp
committed
character(len=16) :: plasticity_ss_type_coh,plasticity_ss_type_phi
character(len=16) :: fvisc_weak_type
Dave Whipp
committed
logical :: is_plastic_temp,vbounded_temp
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
eps=1.d-10
Dave Whipp
committed
select case (matrule)
case (1) ! Assign material properties of volumetric majority material
matel=params%materialn(sum(maxloc(vol_lsf)))
! Check for overriding material transitions
!write (*,*) 'Calling mat_trans for an cut element (majority rule)'
call mat_trans(params,mat,matel,xmean,ymean,zmean)
!write (*,*) 'Finished calling mat_trans for an cut element (majority rule)'
Dave Whipp
committed
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_ss_type_coh=mat(matel)%plasticity_ss_type_coh
plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
Dave Whipp
committed
plasticity_parameters=mat(matel)%plasticity_parameters
fviscosity=mat(matel)%fviscosity
fvisc_weak_type=mat(matel)%fvisc_weak_type
fvisc_weak_onset=mat(matel)%fvisc_weak_onset
fvisc_weak_end=mat(matel)%fvisc_weak_end
fvisc_weak_final=mat(matel)%fvisc_weak_final
Dave Whipp
committed
case(2) ! Assign material properties of volumetric minority material
matel=params%materialn(sum(minloc(vol_lsf)))
! Check for overriding material transitions
!write (*,*) 'Calling mat_trans for an cut element (minority rule)'
call mat_trans(params,mat,matel,xmean,ymean,zmean)
!write (*,*) 'Finished calling mat_trans for an cut element (minority rule)'
Dave Whipp
committed
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_ss_type_coh=mat(matel)%plasticity_ss_type_coh
plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
Dave Whipp
committed
plasticity_parameters=mat(matel)%plasticity_parameters
fviscosity=mat(matel)%fviscosity
fvisc_weak_type=mat(matel)%fvisc_weak_type
fvisc_weak_onset=mat(matel)%fvisc_weak_onset
fvisc_weak_end=mat(matel)%fvisc_weak_end
fvisc_weak_final=mat(matel)%fvisc_weak_final
Dave Whipp
committed
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)
! Check for overriding material transitions
!write (*,*) 'Calling mat_trans for an cut element (divFEM)'
call mat_trans(params,mat,matel,xmean,ymean,zmean)
!write (*,*) 'Finished calling mat_trans for an cut element (divFEM)'
if (vol_lsf(i).ge.0.015625-eps) matnum=matnum+matel
Dave Whipp
committed
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
Dave Whipp
committed
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_ss_type_coh=mat(matel)%plasticity_ss_type_coh
plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
Dave Whipp
committed
plasticity_parameters=mat(matel)%plasticity_parameters
fviscosity=mat(matel)%fviscosity
fvisc_weak_type=mat(matel)%fvisc_weak_type
fvisc_weak_onset=mat(matel)%fvisc_weak_onset
fvisc_weak_end=mat(matel)%fvisc_weak_end
fvisc_weak_final=mat(matel)%fvisc_weak_final
Dave Whipp
committed
end select
if (maxval(vol_lsf).lt.eps) then
matel=params%materialn(0)
plasticity_type=mat(matel)%plasticity_type
plasticity_ss_type_coh=mat(matel)%plasticity_ss_type_coh
plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
Dave Whipp
committed
plasticity_parameters=mat(matel)%plasticity_parameters
fvisc_weak_type=mat(matel)%fvisc_weak_type
fvisc_weak_onset=mat(matel)%fvisc_weak_onset
fvisc_weak_end=mat(matel)%fvisc_weak_end
fvisc_weak_final=mat(matel)%fvisc_weak_final
Dave Whipp
committed
endif
if (params%sstemp_spinup .and. params%first_step .and. params%sstemp_viscosity_spinup>=eps) then
! Use a linear viscous material for the first time step if using a
! steady-state thermal calculation
viscosity = params%sstemp_viscosity_spinup
plasticity_type = 'No'
fviscosity = 1.d0
fvisc_weak_onset = -1.d0
activ = 0.d0
expon = 1.d0
endif
if (params%compaction) density = compaction_density
Dave Whipp
committed
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,viscosity,density, &
penal,expon,activ,expan,diffusivity,heat,plasticity_type, &
plasticity_parameters,fviscosity,fvisc_weak_type, &
fvisc_weak_onset,fvisc_weak_end,fvisc_weak_final,u,v,w,temp, &
pressure,strain,e2dp,is_plastic_temp,nnode,r0,s0,t0,rst, &
ileaves,eviscosityp,vbounded_temp,yield_ratio_temp, &
frict_angle_temp,threadinfo,weight)
is_plastic=(is_plastic.or.is_plastic_temp)
vbounded=(vbounded.or.vbounded_temp)
Dave Whipp
committed
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
eviscosity=eviscosity+eviscosityp/(8.d0**level)
yield_ratio=yield_ratio+yield_ratio_temp/(8.d0**level)
frict_angle=frict_angle+frict_angle_temp/(8.d0**level)
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
Dave Whipp
committed
end subroutine set_elem_props
!-------------------------------------------------------------------------------
subroutine divide_element(params,mat,level,levelmax,ndof,nlsf,nnode,icon,kfix, &
ileaves,leaf_mat_bin,ematnump,matnum,lsf,r0,s0,t0,rst,&
x,y,z,u,v,w,temp,pressure,strain,e2dp, &
compaction_density,eviscosity,weightel,ael,bel, &
yield_ratio,frict_angle,is_plastic,vbounded,threadinfo)
Dave Whipp
committed
use threads
use definitions
implicit none
Dave Whipp
committed
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(in) :: leaf_mat_bin(0:params%nmat)
integer,intent(in) :: ematnump
Dave Whipp
committed
integer,intent(inout) :: matnum
Dave Whipp
committed
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)
Dave Whipp
committed
double precision,intent(in) :: e2dp(nnode)
double precision,intent(in) :: compaction_density
Dave Whipp
committed
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) :: yield_ratio
double precision,intent(inout) :: frict_angle
Dave Whipp
committed
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
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
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
Dave Whipp
committed
call make_cut (levelp,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,u,v,w,&
temp,pressure,strain,e2dp,is_plastic,nnode,lsfp,nlsf, &
Dave Whipp
committed
r0p,s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded, &
yield_ratio,frict_angle,compaction_density,params, &
threadinfo,weightel,matnum,leaf_mat_bin,ematnump)
Dave Whipp
committed
end subroutine divide_element
!-------------------------------------------------------------------------------
subroutine mat_trans (params,mat,matel,xmean,ymean,zmean)
use definitions
implicit none
type(parameters),intent(in) :: params
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(inout) :: matel
double precision,intent(in) :: xmean,ymean,zmean
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
logical :: checkmattrans
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
checkmattrans=.true.
!write (*,*) 'checkmattrans at start of call: ',checkmattrans
do while (checkmattrans)
checkmattrans=.false.
!write (*,*) 'checkmattrans at start of loop: ',checkmattrans
if (mat(matel)%mattrans(1).gt.0.d0 .and. xmean.lt.mat(matel)%mattrans(1)) then
checkmattrans=.true.
!write (*,*) 'checkmattrans for xmean lt mattrans: ',checkmattrans
matel=mat(matel)%transnum(1)
endif
if (mat(matel)%mattrans(2).gt.0.d0 .and. xmean.gt.mat(matel)%mattrans(2)) then
checkmattrans=.true.
!write (*,*) 'checkmattrans for xmean gt mattrans: ',checkmattrans
matel=mat(matel)%transnum(2)
endif
if (mat(matel)%mattrans(3).gt.0.d0 .and. ymean.lt.mat(matel)%mattrans(3)) then
checkmattrans=.true.
!write (*,*) 'checkmattrans for ymean lt mattrans: ',checkmattrans
matel=mat(matel)%transnum(3)
endif
if (mat(matel)%mattrans(4).gt.0.d0 .and. ymean.gt.mat(matel)%mattrans(4)) then
checkmattrans=.true.
!write (*,*) 'checkmattrans for ymean gt mattrans: ',checkmattrans
matel=mat(matel)%transnum(4)
endif
if (mat(matel)%mattrans(5).gt.0.d0 .and. zmean.lt.mat(matel)%mattrans(5)) then
checkmattrans=.true.
!write (*,*) 'checkmattrans for zmean lt mattrans: ',checkmattrans
matel=mat(matel)%transnum(5)
endif
if (mat(matel)%mattrans(6).gt.0.d0 .and. zmean.gt.mat(matel)%mattrans(6)) then
checkmattrans=.true.
!write (*,*) 'checkmattrans for zmean gt mattrans: ',checkmattrans
matel=mat(matel)%transnum(6)
endif
!write (*,*) 'checkmattrans at end of loop: ',checkmattrans
enddo
!write (*,*) 'checkmattrans at end of call: ',checkmattrans
end subroutine mat_trans
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------