resolve_panel_surface_sides Subroutine

  subroutine resolve_panel_surface_sides(mesh, policy, status, message, mesh_id)
    type(mesh_type), intent(inout) :: mesh
    character(len=*), intent(in) :: policy
    integer(i32), intent(out) :: status
    character(len=*), intent(out) :: message
    integer(i32), intent(in), optional :: mesh_id
    character(len=32) :: normalized
    integer(i32), allocatable :: signs(:), neighbor(:, :), queue(:)
    logical, allocatable :: visited(:), active(:)
    integer(i32) :: i, edge, face, next_face, head, tail, component_sign
    real(dp) :: signed_volume, volume_scale

    status = panel_surface_side_ok
    message = ''
    normalized = lower_ascii(trim(policy))
    allocate (signs(mesh%nelem), active(mesh%nelem))
    active = .true.
    if (present(mesh_id)) active = mesh%elem_mesh_id == mesh_id
    if (.not. any(active)) then
      status = panel_surface_side_invalid_policy
      message = 'surface side mesh_id selects no elements'
      return
    end if
    signs = mesh%elem_vacuum_sign
    select case (trim(normalized))
    case ('normal_plus')
      where (active) signs = 1_i32
    case ('normal_minus')
      where (active) signs = -1_i32
    case ('outward_closed')
      allocate (neighbor(3, mesh%nelem), queue(mesh%nelem), visited(mesh%nelem))
      call build_closed_neighbors(mesh, active, neighbor, status, message)
      if (status /= panel_surface_side_ok) return
      where (active) signs = 0_i32
      visited = .false.
      do i = 1, mesh%nelem
        if (.not. active(i)) cycle
        if (visited(i)) cycle
        head = 1_i32
        tail = 1_i32
        queue(1) = i
        visited(i) = .true.
        signed_volume = 0.0_dp
        volume_scale = 0.0_dp
        do while (head <= tail)
          face = queue(head)
          head = head + 1_i32
          signed_volume = signed_volume + dot_product( &
                          mesh%v0(:, face), cross_product(mesh%v1(:, face), mesh%v2(:, face)) &
                          )/6.0_dp
          volume_scale = volume_scale + mesh%panel_area(face)*max(1.0_dp, norm2(mesh%centers(:, face)))
          do edge = 1, 3
            next_face = neighbor(edge, face)
            if (.not. visited(next_face)) then
              tail = tail + 1_i32
              queue(tail) = next_face
              visited(next_face) = .true.
            end if
          end do
        end do
        if (abs(signed_volume) <= 128.0_dp*epsilon(1.0_dp)*max(1.0_dp, volume_scale)) then
          status = panel_surface_side_inconsistent
          message = 'closed component has zero or indeterminate signed volume'
          return
        end if
        component_sign = merge(1_i32, -1_i32, signed_volume > 0.0_dp)
        do face = 1, tail
          signs(queue(face)) = component_sign
        end do
      end do
    case default
      status = panel_surface_side_invalid_policy
      message = 'surface side policy must be normal_plus, normal_minus, or outward_closed'
      return
    end select

    mesh%elem_vacuum_sign = signs
    do i = 1, mesh%nelem
      if (signs(i) /= 0_i32) mesh%vacuum_normals(:, i) = real(signs(i), dp)*mesh%normals(:, i)
    end do
  end subroutine resolve_panel_surface_sides