diff --git a/CMakeLists.txt b/CMakeLists.txt index fd00fcc8..e2304dcf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -209,16 +209,16 @@ if(SIMPLE_ENABLE_CGAL) include(FetchContent) FetchContent_Declare( cgal - URL https://github.com/CGAL/cgal/releases/download/v5.6.1/CGAL-5.6.1.tar.xz + URL https://github.com/CGAL/cgal/releases/download/v6.0.1/CGAL-6.0.1.tar.xz DOWNLOAD_EXTRACT_TIMESTAMP TRUE ) FetchContent_GetProperties(cgal) if(NOT cgal_POPULATED) FetchContent_Populate(cgal) endif() - # CGAL 5.x is header-only, just add include path + # CGAL 6.x is header-only, just add include path set(CGAL_INCLUDE_DIR "${cgal_SOURCE_DIR}/include") - message(STATUS "CGAL enabled via FetchContent: 5.6.1 (header-only)") + message(STATUS "CGAL enabled via FetchContent: 6.0.1 (header-only)") message(STATUS "CGAL include dir: ${CGAL_INCLUDE_DIR}") add_compile_definitions(SIMPLE_ENABLE_CGAL) else() diff --git a/src/netcdf_results_output.f90 b/src/netcdf_results_output.f90 index 32d128ad..b5a8e74d 100644 --- a/src/netcdf_results_output.f90 +++ b/src/netcdf_results_output.f90 @@ -366,13 +366,14 @@ end subroutine write_results_netcdf subroutine compute_cartesian_positions(xstart_cart, xend_cart) !> Convert reference coordinates to Cartesian for all particles. !> - !> For particles that were not traced (zend = 0), xend_cart is set - !> to xstart_cart since they effectively remained at the start. + !> For particles that were not traced (zend = 0), xend_cart is set to zero. !> - !> TODO: Investigate batch coordinate transformation for performance. - !> Current implementation uses O(n) individual calls. + !> For STL wall hits (wall_hit == 1), xend_cart is set directly from + !> wall_hit_cart since that is the authoritative CGAL intersection point. + !> This avoids round-trip coordinate conversion errors when from_cart + !> fails in ill-conditioned chartmap regions. - use params, only: ntestpart, zstart, zend + use params, only: ntestpart, zstart, zend, wall_hit, wall_hit_cart use reference_coordinates, only: ref_coords real(dp), intent(out) :: xstart_cart(:, :) @@ -389,14 +390,19 @@ subroutine compute_cartesian_positions(xstart_cart, xend_cart) do i = 1, ntestpart call ref_coords%evaluate_cart(zstart(1:3, i), xstart_cart(:, i)) - ! Check if particle was traced: zend is set to exactly 0 for untraced - ! particles (see simple_main.f90:493 and classification.f90:94) - zend_is_zero = all(zend(1:3, i) == 0.0_dp) - if (zend_is_zero) then - ! Particle not traced - use start position - xend_cart(:, i) = xstart_cart(:, i) + ! For STL wall hits, use the authoritative CGAL intersection point + if (wall_hit(i) == 1) then + xend_cart(:, i) = wall_hit_cart(:, i) else - call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i)) + ! Check if particle was traced: zend is set to exactly 0 for untraced + ! particles (see simple_main.f90:493 and classification.f90:94) + zend_is_zero = all(zend(1:3, i) == 0.0_dp) + if (zend_is_zero) then + ! Particle not traced - keep xend_cart as zero + xend_cart(:, i) = 0.0_dp + else + call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i)) + end if end if end do diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 12990946..908405a8 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -542,6 +542,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) real(dp) :: x_prev(3), x_cur(3) real(dp) :: x_prev_m(3), x_cur_m(3), x_hit_m(3), x_hit(3) real(dp) :: normal_m(3), vhat(3), vnorm, cos_inc + real(dp) :: segment_length, hit_distance, t_frac integer :: it, ierr_orbit, it_final integer(8) :: kt logical :: passing @@ -626,6 +627,16 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) if (ierr_from_cart == 0) then call ref_to_integ(u_ref_hit, z(1:3)) + else + ! Fallback: linear interpolation of reference coordinates + ! when from_cart fails (ill-conditioned regions of chartmap) + segment_length = sqrt(sum((x_cur_m - x_prev_m)**2)) + if (segment_length > 0.0_dp) then + hit_distance = sqrt(sum((x_hit_m - x_prev_m)**2)) + t_frac = hit_distance / segment_length + u_ref_hit = u_ref_prev + t_frac * (u_ref_cur - u_ref_prev) + call ref_to_integ(u_ref_hit, z(1:3)) + end if end if ierr_orbit = 77