Updated to function similarly to make_cut.f90. Fixed logic for recursive calls

src/pressure_cut.f90

@@ -76,7 +76,6 @@ integer,intent(in) :: matnum
 integer,intent(in) :: leaf_mat_bin(0:params%nmat)
 integer,intent(in) :: ematnump
 
-
 !------------------------------------------------------------------------------|
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
 !------------------------------------------------------------------------------|
@@ -84,12 +83,19 @@ integer,intent(in) :: ematnump
 integer :: matel,matrule,i,j,k,err
 double precision :: prod,vol_lsf0,eps,pressurep,xmean,ymean,zmean,penal
 double precision,allocatable :: vol_lsf(:)
+logical :: cloud_for_mat_props
 
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
 
+
+!write(*,*) 'Hello from pressure_cut'
+
+
+
 eps=1.d-10
 icut=0
+cloud_for_mat_props=.false.
 
 if (level.eq.0) then
   pressure=0.d0
@@ -116,11 +122,24 @@ do i=1,nlsf
      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$')
+
+
+       !write(*,*) 'before get mat vol'
+
+
+
        call get_mat_volume_pc(params,nlsf,lsf,vol_lsf,vol_lsf0)
        ! Added option to use cloud particles for material property assignment
        ! dwhipp/jallen 03/2013
        if (params%materials_on_cloud) then
          if (vol_lsf0.gt.eps .or. level.gt.0) then
+
+
+           !write (*,*) 'vol_lsf0: ',vol_lsf0
+           !write (*,*) 'level: ',level
+
+
+
            matrule=0
          else
            if (params%surface_for_mat_props(i)) then
@@ -131,7 +150,15 @@ do i=1,nlsf
                endif
              endif
            else
+             !write(*,*) 'before get elem mat num'
+
+
              call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
+
+
+             !write(*,*) 'after get elem mat num'
+
+
              cloud_for_mat_props=.true.
            endif
          endif
@@ -159,46 +186,123 @@ do i=1,nlsf
                                    pressure,r0,s0,t0,rst,xmean,ymean,zmean,     &
                                    eviscosity)
            else
-             call divide_element_pc(params,mat,level,levelmax,nlsf,nnode,icon,  &
-                                   ileaves,materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,&
-                                   w,temp,pressure,eviscosity)
+             call divide_element_pc(params,mat,level,levelmax,nlsf,nnode,icon, &
+                                   ileaves,materialn,matnum,leaf_mat_bin,      &
+                                   ematnump,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp, &
+                                   pressure,eviscosity)
            endif
          endif
          deallocate (vol_lsf)
        endif
      endif
    enddo
+
+
+   !write (*,*) 'level before uncut test: ',level
+
+
    if (params%materials_on_cloud .and. level.eq.0) then ! WILL ONLY OCCUR FOR UNCUT ELEMENTS
+
+     !write (*,*) 'Uncut element, right?'
+     !write (*,*) 'lsf ',i,': ',lsf(:,i)
+
+
      if (.not.cloud_for_mat_props .and. icut.eq.0 .and. i.eq.nlsf) then
        ! Assign material properties for given element from cloud
+
+
+       !write(*,*) 'before get elem mat num 2'
+
+
        call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
+
+
+       !write(*,*) 'after get elem mat num 2' 
+
        cloud_for_mat_props=.true.
+
+
      endif
    else
+
+
+     !write (*,*) 'This is a cut element'
+
+
      if (prod.lt.0.d0 .and. icut.eq.0) then
        matel=materialn(i)
+
+
+       !write (*,*) 'Assigning material number matel=',matel
+
+
      endif
    endif
 enddo
 
 ! Check for overriding material transitions
 call mat_trans_pc(params,mat,matel,xmean,ymean,zmean)
-if (matel /= matnum) call stop_run ('matnum disagreement in pressure_cut$')
+
+
+!write (*,*) 'matel: ',matel
+!write (*,*) 'matnum: ',matnum
+
+
+if (level == 0) then
+  if (icut == 0 .or. cloud_for_mat_props) then
+    ! For uncut elements, make sure the assigned material number agrees with osolve%matnum
+    if (matel /= matnum) call stop_run ('matnum disagreement in pressure_cut$')
+  endif
+endif
+
+
+
+!write (*,*) 'matel: ',matel
+!write (*,*) 'matnum: ',matnum
+
+
+
 !matnum=matel
 
 ! If the element is uncut and using divFEM, or the cloud was used for material
 ! property assignment, call make_pressure here
 if (icut.eq.0 .or. cloud_for_mat_props) then
+
+
+  !write (*,*) 'Uncut element'
+
+
+
   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
     penal = params%sstemp_penalty_spinup
   else
+
+
+    !write (*,*) 'Assigning penal'
+
+
     penal = mat(matel)%penalty
+
+
+    !write (*,*) 'penal assigned'
+
+
   endif
+
+
+  !write (*,*) 'before make pressure'
+
+
   call make_pressure(params,icon,x,y,z,penal,u,v,w,temp,pressurep,nnode,r0,  &
                     s0,t0,rst,eviscosity)
 
+
+
+  !write (*,*) 'after make pressure'
+
+
   pressure=pressure+pressurep/(8.d0**level)
 endif
 
@@ -336,8 +440,8 @@ end subroutine set_elem_props_pc
 !-------------------------------------------------------------------------------
 
 subroutine divide_element_pc(params,mat,level,levelmax,nlsf,nnode,icon,ileaves,&
-                            materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp,       &
-                            pressure,eviscosity)
+                            materialn,matnum,leaf_mat_bin,ematnump,lsf,r0,s0,  &
+                            t0,rst,x,y,z,u,v,w,temp,pressure,eviscosity)
 
 use definitions
 
@@ -352,6 +456,9 @@ integer,intent(in) :: nnode
 integer,intent(in) :: icon(params%mpe)
 integer,intent(in) :: ileaves
 integer,intent(in) :: materialn(0:nlsf)
+integer,intent(in) :: matnum
+integer,intent(in) :: leaf_mat_bin(0:params%nmat)
+integer,intent(in) :: ematnump
 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)
@@ -415,7 +522,7 @@ do kk=1,2
 
          call pressure_cut (params,levelp,levelmax,icon,x,y,z,mat,materialn,u, &
                            v,w,temp,pressure,nnode,lsfp,nlsf,r0p,s0p,t0p,rstp, &
-                           jcut,ileaves,eviscosity)
+                           jcut,ileaves,eviscosity,matnum,leaf_mat_bin,ematnump)
       enddo
    enddo
 enddo