From cc0e8354c5091219270590f10c0f0c885ef80dcf Mon Sep 17 00:00:00 2001
From: Dave Whipp <dwhipp@dal.ca>
Date: Mon, 10 Jan 2011 16:23:15 +0000
Subject: [PATCH] Updated to do a proper majority rule for full elements.
Similar to the make_cut.f90 update
---
src/pressure_cut.f90 | 364 +++++++++++++++++++++++++++----------------
1 file changed, 233 insertions(+), 131 deletions(-)
diff --git a/src/pressure_cut.f90 b/src/pressure_cut.f90
index 979eb054..c8a899d9 100644
--- a/src/pressure_cut.f90
+++ b/src/pressure_cut.f90
@@ -20,7 +20,6 @@ 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 ))))))))))))))))))))))))))))))))))))))
@@ -53,166 +52,268 @@ 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
+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 )))))))))))))
!------------------------------------------------------------------------------|
-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)
+integer :: matel,matrule,i,j,k,err
+double precision :: prod,vol_lsf0,eps,pressurep
+double precision,allocatable :: vol_lsf(:)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
eps=1.d-10
+icut=0
if (level.eq.0) then
- pressure=0.d0
+ pressure=0.d0
endif
matel=materialn(0)
+matrule=params%matrule
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
- if (prod*lsf(k,i).le.0.d0) goto 222
+ 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,nlsf,nnode,icon,level, &
+ materialn,vol_lsf,vol_lsf0,x,y,z,u,v,w,temp, &
+ pressure,strain,r0,s0,t0,rst,eviscosity)
+ else ! Use divFEM
+ if (level.eq.levelmax) then
+ call set_elem_props(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)
+ else
+ call divide_element(params,mat,level,levelmax,nlsf,nnode,icon, &
+ ileaves,materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,w, &
+ temp,pressure,strain,eviscosity)
+ endif
+ endif
+ deallocate (vol_lsf)
+ endif
enddo
- if (prod.lt.0.d0) then
- matel=materialn(i)
+ if (prod.lt.0.d0 .and. icut.eq.0) 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)
+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, &
+ t0,rst,mat(matel)%plasticity_type, &
+ mat(matel)%plasticity_parameters,eviscosity)
-pressure=pressure+pressurep/(8.d0**level)
+ pressure=pressure+pressurep/(8.d0**level)
+endif
-icut=0
+end subroutine pressure_cut
-return
+!-------------------------------------------------------------------------------
-222 continue
+subroutine get_mat_volume (params,nlsf,lsf,vol_lsf,vol_lsf0)
-icut=1
+use definitions
-! when we get to the bottom of the division
-! we use an approximate algorithm
+implicit none
-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
+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
- 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
+!------------------------------------------------------------------------------|
+!(((((((((((((((( 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
-! If we are not at the bottom level, we keep dividing
+! 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
+
+!-------------------------------------------------------------------------------
-allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in pressure_cut$')
+subroutine set_elem_props (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)
+
+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
+double precision :: eps,viscosity,penal,expon,pressurep
+double precision :: plasticity_parameters(9)
+character(len=8) :: plasticity_type
+
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+
+eps=1.d-10
+
+select case (matrule)
+case (1) ! Assign material properties of volumetric majority material
+ matel=materialn(sum(maxloc(vol_lsf)))
+ viscosity=mat(matel)%viscosity
+ penal=mat(matel)%penalty
+ expon=mat(matel)%expon
+ plasticity_type=mat(matel)%plasticity_type
+ plasticity_parameters=mat(matel)%plasticity_parameters
+
+case(2) ! Assign material properties of volumetric minority material
+ matel=materialn(sum(minloc(vol_lsf)))
+ viscosity=mat(matel)%viscosity
+ penal=mat(matel)%penalty
+ expon=mat(matel)%expon
+ 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
+ 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
+ matel=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=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)
+
+end subroutine set_elem_props
+
+!-------------------------------------------------------------------------------
+
+subroutine divide_element(params,mat,level,levelmax,nlsf,nnode,icon,ileaves, &
+ materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp,pressure, &
+ strain,eviscosity)
+
+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$')
do kk=1,2
t(1:4)=-1.d0+float(kk-1)
@@ -252,6 +353,7 @@ do kk=1,2
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)
@@ -261,7 +363,7 @@ enddo
deallocate (lsfp)
-return
-end
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
\ No newline at end of file
+end subroutine divide_element
+
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
\ No newline at end of file
--
GitLab