crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Fixed error in mesh_intersects. Closes #75 on github.

Browse source
This commit is contained in:
Paul Romano 2012-03-25 20:55:14 -04:00
parent 93bfdcccdd
commit 5524fde7f0
1 changed files with 36 additions and 28 deletions
Download patch file Download diff file Expand all files Collapse all files

64
src/mesh.F90
View file

@ -218,7 +218,9 @@ contains
end subroutine count_fission_sites
!===============================================================================
! MESH_INTERSECT
! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer
! boundary of the given mesh. This is important for determining whether a track
! will score to a mesh tally.
!===============================================================================
function mesh_intersects(m, xyz0, xyz1) result(intersects)
@ -228,9 +230,11 @@ contains
real(8), intent(in) :: xyz1(3)
logical :: intersects
real(8) :: x0, y0, z0 ! track start point
real(8) :: x1, y1, z1 ! track end point
real(8) :: xi, yi, zi ! track intersection point with mesh
real(8) :: x0, y0, z0 ! track start point
real(8) :: x1, y1, z1 ! track end point
real(8) :: xi, yi, zi ! track intersection point with mesh
real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh
real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
@ -242,62 +246,66 @@ contains
y1 = xyz1(2)
z1 = xyz1(3)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
ym0 = m % lower_left(2)
zm0 = m % lower_left(3)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
ym1 = m % upper_right(2)
zm1 = m % upper_right(3)
! Check if line intersects bottom surface -- calculate the intersection
! point (y,z)
xi = m % lower_left(1)
yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0)
if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects left surface -- calculate the intersection point
! (x,z)
yi = m % lower_left(2)
xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0)
if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects front surface -- calculate the intersection point
! (x,y)
zi = m % lower_left(3)
xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0)
if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
! Check if line intersects top surface -- calculate the intersection
! point (y,z)
xi = m % upper_right(1)
yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0)
if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects right surface -- calculate the intersection point
! (x,z)
yi = m % upper_right(2)
xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0)
if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
! Check if line intersects back surface -- calculate the intersection point
! (x,y)
zi = m % upper_right(3)
xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0)
if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if