ResPart_set_GetCentroidVF Function

    impure function ResPart_set_GetCentroidVF(this,ibVF) result (val)
      !> Computes average solid volume fraction at particle centroids
      implicit none
      class(ResPart_set),  intent(inout) :: this                               !! Collection of Resolved Particles
      type(Eulerian_obj_r),intent(in)    :: ibVF                               !! Solid volume fraction
      real(wp)                           :: val                                !! Mean solid volume fraction at particle centroids
      ! work variable
      real(wp) :: buff

      integer  :: slo(3), shi(3)
      integer  :: n,dir,i,j,k
      integer  :: sx(3)
      real(wp) :: w1D(2,3)
      real(wp) :: w3D(2,2,2)


      ! Compute hydrodynamic forces exerted on particles
      val = 0.0_wp
      select type (particles =>this%p)
      type is (ResPart_obj)
        do n=1,this%count_
          ! Locate centroid w/r to mesh
          ! -------------------------- !
          slo = particles(n)%c

          do dir=1,3
            sx(dir) = -1; if (particles(n)%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)) - particles(n)%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

          val = val + sum(w3D*ibVF%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))
        end do
      end select

      call this%parallel%Sum(val,buff); val = buff/this%count

      return
    end function ResPart_set_GetCentroidVF