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)

一様グリッド + 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

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~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~sample_photo_raycast_particles sample_photo_raycast_particles proc~sample_photo_raycast_particles->proc~find_first_hit proc~sample_photo_species_state sample_photo_species_state proc~sample_photo_species_state->proc~sample_photo_raycast_particles proc~init_particle_batch_from_config init_particle_batch_from_config proc~init_particle_batch_from_config->proc~sample_photo_species_state

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 &
    )
    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

    real(dp), parameter :: eps = 1.0d-12
    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)
    integer(i32) :: axis, nx, ny, cid
    integer(i32) :: cell(3), step(3)
    integer(i32) :: k, start_idx, end_idx, elem_idx
    logical :: ok, hit_grid

    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 (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
    do axis = 1, 3
      cell(axis) = coord_to_cell(mesh, p_entry(axis), int(axis, kind=i32))
      if (abs(d(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 (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
        do while (t_max(axis) < t_cur - eps)
          t_max(axis) = t_max(axis) + t_delta(axis)
        end do
      end if
    end do

    nx = mesh%grid_ncell(1)
    ny = mesh%grid_ncell(2)

    do
      if (t_cur > t_exit + eps) exit
      if (t_cur > best_t + eps) exit

      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 (t_next > t_exit + eps) exit
      if (t_next > best_t + eps) exit

      do axis = 1, 3
        if (t_max(axis) <= t_next + eps) then
          if (step(axis) /= 0_i32) then
            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
      t_cur = t_next
    end do
  end subroutine find_first_hit_base_grid