Method IntersectWithLine

Common/DataModel/vtkHyperTreeGridGeometricLocator.cxx:184–303

------------------------------------------------------------------------------

//------------------------------------------------------------------------------
int vtkHyperTreeGridGeometricLocator::IntersectWithLine(const double p0[3], const double p1[3],
  double tol, double& t, double x[3], double pcoords[3], int& subId, vtkIdType& cellId,
  vtkGenericCell* cell)
{
  // initialize outputs
  cellId = -1;
  t = -1.0;
  subId = 0;
  std::fill(x, x + 3, 0.0);
  std::fill(pcoords, pcoords + 3, 0.0);

  // setup calculation
  unsigned int dim = this->HTG->GetDimension();
  std::vector<double> sizes(dim, 0.0);
  std::vector<double> origin(dim, 0.0);
  this->GetZeroLevelOriginAndSize(origin.data(), sizes.data());

  // find intersection point with entire grid
  bool p0InCell = true;
  for (unsigned int d = 0; d < dim; d++)
  {
    p0InCell &= ((p0[d] - origin[d]) < (sizes[d] + tol));
  }

  if (!p0InCell)
  {
    if (!this->ConstructCell(origin.data(), sizes.data(), cell))
    {
      vtkErrorMacro("Could not construct cell");
      return -1;
    }
    // line does not intersect grid at all
    if (cell->IntersectWithLine(p0, p1, tol, t, x, pcoords, subId) == 0)
    {
      return 0;
    }
    // run FindCell on the intersection point + epsilon so we're sure we're in the cell
    std::vector<double> tangent(3, 0.0);
    vtkMath::Subtract(p1, p0, tangent.data());
    vtkMath::Normalize(tangent.data());
    double epsilon = 0.01 *
      (vtkMath::Norm(sizes.data()) /
        std::pow(this->HTG->GetBranchFactor(), this->HTG->GetNumberOfLevels()));
    epsilon = std::max(epsilon, tol * 2.0);
    vtkMath::MultiplyScalar(tangent.data(), epsilon);
    vtkMath::Add(x, tangent.data(), x);
  }
  else
  {
    std::copy(p0, p0 + dim, x);
  }
  {
    std::vector<double> locWeights(std::pow(2, dim), 0.0);
    std::vector<double> locPCoords(3, 0.0);
    // potential speed-up here by automatically looking for the cell but also could be loss
    cellId = this->FindCell(x, tol, cell, subId, locPCoords.data(), locWeights.data());
  }
  if (cellId >= 0)