Implement z planes and cylinders in C++

CMakeLists.txt

@@ -436,7 +436,10 @@ set(LIBOPENMC_FORTRAN_SRC
 )
 set(LIBOPENMC_CXX_SRC
   src/random_lcg.h
-  src/random_lcg.cpp)
+  src/random_lcg.cpp
+  src/surface_header.C
+  src/pugixml/pugixml.hpp
+  src/pugixml/pugixml.cpp)
 add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
 set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc)
 add_executable(${program} src/main.F90)

src/geometry.F90

@@ -10,8 +10,33 @@ module geometry
   use stl_vector,             only: VectorInt
   use string,                 only: to_str
 
+  use, intrinsic :: ISO_C_BINDING
+
   implicit none
 
+  interface
+    function surface_sense_c(surf_ind, xyz, uvw) &
+         bind(C, name='surface_sense') result(sense)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value :: surf_ind;
+      real(C_DOUBLE) :: xyz(3);
+      real(C_DOUBLE) :: uvw(3);
+      logical(C_BOOL) :: sense;
+    end function surface_sense_c
+
+    function surface_distance_c(surf_ind, xyz, uvw, coincident) &
+         bind(C, name='surface_distance') result(d)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value :: surf_ind;
+      real(C_DOUBLE) :: xyz(3);
+      real(C_DOUBLE) :: uvw(3);
+      logical(C_BOOL), value :: coincident;
+      real(C_DOUBLE) :: d;
+    end function surface_distance_c
+  end interface
+
 contains
 
 !===============================================================================
@@ -28,7 +53,8 @@ contains
 ! expression using a stack, similar to how a RPN calculator would work.
 !===============================================================================
 
-  pure function cell_contains(c, p) result(in_cell)
+  !pure function cell_contains(c, p) result(in_cell)
+  function cell_contains(c, p) result(in_cell)
     type(Cell), intent(in) :: c
     type(Particle), intent(in) :: p
     logical :: in_cell
@@ -40,7 +66,8 @@ contains
     end if
   end function cell_contains
 
-  pure function simple_cell_contains(c, p) result(in_cell)
+  !pure function simple_cell_contains(c, p) result(in_cell)
+  function simple_cell_contains(c, p) result(in_cell)
     type(Cell), intent(in) :: c
     type(Particle), intent(in) :: p
     logical :: in_cell
@@ -48,6 +75,7 @@ contains
     integer :: i
     integer :: token
     logical :: actual_sense    ! sense of particle wrt surface
+    logical :: alt_sense
 
     in_cell = .true.
     do i = 1, size(c%rpn)
@@ -65,6 +93,15 @@ contains
         else
           actual_sense = surfaces(abs(token))%obj%sense(&
                p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
+          alt_sense = surface_sense_c(abs(token)-1, &
+               p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
+          if (actual_sense .neqv. alt_sense) then
+            write(*, *) surfaces(abs(token)) % obj % id
+            write(*, *) p % coord (p % n_coord) % xyz
+            write(*, *) p % coord (p % n_coord) % uvw
+            write(*, *) actual_sense, alt_sense
+            call fatal_error("")
+          end if
           if (actual_sense .neqv. (token > 0)) then
             in_cell = .false.
             exit
@@ -484,6 +521,7 @@ contains
     type(Cell),       pointer :: c
     class(Surface),   pointer :: surf
     class(Lattice),   pointer :: lat
+    real(8) :: d_alt
 
     ! inialize distance to infinity (huge)
     dist = INFINITY
@@ -519,6 +557,11 @@ contains
         ! Calculate distance to surface
         surf => surfaces(index_surf) % obj
         d = surf % distance(p % coord(j) % xyz, p % coord(j) % uvw, coincident)
+        d_alt = surface_distance_c(index_surf-1, p % coord(j) % xyz, &
+             p % coord(j) % uvw, logical(coincident, kind=C_BOOL))
+        if (d /= d_alt) then
+          write(*, *) d, d_alt
+        end if
 
         ! Check if calculated distance is new minimum
         if (d < d_surf) then

src/input_xml.F90

@@ -48,6 +48,14 @@ module input_xml
   implicit none
   save
 
+  interface
+    subroutine read_surfaces(node_ptr) bind(C, name='read_surfaces')
+      use ISO_C_BINDING
+      implicit none
+      type(C_PTR) :: node_ptr
+    end subroutine read_surfaces
+  end interface 
+
 contains
 
 !===============================================================================
@@ -966,6 +974,7 @@ contains
 
     ! get pointer to list of xml <surface>
     call get_node_list(root, "surface", node_surf_list)
+    call read_surfaces(root % ptr)
 
     ! Get number of <surface> tags
     n_surfaces = size(node_surf_list)

src/surface_header.C

@@ -0,0 +1,280 @@
+#include <cstring>  // For strcmp
+#include <iostream>
+#include <limits>  // For numeric_limits
+#include <math.h>  // For fabs
+#include "pugixml/pugixml.hpp"
+
+//==============================================================================
+// Constants
+//==============================================================================
+
+const double FP_COINCIDENT = 1e-12;
+const double INFTY = std::numeric_limits<double>::max();
+
+//==============================================================================
+// Global array of surfaces
+//==============================================================================
+
+class Surface;
+Surface **surfaces_c;
+
+//==============================================================================
+// SURFACE type defines a first- or second-order surface that can be used to
+// construct closed volumes (cells)
+//==============================================================================
+
+class Surface
+{
+  int id;                    // Unique ID
+  int neighbor_pos[],        // List of cells on positive side
+      neighbor_neg[];        // List of cells on negative side
+  int bc;                    // Boundary condition
+  //TODO: swith that zero to a NONE constant.
+  int i_periodic = 0;        // Index of corresponding periodic surface
+  char name[104];            // User-defined name
+
+public:
+  bool sense(double xyz[3], double uvw[3]);
+  void reflect(double xyz[3], double uvw[3]);
+  virtual double evaluate(double xyz[3]) = 0;
+  virtual double distance(double xyz[3], double uvw[3], bool coincident) = 0;
+  virtual void normal(double xyz[3], double uvw[3]) = 0;
+};
+
+bool
+Surface::sense(double xyz[3], double uvw[3])
+{
+  // Evaluate the surface equation at the particle's coordinates to determine
+  // which side the particle is on.
+  const double f = evaluate(xyz);
+
+  // Check which side of surface the point is on.
+  bool s;
+  if (fabs(f) < FP_COINCIDENT)
+  {
+    // Particle may be coincident with this surface. To determine the sense, we
+    // look at the direction of the particle relative to the surface normal (by
+    // default in the positive direction) via their dot product.
+    double norm[3];
+    normal(xyz, norm);
+    s = (uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0);
+  }
+  else
+  {
+    s = (f > 0.0);
+  }
+  return s;
+}
+
+void
+Surface::reflect(double xyz[3], double uvw[3])
+{
+
+  // Determine projection of direction onto normal and squared magnitude of
+  // normal.
+  double norm[3];
+  normal(xyz, norm);
+  const double projection = norm[0]*uvw[0] + norm[1]*uvw[1] + norm[2]*uvw[2];
+  const double magnitude = norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2];
+
+  // Reflect direction according to normal.
+  uvw[0] -= 2.0 * projection / magnitude * norm[0];
+  uvw[1] -= 2.0 * projection / magnitude * norm[1];
+  uvw[2] -= 2.0 * projection / magnitude * norm[2];
+}
+
+//==============================================================================
+
+class SurfaceZPlane : public Surface
+{
+  double z0;
+public:
+  SurfaceZPlane(pugi::xml_node surf_node);
+  double evaluate(double xyz[3]);
+  double distance(double xyz[3], double uvw[3], bool coincident);
+  void normal(double xyz[3], double uvw[3]);
+};
+
+SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node)
+{
+  const char *coeffs = surf_node.attribute("coeffs").value();
+  int stat = sscanf(coeffs, "%lf", &z0);
+  if (stat != 1)
+  {
+    std::cout << "Something went wrong reading surface coeffs!" << std::endl;
+  }
+}
+
+double
+SurfaceZPlane::evaluate(double xyz[3])
+{
+  double f = xyz[2] - z0;
+  return f;
+}
+
+double
+SurfaceZPlane::distance(double xyz[3], double uvw[3], bool coincident)
+{
+  double f = z0 - xyz[2];
+  double d;
+  if (coincident or fabs(f) < FP_COINCIDENT or uvw[2] == 0.0)
+  {
+    d = INFTY;
+  }
+  else
+  {
+    d = f / uvw[2];
+    if (d < 0.0) d = INFTY;
+  }
+  return d;
+}
+
+void
+SurfaceZPlane::normal(double xyz[3], double uvw[3])
+{
+  uvw[0] = 0.0;
+  uvw[1] = 0.0;
+  uvw[2] = 1.0;
+}
+
+//==============================================================================
+
+class SurfaceZCylinder : public Surface
+{
+  double x0, y0, r;
+public:
+  SurfaceZCylinder(pugi::xml_node surf_node);
+  double evaluate(double xyz[3]);
+  double distance(double xyz[3], double uvw[3], bool coincident);
+  void normal(double xyz[3], double uvw[3]);
+};
+
+SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
+{
+  const char *coeffs = surf_node.attribute("coeffs").value();
+  int stat = sscanf(coeffs, "%lf %lf %lf", &x0, &y0, &r);
+  if (stat != 3)
+  {
+    std::cout << "Something went wrong reading surface coeffs!" << std::endl;
+  }
+}
+
+double
+SurfaceZCylinder::evaluate(double xyz[3])
+{
+  const double x = xyz[0] - x0;
+  const double y = xyz[1] - y0;
+  double f = x*x + y*y - r*r;
+  return f;
+}
+
+double
+SurfaceZCylinder::distance(double xyz[3], double uvw[3], bool coincident)
+{
+  double a = 1.0 - uvw[2]*uvw[2];  // u^2 + v^2
+
+  if (a == 0.0) return INFTY;
+
+  double x = xyz[0] - x0;
+  double y = xyz[1] - y0;
+  double k = x*uvw[0] + y*uvw[1];
+  double c = x*x + y*y - r*r;
+  double quad = k*k - a*c;
+
+  if (quad < 0.0)
+  {
+    // No intersection with cylinder.
+    return INFTY;
+  }
+  else if (coincident or fabs(c) < FP_COINCIDENT)
+  {
+    // Particle is on the cylinder, thus one distance is positive/negative
+    // and the other is zero. The sign of k determines if we are facing in or
+    // out.
+    if (k >= 0.0) return INFTY;
+    else return (-k + sqrt(quad)) / a;
+  }
+  else if (c < 0.0)
+  {
+    // Particle is inside the cylinder, thus one distance must be negative
+    // and one must be positive. The positive distance will be the one with
+    // negative sign on sqrt(quad)
+    return (-k + sqrt(quad)) / a;
+  }
+  else
+  {
+    // Particle is outside the cylinder, thus both distances are either
+    // positive or negative. If positive, the smaller distance is the one
+    // with positive sign on sqrt(quad)
+
+    double d = (-k - sqrt(quad)) / a;
+    if (d < 0.0) return INFTY;
+    return d;
+  }
+}
+
+void
+SurfaceZCylinder::normal(double xyz[3], double uvw[3])
+{
+  uvw[0] = 2.0 * (xyz[0] - x0);
+  uvw[1] = 2.0 * (xyz[1] - y0);
+  uvw[2] = 0.0;
+}
+
+//==============================================================================
+
+extern "C" void
+read_surfaces(pugi::xml_node *node)
+{
+  // Count the number of surfaces.
+  int n_surfaces = 0;
+  for (pugi::xml_node surf_node = node->child("surface"); surf_node;
+       surf_node = surf_node.next_sibling("surface"))
+  {
+    n_surfaces += 1;
+  }
+
+  // Allocate the array of Surface pointers.
+  surfaces_c = new Surface* [n_surfaces];
+
+  // Loop over XML surface elements and populate the array.
+  {
+    pugi::xml_node surf_node;
+    int i_surf;
+    for (surf_node = node->child("surface"), i_surf = 0; surf_node;
+         surf_node = surf_node.next_sibling("surface"), i_surf++)
+    {
+      if (surf_node.attribute("type"))
+      {
+        const pugi::char_t *surf_type = surf_node.attribute("type").value();
+        if (strcmp(surf_type, "z-cylinder") == 0)
+        {
+          surfaces_c[i_surf] = new SurfaceZCylinder(surf_node);
+        }
+        else if (strcmp(surf_type, "z-plane") == 0)
+        {
+          surfaces_c[i_surf] = new SurfaceZPlane(surf_node);
+        }
+        else
+        {
+          std::cout << "Call error or handle uppercase here!" << std::endl;
+          std::cout << surf_type << std::endl;
+        }
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+extern "C" bool
+surface_sense(int surf_ind, double xyz[3], double uvw[3])
+{
+  return surfaces_c[surf_ind]->sense(xyz, uvw);
+}
+
+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);
+}