Commit 10354343 authored by kraus's avatar kraus

Merge branch '573-broken-tests-due-to-509' into 'master'

Resolve 'Broken tests due to #509'

Closes #573

See merge request !401
parents 87043626 b84a3eb0
//
// Class IndexMap
//
// This class stores and prints the sequence of elements that the referenc particle passes.
// Each time the reference particle enters or leaves an element an entry is added to the map.
// With help of this map one can determine which element can be found at a given position.
//
// Copyright (c) 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include <map>
#include <limits>
#include <iostream>
......@@ -128,8 +150,7 @@ void IndexMap::tidyUp(double zstop) {
if (rit != mapRange2Element_m.rend() &&
(*rit).second.size() == 0 &&
zstop > (*rit).first.first &&
zstop < (*rit).first.second) {
zstop > (*rit).first.first) {
key_t key((*rit).first.first, zstop);
value_t val;
......
//
// Class IndexMap
//
// This class stores and prints the sequence of elements that the referenc particle passes.
// Each time the reference particle enters or leaves an element an entry is added to the map.
// With help of this map one can determine which element can be found at a given position.
//
// Copyright (c) 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef OPAL_INDEXMAP_H
#define OPAL_INDEXMAP_H
......
......@@ -112,6 +112,7 @@ void OrbitThreader::execute() {
std::set<std::string> visitedElements;
trackBack();
updateBoundingBoxWithCurrentPosition();
Vector_t nextR = r_m / (Physics::c * dt_m);
integrator_m.push(nextR, p_m, dt_m);
......@@ -414,6 +415,7 @@ void OrbitThreader::computeBoundingBox() {
FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
FieldList::iterator it = allElements.begin();
const FieldList::iterator end = allElements.end();
globalBoundingBox_m.lowerLeftCorner = Vector_t(std::numeric_limits<double>::max());
globalBoundingBox_m.upperRightCorner = Vector_t(-std::numeric_limits<double>::max());
......@@ -424,6 +426,15 @@ void OrbitThreader::computeBoundingBox() {
ElementBase::BoundingBox other = it->getBoundingBoxInLabCoords();
globalBoundingBox_m.getCombinedBoundingBox(other);
}
updateBoundingBoxWithCurrentPosition();
}
void OrbitThreader::updateBoundingBoxWithCurrentPosition() {
Vector_t dR = Physics::c * dt_m * p_m / Util::getGamma(p_m);
ElementBase::BoundingBox pos = ElementBase::BoundingBox::getBoundingBox({r_m - 10 * dR, r_m + 10 * dR});
globalBoundingBox_m.getCombinedBoundingBox(pos);
}
double OrbitThreader::computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
......
......@@ -108,7 +108,7 @@ private:
void processElementRegister();
void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements);
void computeBoundingBox();
// double computeMaximalImplicitDrift();
void updateBoundingBoxWithCurrentPosition();
double computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
const Vector_t & position,
const Vector_t & direction) const;
......
......@@ -309,6 +309,17 @@ bool ElementBase::isInsideTransverse(const Vector_t &r) const
}
}
ElementBase::BoundingBox ElementBase::BoundingBox::getBoundingBox(const std::vector<Vector_t> & points) {
const Vector_t & point = points.front();
BoundingBox result{point, point};
for (const Vector_t & point: points) {
BoundingBox tmp{point, point};
result.getCombinedBoundingBox(tmp);
}
return result;
}
bool ElementBase::BoundingBox::isInside(const Vector_t & position) const {
Vector_t relativePosition = position - lowerLeftCorner;
Vector_t diagonal = upperRightCorner - lowerLeftCorner;
......@@ -392,4 +403,26 @@ ElementBase::BoundingBox ElementBase::getBoundingBoxInLabCoords() const {
}
return bb;
}
void ElementBase::BoundingBox::print(std::ostream & out) const {
const Vector_t & ll = lowerLeftCorner;
const Vector_t & ur = upperRightCorner;
Vector_t diagonal = ur - ll;
Vector_t dX(diagonal(0), 0, 0), dY(0, diagonal(1), 0), dZ(0, 0, diagonal(2));
std::vector<Vector_t> corners{ll, ll + dX, ll + dX + dY, ll + dY, ur, ur - dX, ur - dX - dY, ur - dY};
std::vector<std::vector<unsigned int>> paths{{0, 1, 2, 3}, {0, 1, 7, 6}, {1, 2, 4, 7}, {2, 3, 5, 4}, {3, 0, 6, 5}, {4, 5, 6, 7}};
out << std::setprecision(8);
for (const std::vector<unsigned int> path: paths) {
for (unsigned int i : {0, 1, 2, 3, 0}) {
const Vector_t & corner = corners[path[i]];
out << std::setw(16) << corner(0)
<< std::setw(16) << corner(1)
<< std::setw(16) << corner(2)
<< std::endl;
}
out << std::endl;
}
}
\ No newline at end of file
......@@ -345,6 +345,8 @@ public:
Vector_t lowerLeftCorner;
Vector_t upperRightCorner;
static BoundingBox getBoundingBox(const std::vector<Vector_t> & points);
void getCombinedBoundingBox(const BoundingBox & other) {
for (unsigned int d = 0; d < 3; ++ d) {
lowerLeftCorner[d] = std::min(lowerLeftCorner[d], other.lowerLeftCorner[d]);
......@@ -354,6 +356,8 @@ public:
bool isInside(const Vector_t &) const;
void print(std::ostream &) const;
/*! Computes the intersection point between a bounding box and the ray which
* has the direction 'direction' and starts at the position 'position'. If
* the position is inside the box then the algorithm should find an inter-
......@@ -609,4 +613,4 @@ inline
bool ElementBase::isElementPositionSet() const
{ return elemedgeSet_m; }
#endif // CLASSIC_ElementBase_HH
#endif // CLASSIC_ElementBase_HH
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment