Updated to do a proper majority rule for full elements. Similar to the make_cut.f90 update

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
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
+end subroutine divide_element
+
+!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
\ No newline at end of file