Method EliminateXWrap

IO/NetCDF/vtkMPASReader.cxx:2049–2203

//------------------------------------------------------------------------------

int vtkMPASReader::EliminateXWrap()
{
  if (this->NumberOfPoints == 0)
  {
    return 1;
  }

  double xLength;
  double xCenter;
  switch (this->Geometry)
  {
    case vtkMPASReader::Spherical:
      vtkErrorMacro("EliminateXWrap called for spherical geometry.");
      return 0;

    case vtkMPASReader::Projected:
      xLength = 2 * vtkMath::Pi();
      xCenter = this->CenterRad;
      break;

    case vtkMPASReader::Planar:
    {
      // Determine the bounds in the x-dimension
      double xRange[2] = { this->PointX[this->PointOffset], this->PointX[this->PointOffset] };
      for (size_t i = 1; i < this->NumberOfPoints; ++i)
      {
        double x = this->PointX[this->PointOffset + i];
        xRange[0] = std::min(xRange[0], x);
        xRange[1] = std::max(xRange[1], x);
      }

      xLength = xRange[1] - xRange[0];
      xCenter = (xRange[0] + xRange[1]) * 0.5;
    }
    break;

    default:
      vtkErrorMacro("Unrecognized geometry type (" << this->Geometry << ").");
      return 0;
  }

  constexpr double tolerance = 5.5;

  // For each cell, examine vertices
  // Add new points and cells where needed to account for wraparound.
  for (size_t j = this->CellOffset; j < this->NumberOfCells + this->CellOffset; j++)
  {
    int* conns = this->OrigConnections + (j * this->PointsPerCell);
    int* modConns = this->ModConnections + (j * this->PointsPerCell);

    // Determine if we are wrapping in X direction
    size_t lastk = this->PointsPerCell - 1;
    bool xWrap = false;
    for (size_t k = 0; k < this->PointsPerCell; k++)
    {
      if (std::abs(this->PointX[conns[k]] - this->PointX[conns[lastk]]) > tolerance)
      {
        xWrap = true;