make_plate_hole Subroutine

  subroutine make_plate_hole(mesh, size_x, size_y, radius, n_theta, n_r, center)
    type(mesh_type), intent(out) :: mesh
    real(dp), intent(in), optional :: size_x, size_y, radius, center(3)
    integer(i32), intent(in), optional :: n_theta, n_r
    real(dp) :: sx, sy, r_hole, c(3), pi, angle, f0, f1
    real(dp) :: x_min, x_max, y_min, y_max, min_margin
    real(dp) :: p00(3), p01(3), p10(3), p11(3), corners(3, 3)
    integer(i32) :: nt, nr0, it, ir, itri, nelem, n_corner, k
    integer(i32), allocatable :: edge_ids(:)
    real(dp), allocatable :: hole_pts(:, :), outer_pts(:, :)
    real(dp), allocatable :: v0(:, :), v1(:, :), v2(:, :)

    sx = 1.0d0; sy = 1.0d0; r_hole = 0.2d0; nt = 48; nr0 = 4; c = 0.0d0
    if (present(size_x)) sx = size_x
    if (present(size_y)) sy = size_y
    if (present(radius)) r_hole = radius
    if (present(n_theta)) nt = n_theta
    if (present(n_r)) nr0 = n_r
    if (present(center)) c = center
    if (sx <= 0.0d0 .or. sy <= 0.0d0) error stop "size_x and size_y must be positive"
    if (nt < 3 .or. nr0 <= 0) error stop "n_theta >= 3 and n_r > 0"

    x_min = c(1) - 0.5d0*sx
    x_max = c(1) + 0.5d0*sx
    y_min = c(2) - 0.5d0*sy
    y_max = c(2) + 0.5d0*sy
    min_margin = min(min(c(1) - x_min, x_max - c(1)), min(c(2) - y_min, y_max - c(2)))
    if (r_hole <= 0.0d0) error stop "radius must be > 0 for plate_hole"
    if (r_hole >= min_margin) error stop "plate_hole radius must be smaller than half-width/half-height"

    allocate (hole_pts(3, nt + 1), outer_pts(3, nt + 1), edge_ids(nt + 1))
    pi = acos(-1.0d0)
    do it = 0, nt
      angle = 2.0d0*pi*real(it, dp)/real(nt, dp)
      hole_pts(:, it + 1) = [c(1) + r_hole*cos(angle), c(2) + r_hole*sin(angle), c(3)]
      call ray_to_rectangle( &
        c(1), c(2), cos(angle), sin(angle), x_min, x_max, y_min, y_max, c(3), outer_pts(:, it + 1), edge_ids(it + 1) &
        )
    end do

    nelem = 2*nt*nr0
    do it = 1, nt
      nelem = nelem + transition_corner_count(edge_ids(it), edge_ids(it + 1))
    end do
    allocate (v0(3, nelem), v1(3, nelem), v2(3, nelem))

    itri = 0
    do it = 1, nt
      do ir = 0, nr0 - 1
        f0 = real(ir, dp)/real(nr0, dp)
        f1 = real(ir + 1, dp)/real(nr0, dp)
        p00 = hole_pts(:, it) + f0*(outer_pts(:, it) - hole_pts(:, it))
        p10 = hole_pts(:, it) + f1*(outer_pts(:, it) - hole_pts(:, it))
        p01 = hole_pts(:, it + 1) + f0*(outer_pts(:, it + 1) - hole_pts(:, it + 1))
        p11 = hole_pts(:, it + 1) + f1*(outer_pts(:, it + 1) - hole_pts(:, it + 1))
        call push_tri(v0, v1, v2, itri, p00, p10, p11)
        call push_tri(v0, v1, v2, itri, p00, p11, p01)
      end do

      call transition_corners( &
        edge_ids(it), edge_ids(it + 1), x_min, x_max, y_min, y_max, c(3), corners, n_corner &
        )
      if (n_corner > 0) then
        if (n_corner == 1) then
          call push_tri(v0, v1, v2, itri, outer_pts(:, it), corners(:, 1), outer_pts(:, it + 1))
        else
          call push_tri(v0, v1, v2, itri, outer_pts(:, it), corners(:, 1), corners(:, 2))
          do k = 2, n_corner - 1
            call push_tri(v0, v1, v2, itri, outer_pts(:, it), corners(:, k), corners(:, k + 1))
          end do
          call push_tri(v0, v1, v2, itri, outer_pts(:, it), corners(:, n_corner), outer_pts(:, it + 1))
        end if
      end if
    end do

    if (itri /= nelem) error stop "plate_hole triangle bookkeeping mismatch"
    call init_mesh(mesh, v0, v1, v2)
  end subroutine make_plate_hole