Merge branch '585-opal-t-fails-to-write-lossdatasink-if-executed-in-parallel-2' into 'master'

Resolve "Opal-T fails to write LossDataSink if executed in parallel"

Closes #585

See merge request !439

src/Algorithms/IndexMap.cpp

@@ -55,16 +55,16 @@ void IndexMap::print(std::ostream &out) const {
     auto mapIti = mapRange2Element_m.begin();
     auto mapItf = mapRange2Element_m.end();
 
-    double totalLength = (*mapRange2Element_m.rbegin()).first.second;
+    double totalLength = (*mapRange2Element_m.rbegin()).first.end;
     unsigned int numDigits = std::floor(std::max(0.0, log(totalLength) / log(10.0))) + 1;
 
     for (; mapIti != mapItf; mapIti++) {
-        const auto key = (*mapIti).first;
-        const auto val = (*mapIti).second;
+        const key_t key = (*mapIti).first;
+        const value_t val = (*mapIti).second;
         out << "Key: ("
-            << std::setw(numDigits + 7) << std::right << key.first
+            << std::setw(numDigits + 7) << std::right << key.begin
             << " - "
-            << std::setw(numDigits + 7) << std::right << key.second
+            << std::setw(numDigits + 7) << std::right << key.end
             << ") number of overlapping elements " << val.size() << "\n";
 
         for (auto element: val) {
@@ -79,7 +79,7 @@ IndexMap::value_t IndexMap::query(key_t::first_type s, key_t::second_type ds) {
     value_t elementSet;
 
     map_t::reverse_iterator rit = mapRange2Element_m.rbegin();
-    if (rit != mapRange2Element_m.rend() && lowerLimit > (*rit).first.second) {
+    if (rit != mapRange2Element_m.rend() && lowerLimit > (*rit).first.end) {
         throw OutOfBounds("IndexMap::query", "out of bounds");
     }
 
@@ -87,8 +87,8 @@ IndexMap::value_t IndexMap::query(key_t::first_type s, key_t::second_type ds) {
     const map_t::iterator end = mapRange2Element_m.end();
 
     for (; it != end; ++ it) {
-        const double low = (*it).first.first;
-        const double high = (*it).first.second;
+        const double low = (*it).first.begin;
+        const double high = (*it).first.end;
 
         if (lowerLimit < high && upperLimit >= low) break;
     }
@@ -97,7 +97,7 @@ IndexMap::value_t IndexMap::query(key_t::first_type s, key_t::second_type ds) {
 
     map_t::iterator last = std::next(it);
     for (; last != end; ++ last) {
-        const double low = (*last).first.first;
+        const double low = (*last).first.begin;
 
         if (upperLimit < low) break;
     }
@@ -114,10 +114,10 @@ void IndexMap::add(key_t::first_type initialS, key_t::second_type finalS, const
     if (initialS > finalS) {
         std::swap(initialS, finalS);
     }
-    key_t key(initialS, finalS * oneMinusEpsilon_m);
+    key_t key{initialS, finalS * oneMinusEpsilon_m};
 
     mapRange2Element_m.insert(std::pair<key_t, value_t>(key, val));
-    totalPathLength_m = (*mapRange2Element_m.rbegin()).first.second;
+    totalPathLength_m = (*mapRange2Element_m.rbegin()).first.end;
 
     value_t::iterator setIt = val.begin();
     const value_t::iterator setEnd = val.end();
@@ -132,8 +132,8 @@ void IndexMap::add(key_t::first_type initialS, key_t::second_type finalS, const
             for (auto it = itpair.first; it != itpair.second; ++ it) {
                 key_t &currentRange = it->second;
 
-                if (almostEqual(key.first, currentRange.second / oneMinusEpsilon_m)) {
-                    currentRange.second = key.second;
+                if (almostEqual(key.begin, currentRange.end / oneMinusEpsilon_m)) {
+                    currentRange.end = key.end;
                     extendedExisting = true;
                     break;
                 }
@@ -150,9 +150,9 @@ void IndexMap::tidyUp(double zstop) {
 
     if (rit != mapRange2Element_m.rend() &&
         (*rit).second.size() == 0 &&
-        zstop > (*rit).first.first) {
+        zstop > (*rit).first.begin) {
 
-        key_t key((*rit).first.first, zstop);
+        key_t key{(*rit).first.begin, zstop};
         value_t val;
 
         mapRange2Element_m.erase(std::next(rit).base());
@@ -196,8 +196,8 @@ void IndexMap::saveSDDS(double initialPathLength) const {
 
         const auto &sectorRange = (*mapIti).first;
 
-        double sectorBegin = sectorRange.first;
-        double sectorEnd = sectorRange.second;
+        double sectorBegin = sectorRange.begin;
+        double sectorEnd = sectorRange.end;
 
         std::vector<std::tuple<double, std::vector<double>, std::string> > currentSector(4);
         std::get<0>(currentSector[0]) = sectorBegin;
@@ -216,8 +216,8 @@ void IndexMap::saveSDDS(double initialPathLength) const {
             auto end = elementPassages.second;
             for (; passage != end; ++ passage) {
                 const auto &elementRange = (*passage).second;
-                double elementBegin = elementRange.first;
-                double elementEnd = elementRange.second;
+                double elementBegin = elementRange.begin;
+                double elementEnd = elementRange.end;
 
                 if (elementBegin <= sectorBegin &&
                     elementEnd >= sectorEnd) {
@@ -226,7 +226,7 @@ void IndexMap::saveSDDS(double initialPathLength) const {
             }
 
             const auto &elementRange = (*passage).second;
-            if (elementRange.first < sectorBegin) {
+            if (elementRange.begin < sectorBegin) {
                 ::insertFlags(std::get<1>(currentSector[0]), element);
                 std::get<2>(currentSector[0]) += element->getName() + ", ";
             }
@@ -237,7 +237,7 @@ void IndexMap::saveSDDS(double initialPathLength) const {
             ::insertFlags(std::get<1>(currentSector[2]), element);
             std::get<2>(currentSector[2]) += element->getName() + ", ";
 
-            if (elementRange.second > sectorEnd) {
+            if (elementRange.end > sectorEnd) {
                 ::insertFlags(std::get<1>(currentSector[3]), element);
                 std::get<2>(currentSector[3]) += element->getName() + ", ";
             }
@@ -265,7 +265,7 @@ void IndexMap::saveSDDS(double initialPathLength) const {
 
         unsigned int i = 0;
         for (; i < numEntries; ++ i) {
-            if (std::get<0>(sectors[i]) >= range.first) {
+            if (std::get<0>(sectors[i]) >= range.begin) {
                 break;
             }
         }
@@ -273,16 +273,16 @@ void IndexMap::saveSDDS(double initialPathLength) const {
         if (i == numEntries) continue;
 
         unsigned int j = ++ i;
-        while (std::get<0>(sectors[j]) < range.second) {
+        while (std::get<0>(sectors[j]) < range.end) {
             ++ j;
         }
 
-        double length = range.second - range.first;
+        double length = range.end - range.begin;
         for (; i <= j; ++ i) {
             double pos = std::get<0>(sectors[i]);
             auto &items = std::get<1>(sectors[i]);
 
-            items[RFCAVITY] = 1.0 - 2 * (pos - range.first) / length;
+            items[RFCAVITY] = 1.0 - 2 * (pos - range.begin) / length;
         }
     }
 
@@ -372,18 +372,18 @@ namespace {
     }
 }
 
-std::pair<double, double> IndexMap::getRange(const IndexMap::value_t::value_type &element,
-                                             double position) const {
+IndexMap::key_t IndexMap::getRange(const IndexMap::value_t::value_type &element,
+                                   double position) const {
     double minDistance = std::numeric_limits<double>::max();
-    std::pair<double, double> range(0.0, 0.0);
+    key_t range{0.0, 0.0};
     const std::pair<invertedMap_t::const_iterator, invertedMap_t::const_iterator> its = mapElement2Range_m.equal_range(element);
     if (std::distance(its.first, its.second) == 0)
         throw OpalException("IndexMap::getRange()",
                             "Element \"" + element->getName() + "\" not registered");
 
     for (invertedMap_t::const_iterator it = its.first; it != its.second; ++ it) {
-        double distance = std::min(std::abs((*it).second.first - position),
-                                   std::abs((*it).second.second - position));
+        double distance = std::min(std::abs((*it).second.begin - position),
+                                   std::abs((*it).second.end - position));
         if (distance < minDistance) {
             minDistance = distance;
             range = (*it).second;
@@ -393,14 +393,14 @@ std::pair<double, double> IndexMap::getRange(const IndexMap::value_t::value_type
     return range;
 }
 
-IndexMap::value_t IndexMap::getTouchingElements(const std::pair<double, double> &range) const {
+IndexMap::value_t IndexMap::getTouchingElements(const IndexMap::key_t &range) const {
     map_t::const_iterator it = mapRange2Element_m.begin();
     const map_t::const_iterator end = mapRange2Element_m.end();
-    IndexMap::value_t touchingElements;
+    value_t touchingElements;
 
     for (; it != end; ++ it) {
-        if (almostEqual(it->first.first, range.second) ||
-            almostEqual(it->first.second, range.first))
+        if (almostEqual(it->first.begin, range.begin) ||
+            almostEqual(it->first.end, range.end))
             touchingElements.insert((it->second).begin(), (it->second).end());
     }
 

src/Algorithms/IndexMap.h

@@ -36,7 +36,14 @@
 class IndexMap
 {
 public:
-    typedef std::pair<double, double> key_t;
+    struct Range
+    {
+        typedef double first_type;
+        typedef double second_type;
+        first_type begin;
+        second_type end;
+    };
+    typedef Range key_t;
     typedef std::set<std::shared_ptr<Component> > value_t;
 
     IndexMap();
@@ -52,9 +59,8 @@ public:
     size_t size() const;
 
     size_t numElements() const;
-    std::pair<double, double> getRange(const IndexMap::value_t::value_type &element,
-                                       double position) const;
-    IndexMap::value_t getTouchingElements(const std::pair<double, double> &range) const;
+    key_t getRange(const IndexMap::value_t::value_type &element, double position) const;
+    value_t getTouchingElements(const key_t &range) const;
 
     class OutOfBounds: public OpalException {
     public:
@@ -75,10 +81,10 @@ private:
     public:
         bool operator()(const key_t x , const key_t y) const
         {
-            if (x.first < y.first) return true;
+            if (x.begin < y.begin) return true;
 
-            if (x.first == y.first) {
-                if (x.second < y.second) return true;
+            if (x.begin == y.begin) {
+                if (x.end < y.end) return true;
             }
 
             return false;

src/Algorithms/OrbitThreader.h

@@ -50,9 +50,9 @@ public:
     IndexMap::value_t query(IndexMap::key_t::first_type step,
                             IndexMap::key_t::second_type length);
 
-    std::pair<double, double> getRange(const IndexMap::value_t::value_type &element,
-                                       double position) const;
-    IndexMap::value_t getTouchingElements(const std::pair<double, double> &range) const;
+    IndexMap::key_t getRange(const IndexMap::value_t::value_type &element,
+                             double position) const;
+    IndexMap::value_t getTouchingElements(const IndexMap::key_t &range) const;
 
 private:
     /// position of reference particle in lab coordinates
@@ -125,13 +125,13 @@ IndexMap::value_t OrbitThreader::query(IndexMap::key_t::first_type pathLength,
 }
 
 inline
-std::pair<double, double> OrbitThreader::getRange(const IndexMap::value_t::value_type &element,
-                                                  double position) const {
+IndexMap::key_t OrbitThreader::getRange(const IndexMap::value_t::value_type &element,
+                                        double position) const {
     return imap_m.getRange(element, position);
 }
 
 inline
-IndexMap::value_t OrbitThreader::getTouchingElements(const std::pair<double, double> &range) const {
+IndexMap::value_t OrbitThreader::getTouchingElements(const IndexMap::key_t &range) const {
     return imap_m.getTouchingElements(range);
 }
 

src/Algorithms/ParallelTTracker.cpp

@@ -674,7 +674,7 @@ void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elemen
     Inform msg("ParallelTTracker ", *gmsg);
     std::set<IndexMap::value_t::value_type> elementsWithParticleMatterInteraction;
     std::set<ParticleMatterInteractionHandler*> particleMatterinteractionHandlers;
-    std::pair<double, double> currentRange(0.0, 0.0);
+    IndexMap::key_t currentRange{0.0, 0.0};
 
     while (elements.size() > 0) {
         auto it = elements.begin();
@@ -682,9 +682,9 @@ void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elemen
             elementsWithParticleMatterInteraction.insert(*it);
             particleMatterinteractionHandlers.insert((*it)->getParticleMatterInteraction());
 
-            std::pair<double, double> range = oth.getRange(*it, pathLength_m);
-            currentRange.first = std::min(currentRange.first, range.first);
-            currentRange.second = std::max(currentRange.second, range.second);
+            IndexMap::key_t range = oth.getRange(*it, pathLength_m);
+            currentRange.begin = std::min(currentRange.begin, range.begin);
+            currentRange.end = std::max(currentRange.end, range.end);
 
             IndexMap::value_t touching = oth.getTouchingElements(range);
             elements.insert(touching.begin(), touching.end());

src/Classic/Utilities/ClassicField.h

@@ -24,7 +24,8 @@ public:
     void setOff();
 
     static bool SortAsc(const ClassicField &fle1, const ClassicField &fle2) {
-        return (fle1.start_m < fle2.start_m);
+        return (fle1.start_m < fle2.start_m
+                || (fle1.start_m == fle2.start_m && fle1.element_m->getName() < fle2.element_m->getName()));
     }
 
     static bool ZeroLength(const ClassicField &fle) {
@@ -103,4 +104,4 @@ inline
 ElementBase::BoundingBox ClassicField::getBoundingBoxInLabCoords() const {
     return element_m->getBoundingBoxInLabCoords();
 }
-#endif // CLASSIC_FIELD_H
+#endif // CLASSIC_FIELD_H
\ No newline at end of file

src/Elements/OpalBeamline.cpp

@@ -142,23 +142,6 @@ void OpalBeamline::switchElements(const double &min, const double &max, const do
     }
 }
 
-void OpalBeamline::switchElementsOff(const double &min, ElementBase::ElementType eltype) {
-    if(eltype == ElementBase::ANY) {
-        for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) {
-            if((*flit).isOn() && min >= (*flit).getEnd()) {
-                (*flit).setOff();
-            }
-
-        }
-    } else {
-        for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) {
-            if((*flit).isOn() && min >= (*flit).getEnd() && (*flit).getElement()->getType() == eltype) {
-                (*flit).setOff();
-            }
-        }
-    }
-}
-
 void OpalBeamline::switchElementsOff() {
     for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit)
         (*flit).setOff();
@@ -169,6 +152,7 @@ void OpalBeamline::prepareSections() {
         prepared_m = true;
         return;
     }
+    elements_m.sort(ClassicField::SortAsc);
     prepared_m = true;
 }
 
@@ -221,22 +205,11 @@ void OpalBeamline::compute3DLattice() {
     static unsigned int order = 0;
     const FieldList::iterator end = elements_m.end();
 
-
     unsigned int minOrder = order;
     {
         double endPriorPathLength = 0.0;
         CoordinateSystemTrafo currentCoordTrafo = coordTransformationTo_m;
 
-        elements_m.sort([](const ClassicField& a, const ClassicField& b) {
-                double edgeA = 0.0, edgeB = 0.0;
-                if (a.getElement()->isElementPositionSet())
-                    edgeA = a.getElement()->getElementPosition();
-
-                if (b.getElement()->isElementPositionSet())
-                    edgeB = b.getElement()->getElementPosition();
-
-                return edgeA < edgeB;
-            });
         FieldList::iterator it = elements_m.begin();
         for (; it != end; ++ it) {
             if ((*it).isPositioned()) {
@@ -376,8 +349,6 @@ void OpalBeamline::compute3DLattice() {
 
         (*it).fixPosition();
     }
-
-    elements_m.sort(ClassicField::SortAsc);
 }
 
 void OpalBeamline::save3DLattice() {

src/Elements/OpalBeamline.h

@@ -68,8 +68,6 @@ public:
     double getEnd(const Vector_t &) const;
 
     void switchElements(const double &, const double &, const double &kineticEnergy, const bool &nomonitors = false);
-
-    void switchElementsOff(const double &, ElementBase::ElementType eltype = ElementBase::ANY);
     void switchElementsOff();
 
     ParticleMatterInteractionHandler *getParticleMatterInteractionHandler(const unsigned int &);