find_first_hit_base_grid Subroutine

public subroutine find_first_hit_base_grid(mesh, p0, p1, d, seg_min, seg_max, hit, best_t, use_box_filter, box_min, box_max, box_tol, require_elem_inside, status)

一様グリッド + 3D-DDA で候補セルのみ探索し、最初の命中要素を返す。

Arguments

Type IntentOptional Attributes Name
type(mesh_type), intent(in) :: mesh
real(kind=dp), intent(in) :: p0(3)
real(kind=dp), intent(in) :: p1(3)
real(kind=dp), intent(in) :: d(3)
real(kind=dp), intent(in) :: seg_min(3)
real(kind=dp), intent(in) :: seg_max(3)
type(hit_info), intent(inout) :: hit
real(kind=dp), intent(inout) :: best_t
logical, intent(in) :: use_box_filter
real(kind=dp), intent(in) :: box_min(3)
real(kind=dp), intent(in) :: box_max(3)
real(kind=dp), intent(in) :: box_tol
logical, intent(in) :: require_elem_inside
integer(kind=i32), intent(out) :: status

Calls

proc~~find_first_hit_base_grid~~CallsGraph proc~find_first_hit_base_grid find_first_hit_base_grid proc~bbox_inside_box bbox_inside_box proc~find_first_hit_base_grid->proc~bbox_inside_box proc~cell_id cell_id proc~find_first_hit_base_grid->proc~cell_id proc~coord_to_cell coord_to_cell proc~find_first_hit_base_grid->proc~coord_to_cell proc~lower_ascii lower_ascii proc~find_first_hit_base_grid->proc~lower_ascii proc~point_inside_box point_inside_box proc~find_first_hit_base_grid->proc~point_inside_box proc~segment_aabb_intersection_t segment_aabb_intersection_t proc~find_first_hit_base_grid->proc~segment_aabb_intersection_t proc~segment_bbox_overlap_precomputed segment_bbox_overlap_precomputed proc~find_first_hit_base_grid->proc~segment_bbox_overlap_precomputed proc~segment_triangle_intersect segment_triangle_intersect proc~find_first_hit_base_grid->proc~segment_triangle_intersect proc~cross cross proc~segment_triangle_intersect->proc~cross

Called by

proc~~find_first_hit_base_grid~~CalledByGraph proc~find_first_hit_base_grid find_first_hit_base_grid proc~find_first_hit_base find_first_hit_base proc~find_first_hit_base->proc~find_first_hit_base_grid proc~find_first_hit find_first_hit proc~find_first_hit->proc~find_first_hit_base proc~find_first_hit_periodic2 find_first_hit_periodic2 proc~find_first_hit->proc~find_first_hit_periodic2 proc~find_first_hit_periodic2->proc~find_first_hit_base proc~advance_particle_step advance_particle_step proc~advance_particle_step->proc~find_first_hit proc~advance_particle_step_upper_panel_fourier advance_particle_step_upper_panel_fourier proc~advance_particle_step_upper_panel_fourier->proc~find_first_hit proc~resolve_particle_boundary_candidate resolve_particle_boundary_candidate proc~resolve_particle_boundary_candidate->proc~find_first_hit proc~sample_photo_raycast_particles sample_photo_raycast_particles proc~sample_photo_raycast_particles->proc~find_first_hit

Source Code

  subroutine find_first_hit_base_grid( &
    mesh, p0, p1, d, seg_min, seg_max, hit, best_t, use_box_filter, box_min, box_max, box_tol, require_elem_inside, status &
    )
    type(mesh_type), intent(in) :: mesh
    real(dp), intent(in) :: p0(3), p1(3), d(3), seg_min(3), seg_max(3)
    type(hit_info), intent(inout) :: hit
    real(dp), intent(inout) :: best_t
    logical, intent(in) :: use_box_filter
    real(dp), intent(in) :: box_min(3), box_max(3), box_tol
    logical, intent(in) :: require_elem_inside
    integer(i32), intent(out) :: status

    real(dp), parameter :: t_eps = 1.0d-12
    real(dp), parameter :: axis_rel_eps = 64.0d0*epsilon(1.0d0)
    real(dp) :: t_entry, t_exit, t_cur, t_next
    real(dp) :: t_max(3), t_delta(3), cell_size
    real(dp) :: t, h(3), p_entry(3)
    real(dp) :: axis_eps
    integer(i32) :: axis, nx, ny, cid
    integer(i32) :: cell(3), step(3)
    integer(i32) :: k, start_idx, end_idx, elem_idx
    integer(i64) :: align_iterations, traversal_iterations, max_traversal_iterations
    logical :: ok, hit_grid, advanced_cell

    status = collision_query_ok
    if (any(mesh%grid_ncell <= 0_i32) .or. &
        any(.not. ieee_is_finite(mesh%grid_inv_cell)) .or. any(mesh%grid_inv_cell <= 0.0d0) .or. &
        any(.not. ieee_is_finite(mesh%grid_bb_min)) .or. any(.not. ieee_is_finite(mesh%grid_bb_max))) then
      call report_collision_grid_stall('invalid_grid_geometry', mesh, p0, p1)
      status = collision_query_grid_stalled
      return
    end if

    call segment_aabb_intersection_t(p0, d, mesh%grid_bb_min, mesh%grid_bb_max, hit_grid, t_entry, t_exit)
    if (.not. hit_grid) return
    if (.not. ieee_is_finite(t_entry) .or. .not. ieee_is_finite(t_exit)) then
      call report_collision_grid_stall( &
        'nonfinite_aabb_interval', mesh, p0, p1, t_entry=t_entry, t_exit=t_exit &
        )
      status = collision_query_grid_stalled
      return
    end if
    if (t_exit < 0.0d0 .or. t_entry > 1.0d0) return

    t_cur = max(0.0d0, t_entry)
    if (t_cur > t_exit) return

    p_entry = p0 + t_cur*d
    axis_eps = axis_rel_eps*max( &
               maxval(abs(d)), &
               maxval(abs(mesh%grid_bb_max - mesh%grid_bb_min)), &
               tiny(1.0d0) &
               )
    cell = 0_i32
    step = 0_i32
    t_max = huge(1.0d0)
    t_delta = huge(1.0d0)
    do axis = 1, 3
      cell(axis) = coord_to_cell(mesh, p_entry(axis), int(axis, kind=i32))
      if (abs(d(axis)) <= axis_eps) then
        step(axis) = 0_i32
        t_max(axis) = huge(1.0d0)
        t_delta(axis) = huge(1.0d0)
      else
        cell_size = 1.0d0/mesh%grid_inv_cell(axis)
        if (.not. ieee_is_finite(cell_size) .or. cell_size <= 0.0d0) then
          call report_collision_grid_stall( &
            'invalid_cell_size', mesh, p0, p1, cell=cell, step=step, axis=axis, &
            t_entry=t_entry, t_exit=t_exit, t_cur=t_cur &
            )
          status = collision_query_grid_stalled
          return
        end if
        if (d(axis) > 0.0d0) then
          step(axis) = 1_i32
          t_max(axis) = (mesh%grid_bb_min(axis) + real(cell(axis), dp)*cell_size - p0(axis))/d(axis)
          t_delta(axis) = cell_size/d(axis)
        else
          step(axis) = -1_i32
          t_max(axis) = (mesh%grid_bb_min(axis) + real(cell(axis) - 1_i32, dp)*cell_size - p0(axis))/d(axis)
          t_delta(axis) = -cell_size/d(axis)
        end if
        if (.not. ieee_is_finite(t_max(axis)) .or. .not. ieee_is_finite(t_delta(axis)) .or. &
            t_delta(axis) <= 0.0d0) then
          call report_collision_grid_stall( &
            'invalid_axis_progress', mesh, p0, p1, cell=cell, step=step, axis=axis, &
            t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, t_max=t_max, t_delta=t_delta &
            )
          status = collision_query_grid_stalled
          return
        end if
        align_iterations = 0_i64
        do while (t_max(axis) < t_cur - t_eps)
          align_iterations = align_iterations + 1_i64
          if (align_iterations > int(mesh%grid_ncell(axis), i64) + 1_i64) then
            call report_collision_grid_stall( &
              'align_iteration_limit', mesh, p0, p1, cell=cell, step=step, axis=axis, &
              iterations=align_iterations, max_iterations=int(mesh%grid_ncell(axis), i64) + 1_i64, &
              t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, t_max=t_max, t_delta=t_delta &
              )
            status = collision_query_grid_stalled
            return
          end if
          t_max(axis) = t_max(axis) + t_delta(axis)
          if (.not. ieee_is_finite(t_max(axis))) then
            call report_collision_grid_stall( &
              'nonfinite_t_max_after_align', mesh, p0, p1, cell=cell, step=step, axis=axis, &
              iterations=align_iterations, t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, &
              t_max=t_max, t_delta=t_delta &
              )
            status = collision_query_grid_stalled
            return
          end if
        end do
      end if
    end do

    nx = mesh%grid_ncell(1)
    ny = mesh%grid_ncell(2)
    traversal_iterations = 0_i64
    max_traversal_iterations = sum(int(mesh%grid_ncell, i64)) + 3_i64

    do
      traversal_iterations = traversal_iterations + 1_i64
      if (traversal_iterations > max_traversal_iterations) then
        call report_collision_grid_stall( &
          'traversal_iteration_limit', mesh, p0, p1, cell=cell, step=step, &
          iterations=traversal_iterations, max_iterations=max_traversal_iterations, &
          t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, t_max=t_max, t_delta=t_delta &
          )
        status = collision_query_grid_stalled
        return
      end if
      if (.not. ieee_is_finite(t_cur)) then
        call report_collision_grid_stall( &
          'nonfinite_t_cur', mesh, p0, p1, cell=cell, step=step, &
          iterations=traversal_iterations, t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, &
          t_max=t_max, t_delta=t_delta &
          )
        status = collision_query_grid_stalled
        return
      end if
      if (t_cur > t_exit + t_eps) exit
      if (t_cur > best_t + t_eps) exit

      if (any(cell < 1_i32) .or. any(cell > mesh%grid_ncell)) then
        call report_collision_grid_stall( &
          'cell_out_of_range', mesh, p0, p1, cell=cell, step=step, &
          iterations=traversal_iterations, t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, &
          t_max=t_max, t_delta=t_delta &
          )
        status = collision_query_grid_stalled
        return
      end if

      cid = cell_id(cell(1), cell(2), cell(3), nx, ny)
      start_idx = mesh%grid_cell_start(cid)
      end_idx = mesh%grid_cell_start(cid + 1_i32) - 1_i32
      do k = start_idx, end_idx
        elem_idx = mesh%grid_cell_elem(k)
        if (use_box_filter) then
          if (.not. segment_bbox_overlap_precomputed( &
              mesh%bb_min(:, elem_idx), mesh%bb_max(:, elem_idx), box_min, box_max)) cycle
          if (require_elem_inside) then
            if (.not. bbox_inside_box( &
                mesh%bb_min(:, elem_idx), mesh%bb_max(:, elem_idx), box_min, box_max, box_tol)) cycle
          end if
        end if
        if (.not. segment_bbox_overlap_precomputed( &
            seg_min, seg_max, mesh%bb_min(:, elem_idx), mesh%bb_max(:, elem_idx))) cycle
        call segment_triangle_intersect( &
          p0, p1, mesh%v0(:, elem_idx), mesh%v1(:, elem_idx), mesh%v2(:, elem_idx), ok, t, h &
          )
        if (.not. ok) cycle
        if (use_box_filter) then
          if (.not. point_inside_box(h, box_min, box_max, box_tol)) cycle
        end if
        if (t < best_t) then
          best_t = t
          hit%has_hit = .true.
          hit%elem_idx = elem_idx
          hit%t = t
          hit%pos = h
          hit%pos_wrapped = h
          hit%image_shift = 0_i32
        end if
      end do

      t_next = min(t_max(1), min(t_max(2), t_max(3)))
      if (.not. ieee_is_finite(t_next) .or. t_next < t_cur - t_eps) then
        call report_collision_grid_stall( &
          'invalid_t_next', mesh, p0, p1, cell=cell, step=step, &
          iterations=traversal_iterations, t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, t_next=t_next, &
          t_max=t_max, t_delta=t_delta &
          )
        status = collision_query_grid_stalled
        return
      end if
      if (t_next > t_exit + t_eps) exit
      if (t_next > best_t + t_eps) exit

      advanced_cell = .false.
      do axis = 1, 3
        if (t_max(axis) <= t_next + t_eps) then
          if (step(axis) /= 0_i32) then
            advanced_cell = .true.
            cell(axis) = cell(axis) + step(axis)
            if (cell(axis) < 1_i32 .or. cell(axis) > mesh%grid_ncell(axis)) return
            t_max(axis) = t_max(axis) + t_delta(axis)
          end if
        end if
      end do
      if (.not. advanced_cell) then
        call report_collision_grid_stall( &
          'no_cell_advanced', mesh, p0, p1, cell=cell, step=step, &
          iterations=traversal_iterations, t_entry=t_entry, t_exit=t_exit, t_cur=t_cur, t_next=t_next, &
          t_max=t_max, t_delta=t_delta &
          )
        status = collision_query_grid_stalled
        return
      end if
      t_cur = t_next
    end do
  end subroutine find_first_hit_base_grid