Merge #462: cd, xyz and box geometry errors

Fix #388, #383: concave CD, XYZ inside CD, conestack, and box geometry errors
nrc-cnrc · Mar 7, 2019 · 6f90abd · 6f90abd
2 parents c65b9cc + 116dfd3
commit 6f90abd
Show file tree

Hide file tree

Showing 8 changed files with 63 additions and 22 deletions.
diff --git a/HEN_HOUSE/egs++/egs_base_geometry.cpp b/HEN_HOUSE/egs++/egs_base_geometry.cpp
@@ -431,6 +431,8 @@ EGS_BaseGeometry::EGS_BaseGeometry(const string &Name) : nreg(0), name(Name),
  has_B_scaling(false), has_Ref_rho(false), bfactor(0), rhoRef(1.0),
  nref(0), debug(false), is_convex(true), bproperty(0), bp_array(0),
  boundaryTolerance(epsilon) {
+
+ halfBoundaryTolerance = boundaryTolerance/2.;
  if (!egs_geometries.size()) {
  egs_geometries.addList(new EGS_GeometryPrivate);
  }
@@ -716,7 +718,9 @@ void EGS_BaseGeometry::setBoundaryTolerance(EGS_Input *i) {
  int err = i->getInput("boundary tolerance", boundaryTolerance);
  if (err > 0) {
  egsWarning("EGS_BaseGeometry::setBoundaryTolerance(): error while reading 'boundary tolerance' input\n");
+ return;
  }
+ halfBoundaryTolerance = boundaryTolerance/2.;
 }
 
 void EGS_BaseGeometry::printInfo() const {

diff --git a/HEN_HOUSE/egs++/egs_base_geometry.h b/HEN_HOUSE/egs++/egs_base_geometry.h
@@ -604,6 +604,7 @@ class EGS_EXPORT EGS_BaseGeometry {
  */
  void setBoundaryTolerance(EGS_Float tol) {
  boundaryTolerance = tol;
+ halfBoundaryTolerance = tol/2.;
  }
 
  /*! \brief Is the boolean property \a prop set for region \a ireg ?
@@ -833,7 +834,7 @@ class EGS_EXPORT EGS_BaseGeometry {
  EGS_BPType *bp_array;
 
  /*! \brief Boundary tolerance for geometries that need it */
- EGS_Float boundaryTolerance;
+ EGS_Float boundaryTolerance, halfBoundaryTolerance;
 
  /*! \brief Set to non-zero status if a geometry problem is encountered */
  static int error_flag;

diff --git a/HEN_HOUSE/egs++/geometry/egs_box/egs_box.h b/HEN_HOUSE/egs++/geometry/egs_box/egs_box.h
@@ -236,8 +236,8 @@ class EGS_BOX_EXPORT EGS_Box : public EGS_BaseGeometry {
  }
  if (t1 < t) {
  EGS_Float y1 = xp.y + up.y*t1, z1 = xp.z + up.z*t1;
- if (2*y1 + ay >= 0 && 2*y1 - ay <= 0 &&
- 2*z1 + az >= 0 && 2*z1 - az <= 0) {
+ if (2*y1 + ay > 0 && 2*y1 - ay < 0 &&
+ 2*z1 + az > 0 && 2*z1 - az < 0) {
  t = t1;
  if (newmed) {
  *newmed = med;
@@ -273,8 +273,8 @@ class EGS_BOX_EXPORT EGS_Box : public EGS_BaseGeometry {
  }
  if (t1 < t) {
  EGS_Float x1 = xp.x + up.x*t1, z1 = xp.z + up.z*t1;
- if (2*x1 + ax >= 0 && 2*x1 - ax <= 0 &&
- 2*z1 + az >= 0 && 2*z1 - az <= 0) {
+ if (2*x1 + ax > 0 && 2*x1 - ax < 0 &&
+ 2*z1 + az > 0 && 2*z1 - az < 0) {
  t = t1;
  if (newmed) {
  *newmed = med;
@@ -310,8 +310,8 @@ class EGS_BOX_EXPORT EGS_Box : public EGS_BaseGeometry {
  }
  if (t1 < t) {
  EGS_Float x1 = xp.x + up.x*t1, y1 = xp.y + up.y*t1;
- if (2*x1 + ax >= 0 && 2*x1 - ax <= 0 &&
- 2*y1 + ay >= 0 && 2*y1 - ay <= 0) {
+ if (2*x1 + ax > 0 && 2*x1 - ax < 0 &&
+ 2*y1 + ay > 0 && 2*y1 - ay < 0) {
  t = t1;
  if (newmed) {
  *newmed = med;

diff --git a/HEN_HOUSE/egs++/geometry/egs_cd_geometry/egs_cd_geometry.h b/HEN_HOUSE/egs++/geometry/egs_cd_geometry/egs_cd_geometry.h
@@ -565,11 +565,15 @@ class EGS_CDGEOMETRY_EXPORT EGS_CDGeometry : public EGS_BaseGeometry {
  }
  // Skip for concave geometries: particles can reenter the CD geometry after leaving.
  // Check 2. below will catch these cases.
- else if ((ixnew < 0 && tb <= boundaryTolerance) && (bg->isConvex() &&
- (g[ibase] && g[ibase]->isConvex()))) {
- // (a) is true
- return ixnew;
- }
+ // This follow has been commented out because it fails in cases
+ // where the inscribe geometries are themselves convex but
+ // form a concave shape. Instead, just continue to do more checks.
+// else if ((ixnew < 0 && tb <= boundaryTolerance) && (bg->isConvex() &&
+// (g[ibase] && g[ibase]->isConvex()))) {
+// // (a) is true
+// return ixnew;
+// }
+
  // If a particle approaching the geometry sits on a boundary, we look back to see
  // if we just entered the geometry (the previous checks fail to catch this case).
  EGS_Float tb_neg = veryFar;

diff --git a/HEN_HOUSE/egs++/geometry/egs_cones/egs_cones.h b/HEN_HOUSE/egs++/geometry/egs_cones/egs_cones.h
@@ -1554,6 +1554,12 @@ class EGS_ConeStack : public EGS_BaseGeometry {
  t = tp; // set t = distance to plane
  }
 
+ // If we are right on a plane, return -1
+ if(tp < boundaryTolerance) {
+ t = 0;
+ return -1;
+ }
+
  // distance to outer cone
  int ir = ireg - il*nmax; // cone index in current layer
  bool hitc = false; // assume we don't hit the cone

diff --git a/HEN_HOUSE/egs++/geometry/egs_cylinders/egs_cylinders.h b/HEN_HOUSE/egs++/geometry/egs_cylinders/egs_cylinders.h
@@ -333,7 +333,7 @@ class EGS_CYLINDERS_EXPORT EGS_CylindersT : public EGS_BaseGeometry {
  egsWarning("t=%g d=%g\n",last_t,last_d);
 #endif
  }
- d = boundaryTolerance;
+ d = halfBoundaryTolerance;
  }
  }
 
@@ -377,7 +377,7 @@ class EGS_CYLINDERS_EXPORT EGS_CylindersT : public EGS_BaseGeometry {
  if (C < -boundaryTolerance) egsWarning("EGS_CylindersT::howfar(): "
  "the particle may not be in the region we think it "
  "is as Cin = %g\n",C);
- d = boundaryTolerance;
+ d = halfBoundaryTolerance;
  }
  }
  }
@@ -401,7 +401,7 @@ class EGS_CYLINDERS_EXPORT EGS_CylindersT : public EGS_BaseGeometry {
  egsWarning(" ireg=%d x=(%g,%g,%g) u=(%g,%g,%g)\n",
  ireg,x.x,x.y,x.z,u.x,u.y,u.z);
  }
- d = boundaryTolerance;
+ d = halfBoundaryTolerance;
  }
  }
  }

diff --git a/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h b/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h
@@ -901,7 +901,13 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
  if (u.x > 0) {
  EGS_Float d = (xpos[ix+1]-x.x)/u.x;
  if (d <= t) {
- t = d;
+ if(d <= boundaryTolerance) {
+ // t=0 works on most cases but can result in
+ // getting stuck on an edge, so use boundaryTolerance
+ t = boundaryTolerance;
+ } else {
+ t = d;
+ }
  if ((++ix) < nx) {
  inew = ireg + 1;
  }
@@ -916,7 +922,11 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
  else if (u.x < 0) {
  EGS_Float d = (xpos[ix]-x.x)/u.x;
  if (d <= t) {
- t = d;
+ if(d <= boundaryTolerance) {
+ t = boundaryTolerance;
+ } else {
+ t = d;
+ }
  if ((--ix) >= 0) {
  inew = ireg - 1;
  }
@@ -931,7 +941,11 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
  if (u.y > 0) {
  EGS_Float d = (ypos[iy+1]-x.y)/u.y;
  if (d <= t) {
- t = d;
+ if(d <= boundaryTolerance) {
+ t = boundaryTolerance;
+ } else {
+ t = d;
+ }
  if ((++iy) < ny) {
  inew = ireg + nx;
  }
@@ -946,7 +960,11 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
  else if (u.y < 0) {
  EGS_Float d = (ypos[iy]-x.y)/u.y;
  if (d <= t) {
- t = d;
+ if(d <= boundaryTolerance) {
+ t = boundaryTolerance;
+ } else {
+ t = d;
+ }
  if ((--iy) >= 0) {
  inew = ireg - nx;
  }
@@ -961,7 +979,11 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
  if (u.z > 0) {
  EGS_Float d = (zpos[iz+1]-x.z)/u.z;
  if (d <= t) {
- t = d;
+ if(d <= boundaryTolerance) {
+ t = boundaryTolerance;
+ } else {
+ t = d;
+ }
  if ((++iz) < nz) {
  inew = ireg+nxy;
  }
@@ -976,7 +998,11 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
  else if (u.z < 0) {
  EGS_Float d = (zpos[iz]-x.z)/u.z;
  if (d <= t) {
- t = d;
+ if(d <= boundaryTolerance) {
+ t = boundaryTolerance;
+ } else {
+ t = d;
+ }
  if ((--iz) >= 0) {
  inew = ireg-nxy;
  }

diff --git a/HEN_HOUSE/egs++/geometry/egs_spheres/egs_spheres.cpp b/HEN_HOUSE/egs++/geometry/egs_spheres/egs_spheres.cpp
@@ -173,7 +173,7 @@ int EGS_cSpheres::howfar(int ireg, const EGS_Vector &x,
  */
  R2b2 = R2[ireg] - bb2;
  if (R2b2 <= 0 && aa > 0) {
- d = boundaryTolerance; // hopefully a truncation problem
+ d = halfBoundaryTolerance; // hopefully a truncation problem
  }
  else {
  tmp = aa2 + R2b2;