Merge pull request #1816 from paulromano/mesh-nan-fix

Prevent divide-by-zero in bins_crossed methods for meshes

.github/actions/fix-etc-hosts/action.yml

@@ -0,0 +1,18 @@
+name: Fix /etc/hosts
+description: |
+  Workaround for
+  "Reverse name lookup is broken for current hostname in ubuntu-latest VMs",
+  reported as https://github.com/actions/virtual-environments/issues/3185
+runs:
+  using: composite
+  steps:
+    - run: |
+        if [[ "$OSTYPE" == "linux-gnu"* ]]; then
+            # Ensure that reverse lookups for current hostname are handled properly
+            # Add the current IP address, long hostname and short hostname record to /etc/hosts file
+            eth0_ip_addr=$(ip addr show eth0 | grep "inet\b" | awk '{print $2}' | cut -d/ -f1)
+            hostname_fqdn=$(hostname -f)
+            hostname_short=$(hostname -s)
+            echo -e "${eth0_ip_addr}\t${hostname_fqdn} ${hostname_short}" | sudo tee -a /etc/hosts
+        fi
+      shell: bash

.github/workflows/ci.yml

@@ -76,8 +76,9 @@ jobs:
       LIBMESH: ${{ matrix.libmesh }}
 
     steps:
-      -
-        uses: actions/checkout@v2
+      - uses: actions/checkout@v2
+
+      - uses: ./.github/actions/fix-etc-hosts
 
       -
         name: Set up Python ${{ matrix.python-version }}

openmc/model/funcs.py

@@ -175,10 +175,17 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.),
     max_x1 = plane(x1, 'maximum', width/2 + origin[0])
     min_x2 = plane(x2, 'minimum', -height/2 + origin[1])
     max_x2 = plane(x2, 'maximum', height/2 + origin[1])
+    if boundary_type == 'periodic':
+        min_x1.periodic_surface = max_x1
+        min_x2.periodic_surface = max_x2
     prism = +min_x1 & -max_x1 & +min_x2 & -max_x2
 
     # Handle rounded corners if given
     if corner_radius > 0.:
+        if boundary_type == 'periodic':
+            raise ValueError('Periodic boundary conditions not permitted when '
+                             'rounded corners are used.')
+
         args = {'R': corner_radius, 'boundary_type': boundary_type}
 
         args[x1 + '0'] = origin[0] - width/2 + corner_radius

src/mesh.cpp

@@ -594,6 +594,7 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
 
   // Compute the length of the entire track.
   double total_distance = (r - last_r).norm();
+  if (total_distance == 0.0) return;
 
   // While determining if this track intersects the mesh, offset the starting
   // and ending coords by a bit.  This avoid finite-precision errors that can
@@ -1713,6 +1714,7 @@ MOABMesh::bins_crossed(const Particle& p,
   moab::CartVect dir(u.x, u.y, u.z);
 
   double track_len = (r1 - r0).length();
+  if (track_len == 0.0) return;
 
   r0 -= TINY_BIT * dir;
   r1 += TINY_BIT * dir;