Merge pull request #357 from wuwenbin2006/some-fixes-for-explict3D

Some fixes for explicit3D ray tracing

src/Cell.cpp

@@ -1003,8 +1003,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;
 }
 
 
@@ -1062,7 +1065,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());
   }
@@ -1332,7 +1335,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 (fabs(x1 - x2) > FLT_EPSILON && num_zcylinders == 2)
@@ -1368,8 +1371,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 && fabs(radius1 - max_radius) < FLT_EPSILON
-        && _num_rings < 3)
+    if (fabs(radius1 - max_radius) < FLT_EPSILON && _num_rings < 3)
       increment = 1.5 * (radius1 - radius2) / _num_rings;
   }
 
@@ -1409,6 +1411,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) */
@@ -1420,8 +1423,10 @@ void Cell::ringify(std::vector<Cell*>& subcells, double max_radius) {
           rings.push_back(ring);
           continue;
         }
-        else
+        else {
           ring->addSurface(+1, *(iter2+1));
+          ring->removeSurface(zcylinder2);
+        }
 
         /* Store the clone in the parent Cell's container of ring Cells */
         rings.push_back(ring);
@@ -1436,6 +1441,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) */
@@ -1447,8 +1453,10 @@ void Cell::ringify(std::vector<Cell*>& subcells, double max_radius) {
         rings.push_back(ring);
         break;
       }
-      else
+      else {
         ring->addSurface(+1, *(iter2+1));
+        ring->removeSurface(zcylinder2);
+      }
 
       /* Store the clone in the parent Cell's container of ring Cells */
       rings.push_back(ring);
@@ -1559,7 +1567,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,15 @@ 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;
-
+  if (_loaded_from_file) {
+    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();
@@ -2070,12 +2072,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;

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(DEBUG, "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

@@ -617,7 +617,8 @@ std::string YPlane::toString() {
          << ", name = " << _name
          << ", type = YPLANE "
          << ", A = " << _A << ", B = " << _B
-         << ", C = " << _C << ", y = " << _y;
+         << ", C = " << _C << ", D = " << _D
+         << ", y = " << _y;
 
   return string.str();
 }
@@ -698,7 +699,8 @@ std::string ZPlane::toString() {
          << ", name = " << _name
          << ", type = ZPLANE "
          << ", A = " << _A << ", B = " << _B
-         << ", C = " << _C << ", z = " << _z;
+         << ", C = " << _C << ", D = " << _D
+         << ", z = " << _z;
 
   return string.str();
 }
@@ -991,11 +993,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 +1021,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 +1061,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/TrackGenerator3D.cpp

@@ -1133,7 +1133,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;
@@ -1263,8 +1263,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;