marker_set_AddCylinder Subroutine

    pure subroutine marker_set_AddCylinder(this,base,L,radius,vel,dl,tag)
      !> Adds an IB cylinder (with open faces).
      class(marker_set),   intent(inout) :: this                               !! A collection of tessellation elements
      real(wp),            intent(in)    :: base(3)                            !! Base point for extrusion
      real(wp),            intent(in)    :: L(3)                               !! Extrusion length
      real(wp),            intent(in)    :: radius                             !! Cylinder radius
      real(wp),            intent(in)    :: vel(3)                             !! Cylinder translational velocity
      real(wp),            intent(in)    :: dl                                 !! Element size
      integer(kind=8),     intent(in),    &
                                optional :: tag                                !! Tag
      ! Work variables
      integer         :: n_axis
      real(wp)        :: d_axis
      integer         :: n_theta
      real(wp)        :: d_theta
      real(wp)        :: theta
      real(wp)        :: u1(3), u2(3)
      integer         :: i,j
      integer(kind=8) :: id
      real(wp),        &
            parameter :: Pi = 4.0_wp*atan(1.0_wp)
      integer         :: tag_

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

      ! Discretize axis
      n_axis = max(1, floor(norm2(L)/dl))
      d_axis = norm2(L)/real(n_axis,wp)

      ! Discretize the theta angle
      n_theta = ceiling(Pi*2.0_wp*radius/dl)
      d_theta = 2.0_wp*Pi/real(n_theta,wp)

      call this%Resize(this%count_ + n_theta*n_axis)

      ! Chose a vector that's not parallel to the axis
      select case (maxloc([abs(dot_product([1,0,0], L/norm2(L))), &
                           abs(dot_product([0,1,0], L/norm2(L))), &
                           abs(dot_product([0,0,1], L/norm2(L)))],dim=1))
      case (1)
        u1 = [1,0,0]
      case (2)
        u1 = [0,1,0]
      case (3)
        u1 = [0,0,1]
      end select

      ! Build two vectors orthogonal to the axis
      u1 = cross_product(L,u1)
      u1 = u1/norm2(u1)

      u2 = cross_product(L,u1)
      u2 = u2/norm2(u2)

      ! Place markers at center of facets
      id = this%count_ - n_theta*n_axis
      select type(markers =>this%p)
      type is (marker_obj)
        do j=1,n_axis
          do i=1,n_theta

            ! Polar angle
            theta = (i-1)*d_theta + d_theta/2.0_wp

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

            ! Position
            markers(id)%p  = base + real(j-1,wp)*d_axis*L/norm2(L) &
                                  + radius*cos(theta)*u1           &
                                  + radius*sin(theta)*u2

            ! Surface area
            markers(id)%SA = radius*d_theta*d_axis

            ! Normal vector
            markers(id)%n  = cos(theta)*u1+sin(theta)*u2

            ! Velocity
            markers(id)%v  = vel

            ! Tag
            markers(id)%s  = tag_

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

      end select

      return
    end subroutine marker_set_AddCylinder