marker_set_AddPlane Subroutine

    impure subroutine marker_set_AddPlane(this,center,normal,width,vel,dl,tag)
      !> Adds an IB plane.
      class(marker_set),   intent(inout) :: this                               !! A collection of tessellation elements
      real(wp),            intent(in)    :: center(3)                          !! Plane center
      character(len=*),    intent(in)    :: normal                             !! Plane normal
      real(wp),            intent(in)    :: width(3)                           !! Plane extents
      real(wp),            intent(in)    :: vel(3)                             !! Plane velocity
      real(wp),            intent(in)    :: dl                                 !! Element size
      integer(kind=8),     intent(in),    &
                                optional :: tag                                !! Tag
      ! Work variables
      integer          :: ns(3)
      real(wp)         :: ds(3)
      integer          :: i,j,k,m
      real(wp)         :: vec(3)
      real(wp)         :: dA
      integer          :: tag_

      if (present(tag)) then
        tag_ = int(tag,kind=4)
      else
        tag_ = 0
      end if

      ! Count the number of markers
      ns(1) = 1 + floor(width(1)/dl)
      ds(1) = width(1)/real(ns(1),wp)

      ns(2) = 1 + floor(width(2)/dl)
      ds(2) = width(2)/real(ns(2),wp)

      ns(3) = 1 + floor(width(3)/dl)
      ds(3) = width(3)/real(ns(3),wp)

      call this%Resize(this%count_ + ns(1)*ns(2)*ns(3))

      vec = 0.0_wp
      select case (trim(adjustl(normal)))
      case ('+x1','+X1','+x','+X')
        vec(1) = 1.0_wp
        dA     = ds(2)*ds(3)
      case ('-x1','-X1','-x','-X')
        vec(1) = -1.0_wp
        dA     = ds(2)*ds(3)
      case ('+x2','+X2','+y','+Y')
        vec(2) = 1.0_wp
        dA     = ds(1)*ds(3)
      case ('-x2','-X2','-y','-Y')
        vec(2) = -1.0_wp
        dA     = ds(1)*ds(3)
      case ('+x3','+X3','+z','+Z')
        vec(3) = 1.0_wp
        dA     = ds(1)*ds(2)
      case ('-x3','-X3','-z','-Z')
        vec(3) = -1.0_wp
        dA     = ds(1)*ds(2)
      case default
        call this%parallel%Stop("Invalid plane normal")
      end select

      if (width(1)*width(2)*width(3).ne.0.0_wp) then
        call this%parallel%Stop("Invalid plane: must have zero width")
      end if

      select type (markers=>this%p)
      type is (marker_obj)
        ! Add box walls
        m = this%count_ - ns(1)*ns(2)*ns(3)
        do k=1,ns(3)
          do j=1,ns(2)
            do i=1,ns(1)
              m = m + 1

              ! Facet ID
              markers(m)%id = int(m,kind=8)

              ! Position
              markers(m)%p  = center - 0.5_wp*width &
                            + [(real(i,wp)+0.5_wp)*ds(1), (real(j,wp)+0.5_wp)*ds(2),(real(k,wp)+0.5_wp)*ds(3)]

              ! Surface area
              markers(m)%SA = dA

              ! Normal vector
              markers(m)%n  = vec

              ! Velocity
              markers(m)%v  = vel

              ! Tag
              markers(m)%s  = tag_

              ! Force
              markers(m)%f  = 0.0_wp
            end do
          end do
        end do

      end select

      return
    end subroutine marker_set_AddPlane