Some fixes for explicit3D ray tracing

src/Cell.cpp

@@ -299,7 +299,7 @@ double* Cell::getRotationMatrix() {
  * @param units the angular units in "radians" or "degrees" (default)
  */
 void Cell::retrieveRotation(double* rotations, int num_axes,
-			    std::string units) {
+          std::string units) {
   if (num_axes != 3)
     log_printf(ERROR, "Unable to get rotation with %d axes for Cell %d. "
                "The rotation array should be length 3.", num_axes, _id);
@@ -1002,8 +1002,11 @@ void Cell::setNumRings(int num_rings) {
   if (num_rings < 0)
     log_printf(ERROR, "Unable to give %d rings to Cell %d since this is "
                "a negative number", num_rings, _id);
-
-  _num_rings = num_rings;
+
+  if (num_rings == 1)
+    _num_rings = 0;
+  else
+    _num_rings = num_rings;
 }
 
 
@@ -1061,7 +1064,7 @@ void Cell::addSurface(int halfspace, Surface* surface) {
  */
 void Cell::removeSurface(Surface* surface) {
 
-  if (_surfaces.find(surface->getId()) != _surfaces.end()) {
+  if (surface != NULL && _surfaces.find(surface->getId()) != _surfaces.end()) {
     delete _surfaces[surface->getId()];
     _surfaces.erase(surface->getId());
   }
@@ -1331,7 +1334,7 @@ void Cell::ringify(std::vector<Cell*>& subcells, double max_radius) {
   }
 
   if (num_zcylinders > 2)
-    log_printf(NORMAL, "Unable to ringify Cell %d since it "
+    log_printf(ERROR, "Unable to ringify Cell %d since it "
                "contains more than 2 ZCYLINDER Surfaces", _id);
 
   if (x1 != x2 && num_zcylinders == 2)
@@ -1367,7 +1370,7 @@ void Cell::ringify(std::vector<Cell*>& subcells, double max_radius) {
     increment = fabs(radius1 - radius2) / _num_rings;
 
     /* Heuristic to improve area-balancing for low number of rings */
-    if (halfspace1 == 0 && radius1 == max_radius && _num_rings < 3)
+    if (radius1 == max_radius && _num_rings < 3)
       increment = 1.5 * (radius1 - radius2) / _num_rings;
   }
 
@@ -1407,6 +1410,7 @@ void Cell::ringify(std::vector<Cell*>& subcells, double max_radius) {
         Cell* ring = (*iter3)->clone();
         ring->setNumSectors(0);
         ring->setNumRings(0);
+        ring->removeSurface(zcylinder1);
 
         /* Add ZCylinder only if this is not the outermost ring in an
          * unbounded Cell (i.e. the moderator in a fuel pin cell) */
@@ -1434,6 +1438,7 @@ void Cell::ringify(std::vector<Cell*>& subcells, double max_radius) {
       Cell* ring = clone();
       ring->setNumSectors(0);
       ring->setNumRings(0);
+      ring->removeSurface(zcylinder1);
 
       /* Add ZCylinder only if this is not the outermost ring in an
        * unbounded Cell (i.e. the moderator in a fuel pin cell) */
@@ -1557,7 +1562,8 @@ std::string Cell::toString() {
   std::map<int, surface_halfspace*>::iterator iter;
   string << ", Surfaces: ";
   for (iter = _surfaces.begin(); iter != _surfaces.end(); ++iter)
-    string <<  iter->second->_surface->toString() << ", ";
+    string << "\nhalfspace = " << iter->second->_halfspace << ", " <<
+           iter->second->_surface->toString();
 
   return string.str();
 }

src/Geometry.cpp

@@ -42,13 +42,6 @@ Geometry::Geometry() {
  * @brief Destructor clears FSR to Cells and Materials maps.
  */
 Geometry::~Geometry() {
-
-  /* Free all materials */
-  std::map<int, Material*> materials = _root_universe->getAllMaterials();
-  std::map<int, Material*>::iterator iter;
-  for (iter = materials.begin(); iter != materials.end(); ++iter)
-    delete iter->second;
-
   /* Free all surfaces */
   if (_loaded_from_file) {
     std::map<int, Surface*> surfaces = getAllSurfaces();
@@ -422,8 +415,8 @@ std::map<int, Surface*> Geometry::getAllSurfaces() {
       surfs = cell->getSurfaces();
 
       for (s_iter = surfs.begin(); s_iter != surfs.end(); ++s_iter) {
-	surf = (*s_iter).second->_surface;
-	all_surfs[surf->getId()] = surf;
+  surf = (*s_iter).second->_surface;
+  all_surfs[surf->getId()] = surf;
       }
     }
   }
@@ -2052,12 +2045,12 @@ void Geometry::segmentize3D(Track3D* track, bool setup) {
     }
 
     /* Calculate the local centroid of the segment if available */
-    if (_contains_FSR_centroids && !setup) {
-      Point* centroid = getFSRCentroid(fsr_id);
+    if (!setup) {
+      //Point* centroid = getFSRCentroid(fsr_id);
       Point* starting_point = start.getHighestLevel()->getPoint();
-      double x_start = starting_point->getX() - centroid->getX();
-      double y_start = starting_point->getY() - centroid->getY();
-      double z_start = starting_point->getZ() - centroid->getZ();
+      double x_start = starting_point->getX();
+      double y_start = starting_point->getY();
+      double z_start = starting_point->getZ();
       new_segment->_starting_position[0] = x_start;
       new_segment->_starting_position[1] = y_start;
       new_segment->_starting_position[2] = z_start;
@@ -2625,11 +2618,11 @@ void Geometry::computeFissionability(Universe* univ) {
 
   Material* material;
   std::map<int, Material*> materials;
-	std::map<int, Material*>::iterator mat_iter;
+  std::map<int, Material*>::iterator mat_iter;
 
   Universe* universe;
   std::map<int, Universe*> universes;
-	std::map<int, Universe*>::iterator univ_iter;
+  std::map<int, Universe*>::iterator univ_iter;
 
   /* If no Universe was passed in as an argument, then this is the first
    * recursive call from a user via Python, so get the base Universe */

src/LocalCoords.cpp

@@ -203,12 +203,14 @@ void LocalCoords::detectLoop() {
   int n = 0;
 
   LocalCoords* iter = _next;
-  while (_next != NULL) {
+  while (iter != NULL) {
       iter = iter->getNext();
       n++;
       if (n > 1000)
         log_printf(ERROR, "Infinite loop of coords");
   }
+  log_printf(DEBUG, "The LocalCoords is: %s\n", toString().c_str());
+  log_printf(NORMAL, "The depth of the chain is %d \n", n);
 }
 
 
@@ -413,7 +415,7 @@ LocalCoords* LocalCoords::getHighestLevel() {
 
 
 /**
- * @brief Translate all of the x,y coordinates for each LocalCoords object in
+ * @brief Translate all of the x,y,z coordinates for each LocalCoords object in
  *        the linked list.
  * @details This method will traverse the entire linked list and apply the
  *          translation to each element.
@@ -475,7 +477,7 @@ void LocalCoords::prune() {
   LocalCoords* next = curr->getPrev();
 
   /* Iterate over LocalCoords beneath this one in the linked list */
-  while (curr != this) {
+  while (curr != this) {// it seems nothing has been down in this while loop
     next = curr->getPrev();
     if (curr->getPosition() == -1)
       curr->deleteArray();

src/Surface.cpp

@@ -991,11 +991,9 @@ int ZCylinder::intersection(Point* point, double azim, double polar, Point* poin
     else {
 
       /* Determine first point of intersection */
-      double interior = pow(ycurr - y0, 2.0) + pow(xcurr - x0, 2.0);
-      if (interior < 0.0)
-        interior = 0.0;
       xcurr = (-b + sqrt(discr)) / (2*a);
       ycurr = y0 + m * (xcurr - x0);
+      double interior = pow(ycurr - y0, 2.0) + pow(xcurr - x0, 2.0);
       zcurr = z0 + sqrt(interior) * tan(M_PI_2 - polar);
       points[num].setCoords(xcurr, ycurr, zcurr);
 
@@ -1021,6 +1019,7 @@ int ZCylinder::intersection(Point* point, double azim, double polar, Point* poin
       /* Determine second point of intersection */
       xcurr = (-b - sqrt(discr)) / (2*a);
       ycurr = y0 + m * (xcurr - x0);
+      interior = pow(ycurr - y0, 2.0) + pow(xcurr - x0, 2.0);
       zcurr = z0 + sqrt(interior) * tan(M_PI_2 - polar);
       points[num].setCoords(xcurr, ycurr, zcurr);
 
@@ -1060,7 +1059,7 @@ std::string ZCylinder::toString() {
   std::stringstream string;
 
   string << "Surface ID = " << _id
-         << ", name " << _name
+         << ", name = " << _name
          << ", type = ZCYLINDER "
          << ", A = " << _A << ", B = " << _B
          << ", C = " << _C << ", D = " << _D << ", E = " << _E

src/TrackGenerator.cpp

@@ -1625,7 +1625,15 @@ void TrackGenerator::resetStatus() {
 void TrackGenerator::allocateTemporarySegments() {}
 
 
-//FIXME: description
+/**
+ * @brief Compute and return the track UID number
+ * @details Compute and return the track UID number based on the azimuthal 
+            number and the xy_index. Track UID number is also the index in the 
+            Track** TrackGenerator::_tracks_2D_array
+ * @param a the azimuthal number of a track
+ * @param x the xy_index of a track
+ * @return the track UID number of a track 
+ */
 int TrackGenerator::get2DTrackID(int a, int x) {
 
   int uid = 0;

src/TrackGenerator3D.cpp

@@ -908,7 +908,13 @@ void TrackGenerator3D::initialize2DTrackChains() {
 }
 
 
-//FIXME: description
+/**
+ * @brief Traverses all the 3D chain tracks, and split each 3D chain track into 
+          3D tracks belonging to the z-stacks
+ * @param tcis The array of track chain indices
+ * @param num_chains The size of the array of track chain indices
+ * @param save_tracks Whether or not to save the 3D tracks
+ */
 void TrackGenerator3D::getCycleTrackData(TrackChainIndexes* tcis,
                                          int num_chains, bool save_tracks) {
 
@@ -942,7 +948,13 @@ void TrackGenerator3D::getCycleTrackData(TrackChainIndexes* tcis,
 }
 
 
-//FIXME: description
+/**
+ * @brief Compute the coordinates of a 3D chain track, and the link number of 
+          the 2D track where its start point reseide on
+ * @param tci the track chain index of a chain track
+ * @param track_3D the chain track
+ * @return the link number of the 2D track where its start point reseide on
+ */
 int TrackGenerator3D::getFirst2DTrackLinkIndex(TrackChainIndexes* tci,
                                                Track3D* track_3D) {
 
@@ -1035,7 +1047,7 @@ int TrackGenerator3D::getFirst2DTrackLinkIndex(TrackChainIndexes* tci,
  *          boundaries is split into multiple tracks by finding those
  *          x and y intersections. If create_tracks is set to true, the Tracks
  *          will be altered to represent the correct 3D Tracks. If not, the
- *          number of tracks in the train will simply be counted. Whenever this
+ *          number of tracks in the chain will simply be counted. Whenever this
  *          function is called, the number of tracks in the associated z-stacks
  *          are incremented.
  * @param track The 3D track to be decomposed
@@ -1128,7 +1140,7 @@ void TrackGenerator3D::set3DTrackData(TrackChainIndexes* tci,
       break;
 
     /* If the Z boundary or last link or the desired link has been reached,
-     * exit */
+     * exit. tci->_link appoints the desired track*/
     if (dl_z < dl_xy || track_2D->getXYIndex() >= _num_y[tci->_azim] ||
         tci->_link == link - first_link)
       end_of_chain = true;
@@ -1258,8 +1270,15 @@ void TrackGenerator3D::segmentize() {
 #pragma omp parallel for
     for (int i=0; i < _num_x[a] + _num_y[a]; i++) {
       for (int p=0; p < _num_polar; p++) {
-        for (int z=0; z < _tracks_per_stack[a][i][p]; z++)
+        for (int z=0; z < _tracks_per_stack[a][i][p]; z++){
           _geometry->segmentize3D(&_tracks_3D[a][i][p][z]);
+          TrackStackIndexes tsi;
+          tsi._azim = a;
+          tsi._xy = i;
+          tsi._polar = p;
+          tsi._z = z;
+          _tracks_3D[a][i][p][z].setUid(get3DTrackID(&tsi));    
+        }
       }
     }
 

src/TraverseSegments.cpp

@@ -572,7 +572,7 @@ void TraverseSegments::traceStackOTF(Track* flattened_track, int polar_index,
 
   /* Set the current x and y coordinates */
   double x_curr = x_start_2D;
-  double y_curr = flattened_track->getStart()->getY();
+  double y_curr = y_start_2D;
 
   /* Loop over 2D segments */
   double first_start_z = start_z;