diff --git a/src/Cell.cpp b/src/Cell.cpp index 2f9b5e33b..08229ffd6 100644 --- a/src/Cell.cpp +++ b/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& 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& 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& 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& 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& 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& 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::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(); } diff --git a/src/Geometry.cpp b/src/Geometry.cpp index 3f4a12bd6..d4aefd329 100644 --- a/src/Geometry.cpp +++ b/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 materials = _root_universe->getAllMaterials(); - std::map::iterator iter; - for (iter = materials.begin(); iter != materials.end(); ++iter) - delete iter->second; - + if (_loaded_from_file) { + std::map materials = _root_universe->getAllMaterials(); + std::map::iterator iter; + for (iter = materials.begin(); iter != materials.end(); ++iter) + delete iter->second; + } + /* Free all surfaces */ if (_loaded_from_file) { std::map 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; diff --git a/src/LocalCoords.cpp b/src/LocalCoords.cpp index c072c16ed..6c8f9196c 100644 --- a/src/LocalCoords.cpp +++ b/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(); diff --git a/src/Surface.cpp b/src/Surface.cpp index 85a662455..fa815c948 100644 --- a/src/Surface.cpp +++ b/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 diff --git a/src/TrackGenerator3D.cpp b/src/TrackGenerator3D.cpp index e851e7ce2..1b56d891b 100644 --- a/src/TrackGenerator3D.cpp +++ b/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)); + } } } diff --git a/src/TraverseSegments.cpp b/src/TraverseSegments.cpp index 278451a51..d0705a103 100644 --- a/src/TraverseSegments.cpp +++ b/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;