ResPart_set_GetSurfaceStresses Subroutine

    impure subroutine ResPart_set_GetSurfaceStresses(this,P,U,V,W,ibVF,visc)
      !> Computes hydrodynamic stresses on markers
      implicit none
      class(ResPart_set),   intent(inout) :: this                              !! Collection of Resolved Particles
      type(Eulerian_obj_r), intent(in)    :: P                                 !! Fluid pressure field
      type(Eulerian_obj_r), intent(in)    :: U                                 !! Fluid velocity field in 1-dir
      type(Eulerian_obj_r), intent(in)    :: V                                 !! Fluid velocity field in 2-dir
      type(Eulerian_obj_r), intent(in)    :: W                                 !! Fluid velocity field in 3-dir
      type(Eulerian_obj_r), intent(in)    :: ibVF                              !! Solid volume fraction
      real(wp),             intent(in)    :: visc                              !! Fluid viscosity
      ! work variable
      type(Eulerian_obj_r) :: sig(6)                                           !! Shear Stress tensor
      integer  :: slo(3), shi(3)                                               !! Stencil bounds
      real(wp) :: pm(3)                                                        !! Resolved pressure at marker location
      real(wp) :: sigm(6)                                                      !! Resolved shear stress at marker location
      real(wp) :: vf                                                           !! Fluid volume fraction at marker location
      integer  :: n,m,dir,i,j,k
      integer  :: sx(3)
      real(wp) :: w1D(2,3)
      real(wp) :: w3D(2,2,2)
      logical  :: is_active

      ! Compute viscous stress tensor
      sig = this%op%StrainRate(U,V,W)

      do n=1,6
       sig(n) = sig(n)*(2.0_wp*visc)
      end do

      ! Get one-sided stresses at each marker
      select type (markers=>this%ib%p)
      type is (marker_obj)
        do m=1,this%ib%count_

          is_active =.true.
          do dir=1,3
            if (markers(m)%p(dir).lt.this%block%pmin(dir).or. &
                markers(m)%p(dir).gt.this%block%pmax(dir)) then
                is_active=.false.
            end if
          end do

          if (is_active) then
            ! Locate marker w/r to mesh
            ! -------------------------- !
            slo = markers(m)%c

            do dir=1,3
              sx(dir) = -1; if (markers(m)%p(dir).ge.this%block%axis(dir)%xm(slo(dir))) sx(dir) = 0
            end do
            slo = slo + sx
            shi = slo + 1

            ! - Build interpolator
            do dir=1,3
              w1D(1,dir) = (this%block%axis(dir)%xm(shi(dir)) - markers(m)%p(dir))/(this%block%axis(dir)%xm(shi(dir))-this%block%axis(dir)%xm(slo(dir)))
              w1D(2,dir) = 1.0_wp - w1D(1,dir)
            end do

            do k=1,2
              do j=1,2
                do i=1,2
                  w3D(i,j,k) = w1D(i,1)*w1D(j,2)*w1D(k,3)
                end do
              end do
            end do
            ! - Interpolate
            ! --- Pressure
            pm       = sum(w3D*P%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))
            ! --- Shear stress tensor
            do n=1,6
              sigm(n) = sum(w3D*Sig(n)%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))
            end do
            ! --- Fluid volume fraction
            vf       = 1.0_wp - sum(w3D*ibVF%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))

            markers(m)%f = - max(0.0_wp,min(1.0_wp,vf))*markers(m)%f               &
                           - Pm*markers(m)%n                                       &
                           + [dot_product(markers(m)%n,[sigm(1),sigm(4),sigm(6)]), &
                              dot_product(markers(m)%n,[sigm(4),sigm(2),sigm(5)]), &
                              dot_product(markers(m)%n,[sigm(6),sigm(5),sigm(3)])]

          else
            markers(m)%f = 0.0_wp
          end if
        end do
      end select

      return
    end subroutine ResPart_set_GetSurfaceStresses