marker_set_GetIBForcing Subroutine

    impure subroutine marker_set_GetIBForcing(this,Um,Vm,Wm,rhof,SDF,ibF,dt)
      !> Computes the IB forcing
      ! Interpolation are carried out by trilinear interpolations.
      implicit none
      class(marker_set),    intent(inout) :: this                              !! A collection of tessellation elements
      type(Eulerian_obj_r), intent(in)    :: Um                                !! Velocity in 1-dir
      type(Eulerian_obj_r), intent(in)    :: Vm                                !! Velocity in 2-dir
      type(Eulerian_obj_r), intent(in)    :: Wm                                !! Velocity in 3-dir
      real(wp),             intent(in)    :: rhof                              !! Fluid density
      type(Eulerian_obj_r), intent(in)    :: SDF                               !! Surface density function
      type(Eulerian_obj_r), intent(inout) :: ibF(3)                            !! IB forcing
      real(wp),             intent(in)    :: dt                                !! Timestep
      ! work variable
      type(eulerian_obj_r):: tmp
      real(wp):: uf(3)
      integer :: shi(3), slo(3), sx(3)                                         !! Stencil extents
      real(wp):: SDFm
      real(wp):: w1D(2,3)
      real(wp):: w3D(2,2,2)
      integer :: m,dir,i,j,k
      logical :: is_active

      select type (markers=>this%p)
      type is (marker_obj)
        do m=1,this%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 centroid 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
            uf(1)    = sum(w3D*Um%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))
            uf(2)    = sum(w3D*Vm%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))
            uf(3)    = sum(w3D*Wm%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))
            SDFm     = sum(w3D*SDF%cell(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)))

            markers(m)%f = (1.0_wp/SDFm)*rhof*(markers(m)%v-uf)/dt

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

      ! Compute forcing on mid points
      call tmp%Initialize('tmp',this%block,this%parallel,0)

      call this%Filter('F1',tmp); ibF(1) = tmp + ibF(1)
      call this%Filter('F2',tmp); ibF(2) = tmp + ibF(2)
      call this%Filter('F3',tmp); ibF(3) = tmp + ibF(3)

      call tmp%Finalize()

      return
    end subroutine marker_set_GetIBForcing