Bug report: particles not hitting surface for very simple case

[jguterl]

Setup: as simple as possible. 3 particles start from slightly above the surface at z=0 of a rectangular box. Bfield along z direction. No ionization, no recombination, no collision, no diff, no Efield. Only boris and check_geom. vx,vy=0, vz=-10000, x and y close to center of the box at x,y=0,0.

Problem: not all particles are collected at the bottom surface.
Input Files available here:
https://www.dropbox.com/sh/poi1vy6fbg3ldyd/AAApz5DE3mbKRY7Ys3E3wNAZa?dl=0

[jguterl]

To analyze output:

import numpy as np
FileName='/home/jguterl/Dropbox/python/pyGITR/examples/simple_cube/output/surface.nc'

import netCDF4
from netCDF4 import Dataset
SurfaceData = Dataset(FileName, "r", format="NETCDF4")
Distrib = np.array(SurfaceData.variables['surfEDist'])
print('Total particle={}'.format(np.sum(Distrib,(0,1,2))))

[jguterl]

Solution:
In geometryCheck.h,
replace

if (vectorNorm(crossABAp) == 0 || vectorNorm(crossBCBp) == 0 ||
              vectorNorm(crossCACp) == 0) {
            hitSurface = 1;
          }

by

if (vectorNorm(crossABAp) < 1e-13 || vectorNorm(crossBCBp) < 1e-13 ||
              vectorNorm(crossCACp) < 1e-13) {
            hitSurface = 1;
          }

[jguterl]

Faster and more robust algo to find intersection:

        if (signPoint0 != signPoint1) {
          t = -(a * p0[0] + b * p0[1] + c * p0[2] + d) /
              (a * (p1[0] - p0[0]) + b * (p1[1] - p0[1]) + c * (p1[2] - p0[2]));
          vectorAssign(boundaryVector[i].x1-p0[0] - t * (p1[0] - p0[0]), boundaryVector[i].y1-p0[1] - t * (p1[1] - p0[1]),
                       boundaryVector[i].z1- p0[2] - t * (p1[2] - p0[2]), A);
          vectorAssign(boundaryVector[i].x2-p0[0] + t * (p1[0] - p0[0]), boundaryVector[i].y2-p0[1] - t * (p1[1] - p0[1]),
                       boundaryVector[i].z2- p0[2] - t * (p1[2] - p0[2]), B);
          vectorAssign(boundaryVector[i].x3-p0[0] + t * (p1[0] - p0[0]), boundaryVector[i].y3-p0[1] - t * (p1[1] - p0[1]),
                       boundaryVector[i].z3- p0[2] - t * (p1[2] - p0[2]), C);
          vectorCrossProduct(B, C, u);
          vectorCrossProduct(C, A, v);
          vectorCrossProduct(A, B, w);

          if (vectorDotProduct(u,v) >= 0.0f && vectorDotProduct(u,w) >= 0.0f)
          {hitSurface = 1;}