Merge branch 'fix-301-premature-termination-of-integration' into 'master'

Fix for problems caused by #301

See merge request !95

src/Algorithms/IndexMap.cpp

@@ -117,11 +117,19 @@ void IndexMap::add(key_t::first_type initialS, key_t::second_type finalS, const
     }
 }
 
-void IndexMap::tidyUp() {
+void IndexMap::tidyUp(double zstop) {
     map_t::reverse_iterator rit = mapRange2Element_m.rbegin();
 
-    if (rit != mapRange2Element_m.rend() && (*rit).second.size() == 0) {
+    if (rit != mapRange2Element_m.rend() &&
+        (*rit).second.size() == 0 &&
+        zstop > (*rit).first.first &&
+        zstop < (*rit).first.second) {
+
+        key_t key((*rit).first.first, zstop);
+        value_t val;
+
         mapRange2Element_m.erase(std::next(rit).base());
+        mapRange2Element_m.insert(std::pair<key_t, value_t>(key, val));
     }
 }
 
@@ -135,17 +143,18 @@ enum elements {
     SOLENOID,
     RFCAVITY,
     MONITOR,
+    OTHER,
     SIZE
 };
 
-void IndexMap::saveSDDS(double startS) const {
+void IndexMap::saveSDDS(double initialPathLength) const {
     auto opal = OpalData::getInstance();
     if (opal->isOptimizerRun()) return;
 
     std::string fileName("data/" + OpalData::getInstance()->getInputBasename() + "_ElementPositions.sdds");
     std::ofstream sdds;
     if (OpalData::getInstance()->hasPriorTrack() && boost::filesystem::exists(fileName)) {
-        Util::rewindLinesSDDS(fileName, startS, false);
+        Util::rewindLinesSDDS(fileName, initialPathLength, false);
         sdds.open(fileName, std::ios::app);
     } else {
         sdds.open(fileName);
@@ -222,6 +231,12 @@ void IndexMap::saveSDDS(double startS) const {
              << indent << "units=1, \n"
              << indent << "description=\"10 monitor present\" \n"
              << "&end\n";
+        sdds << "&column \n"
+             << indent << "name=other, \n"
+             << indent << "type=float, \n"
+             << indent << "units=1, \n"
+             << indent << "description=\"11 other element present\" \n"
+             << "&end\n";
         sdds << "&column \n"
              << indent << "name=element_names, \n"
              << indent << "type=string, \n"
@@ -238,7 +253,7 @@ void IndexMap::saveSDDS(double startS) const {
     std::vector<std::vector<int> > allItems(SIZE);
     std::vector<double> allPositions;
     std::vector<std::string> allNames;
-    std::vector<double> typeMultipliers = {3.3333e-1, 1.0, 0.5, 0.25, 1.0, 1.0, 1.0, 1.0, 1.0};
+    std::vector<double> typeMultipliers = {3.3333e-1, 1.0, 0.5, 0.25, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
     unsigned int counter = 0;
 
     auto mapIti = mapRange2Element_m.begin();
@@ -298,7 +313,8 @@ void IndexMap::saveSDDS(double startS) const {
                 items[MONITOR] = 1;
                 break;
             default:
-                continue;
+                items[OTHER] = 1;
+                break;
             }
 
             names += element->getName() + ", ";
@@ -328,6 +344,26 @@ void IndexMap::saveSDDS(double startS) const {
 
     if (counter == 0) return;
 
+    if (allPositions.front() > initialPathLength) {
+        {
+            auto tmp = allPositions;
+            allPositions = std::vector<double>(4, initialPathLength);
+            allPositions.insert(allPositions.end(), tmp.begin(), tmp.end());
+        }
+
+        {
+            auto tmp = allNames;
+            allNames = std::vector<std::string>(4, "");
+            allNames.insert(allNames.end(), tmp.begin(), tmp.end());
+        }
+
+        for (unsigned int i = 0; i < SIZE; ++ i) {
+            auto tmp = allItems[i];
+            allItems[i] = std::vector<int>(4, 0);
+            allItems[i].insert(allItems[i].end(), tmp.begin(), tmp.end());
+        }
+    }
+
     const unsigned int totalSize = counter;
     for (unsigned int i = 0; i < SIZE; ++ i) {
         for (unsigned int j = 0; j < totalSize - 1; ++ j) {

src/Algorithms/IndexMap.h

@@ -23,7 +23,7 @@ public:
 
     value_t query(key_t::first_type s, key_t::second_type ds);
 
-    void tidyUp();
+    void tidyUp(double zstop);
 
     void print(std::ostream&) const;
     void saveSDDS(double startS) const;
@@ -92,4 +92,4 @@ Inform& operator<< (Inform &out, const IndexMap &im) {
     return out;
 }
 
-#endif
+#endif
\ No newline at end of file

src/Algorithms/OrbitThreader.cpp

@@ -4,6 +4,7 @@
 #include "AbsBeamline/RFCavity.h"
 #include "AbsBeamline/BendBase.h"
 #include "AbsBeamline/TravelingWave.h"
+#include "BeamlineCore/MarkerRep.h"
 
 #include "AbstractObjects/OpalData.h"
 #include "BasicActions/Option.h"
@@ -132,7 +133,7 @@ void OrbitThreader::execute() {
     } while (errorFlag_m != HITMATERIAL &&
              errorFlag_m != EOL);
 
-    imap_m.tidyUp();
+    imap_m.tidyUp(zstop_m);
     *gmsg << level1 << "\n" << imap_m << endl;
     imap_m.saveSDDS(initialPathLength);
 
@@ -387,6 +388,14 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
 double OrbitThreader::computeMaximalImplicitDrift() {
     FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
     double maxDrift = 0.0;
+
+    MarkerRep start("#S");
+    CoordinateSystemTrafo toEdge(r_m, getQuaternion(p_m, Vector_t(0, 0, 1)));
+    start.setElementLength(0.0);
+    start.setCSTrafoGlobal2Local(toEdge);
+    std::shared_ptr<Component> startPtr(static_cast<Marker*>(start.clone()));
+    allElements.push_front(ClassicField(startPtr, 0.0, 0.0));
+
     FieldList::iterator it = allElements.begin();
     const FieldList::iterator end = allElements.end();
 
@@ -414,7 +423,9 @@ double OrbitThreader::computeMaximalImplicitDrift() {
             auto element2 = it2->getElement();
             const auto &toEdge = element2->getCSTrafoGlobal2Local();
             auto toBegin = element2->getEdgeToBegin() * toEdge;
+            auto toEnd = element2->getEdgeToEnd() * toEdge;
             Vector_t begin2 = toBegin.transformFrom(Vector_t(0, 0, 0));
+            Vector_t end2 = toEnd.transformFrom(Vector_t(0, 0, 0));
             Vector_t directionBegin = toBegin.rotateFrom(Vector_t(0, 0, 1));
             if (element2->getType() == ElementBase::RBEND ||
                 element2->getType() == ElementBase::SBEND ||
@@ -427,15 +438,22 @@ double OrbitThreader::computeMaximalImplicitDrift() {
 
             double distance = euclidean_norm(begin2 - end1);
             double directionProjection = dot(directionEnd, directionBegin);
-            if (directionProjection > 0.999 && minDistanceLocal > distance) {
+            bool overlapping = dot(begin2 - end1, directionBegin) < 0.0? true: false;
+
+            if (!overlapping &&
+                directionProjection > 0.999 &&
+                minDistanceLocal > distance) {
                 minDistanceLocal = distance;
             }
         }
 
-        if (maxDrift < minDistanceLocal) {
+        if (maxDrift < minDistanceLocal &&
+            minDistanceLocal != std::numeric_limits<double>::max()) {
             maxDrift = minDistanceLocal;
         }
     }
 
+    maxDrift = std::min(maxIntegSteps_m * dt_m * Physics::c, maxDrift);
+
     return maxDrift;
 }
\ No newline at end of file