diff --git a/src/geometry.F90 b/src/geometry.F90 index bc1eca2ae9..8a449a284b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -35,6 +35,15 @@ module geometry logical(C_BOOL), value :: coincident; real(C_DOUBLE) :: d; end function surface_distance_c + + subroutine surface_normal_c(surf_ind, xyz, uvw) & + bind(C, name='surface_normal') + use ISO_C_BINDING + implicit none + integer(C_INT), value :: surf_ind; + real(C_DOUBLE) :: xyz(3); + real(C_DOUBLE) :: uvw(3); + end subroutine surface_normal_c end interface contains @@ -76,6 +85,7 @@ contains integer :: token logical :: actual_sense ! sense of particle wrt surface logical :: alt_sense + real(8) :: uvw1(3), uvw2(3) in_cell = .true. do i = 1, size(c%rpn) @@ -102,6 +112,14 @@ contains write(*, *) actual_sense, alt_sense call fatal_error("") end if + uvw1 = surfaces(abs(token))%obj%normal(p%coord(p%n_coord)%xyz) + call surface_normal_c(abs(token)-1, p%coord(p%n_coord)%xyz, uvw2) + if (.not. all(uvw1 - uvw2 == ZERO)) then + write(*, *) "Surface normals don't match" + write(*, *) uvw1 + write(*, *) uvw2 + call fatal_error("") + end if if (actual_sense .neqv. (token > 0)) then in_cell = .false. exit diff --git a/src/surface_header.C b/src/surface_header.C index 3ac3d3c532..43209fc803 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -903,3 +903,8 @@ extern "C" double surface_distance(int surf_ind, double xyz[3], double uvw[3], bool coincident) { return surfaces_c[surf_ind]->distance(xyz, uvw, coincident); } + +extern "C" void +surface_normal(int surf_ind, double xyz[3], double uvw[3]) { + return surfaces_c[surf_ind]->normal(xyz, uvw); +}