Yes on 1 and 2. I think 3 should always be beta * lambda length, where lambda is the wavelength of the RF. If there is acceleration (RFQ), then the grid length will have to change from cell to cell.

EVoverC, EV_over_C?

@gsell Ok, that is "good" behaviour for the SAAMG/ArbitraryDomain where P0 is inside. Then the distance |I - P| will never be 0 in the stencil calculation.

The reference point should always be inside or always be outside, we should make either one the standard, otherwise we set ourselves up for more errors. I vote for P0 always inside (for the above-mentioned reason).

@gsell Quick question: If P lies directly on a boundary (i.e. |I - P| <= 1e-20 or whatever the number ends up being), is it counted as inside or outside? Asking because if it is counted as inside, I need to catch this special case during linear interpolation for the stencil creation, otherwise we have division by 0 there.

@gsell Were you able to reproduce the failing inside and intersection tests in the latest version (after merge)?

@gsell Happy to move, could you do it? As far as I know, I am using the exact same input files that I sent you. Did you try different mesh sizes?

But here is an OPAL-T .in file that I just used a few minutes ago :)

@gsell I pulled the latest changes and commented out my hardcoded beam pipe. Sadly, I still get a few points that are not correctly recognized as *outside*. In the following plot projecting points into the three planes, the upper row are all points that register as outside the pipe, the lower row are points registered as inside, which then consequently fail the boundary intersection tests. I manually checked, these red dots do not appear in the blue set, it just looks like it because of the regular grid and because I'm projecting them into planes.

@gsell are all changes merged into this branch?

With ArbitraryDomain, it should be the Dirichlet boundaries that are determined by the h5 file. In OPAL-T we could discuss whether to put in a flag that makes them open, but it should be the user's choice. In OPAL-cycl, I don't see a use case for open BC. We are only using it for the central region and there is no real "longitudinal direction", gamma is set to 1 (relativity = off) and the whole complicated geometry is meshed up.

Good point, I will look at that when the solver actually works again ^^

@frey_m I do use a beam pipe from h5 file, but I hardcoded the intersection points in the compute() function to circumvent ray boundary intersection tests from the geometry class. (The h5 geometry is still used in other places of the code, like finding the limits of the mesh) Here is the snippet:

```
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);
setHr(hr);
int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr_m[2] - 1) ? 0 : 1;
int yGhostOffsetLeft = (localId[1].first()== 0) ? 0 : 1;
int yGhostOffsetRight = (localId[1].last() == nr_m[1] - 1) ? 0 : 1;
int xGhostOffsetLeft = (localId[0].first()== 0) ? 0 : 1;
int xGhostOffsetRight = (localId[0].last() == nr_m[0] - 1) ? 0 : 1;
hasGeometryChanged_m = true;
intersectLoX_m.clear();
intersectHiX_m.clear();
intersectLoY_m.clear();
intersectHiY_m.clear();
intersectLoZ_m.clear();
intersectHiZ_m.clear();
// Calculate intersection
Vector_t P, dir, I;
Vector_t P0 = globalInsideP0_m;
isInsideMap_m.clear();
// TODO -DW 2020: I am hardcoding the beam pipe for now, to see if finding the correct intersection points
// is the issue.
double rp_temp = 0.007; // pipe radius in m
double lp_temp = 0.024; // pipe length in m
double r_temp_sq = 0.0; // P in cylinder coordinate r
double x_inters_temp = 0.0;
double y_inters_temp = 0.0;
// TODO -DW 2020: Remove the following debug output
*gmsg << "*** Some Domain shape/size DEBUG output in ArbitraryDomain.cpp: " << endl;
*gmsg << "get<i>RangeMin(): " << getXRangeMin() << ", " << getYRangeMin() << ", " << getZRangeMin() << endl;
*gmsg << "get<i>RangeMax(): " << getXRangeMax() << ", " << getYRangeMax() << ", " << getZRangeMax() << endl;
*gmsg << "Minima in Loop: " << getXRangeMin() + (localId[0].first() + 0.5) * hr[0] << ", "
<< getYRangeMin() + (localId[1].first() + 0.5) * hr[1] << ", "
<< getZRangeMin() + (localId[2].first() + 0.5) * hr[2] << endl;
*gmsg << "Maxima in Loop: " << getXRangeMin() + (localId[0].last() + 0.5) * hr[0] << ", "
<< getYRangeMin() + (localId[1].last() + 0.5) * hr[1] << ", "
<< getZRangeMin() + (localId[2].last() + 0.5) * hr[2] << endl;
*gmsg << "Indices: " << localId[0].first() << ", "
<< localId[1].first() << ", "
<< localId[2].first() << ", "
<< localId[0].last() << ", "
<< localId[1].last() << ", "
<< localId[2].last() << endl;
*gmsg << "*** " << endl;
// END TODO -DW 2020
// Note that the ArbitraryDomain, by definition, has no symmetries -DW 2020
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
P[2] = getZRangeMin() + (idz + 0.5) * hr[2];
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
P[1] = getYRangeMin() + (idy + 0.5) * hr[1];
for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
P[0] = getXRangeMin() + (idx + 0.5) * hr[0];
r_temp_sq = P[0] * P[0] + P[1] * P[1]; // TODO -DW 2020 to be removed
if (r_temp_sq < rp_temp * rp_temp && 0.0 < P[2] && P[2] < lp_temp) { // TODO -DW 2020: to be removed
// if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
// Fill the map with true or false values for very fast isInside tests
// during the rest of the fieldsolve.
isInsideMap_m[toCoordIdx(idx, idy, idz)] = true;
// Replace the old reference point with the new point (which we know is
// inside because we just tested for it. This makes the algorithm faster
// because fastIsInside() creates a number of segments that depends on the
// distance between P and P0. Using the previous P as the new P0
// assures the smallest possible distance in most cases. -DW
//P0 = P;
std::tuple<int, int, int> pos(idx, idy, idz);
// TODO -DW 2020: Analytically calculate the intersection points (to be removed)
// z coordinate is easy
intersectLoZ_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, 0.0));
intersectHiZ_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, lp_temp));
// x and y have to be calculated for the respective position
// x_intersections
x_inters_temp = sqrt(rp_temp*rp_temp - P[1]*P[1]);
intersectLoX_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, -x_inters_temp));
intersectHiX_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, x_inters_temp));
// y
y_inters_temp = sqrt(rp_temp*rp_temp - P[0]*P[0]);
intersectLoY_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, -y_inters_temp));
intersectHiY_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, y_inters_temp));
...
```

Another issue that I now see is that there is a longitudinal offset in the charge density. In this image, I create a flattop beam centered at 0.0 with length 20 mm. Clearly rho it is not centered in z. As far as I can see, the saving of rho in the lab frame works correctly, so it has to be the scatter that is off, but maybe you could double-check @frey_m?

@gsell Ok. in the meantime, I have analyzed the Problems with ArbitraryDomain a little more. This should also be interesting for @adelmann and @frey_m. There are four major issues:

- The aforementioned failures to find inside points and intersections. @gsell is working on that. Everything I do is with a hardcoded cylindrical beam pipe for now. Until the geometry class is fixed, we can't activate the ArbitraryDomain.
- Relativity. The mesh size in z direction is scaled with gammaz, which then leads to incorrect comparison of the bunch rest frame coordinates with the geometry, which is still in the lab frame. I am unsure how to really handle this. In OPAL-cycl, in the SPIRAL mode, I set gamma to 1 and force the whole calculation to be done non-relativistically. In OPAL-T, we could do the same, or we could pass the scaling factor from PartBunch to the MGSolver to the ArbitraryDomain and do the intersection checks in the bunch rest frame. Or the geometry class could handle Lorentz-contractions by passing gammaz into the intersection test?
- Linear interpolation isn't working atm. Everything I'm talking about so far is done in constant interpolation. This is likely a matter of building up the maps wrong. I will keep working on this.
- Parallelization. I haven't touched this one at all yet, but looking over the code, I think it shouldn't be too hard to extend it from Z to three dimensions.

@gsell Thanks for working on it! When you say parts of the code, you mean the geometry class or also the fieldsolver?

Apologies, I just ran a test case (different from the one where I noticed the problem) and opal 2.0, 2.2 and 2.4 are producing identical (correct looking) results. Weird. Let's not close the issue just yet, though, I'll do some more testing.

h5hut field maps (.h5part), loaded as part of the BANDRF cyclotron type that produce the desired effect in OPAL 2.0, don't seem to have any effect in OPAL 2.2 and up. Cyclotron units used to be mm and kV/mm (same as MV/m). Have these input units been changed somehow?

My current example is that of an electrostatic extraction septum that correctly pushes the final turn out by ~2 cm in OPAL 2.0, but does nothing in OPAL 2.2 and up.

This one should reproduce it, I think. You have to #define DEBUG_INTERSECT_RAY_BOUNDARY to see the debug output in ArbitraryDomain.cpp. Any line in the .o file (or stdout) that has one of "xdir=+1", "xdir=-1", "ydir=+1", "ydir=-1", "zdir=+1", "zdir=-1" in it basically means a mesh cell center was identified as inside point, but the line intersection test failed.