XmUGrid.cpp

Go to the documentation of this file.

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

//----- Included files ---------------------------------------------------------

// 1. Precompiled header

// 2. My own header
#include <xmsgrid/ugrid/XmUGrid.h>

// 3. Standard library headers
#include <cmath>

// 4. External library headers
#include <boost/container/flat_set.hpp>

// 5. Shared code headers
#include <xmscore/misc/XmConst.h>
#include <xmscore/misc/XmError.h>
#include <xmscore/misc/XmLog.h>
#include <xmscore/misc/xmstype.h>
#include <xmscore/stl/set.h>

// 6. Non-shared code headers
#include <xmsgrid/geometry/geoms.h>
#include <xmsgrid/ugrid/detail/XmGeometry.h>
#include <xmsgrid/ugrid/XmEdge.h>

//----- Forward declarations ---------------------------------------------------

//----- External globals -------------------------------------------------------

//----- Namespace declaration --------------------------------------------------

//----- Constants / Enumerations -----------------------------------------------

//----- Classes / Structs ------------------------------------------------------

//----- Class / Function definitions -------------------------------------------

namespace xms
{
class XmUGrid::Impl
{
public:
  Impl() = default;
  Impl(const Impl&) = default;
  Impl& operator=(const Impl& a_xmUGrid) = default;

  // Misc
  bool GetModified() const;
  void SetUnmodified();
  void SetUseCache(bool a_useCache);

  // Points

  int GetPointCount() const;

  const VecPt3d& GetLocations() const;
  void SetLocations(const VecPt3d& a_locations);

  Pt3d GetPointLocation(int a_pointIdx) const;
  bool SetPointLocation(int a_pointIdx, const Pt3d& a_location);

  Pt3d GetPointXy0(int a_pointIdx) const;
  VecPt3d GetPointsLocations(const VecInt& a_points) const;

  void GetExtents(Pt3d& a_min, Pt3d& a_max) const;

  int GetPointAdjacentCellCount(int a_pointIdx) const;
  VecInt GetPointAdjacentCells(int a_pointIdx) const; // cells associated with point
  void GetPointAdjacentCells(int a_pointIdx, VecInt& a_adjacentCells) const;

  VecInt GetPointsAdjacentCells(const VecInt& a_points) const;
  void GetPointsAdjacentCells(const int* a_pointIdxs,
                              int a_numpointIdxs,
                              VecInt& a_commonCellIdxs) const;
  void GetPointsAdjacentCells(const VecInt& a_pointIdxs, VecInt& a_commonCellIdxs) const;
  void GetPointsAdjacentCells(int a_pointIdx1, int a_pointIdx2, VecInt& a_adjacentCellIdxs) const;
  bool IsValidPointChange(int a_changedPtIdx, const Pt3d& a_newPosition) const;

  // Cells
  int GetCellCount() const;

  int GetCellPointCount(int a_cellIdx) const;
  VecInt GetCellPoints(int a_cellIdx) const;
  bool GetCellPoints(int a_cellIdx, VecInt& a_cellPoints) const;
  void GetCellLocations(int a_cellIdx, VecPt3d& a_cellLocations) const;

  XmUGridCellType GetCellType(int a_cellIdx) const;
  std::vector<int> GetDimensionCounts() const;
  int GetCellDimension(int a_cellIdx) const;
  void GetCellExtents(int a_cellIdx, Pt3d& a_min, Pt3d& a_max) const;

  const VecInt& GetCellstream() const;
  bool SetCellstream(const VecInt& a_cellstream);
  bool GetCellCellstream(int a_cellIdx, VecInt& a_cellstream) const;

  VecInt GetCellAdjacentCells(int a_cellIdx) const;
  void GetCellAdjacentCells(int a_cellIdx, VecInt& a_cellNeighbors) const;
  bool GetCellPlanViewPolygon(int a_cellIdx, VecPt3d& a_polygon) const;
  bool GetCellCentroid(int a_cellIdx, Pt3d& a_centroid) const;

  // Edges
  int GetCellEdgeCount(int a_cellIdx) const;
  XmEdge GetCellEdge(int a_cellIdx, int a_edgeIdx) const;
  VecInt GetCellEdgeAdjacentCells(int a_cellIdx, int a_edgeIdx) const;
  void GetCellEdgeAdjacentCells(int a_cellIdx, int a_edgeIdx, VecInt& a_adjacentCellIdxs) const;
  int GetCell2dEdgeAdjacentCell(int a_cellIdx, int a_edgeIdx) const;
  void GetEdgeAdjacentCells(const XmEdge& a_edge, VecInt& a_adjacentCellIdxs) const;

  VecInt GetEdgeAdjacentCells(const XmEdge& a_edge) const;

  std::vector<XmEdge> GetCellEdges(int a_cellIdx) const;
  void GetCellEdges(int a_cellIdx, std::vector<XmEdge>& a_edges) const;
  void GetPointAdjacentPoints(int a_pointIdx, VecInt& a_edgePoints) const;
  void GetPointAdjacentLocations(int a_pointIdx, VecPt3d& a_edgePoints) const;

  // Faces
  int GetCell3dFaceCount(int a_cellIdx) const;
  int GetCell3dFacePointCount(int a_cellIdx, int a_faceIdx) const;

  VecInt GetCell3dFacePoints(int a_cellIdx, int a_faceIdx) const;
  void GetCell3dFacePoints(int a_cellIdx, int a_faceIdx, VecInt& a_facePtIdxs) const;
  VecInt2d GetCell3dFacesPoints(int a_cellIdx) const;

  int GetCell3dFaceAdjacentCell(int a_cellIdx, int a_faceIdx) const;
  bool GetCell3dFaceAdjacentCell(int a_cellIdx,
                                 int a_faceIdx,
                                 int& a_neighborCell,
                                 int& a_neighborFace) const;

  XmUGridFaceOrientation GetCell3dFaceOrientation(int a_cellIdx, int a_faceIdx) const;
  XmUGridFaceOrientation FaceOrientation(int a_cellIdx, int a_faceIdx) const;
  XmUGridFaceOrientation VerticalOrientationFromOpposing(int a_cellIdx, int a_faceIdx) const;
  bool IsSideFace(int a_cellIdx, int a_faceIdx) const;
  XmUGridFaceOrientation ConnectedTopOrBottom(int a_cellIdx, int a_faceIdx) const;
  XmUGridFaceOrientation GetOrientationFromArea(int a_cellIdx, int a_faceIdx) const;

private:
  void UpdateLinks(); // Calls UpdateCellLinks & UpdatePointLinks
  void UpdateCellLinks();
  void UpdatePointLinks();

  void SetModified();

  bool IsCellValidWithPointChange(int a_cellIdx,
                                  int a_changedPtIdx,
                                  const Pt3d& a_newPosition) const;
  bool IsValidCellIdx(int a_cellIdx) const;

  static int DimensionFromCellType(XmUGridCellType a_cellType);

  int GetNumberOfItemsForCell(int a_cellIdx) const;

  // Optimized for efficiency
  void GetCellCellstream(int a_cellIdx, const int** a_start, int& a_length) const;

  int GetNumberOfPolyhedronEdges(int a_cellIdx) const;
  static void GetUniquePointsFromPolyhedronCellstream(const VecInt& a_cellstream,
                                                      int a_numCellItems,
                                                      int& a_currIdx,
                                                      VecInt& a_uniqueGetCellPoints,
                                                      VecInt& a_pointLastUsedIdx);
  static bool GetUniquePointsFromPolyhedronSingleCellstream(const VecInt& a_cellstream,
                                                            VecInt& a_cellPoints);
  static void GetUniqueEdgesFromPolyhedronCellstream(
    const int* a_start,
    int& a_length,
    boost::container::flat_set<XmEdge>& a_cellEdges,
    int& a_currIdx);

  bool GetPlanViewPolygon2d(int a_cellIdx, VecPt3d& a_polygon) const;
  bool GetPlanViewPolygon3d(int a_cellIdx, VecPt3d& a_polygon) const;

  bool IsFaceSide(const VecInt& a_facePts) const;                 // plan view
  bool GetCellXySegments(int cellIdx, VecPt3d& a_segments) const; // plan view
  void GetEdgesOfFace(const VecInt& a_face, std::vector<XmEdge>& a_edges) const;
  bool DoEdgesCrossWithPointChange(int a_changedPtIdx,
                                   const Pt3d& a_newPosition,
                                   const std::vector<XmEdge>& a_edges) const;
  void GetExtentsFromPoints(const VecPt3d& a_locations, Pt3d& a_min, Pt3d& a_max) const;
  bool GetFaceXySegments(int a_cellIdx, int a_faceIdx, VecPt3d& a_segments) const; // plan view

  void CalculateCacheValues() const;
  void ClearCacheValues();
  int GetCell3dFaceCountNoCache(int a_cellIdx) const;
  int GetCell3dFaceAdjacentCellNoCache(int a_cellIdx, int a_faceIdx) const;
  XmUGridFaceOrientation GetCell3dFaceOrientationNoCache(int a_cellIdx, int a_faceIdx) const;
  bool GetNextFaceColumn(const VecInt& a_facePoints,
                         int a_starti,
                         int& a_columnBegin,
                         int& a_columnEnd) const;
  void GetFacePointSegments(const VecInt& a_facePts,
                            int a_columnBegin,
                            int a_columnEnd,
                            VecPt3d& a_segments) const;
  XmUGridFaceOrientation FaceOrientationWithFail(int a_cellIdx, int a_faceIdx) const;

  enum XmUGridCacheHolder {
    NEEDS_CALCULATION = -2 
  };

  VecPt3d m_locations;              
  VecInt m_cellstream;              
  VecInt m_cellIdxToStreamIdx;      
  VecInt m_pointsToCells;           
  VecInt m_pointIdxToPointsToCells; 
  bool m_modified = false;          
  bool m_useCache = true;           
  mutable VecInt m_numberOfFaces;   
  mutable VecInt m_cellFaceOffset;  
  mutable VecInt m_faceOrientation; 
  mutable VecInt m_faceNeighbor;    
  mutable VecInt m_cellDimensionCounts; 
};

//----- Internal functions -----------------------------------------------------
namespace
{
const char* CELLSTREAM_OFFSETS = "CELLSTREAM_OFFSETS";
const char* POINTS_TO_CELLS = "POINTS_TO_CELLS";
const char* POINTS_TO_CELLS_OFFSETS = "POINTS_TO_CELLS_OFFSETS";

typedef std::vector<XmEdge> VecEdge; 

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecEdge& iGetEdgeOffsetTable(int a_cellType)
{
  static const VecEdge fg_empty;
  static const VecEdge fg_line = {{0, 1}};
  static const VecEdge fg_triangle = {{0, 1}, {1, 2}, {2, 0}};
  static const VecEdge fg_pixel = {{0, 1}, {1, 3}, {3, 2}, {2, 0}};
  static const VecEdge fg_quad = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
  static const VecEdge fg_tetra = {{0, 1}, {1, 2}, {2, 0}, {0, 3}, {1, 3}, {2, 3}};
  static const VecEdge fg_wedge = {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5},
                                   {5, 3}, {0, 3}, {1, 4}, {2, 5}};
  static const VecEdge fg_voxel = {{0, 1}, {1, 3}, {2, 3}, {0, 2}, {4, 5}, {5, 7},
                                   {6, 7}, {4, 6}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
  static const VecEdge fg_hexahedron = {{0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6},
                                        {7, 6}, {4, 7}, {0, 4}, {1, 5}, {3, 7}, {2, 6}};
  static const VecEdge fg_pyramid = {{0, 1}, {1, 2}, {2, 3}, {3, 0},
                                     {0, 4}, {1, 4}, {2, 4}, {3, 4}};
  static const VecEdge fg_pentagonalPrism = {{0, 1}, {1, 2}, {2, 3}, {3, 4}, {4, 0},
                                             {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 5},
                                             {0, 5}, {1, 6}, {2, 7}, {3, 8}, {4, 9}};
  static const VecEdge fg_hexagonalPrism = {{0, 1}, {1, 2}, {2, 3}, {3, 4},  {4, 5},   {5, 0},
                                            {6, 7}, {7, 8}, {8, 9}, {9, 10}, {10, 11}, {11, 6},
                                            {0, 6}, {1, 7}, {2, 8}, {3, 9},  {4, 10},  {5, 11}};

  switch (a_cellType)
  {
  // invalid
  case XMU_INVALID_CELL_TYPE:
    return fg_empty;

  // 0D
  case XMU_EMPTY_CELL:
  case XMU_VERTEX:
  case XMU_POLY_VERTEX:
  case XMU_CONVEX_POINT_SET: // Special class of cells formed by convex group of
                             // points
    return fg_empty;

  // 1D
  case XMU_LINE:
    return fg_line;

  // 2D
  case XMU_TRIANGLE:
    return fg_triangle;
  case XMU_PIXEL:
    return fg_pixel;
  case XMU_QUAD:
    return fg_quad;

  case XMU_POLYGON:
    // should be using cell stream
    XM_ASSERT(0);
    return fg_empty;

  // 2D Not yet supported
  case XMU_QUADRATIC_EDGE:
  case XMU_PARAMETRIC_CURVE:
  case XMU_HIGHER_ORDER_EDGE:
  case XMU_CUBIC_LINE: // Cubic, isoparametric cell
  case XMU_TRIANGLE_STRIP:
  case XMU_QUADRATIC_TRIANGLE:
  case XMU_BIQUADRATIC_TRIANGLE:
  case XMU_HIGHER_ORDER_TRIANGLE:
  case XMU_HIGHER_ORDER_POLYGON:
  case XMU_QUADRATIC_POLYGON:
  case XMU_PARAMETRIC_SURFACE:
  case XMU_PARAMETRIC_TRI_SURFACE:
  case XMU_PARAMETRIC_QUAD_SURFACE:
    XM_ASSERT(0);
    return fg_empty;

  // 3D dimensions
  case XMU_TETRA:
    return fg_tetra;
  case XMU_WEDGE:
    return fg_wedge;
  case XMU_VOXEL:
    return fg_voxel;
  case XMU_HEXAHEDRON:
    return fg_hexahedron;
  case XMU_PYRAMID:
    return fg_pyramid;
  case XMU_PENTAGONAL_PRISM:
    return fg_pentagonalPrism;
  case XMU_HEXAGONAL_PRISM:
    return fg_hexagonalPrism;
  case XMU_POLYHEDRON:
    // should be using cell stream
    XM_ASSERT(0);
    return fg_empty;

  // 3D Not yet supported
  case XMU_QUADRATIC_TETRA:
  case XMU_HIGHER_ORDER_TETRAHEDRON:
  case XMU_QUADRATIC_WEDGE:
  case XMU_QUADRATIC_LINEAR_WEDGE:
  case XMU_BIQUADRATIC_QUADRATIC_WEDGE:
  case XMU_HIGHER_ORDER_WEDGE:
  case XMU_QUADRATIC_HEXAHEDRON:
  case XMU_TRIQUADRATIC_HEXAHEDRON:
  case XMU_BIQUADRATIC_QUADRATIC_HEXAHEDRON:
  case XMU_HIGHER_ORDER_HEXAHEDRON:
  case XMU_QUADRATIC_PYRAMID:
  case XMU_HIGHER_ORDER_PYRAMID:
  case XMU_PARAMETRIC_TETRA_REGION:
  case XMU_PARAMETRIC_HEX_REGION:
    XM_ASSERT(0);
    return fg_empty;
  default:
    XM_ASSERT(0);
    return fg_empty;
  }

  XM_ASSERT(0);
  return fg_empty;
} // iGetEdgeOffsetTable
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecInt2d& iGetFaceOffsetTable(int a_cellType)
{
  static const VecInt2d fg_empty;
  static const VecInt2d fg_tetra = {{0, 1, 3}, {1, 2, 3}, {2, 0, 3}, {0, 2, 1}};
  static const VecInt2d fg_wedge = {{0, 1, 2}, {3, 5, 4}, {0, 3, 4, 1}, {1, 4, 5, 2}, {2, 5, 3, 0}};
  static const VecInt2d fg_voxel = {{0, 4, 6, 2}, {1, 3, 7, 5}, {0, 1, 5, 4},
                                    {2, 6, 7, 3}, {0, 2, 3, 1}, {4, 5, 7, 6}};
  static const VecInt2d fg_hexahedron = {{0, 4, 7, 3}, {1, 2, 6, 5}, {0, 1, 5, 4},
                                         {3, 7, 6, 2}, {0, 3, 2, 1}, {4, 5, 6, 7}};
  static const VecInt2d fg_pyramid = {{0, 3, 2, 1}, {0, 1, 4}, {1, 2, 4}, {2, 3, 4}, {3, 0, 4}};
  static const VecInt2d fg_pentagonalPrism = {{0, 4, 3, 2, 1}, {5, 6, 7, 8, 9}, {0, 1, 6, 5},
                                              {1, 2, 7, 6},    {2, 3, 8, 7},    {3, 4, 9, 8},
                                              {4, 0, 5, 9}};
  static const VecInt2d fg_hexagonalPrism = {
    {0, 5, 4, 3, 2, 1}, {6, 7, 8, 9, 10, 11}, {0, 1, 7, 6},   {1, 2, 8, 7},
    {2, 3, 9, 8},       {3, 4, 10, 9},        {4, 5, 11, 10}, {5, 0, 6, 11}};

  switch (a_cellType)
  {
  // invalid
  // 0D
  // 1D
  case XMU_INVALID_CELL_TYPE:
  case XMU_EMPTY_CELL:
  case XMU_VERTEX:
  case XMU_POLY_VERTEX:
  case XMU_CONVEX_POINT_SET: // Special class of cells formed by convex group of
                             // points
  case XMU_LINE:
  case XMU_TRIANGLE:
  case XMU_PIXEL:
  case XMU_QUAD:
  case XMU_POLYGON:
  case XMU_POLY_LINE:
  case XMU_QUADRATIC_EDGE:
  case XMU_PARAMETRIC_CURVE:
  case XMU_HIGHER_ORDER_EDGE:
  case XMU_CUBIC_LINE: // Cubic, isoparametric cell
  case XMU_TRIANGLE_STRIP:
  case XMU_QUADRATIC_TRIANGLE:
  case XMU_BIQUADRATIC_TRIANGLE:
  case XMU_HIGHER_ORDER_TRIANGLE:
  case XMU_HIGHER_ORDER_POLYGON:
  case XMU_QUADRATIC_POLYGON:
  case XMU_PARAMETRIC_SURFACE:
  case XMU_PARAMETRIC_TRI_SURFACE:
  case XMU_PARAMETRIC_QUAD_SURFACE:
    return fg_empty;

  // 3D dimensions
  case XMU_TETRA:
    return fg_tetra;
  case XMU_WEDGE:
    return fg_wedge;
  case XMU_VOXEL:
    return fg_voxel;
  case XMU_HEXAHEDRON:
    return fg_hexahedron;
  case XMU_PYRAMID:
    return fg_pyramid;
  case XMU_PENTAGONAL_PRISM:
    return fg_pentagonalPrism;
  case XMU_HEXAGONAL_PRISM:
    return fg_hexagonalPrism;
  case XMU_POLYHEDRON:
    // should be using cell stream
    XM_ASSERT(0);
    return fg_empty;

  // 3D Not yet supported
  case XMU_QUADRATIC_TETRA:
  case XMU_HIGHER_ORDER_TETRAHEDRON:
  case XMU_QUADRATIC_WEDGE:
  case XMU_QUADRATIC_LINEAR_WEDGE:
  case XMU_BIQUADRATIC_QUADRATIC_WEDGE:
  case XMU_HIGHER_ORDER_WEDGE:
  case XMU_QUADRATIC_HEXAHEDRON:
  case XMU_TRIQUADRATIC_HEXAHEDRON:
  case XMU_BIQUADRATIC_QUADRATIC_HEXAHEDRON:
  case XMU_HIGHER_ORDER_HEXAHEDRON:
  case XMU_QUADRATIC_PYRAMID:
  case XMU_HIGHER_ORDER_PYRAMID:
  case XMU_PARAMETRIC_TETRA_REGION:
  case XMU_PARAMETRIC_HEX_REGION:
    XM_ASSERT(0);
    return fg_empty;
  default:
    XM_ASSERT(0);
    return fg_empty;
  }

  XM_ASSERT(0);
  return fg_empty;
} // iGetFaceOffsetTable
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool iPointInSegmentBounds(const Pt3d& a_location, const Pt3d& a_first, const Pt3d& a_second)
{
  bool inBounds = a_location.x >= std::min(a_first.x, a_second.x) &&
                  a_location.y >= std::min(a_first.y, a_second.y) &&
                  a_location.x <= std::max(a_first.x, a_second.x) &&
                  a_location.y <= std::max(a_first.y, a_second.y);
  return inBounds;
} // iPointInSegmentBounds
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <typename T>
std::vector<T> iGetUniquePoints(const std::vector<T>& a_segments)
{
  std::vector<T> points = a_segments;
  std::sort(points.begin(), points.end());
  auto pIt = std::unique(points.begin(), points.end());
  points.resize(pIt - points.begin());
  return points;
} // iGetUniquePoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <typename T>
void iBuildPolygon(std::vector<std::pair<T, T>>& a_segs, std::vector<T>& a_polygon)
{
  xms::VecChar placed(a_segs.size(), false);
  a_polygon.push_back(a_segs[0].first);
  a_polygon.push_back(a_segs[0].second);
  placed[0] = true;
  for (size_t i = 1; i != a_segs.size(); ++i)
  {
    T& toMatch = a_polygon.back();
    for (size_t j = 1; j != a_segs.size(); ++j)
    {
      if (!placed[j])
      {
        if (a_segs[j].first == toMatch)
        {
          a_polygon.push_back(a_segs[j].second);
          placed[j] = true;
        }
        else if (a_segs[j].second == toMatch)
        {
          a_polygon.push_back(a_segs[j].first);
          placed[j] = true;
        }
      }
    }
  }
} // iBuildPolygon
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <typename T>
std::vector<std::pair<T, T>> iGetUniqueSegments(const std::vector<T>& a_segments)
{
  std::vector<std::pair<T, T>> segs;
  for (size_t i = 0; i < a_segments.size(); i += 2)
  {
    T a = a_segments[i];
    T b = a_segments[i + 1];
    if (a < b)
      segs.push_back(std::pair<T, T>(a, b));
    else
      segs.push_back(std::pair<T, T>(b, a));
  }
  std::sort(segs.begin(), segs.end());
  auto segIt = std::unique(segs.begin(), segs.end());
  segs.resize(segIt - segs.begin());
  return segs;
} // iGetUniqueSegments
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void iMergeSegmentsToPoly(const VecPt3d& a_segments, VecPt3d& a_polygon)
{
  if (a_segments.empty())
    return;

  VecPt3d points = iGetUniquePoints(a_segments);
  std::vector<std::pair<Pt3d, Pt3d>> segs = iGetUniqueSegments(a_segments);

  // if any point is on any other segment split up the segment
  for (auto& point : points)
  {
    size_t segsSize = segs.size();
    for (size_t i = 0; i < segsSize; ++i)
    {
      Pt3d pt1 = segs[i].first;
      Pt3d pt2 = segs[i].second;
      if (point != pt1 && point != pt2 && iPointInSegmentBounds(point, pt1, pt2))
      {
        Pt3d diff = pt1 - pt2;
        double distance = std::max(fabs(diff.x), fabs(diff.y));
        double tolerance = distance / 1.0e6;
        if (gmOnLineWithTol(pt1, pt2, point.x, point.y, tolerance))
        {
          segs.emplace_back(point, pt2);
          segs[i].second = point;
        }
      }
    }
  }
  std::sort(segs.begin(), segs.end());
  auto segIt = std::unique(segs.begin(), segs.end());
  segs.resize(segIt - segs.begin());

  iBuildPolygon(segs, a_polygon);

  if (a_polygon.size() > 3 && a_polygon.front() == a_polygon.back())
  {
    a_polygon.pop_back();
    double area = gmPolygonArea(&a_polygon[0], (int)a_polygon.size());
    if (area < 0)
      std::reverse(a_polygon.begin(), a_polygon.end());
  }
  else
  {
    XM_ASSERT(0);
    a_polygon.clear();
  }
} // iMergeSegmentsToPoly

} // namespace

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetModified() const
{
  return m_modified;
} // XmUGrid::Impl::GetModified
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::SetUnmodified()
{
  m_modified = false;
} // XmUGrid::Impl::SetUnmodified
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::SetModified()
{
  ClearCacheValues();
  m_modified = true;
} // XmUGrid::Impl::SetModified
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::SetUseCache(bool a_useCache)
{
  m_useCache = a_useCache;
} // XmUGrid::Impl::SetUseCache
// Points
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetPointCount() const
{
  return (int)m_locations.size();
} // XmUGrid::Impl::GetPointCount
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecPt3d& XmUGrid::Impl::GetLocations() const
{
  return m_locations;
} // XmUGrid::Impl::GetLocations
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::SetLocations(const VecPt3d& a_locations)
{
  m_locations = a_locations;
  SetModified();
} // XmUGrid::Impl::SetLocations
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
Pt3d XmUGrid::Impl::GetPointLocation(int a_pointIdx) const
{
  if (a_pointIdx >= 0 && a_pointIdx < m_locations.size())
    return m_locations[a_pointIdx];
  return Pt3d();
} // XmUGrid::Impl::GetPointLocation
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::SetPointLocation(int a_pointIdx, const Pt3d& a_location)
{
  if (a_pointIdx >= 0 && a_pointIdx < m_locations.size())
  {
    m_locations[a_pointIdx] = a_location;
    SetModified();
    return true;
  }
  return false;
} // XmUGrid::Impl::SetPointLocation

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
Pt3d XmUGrid::Impl::GetPointXy0(int a_pointIdx) const
{
  Pt3d pt = GetPointLocation(a_pointIdx);
  pt.z = 0.0;
  return pt;
} // XmUGrid::Impl::GetPointXy0

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecPt3d XmUGrid::Impl::GetPointsLocations(const VecInt& a_points) const
{
  VecPt3d point3d;
  for (auto point : a_points)
  {
    point3d.push_back(GetPointLocation(point));
  }
  return point3d;
} // XmUGrid::Impl::GetPointsLocations

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetExtents(Pt3d& a_min, Pt3d& a_max) const
{
  GetExtentsFromPoints(m_locations, a_min, a_max);
} // XmUGrid::Impl::GetExtents

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetPointAdjacentCellCount(int a_pointIdx) const
{
  if (a_pointIdx < 0 || a_pointIdx >= m_locations.size())
    return 0;
  int numCells = m_pointsToCells[m_pointIdxToPointsToCells[a_pointIdx]];
  return numCells;
} // XmUGrid::Impl::GetPointAdjacentCellCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetPointAdjacentCells(int a_pointIdx) const
{
  VecInt cellsOfPoint;
  GetPointAdjacentCells(a_pointIdx, cellsOfPoint);
  return cellsOfPoint;
} // XmUGrid::Impl::GetPointAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetPointAdjacentCells(int a_pointIdx, VecInt& a_adjacentCells) const
{
  a_adjacentCells.clear();
  int numCells = GetPointAdjacentCellCount(a_pointIdx);
  for (int cellIdx = 0; cellIdx < numCells; cellIdx++)
  {
    a_adjacentCells.push_back(m_pointsToCells[m_pointIdxToPointsToCells[a_pointIdx] + cellIdx + 1]);
  }
} // XmUGrid::Impl::GetPointAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetPointsAdjacentCells(const VecInt& a_points) const
{
  VecInt commonCells;
  if (!a_points.empty())
    GetPointsAdjacentCells(&a_points[0], (int)a_points.size(), commonCells);
  return commonCells;
} // XmUGrid::Impl::GetPointsAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetPointsAdjacentCells(const int* a_pointIdxs,
                                           int a_numPointIdxs,
                                           VecInt& a_commonCellIdxs) const
{
  a_commonCellIdxs.clear();
  if (a_numPointIdxs == 0)
    return;
  a_commonCellIdxs = GetPointAdjacentCells(a_pointIdxs[0]);
  for (int i = 1; i < a_numPointIdxs; ++i)
  {
    VecInt tempAssociatedCells = GetPointAdjacentCells(a_pointIdxs[i]);
    VecInt remove;
    for (int j = 0; j < a_commonCellIdxs.size(); ++j)
    {
      bool found = false;
      for (int tempAssociatedCell : tempAssociatedCells)
      {
        if (a_commonCellIdxs[j] == tempAssociatedCell)
        {
          found = true;
          break;
        }
      }
      if (!found)
        remove.push_back(j);
    }
    for (int j = (int)remove.size() - 1; j >= 0; --j)
    {
      a_commonCellIdxs.erase(a_commonCellIdxs.begin() + remove[j]);
    }
    if (a_commonCellIdxs.empty())
      break;
  }
} // XmUGrid::Impl::GetPointsAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetPointsAdjacentCells(const VecInt& a_pointIdxs,
                                           VecInt& a_adjacentCellIdxs) const
{
  if (!a_pointIdxs.empty())
    GetPointsAdjacentCells(&a_pointIdxs[0], (int)a_pointIdxs.size(), a_adjacentCellIdxs);
} // XmUGrid::Impl::GetPointsAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetPointsAdjacentCells(int a_pointIdx1,
                                           int a_pointIdx2,
                                           VecInt& a_adjacentCellIdxs) const
{
  int points[] = {a_pointIdx1, a_pointIdx2};
  GetPointsAdjacentCells(points, 2, a_adjacentCellIdxs);
} // XmUGrid::Impl::GetPointsAdjacentCells
// Cells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCellCount() const
{
  if (m_cellIdxToStreamIdx.empty())
    return 0;
  else
    return (int)m_cellIdxToStreamIdx.size() - 1;
} // XmUGrid::Impl::GetCellCount
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCellPointCount(int a_cellIdx) const
{
  int pointCount = 0;
  if (GetCellType(a_cellIdx) == XMU_POLYHEDRON)
  {
    boost::container::flat_set<int> uniqueIdxs;
    const int* cellstream = nullptr;
    int streamLength;
    GetCellCellstream(a_cellIdx, &cellstream, streamLength);
    if (cellstream)
    {
      int currIdx = 1;
      int numFaces = cellstream[currIdx++];
      for (int faceIdx = 0; faceIdx < numFaces; faceIdx++)
      {
        int numPoints = cellstream[currIdx++];
        for (int pointIdx = 0; pointIdx < numPoints; ++pointIdx)
        {
          int point = cellstream[currIdx++];
          uniqueIdxs.insert(point);
        }
      }
    }
    pointCount = (int)uniqueIdxs.size();
  }
  else
  {
    pointCount = GetNumberOfItemsForCell(a_cellIdx);
  }

  return pointCount;
} // XmUGrid::Impl::GetCellPointCount
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetCellPoints(int a_cellIdx) const
{
  VecInt pointsOfCell;
  GetCellPoints(a_cellIdx, pointsOfCell);
  return pointsOfCell;
} // XmUGrid::Impl::GetCellPoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetCellPoints(int a_cellIdx, VecInt& a_cellPoints) const
{
  a_cellPoints.clear();
  VecInt cellstream;
  if (GetCellCellstream(a_cellIdx, cellstream))
  {
    XM_ENSURE_TRUE_NO_ASSERT(!cellstream.empty(), false);
    int cellType = cellstream[0];
    if (cellType == XMU_POLYHEDRON)
    {
      GetUniquePointsFromPolyhedronSingleCellstream(cellstream, a_cellPoints);
      return true;
    }
    else if (cellType == XMU_PIXEL)
    {
      a_cellPoints.push_back(cellstream[2]);
      a_cellPoints.push_back(cellstream[3]);
      a_cellPoints.push_back(cellstream[5]);
      a_cellPoints.push_back(cellstream[4]);
      return true;
    }
    else
    {
      a_cellPoints.assign(cellstream.begin() + 2, cellstream.end());
      return true;
    }
  }
  return false;
} // XmUGrid::Impl::GetCellPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCellLocations(int a_cellIdx, VecPt3d& a_cellLocations) const
{
  VecInt ptIdxs = GetCellPoints(a_cellIdx);
  a_cellLocations = GetPointsLocations(ptIdxs);
} // XmUGrid::Impl::GetCellPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridCellType XmUGrid::Impl::GetCellType(int a_cellIdx) const
{
  if (a_cellIdx < 0 || m_cellIdxToStreamIdx.size() < 2 ||
      a_cellIdx > m_cellIdxToStreamIdx.size() - 2)
  {
    return XMU_INVALID_CELL_TYPE;
  }
  else
  {
    int cellStart = m_cellIdxToStreamIdx[a_cellIdx];
#if _DEBUG
    if ((m_cellstream[cellStart] > XMU_PYRAMID) && (m_cellstream[cellStart] != XMU_POLYHEDRON) ||
        (m_cellstream[cellStart] == XMU_TRIANGLE_STRIP))
    {
      assert("UNSUPPORTED TYPE!");
    }
#endif
    return (XmUGridCellType)m_cellstream[cellStart];
  }
} // XmUGrid::Impl::GetCellType
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::vector<int> XmUGrid::Impl::GetDimensionCounts() const
{
  if (!m_cellDimensionCounts.empty())
  {
    return m_cellDimensionCounts;
  }
  else
  {
    m_cellDimensionCounts.clear();
    m_cellDimensionCounts.resize(4, 0);
    int itemp = 0;
    int cellCount = GetCellCount();
    for (int i = 0; i < cellCount; i++)
    {
      itemp = GetCellDimension(i);
      if (itemp >= 0)
      {
        m_cellDimensionCounts[itemp]++;
      }
    }
    return m_cellDimensionCounts;
  }
} // XmUGrid::Impl::GetDimensionCounts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCellDimension(int a_cellIdx) const
{
  return DimensionFromCellType(GetCellType(a_cellIdx));
} // XmUGrid::Impl::GetCellDimension

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCellExtents(int a_cellIdx, Pt3d& a_min, Pt3d& a_max) const
{
  VecPt3d pts;
  GetCellLocations(a_cellIdx, pts);
  GetExtentsFromPoints(pts, a_min, a_max);
} // XmUGrid::Impl::GetCellExtents

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecInt& XmUGrid::Impl::GetCellstream() const
{
  return m_cellstream;
} // XmUGrid::Impl::GetCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::SetCellstream(const VecInt& a_cellstream)
{
  if (IsValidCellstream(a_cellstream))
  {
    m_cellstream = a_cellstream;
    UpdateLinks();
    SetModified();
    return true;
  }
  else
  {
    XM_LOG(xmlog::error, "Invalid cell stream data.");
    return false;
  }
} // XmUGrid::Impl::SetCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetCellCellstream(int a_cellIdx, VecInt& a_cellstream) const
{
  a_cellstream.clear();
  if (a_cellIdx < 0 || m_cellIdxToStreamIdx.size() < 2 ||
      a_cellIdx > m_cellIdxToStreamIdx.size() - 2)
  {
    return false;
  }
  int startIndex(m_cellIdxToStreamIdx[a_cellIdx]), endIndex(m_cellIdxToStreamIdx[a_cellIdx + 1]);
  a_cellstream.assign(m_cellstream.begin() + startIndex, m_cellstream.begin() + endIndex);
  return true;
} // XmUGrid::Impl::GetCellCellstream

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetCellAdjacentCells(int a_cellIdx) const
{
  VecInt neighbors;
  GetCellAdjacentCells(a_cellIdx, neighbors);
  return neighbors;
} // XmUGrid::Impl::GetCellAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCellAdjacentCells(int a_cellIdx, VecInt& a_cellNeighbors) const
{
  a_cellNeighbors.clear();
  VecInt pointsOfCell;
  pointsOfCell = GetCellPoints(a_cellIdx);
  for (int pointIdx : pointsOfCell)
  {
    VecInt associatedCells = GetPointAdjacentCells(pointIdx);
    for (int associatedCell : associatedCells)
    {
      if (associatedCell != a_cellIdx)
      {
        bool found(false);
        for (int a_cellNeighbor : a_cellNeighbors)
        {
          if (a_cellNeighbor == associatedCell)
          {
            found = true;
          }
        }
        if (!found)
        {
          a_cellNeighbors.push_back(associatedCell);
        }
      }
    }
  }
} // XmUGrid::Impl::GetCellAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetCellPlanViewPolygon(int a_cellIdx, VecPt3d& a_polygon) const
{
  a_polygon.clear();
  if (!IsValidCellIdx(a_cellIdx))
    return false;
  int dimension = GetCellDimension(a_cellIdx);
  if (dimension == 3)
    return GetPlanViewPolygon3d(a_cellIdx, a_polygon);
  else if (dimension == 2)
    return GetPlanViewPolygon2d(a_cellIdx, a_polygon);
  else
    return false;
} // XmUGrid::Impl::GetCellPlanViewPolygon

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetCellCentroid(int a_cellIdx, Pt3d& a_centroid) const
{
  VecPt3d polyPoints;
  bool retVal = true;
  Pt3d centroid = Pt3d(0.0, 0.0, 0.0);
  if (GetCellPlanViewPolygon(a_cellIdx, polyPoints))
  {
    centroid = gmComputePolygonCentroid(polyPoints);
  }
  else if (!IsValidCellIdx(a_cellIdx))
  {
    retVal = false;
  }
  else
  {
    VecPt3d cellPoints;
    GetCellLocations(a_cellIdx, cellPoints);
    for (Pt3d& pt : cellPoints)
    {
      centroid += pt;
    }
    centroid /= (double)cellPoints.size();
  }
  a_centroid = centroid;
  return retVal;
} // XmUGrid::Impl::GetCellCentroid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::IsCellValidWithPointChange(int a_cellIdx,
                                               int a_changedPtIdx,
                                               const Pt3d& a_newPosition) const
{
  if (GetCellDimension(a_cellIdx) == 2)
  {
    std::vector<XmEdge> edges = GetCellEdges(a_cellIdx);
    return !DoEdgesCrossWithPointChange(a_changedPtIdx, a_newPosition, edges);
  }
  else if (GetCellDimension(a_cellIdx) == 3)
  {
    //  Go through the affected faces and if they are not a side face
    // check if the edges cross each other
    VecInt2d faces = GetCell3dFacesPoints(a_cellIdx);
    for (int faceIdx = 0; faceIdx < faces.size(); faceIdx++)
    {
      bool faceIsAffected = false;
      for (int ptIdx = 0; ptIdx < faces[faceIdx].size(); ptIdx++)
      {
        if (faces[faceIdx][ptIdx] == a_changedPtIdx)
        {
          faceIsAffected = true;
        }
      }
      if (faceIsAffected)
      {
        std::vector<XmEdge> edges;
        GetEdgesOfFace(faces[faceIdx], edges);
        if (DoEdgesCrossWithPointChange(a_changedPtIdx, a_newPosition, edges))
          return false;
      }
    }
  }
  return true;
} // XmUGrid::Impl::IsCellValidWithPointChange
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::IsValidCellIdx(int a_cellIdx) const
{
  return a_cellIdx >= 0 && a_cellIdx < GetCellCount();
} // XmUGrid::Impl::IsValidCellIdx
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::IsValidPointChange(int a_changedPtIdx, const Pt3d& a_newPosition) const
{
  bool validChange = false;
  if (a_changedPtIdx >= 0 && a_changedPtIdx < m_locations.size())
  {
    validChange = true;
    VecInt affectedCells = GetPointAdjacentCells(a_changedPtIdx);
    for (int affectedCell : affectedCells)
    {
      if (!IsCellValidWithPointChange(affectedCell, a_changedPtIdx, a_newPosition))
      {
        validChange = false;
        break;
      }
    }
  }
  return validChange;
} // XmUGrid::Impl::IsValidPointChange
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetEdgesOfFace(const VecInt& a_face, std::vector<XmEdge>& a_edges) const
{
  a_edges.reserve(a_face.size());
  for (int i = 1; i < a_face.size(); ++i)
  {
    a_edges.emplace_back(a_face[i - 1], a_face[i]);
  }
  if (a_edges.size() > 1)
    a_edges.emplace_back(a_face[a_face.size() - 1], a_face[0]);
} // XmUGrid::Impl::GetEdgesOfFace

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::DoEdgesCrossWithPointChange(int a_changedPtIdx,
                                                const Pt3d& a_newPosition,
                                                const std::vector<XmEdge>& a_edges) const
{
  std::vector<std::pair<Pt3d, Pt3d>> changedEdges;
  std::vector<std::pair<Pt3d, Pt3d>> unChangedEdges;
  for (const auto& edge : a_edges)
  {
    if (edge.GetFirst() == a_changedPtIdx)
    {
      changedEdges.emplace_back(a_newPosition, GetPointLocation(edge.GetSecond()));
    }
    else if (edge.GetSecond() == a_changedPtIdx)
    {
      changedEdges.emplace_back(GetPointLocation(edge.GetFirst()), a_newPosition);
    }
    else
    {
      unChangedEdges.emplace_back(GetPointLocation(edge.GetFirst()),
                                  GetPointLocation(edge.GetSecond()));
    }
  }
  for (const auto& changedEdge : changedEdges)
  {
    for (const auto& unChangedEdge : unChangedEdges)
    {
      if (gmLinesCross(changedEdge.first, changedEdge.second, unChangedEdge.first,
                       unChangedEdge.second))
        return true;
    }
  }
  return false;
} // XmUGrid::Impl::DoEdgesCrossWithPointChange
// Edges
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCellEdgeCount(int a_cellIdx) const
{
  int cellType(GetCellType(a_cellIdx));
  switch (cellType)
  {
  // invalid
  case XMU_INVALID_CELL_TYPE:
    return -1;
    break;

  case XMU_POLY_LINE:
    return GetNumberOfItemsForCell(a_cellIdx) - 1;
    break;

  case XMU_POLYGON:
    return GetNumberOfItemsForCell(a_cellIdx);
    break;

  case XMU_POLYHEDRON:
    return GetNumberOfPolyhedronEdges(a_cellIdx);
    break;

  default:
  {
    const VecEdge& edgeTable = iGetEdgeOffsetTable(cellType);
    if (!edgeTable.empty())
      return (int)edgeTable.size();
    break;
  }
  }

  XM_ASSERT(0);
  return -1;
} // XmUGrid::Impl::GetCellEdgeCount
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmEdge XmUGrid::Impl::GetCellEdge(int a_cellIdx, int a_edgeIdx) const
{
  XmEdge edge;
  if (a_edgeIdx < 0)
    return edge;

  const int* cellstream = nullptr;
  int streamLength;
  GetCellCellstream(a_cellIdx, &cellstream, streamLength);
  if (cellstream)
  {
    if (streamLength < 4)
      return edge;

    int cellType = cellstream[0];
    switch (cellType)
    {
    case XMU_POLY_LINE:
    case XMU_POLYGON:
    {
      int numPoints = cellstream[1];
      const int* cellPoints = cellstream + 2;
      int idx1 = cellPoints[a_edgeIdx];
      int idx2 = cellPoints[(a_edgeIdx + 1) % numPoints];
      edge = XmEdge(idx1, idx2);
      break;
    }

    case XMU_POLYHEDRON:
    {
      boost::container::flat_set<XmEdge> cellEdges;
      int currIdx = 2;
      GetUniqueEdgesFromPolyhedronCellstream(cellstream, streamLength, cellEdges, currIdx);

      if (a_edgeIdx < (int)cellEdges.size())
        edge = *(cellEdges.begin() + a_edgeIdx);
      break;
    }

    default:
    {
      const VecEdge& edgeTable = iGetEdgeOffsetTable(cellType);
      if (a_edgeIdx < (int)edgeTable.size())
      {
        const XmEdge& edgeOffset = edgeTable[a_edgeIdx];
        const int* cellPoints = cellstream + 2;
        int idx1 = cellPoints[edgeOffset.GetFirst()];
        int idx2 = cellPoints[edgeOffset.GetSecond()];
        edge = XmEdge(idx1, idx2);
      }
      break;
    }
    }
  }

  return edge;
} // XmUGrid::Impl::GetCellEdge
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetCellEdgeAdjacentCells(int a_cellIdx, int a_edgeIdx) const
{
  VecInt adjacentCells;
  GetCellEdgeAdjacentCells(a_cellIdx, a_edgeIdx, adjacentCells);
  return adjacentCells;
} // XmUGrid::Impl::GetCellEdgeAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCellEdgeAdjacentCells(int a_cellIdx,
                                             int a_edgeIdx,
                                             VecInt& a_adjacentCellIdxs) const
{
  a_adjacentCellIdxs.clear();
  VecInt adjacentCells;
  VecInt cellNeighbors = GetCellAdjacentCells(a_cellIdx);
  if (cellNeighbors.empty())
  {
    return;
  }
  if (a_edgeIdx < 0 || a_edgeIdx >= GetCellEdgeCount(a_cellIdx))
  {
    return;
  }

  XmEdge currEdge, neighborEdge;
  currEdge = GetCellEdge(a_cellIdx, a_edgeIdx);
  for (int cellNeighbor : cellNeighbors)
  {
    for (int k(0); k < GetCellEdgeCount(cellNeighbor); k++)
    {
      neighborEdge = GetCellEdge(cellNeighbor, k);
      if (currEdge.IsEquivalent(neighborEdge))
      {
        a_adjacentCellIdxs.push_back(cellNeighbor);
      }
    }
  }
} // XmUGrid::Impl::GetCellEdgeAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCell2dEdgeAdjacentCell(int a_cellIdx, int a_edgeIdx) const
{
  if (GetCellDimension(a_cellIdx) != 2)
    return -1;

  XmEdge edge = GetCellEdge(a_cellIdx, a_edgeIdx);
  int pointIdx = edge.GetFirst();

  int adjCellCount = GetPointAdjacentCellCount(pointIdx);
  for (int cellIdx = 0; cellIdx < adjCellCount; cellIdx++)
  {
    int adjacentCell = m_pointsToCells[m_pointIdxToPointsToCells[pointIdx] + cellIdx + 1];
    if (adjacentCell == a_cellIdx)
      continue;

    int edgeCount = GetCellEdgeCount(adjacentCell);
    for (int adjacentEdgeIdx = 0; adjacentEdgeIdx < edgeCount; ++adjacentEdgeIdx)
    {
      XmEdge adjacentEdge = GetCellEdge(adjacentCell, adjacentEdgeIdx);
      if (adjacentEdge.IsEquivalent(edge))
        return adjacentCell;
    }
  }

  return -1;
} // XmUGrid::Impl::GetCell2dEdgeAdjacentCell
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetEdgeAdjacentCells(const XmEdge& a_edge, VecInt& a_adjacentCellIdxs) const
{
  a_adjacentCellIdxs.clear();
  GetPointsAdjacentCells(a_edge.GetFirst(), a_edge.GetSecond(), a_adjacentCellIdxs);
} // XmUGrid::Impl::GetEdgeAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetEdgeAdjacentCells(const XmEdge& a_edge) const
{
  VecInt adjacentCellIdxs;
  GetEdgeAdjacentCells(a_edge, adjacentCellIdxs);
  return adjacentCellIdxs;
} // XmUGrid::Impl::GetEdgeAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::vector<XmEdge> XmUGrid::Impl::GetCellEdges(int a_cellIdx) const
{
  std::vector<XmEdge> edges;
  GetCellEdges(a_cellIdx, edges);
  return edges;
} // XmUGrid::Impl::GetCellEdges
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCellEdges(int a_cellIdx, std::vector<XmEdge>& a_edges) const
{
  a_edges.clear();
  int numEdges = GetCellEdgeCount(a_cellIdx);
  for (int i = 0; i < numEdges; ++i)
  {
    a_edges.push_back(GetCellEdge(a_cellIdx, i));
  }
} // XmUGrid::Impl::GetCellEdges

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetPointAdjacentPoints(int a_pointIdx, VecInt& a_edgePoints) const
{
  a_edgePoints.clear();
  VecInt associatedCells = GetPointAdjacentCells(a_pointIdx);
  if (associatedCells.empty())
  {
    return;
  }
  for (int associatedCell : associatedCells)
  {
    for (int j = 0; j < GetCellEdgeCount(associatedCell); ++j)
    {
      XmEdge temp = GetCellEdge(associatedCell, j);
      if (temp.GetFirst() == a_pointIdx)
      {
        a_edgePoints.push_back(temp.GetSecond());
      }
      else if (temp.GetSecond() == a_pointIdx)
      {
        a_edgePoints.push_back(temp.GetFirst());
      }
    }
  }
  std::sort(a_edgePoints.begin(), a_edgePoints.end());
  auto it = std::unique(a_edgePoints.begin(), a_edgePoints.end());
  a_edgePoints.erase(it, a_edgePoints.end());
} // XmUGrid::Impl::GetPointAdjacentPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetPointAdjacentLocations(int a_pointIdx, VecPt3d& a_edgePoints) const
{
  VecInt edgePtIdxs;
  GetPointAdjacentPoints(a_pointIdx, edgePtIdxs);
  a_edgePoints = GetPointsLocations(edgePtIdxs);
} // XmUGrid::Impl::GetPointAdjacentLocations

// Faces
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCell3dFaceCount(int a_cellIdx) const
{
  if (m_useCache)
  {
    CalculateCacheValues();
    int faceCount = IsValidCellIdx(a_cellIdx) ? m_numberOfFaces[a_cellIdx] : -1;
    return faceCount;
  }
  else
  {
    return GetCell3dFaceCountNoCache(a_cellIdx);
  }
} // XmUGrid::Impl::GetCell3dFaceCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCell3dFacePointCount(int a_cellIdx, int a_faceIdx) const
{
  if (!IsValidCellIdx(a_cellIdx))
  {
    return -1;
  }

  int cellType = GetCellType(a_cellIdx);
  if (cellType == XMU_POLYHEDRON)
  {
    const int* cellstream;
    int length;
    GetCellCellstream(a_cellIdx, &cellstream, length);
    if (cellstream != nullptr)
    {
      auto currItem = cellstream;
      currItem++; // skip cell type
      int numFaces = *currItem++;
      for (int faceIdx = 0; faceIdx < numFaces; ++faceIdx)
      {
        int numFacePoints = *currItem++;
        if (a_faceIdx == faceIdx)
        {
          return numFacePoints;
        }
        currItem += numFacePoints;
      }
    }
  }
  else
  {
    const VecInt2d& faceTable = iGetFaceOffsetTable(cellType);
    if (a_faceIdx >= 0 && a_faceIdx < (int)faceTable.size())
    {
      return (int)faceTable[a_faceIdx].size();
    }
  }

  return 0;
} // XmUGrid::Impl::GetCell3dFacePointCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::Impl::GetCell3dFacePoints(int a_cellIdx, int a_faceIdx) const
{
  VecInt facePoints;
  GetCell3dFacePoints(a_cellIdx, a_faceIdx, facePoints);
  return facePoints;
} // XmUGrid::Impl::GetCell3dFacePoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCell3dFacePoints(int a_cellIdx, int a_faceIdx, VecInt& a_facePtIdxs) const
{
  a_facePtIdxs.clear();
  if (a_faceIdx < 0)
  {
    return;
  }

  const int* cellstream = nullptr;
  int streamLength;
  GetCellCellstream(a_cellIdx, &cellstream, streamLength);
  if (!cellstream)
    return;

  int cellType = *cellstream++;
  if (cellType == XMU_POLYHEDRON)
  {
    int numFaces = *cellstream++;
    for (int faceIdx = 0; faceIdx < numFaces; ++faceIdx)
    {
      int numFacePoints = *cellstream++;
      if (faceIdx == a_faceIdx)
      {
        a_facePtIdxs.assign(cellstream, cellstream + numFacePoints);
        break;
      }
      cellstream += numFacePoints;
    }
  }
  else
  {
    const VecInt2d& faceTable = iGetFaceOffsetTable(cellType);
    if (a_faceIdx < (int)faceTable.size())
    {
      if (cellstream)
      {
        const int* cellPoints = ++cellstream;
        const VecInt& faceOffsets = faceTable[a_faceIdx];
        for (auto offset : faceOffsets)
        {
          int idx = cellPoints[offset];
          a_facePtIdxs.push_back(idx);
        }
      }
    }
  }
} // XmUGrid::Impl::GetCell3dFacePoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt2d XmUGrid::Impl::GetCell3dFacesPoints(int a_cellIdx) const
{
  int numFaces = GetCell3dFaceCount(a_cellIdx);

  VecInt2d faces(numFaces);
  for (int i(0); i < numFaces; i++)
  {
    faces[i] = GetCell3dFacePoints(a_cellIdx, i);
  }
  return faces;
} // XmUGrid::Impl::GetCell3dFacesPoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCell3dFaceAdjacentCell(int a_cellIdx, int a_faceIdx) const
{
  if (m_useCache)
  {
    CalculateCacheValues();
    int faceNeighbor = -1;
    int numFaces = GetCell3dFaceCount(a_cellIdx);
    if (a_faceIdx >= 0 && a_faceIdx < numFaces)
    {
      int faceOffset = m_cellFaceOffset[a_cellIdx];
      if (faceOffset >= 0)
      {
        faceOffset += a_faceIdx;
        faceNeighbor = m_faceNeighbor[faceOffset];
        if (faceNeighbor == NEEDS_CALCULATION)
        {
          faceNeighbor = GetCell3dFaceAdjacentCellNoCache(a_cellIdx, a_faceIdx);
          m_faceNeighbor[faceOffset] = faceNeighbor;
        }
        else
        {
          faceNeighbor = m_faceNeighbor[faceOffset];
        }
      }
    }
    return faceNeighbor;
  }
  else
  {
    return GetCell3dFaceAdjacentCellNoCache(a_cellIdx, a_faceIdx);
  }
} // XmUGrid::Impl::GetCell3dFaceAdjacentCell
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetCell3dFaceAdjacentCell(int a_cellIdx,
                                              int a_faceIdx,
                                              int& a_neighborCell,
                                              int& a_neighborFace) const
{
  a_neighborCell = GetCell3dFaceAdjacentCell(a_cellIdx, a_faceIdx);
  if (a_neighborCell < 0)
  {
    a_neighborCell = a_neighborFace = -1;
    return false;
  }
  VecInt cellFacePts = GetCell3dFacePoints(a_cellIdx, a_faceIdx);
  for (int faceIdx(0); faceIdx < GetCell3dFaceCount(a_neighborCell); faceIdx++)
  {
    VecInt curCellFacePts = GetCell3dFacePoints(a_neighborCell, faceIdx);
    if (cellFacePts.size() == curCellFacePts.size())
    {
      int cellFacePtsIdx(0);
      for (; cellFacePtsIdx < (int)cellFacePts.size(); ++cellFacePtsIdx)
      {
        bool found = false;
        for (int curCellFacePt : curCellFacePts)
        {
          if (cellFacePts[cellFacePtsIdx] == curCellFacePt)
          {
            found = true;
            break;
          }
        }
        if (!found)
          break; // This face is not the face we are l.
      }
      if (cellFacePtsIdx == (int)cellFacePts.size()) // Every point is found
      {
        a_neighborFace = faceIdx;
        a_neighborFace = faceIdx;
        return true;
      }
    }
  }
  a_neighborCell = a_neighborFace = -1;
  return false;
} // XmUGrid::Impl::GetCell3dFaceAdjacentCell
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::GetCell3dFaceOrientation(int a_cellIdx, int a_faceIdx) const
{
  if (m_useCache)
  {
    CalculateCacheValues();
    int faceOrientation = XMU_ORIENTATION_UNKNOWN;
    int numFaces = GetCell3dFaceCount(a_cellIdx);
    if (a_faceIdx >= 0 && a_faceIdx < numFaces)
    {
      int faceOffset = m_cellFaceOffset[a_cellIdx];
      if (faceOffset >= 0)
      {
        faceOffset += a_faceIdx;
        faceOrientation = m_faceOrientation[faceOffset];
        if (faceOrientation == NEEDS_CALCULATION)
        {
          faceOrientation = (int)GetCell3dFaceOrientationNoCache(a_cellIdx, a_faceIdx);
          m_faceOrientation[faceOffset] = faceOrientation;
        }
        else
        {
          faceOrientation = m_faceOrientation[faceOffset];
        }
      }
    }
    return (XmUGridFaceOrientation)faceOrientation;
  }
  else
  {
    return GetCell3dFaceOrientationNoCache(a_cellIdx, a_faceIdx);
  }
} // XmUGrid::Impl::GetCell3dFaceOrientation
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::FaceOrientation(int a_cellIdx, int a_faceIdx) const
{
  XmUGridFaceOrientation orientation = FaceOrientationWithFail(a_cellIdx, a_faceIdx);
  if (orientation == XMU_ORIENTATION_UNKNOWN)
  {
    // Is there any way to know if we reach this point? First face top and next
    // bottom? What about multi-panel top with no cell above?
    // XM_ASSERT(0);
    orientation = VerticalOrientationFromOpposing(a_cellIdx, a_faceIdx);
  }

  XM_ASSERT(orientation != XMU_ORIENTATION_UNKNOWN);
  return orientation;
} // XmUGrid::Impl::FaceOrientation
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::UpdateLinks()
{
  UpdateCellLinks();
  UpdatePointLinks();
} // XmUGrid::Impl::UpdateLinks
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::UpdateCellLinks()
{
  m_cellIdxToStreamIdx.clear();

  int currIdx = 0;
  if (m_cellstream.empty())
  {
    m_cellIdxToStreamIdx.push_back(currIdx);
    return;
  }

  int numItems = (int)m_cellstream.size();
  while (currIdx < numItems)
  {
    m_cellIdxToStreamIdx.push_back(currIdx);

    // get cell type
    int cellType = m_cellstream[currIdx++];
    if (currIdx >= numItems)
      return;

    // get the number of items
    int numPoints = m_cellstream[currIdx++];
    if (currIdx >= numItems)
      return;

    if (cellType == XMU_POLYHEDRON)
    {
      int numFaces = numPoints;
      for (int faceIdx = 0; faceIdx < numFaces; ++faceIdx)
      {
        int numFacePoints = m_cellstream[currIdx++];
        currIdx += numFacePoints;
      }
    }
    else
    {
      currIdx += numPoints;
    }
  }

  m_cellIdxToStreamIdx.push_back(currIdx);
} // XmUGrid::Impl::UpdateCellLinks
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::UpdatePointLinks()
{
  m_pointIdxToPointsToCells.clear();
  m_pointIdxToPointsToCells.resize(m_locations.size(), 0);
  m_pointsToCells.clear();

  // get number of cells for each point
  int numStreamItems = (int)m_cellstream.size();
  int cellIdx = 0;
  int currIdx = 0;
  VecInt pointLastUsedIdx(m_locations.size(), -1);
  VecInt cellPoints;
  while (currIdx < numStreamItems)
  {
    // get cell type
    int cellType = m_cellstream[currIdx++];
    if (currIdx >= numStreamItems)
      return;

    // get the number of items (points or faces depending on cell type)
    int numCellItems = m_cellstream[currIdx++];
    if (currIdx >= numStreamItems)
      return;

    if (cellType == XMU_POLYHEDRON)
    {
      GetUniquePointsFromPolyhedronCellstream(m_cellstream, numCellItems, currIdx, cellPoints,
                                              pointLastUsedIdx);

      // Deemed to be slower than flat set-- Left only as a warning to others!
      // std::stable_sort(cellPoints.begin(), cellPoints.end());
      // auto uniqueEnd = std::unique(cellPoints.begin(), cellPoints.end());

      // flat set deemed to be slower than marking when a point was last used (pointLastUsedIdx)

      // for (auto pt = cellPoints.begin(); pt != uniqueEnd; ++pt)
      for (auto cellPoint : cellPoints)
      {
        m_pointIdxToPointsToCells[cellPoint] += 1;
      }
    }
    else
    {
      // iterate on points
      int numPoints = numCellItems;
      for (int ptIdx = 0; ptIdx < numPoints; ++ptIdx)
      {
        int pt = m_cellstream[currIdx++];
        m_pointIdxToPointsToCells[pt] += 1;
      }
    }

    ++cellIdx;
  }
  m_pointIdxToPointsToCells.push_back(0);

  // change array of counts to array of offsets
  int currCount = 0;
  for (auto& item : m_pointIdxToPointsToCells)
  {
    int count = item + 1;
    item = currCount;
    currCount += count;
  }

  cellIdx = 0;
  currIdx = 0;
  m_pointsToCells.resize(currCount, 0);
  std::fill(pointLastUsedIdx.begin(), pointLastUsedIdx.end(), -1);
  while (currIdx < numStreamItems)
  {
    // get cell type
    int cellType = m_cellstream[currIdx++];
    if (currIdx >= numStreamItems)
      return;

    // get the number of items (points or faces depending on cell type)
    int numCellItems = m_cellstream[currIdx++];
    if (currIdx >= numStreamItems)
      return;

    if (cellType == XMU_POLYHEDRON)
    {
      GetUniquePointsFromPolyhedronCellstream(m_cellstream, numCellItems, currIdx, cellPoints,
                                              pointLastUsedIdx);

      // Deemed to be slower than flat set-- Left only as a warning to others!
      // std::stable_sort(cellPoints.begin(), cellPoints.end());
      // auto uniqueEnd = std::unique(cellPoints.begin(), cellPoints.end());

      // for (auto pt = cellPoints.begin(); pt != uniqueEnd; ++pt)
      for (auto cellPoint : cellPoints)
      {
        int countIdx = m_pointIdxToPointsToCells[cellPoint];
        int& count = m_pointsToCells[countIdx]; // point's cell count
        ++count;                                // incrementing point's cell count
        int updateIdx = countIdx + count;
        m_pointsToCells[updateIdx] = cellIdx;
      }
    }
    else
    {
      // iterate on points
      int numPoints = numCellItems;
      for (int ptIdx = 0; ptIdx < numPoints; ++ptIdx)
      {
        int pt = m_cellstream[currIdx++];
        int countIdx = m_pointIdxToPointsToCells[pt];
        int& count = m_pointsToCells[countIdx]; // point's cell count
        ++count;                                // incrementing point's cell count
        int updateIdx = countIdx + count;
        m_pointsToCells[updateIdx] = cellIdx;
      }
    }

    ++cellIdx;
  }
} // XmUGrid::Impl::UpdatePointLinks
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::DimensionFromCellType(XmUGridCellType a_cellType)
{
  switch (a_cellType)
  {
  // invalid
  case XMU_INVALID_CELL_TYPE:
    return -1;
    break;

  // 0D
  case XMU_EMPTY_CELL:
  case XMU_VERTEX:
  case XMU_POLY_VERTEX:
  case XMU_CONVEX_POINT_SET: // Special class of cells formed by convex group of
                             // points
    return 0;
    break;

  // 1D
  case XMU_LINE:
  case XMU_POLY_LINE:

  case XMU_QUADRATIC_EDGE:
  case XMU_PARAMETRIC_CURVE:
  case XMU_HIGHER_ORDER_EDGE:
  case XMU_CUBIC_LINE: // Cubic, isoparametric cell
    return 1;
    break;
  // 2D
  case XMU_TRIANGLE:
  case XMU_TRIANGLE_STRIP:
  case XMU_POLYGON:
  case XMU_PIXEL:
  case XMU_QUAD:

  case XMU_QUADRATIC_TRIANGLE:
  case XMU_QUADRATIC_QUAD:
  case XMU_BIQUADRATIC_QUAD:
  case XMU_QUADRATIC_LINEAR_QUAD:
  case XMU_BIQUADRATIC_TRIANGLE:
  case XMU_QUADRATIC_POLYGON:

  case XMU_PARAMETRIC_SURFACE:
  case XMU_PARAMETRIC_TRI_SURFACE:
  case XMU_PARAMETRIC_QUAD_SURFACE:

  case XMU_HIGHER_ORDER_TRIANGLE:
  case XMU_HIGHER_ORDER_QUAD:
  case XMU_HIGHER_ORDER_POLYGON:

    return 2;
    break;

  // assuming the rest are 3D
  default:
    return 3;
    break;
#if 0 // Remaining definitions
    XMU_TETRA = 10,
    XMU_VOXEL = 11,
    XMU_HEXAHEDRON = 12,
    XMU_WEDGE = 13,
    XMU_PYRAMID = 14,
    XMU_PENTAGONAL_PRISM = 15,
    XMU_HEXAGONAL_PRISM = 16,

    // Quadratic, isoparametric cells
    XMU_QUADRATIC_TETRA = 24,
    XMU_QUADRATIC_HEXAHEDRON = 25,
    XMU_QUADRATIC_WEDGE = 26,
    XMU_QUADRATIC_PYRAMID = 27,
    XMU_TRIQUADRATIC_HEXAHEDRON = 29,
    XMU_QUADRATIC_LINEAR_WEDGE = 31,
    XMU_BIQUADRATIC_QUADRATIC_WEDGE = 32,
    XMU_BIQUADRATIC_QUADRATIC_HEXAHEDRON = 33,

    // Polyhedron cell (consisting of polygonal faces)
    XMU_POLYHEDRON = 42,

    // Higher order cells in parametric form
    XMU_PARAMETRIC_TETRA_REGION = 55,
    XMU_PARAMETRIC_HEX_REGION = 56,

    // Higher order cells
    XMU_HIGHER_ORDER_TETRAHEDRON = 64,
    XMU_HIGHER_ORDER_WEDGE = 65,
    XMU_HIGHER_ORDER_PYRAMID = 66,
    XMU_HIGHER_ORDER_HEXAHEDRON = 67,
#endif
  }
  return -1;

} // XmUGrid::Impl::DimensionFromCellType

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetNumberOfItemsForCell(int a_cellIdx) const
{
  int cellType = GetCellType(a_cellIdx);
  if (cellType == XMU_INVALID_CELL_TYPE)
  {
    return 0;
  }
  else if (cellType != XMU_POLYHEDRON)
  {
    // Number of points is right after cell type
    int startIndex = m_cellIdxToStreamIdx[a_cellIdx];
    return m_cellstream[startIndex + 1];
  }
  else
  { // Polyhedron case, same value is number of faces
    int startIndex = m_cellIdxToStreamIdx[a_cellIdx];
    return m_cellstream[startIndex + 1];
  }
} // XmUGrid::Impl::GetNumberOfItemsForCell
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetCellCellstream(int a_cellIdx, const int** a_start, int& a_length) const
{
  int cellType = GetCellType(a_cellIdx);
  if (cellType == XMU_INVALID_CELL_TYPE)
  {
    *a_start = nullptr;
    a_length = 0;
  }
  else
  {
    int startIndex = m_cellIdxToStreamIdx[a_cellIdx];
    int nextCellIndex = m_cellIdxToStreamIdx[a_cellIdx + 1];
    *a_start = &m_cellstream[startIndex];
    a_length = nextCellIndex - startIndex;
  }
} // XmUGrid::Impl::GetCellCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetNumberOfPolyhedronEdges(int a_cellIdx) const
{
  const int* cellstream;
  int streamLength;
  GetCellCellstream(a_cellIdx, &cellstream, streamLength);
  if (cellstream && streamLength > 0 && cellstream[0] == XMU_POLYHEDRON)
  {
    boost::container::flat_set<XmEdge> edges;
    int currItem = 2;
    while (currItem < streamLength)
    {
      GetUniqueEdgesFromPolyhedronCellstream(cellstream, streamLength, edges, currItem);
    }
    return (int)edges.size();
  }
  return 0;
} // XmUGrid::Impl::GetNumberOfPolyhedronEdges
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetUniquePointsFromPolyhedronCellstream(const VecInt& a_cellstream,
                                                            int a_numCellItems,
                                                            int& a_currIdx,
                                                            VecInt& a_uniqueGetCellPoints,
                                                            VecInt& a_pointLastUsedIdx)
{
  int stable = a_currIdx;
  a_uniqueGetCellPoints.clear();
  int numFaces = a_numCellItems;
  for (int faceIdx = 0; faceIdx < numFaces; ++faceIdx)
  {
    int numFacePoints = a_cellstream[a_currIdx++];
    for (int ptIdx = 0; ptIdx < numFacePoints; ++ptIdx)
    {
      int pt = a_cellstream[a_currIdx++];
      if (a_pointLastUsedIdx[pt] < stable)
      {
        a_pointLastUsedIdx[pt] = stable;
        a_uniqueGetCellPoints.push_back(pt);
      }
    }
  }
} // XmUGrid::Impl::GetUniquePointsFromPolyhedronCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetUniquePointsFromPolyhedronSingleCellstream(const VecInt& a_cellstream,
                                                                  VecInt& a_cellPoints)
{
  a_cellPoints.clear();
  if (a_cellstream.size() < 2)
  {
    return false;
  }

  auto currItem = a_cellstream.begin();
  int cellType = *currItem++;
  if (cellType != XMU_POLYHEDRON)
  {
    return false;
  }

  int numFaces = *currItem++;
  for (int faceIdx = 0; faceIdx < numFaces; ++faceIdx)
  {
    int numFacePoints = *currItem++;
    for (int ptIdx = 0; ptIdx < numFacePoints; ++ptIdx)
    {
      int pt = *currItem++;
      if (std::find(a_cellPoints.begin(), a_cellPoints.end(), pt) == a_cellPoints.end())
      {
        a_cellPoints.push_back(pt);
      }
    }
  }

  return true;
} // XmUGrid::Impl::GetUniquePointsFromPolyhedronSingleCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetUniqueEdgesFromPolyhedronCellstream(
  const int* a_start,
  int& a_length,
  boost::container::flat_set<XmEdge>& a_cellEdges,
  int& a_currIdx)
{
  int numFaces = a_start[1];
  for (int i(0); i < numFaces; i++)
  {
    int numPoints = a_start[a_currIdx++];
    for (int pointIdx = 0; pointIdx < numPoints; ++pointIdx)
    {
      int pt1Idx = pointIdx;
      int pt2Idx = (pointIdx + 1) % numPoints;
      // The % operator will reset the index back 0 so the final loop will
      // provide the last and first point
      int pt1 = a_start[a_currIdx + pt1Idx];
      int pt2 = a_start[a_currIdx + pt2Idx];
      // We want unique edges, so we add the lower point index first
      if (pt1 < pt2)
        a_cellEdges.insert(XmEdge(pt1, pt2));
      else
        a_cellEdges.insert(XmEdge(pt2, pt1));
    }
    a_currIdx += numPoints;
  }
} // XmUGrid::Impl::GetUniqueEdgesFromPolyhedronCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetPlanViewPolygon2d(int a_cellIdx, VecPt3d& a_polygon) const
{
  GetCellLocations(a_cellIdx, a_polygon);
  if (!a_polygon.empty())
  {
    return true;
  }
  return false;
} // XmUGrid::Impl::GetPlanViewPolygon2d
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetPlanViewPolygon3d(int a_cellIdx, VecPt3d& a_polygon) const
{
  VecPt3d segments;
  if (GetCellXySegments(a_cellIdx, segments))
  {
    // Prismatic cell
    iMergeSegmentsToPoly(segments, a_polygon);
    return true;
  }
  else
  {
    // Non-prismatic cell
    VecInt uniquePoints = GetCellPoints(a_cellIdx);
    VecPt3d cellPoints(uniquePoints.size());
    for (int i(0); i < uniquePoints.size(); i++)
    {
      cellPoints[i] = GetPointLocation(uniquePoints[i]);
      cellPoints[i].z = 0.0; // Insist that our z values are 0.0 for plan view
    }
    gmGetConvexHull(cellPoints, a_polygon, false);
    return true;
  }
  return false;
} // XmUGrid::Impl::GetPlanViewPolygon3d
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::IsFaceSide(const VecInt& a_facePts) const
{
  // if any face point has same x and y as any other face point then
  // the face is a side face
  if (a_facePts.size() > 2)
  {
    Pt3d ptLast = GetPointLocation(a_facePts[0]);
    for (size_t facePtIdx = 1; facePtIdx < a_facePts.size(); ++facePtIdx)
    {
      Pt3d ptCurr = GetPointLocation(a_facePts[facePtIdx]);
      if ((ptLast.x == ptCurr.x) && (ptLast.y == ptCurr.y))
        return true;

      ptLast = ptCurr;
    }
  }
  return false;
} // XmUGrid::Impl::IsFaceSide
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetCellXySegments(int a_cellIdx, VecPt3d& a_segments) const
{
  bool foundSideFace = false;
  for (int faceIdx = 0; faceIdx < GetCell3dFaceCount(a_cellIdx); ++faceIdx)
  {
    VecInt ptIdxs = GetCell3dFacePoints(a_cellIdx, faceIdx);
    if (IsFaceSide(ptIdxs))
    {
      foundSideFace = true;
      if (!GetFaceXySegments(a_cellIdx, faceIdx, a_segments))
      {
        return false;
      }
    }
  }
  return foundSideFace;
} // XmUGrid::Impl::GetCellXySegments

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetExtentsFromPoints(const VecPt3d& a_locations, Pt3d& a_min, Pt3d& a_max) const
{
  a_min.x = a_min.y = a_min.z = xms::XM_DBL_HIGHEST;
  a_max.x = a_max.y = a_max.z = xms::XM_DBL_LOWEST;
  for (const Pt3d& pt : a_locations)
  {
    gmAddToExtents(pt, a_min, a_max);
  }
} // XmUGrid::Impl::GetExtentsFromPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetFaceXySegments(int a_cellIdx, int a_faceIdx, VecPt3d& a_segments) const
{
  VecInt facePts;
  GetCell3dFacePoints(a_cellIdx, a_faceIdx, facePts);
  if (facePts.empty())
    return false;

  int column1Begin, column1End;
  if (!GetNextFaceColumn(facePts, 0, column1Begin, column1End))
    return false;

  int column2Begin, column2End;
  if (!GetNextFaceColumn(facePts, column1End, column2Begin, column2End))
    return false;

  GetFacePointSegments(facePts, column1End, column2Begin, a_segments);
  GetFacePointSegments(facePts, column2End, column1Begin, a_segments);
  return true;
} // XmUGrid::Impl::GetFaceXySegments
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::CalculateCacheValues() const
{
  if (m_numberOfFaces.empty() && GetCellCount() != 0)
  {
    int cellCount = GetCellCount();
    m_numberOfFaces.assign(cellCount, 0);
    m_cellFaceOffset.assign(cellCount + 1, 0);
    int faceCount = 0;
    for (int cellIdx = 0; cellIdx < cellCount; ++cellIdx)
    {
      int numberOfFaces = GetCell3dFaceCountNoCache(cellIdx);
      m_numberOfFaces[cellIdx] = numberOfFaces;
      faceCount += numberOfFaces;
      m_cellFaceOffset[cellIdx + 1] = faceCount;
    }

    m_faceOrientation.assign(faceCount, NEEDS_CALCULATION);
    m_faceNeighbor.assign(faceCount, NEEDS_CALCULATION);
  }
} // XmUGrid::Impl::CalculateCacheValues
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::ClearCacheValues()
{
  m_numberOfFaces.clear();
  m_cellFaceOffset.clear();
  m_faceOrientation.clear();
  m_faceNeighbor.clear();
  m_cellDimensionCounts.clear();
} // XmUGrid::Impl::ClearCacheValues
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCell3dFaceCountNoCache(int a_cellIdx) const
{
  if (!IsValidCellIdx(a_cellIdx))
    return -1;

  int cellType(GetCellType(a_cellIdx));
  switch (cellType)
  {
  case XMU_POLYHEDRON:
    return GetNumberOfItemsForCell(a_cellIdx);
    break;

  default:
  {
    const VecInt2d& faceTable = iGetFaceOffsetTable(cellType);
    return (int)faceTable.size();
    break;
  }
  }

  return 0;
} // XmUGrid::Impl::GetCell3dFaceCountNoCache
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::Impl::GetCell3dFaceAdjacentCellNoCache(int a_cellIdx, int a_faceIdx) const
{
  VecInt cellFace = GetCell3dFacePoints(a_cellIdx, a_faceIdx);
  if (cellFace.empty())
    return -1;
  VecInt neighborCellFace = GetPointsAdjacentCells(cellFace);
  if (neighborCellFace.size() <= 1)
    return -1;
  if (neighborCellFace.size() > 2)
  {
    assert("Cell definitions are invalid; more than 2 cells found sharing the same face.");
  }
  for (int i : neighborCellFace)
  {
    if (i != a_cellIdx)
    {
      return i;
    }
  }
  return -1;
} // XmUGrid::Impl::GetCell3dFaceAdjacentCellNoCache
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::GetCell3dFaceOrientationNoCache(int a_cellIdx,
                                                                      int a_faceIdx) const
{
  XmUGridFaceOrientation faceOrientation = FaceOrientation(a_cellIdx, a_faceIdx);
  return faceOrientation;
} // XmUGrid::Impl::GetCell3dFaceOrientationNoCache
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::GetNextFaceColumn(const VecInt& a_facePoints,
                                      int a_starti,
                                      int& a_columnBegin,
                                      int& a_columnEnd) const
{
  int facePointsSize = (int)a_facePoints.size();

  // find next start of column of points (matching x/y values)
  Pt3d lastPt = GetPointXy0(a_facePoints[a_starti]);
  int lasti = a_starti;
  int i = (a_starti + 1) % facePointsSize;
  Pt3d pt;
  bool found = false;
  while (!found && i != a_starti)
  {
    pt = GetPointXy0(a_facePoints[i]);
    if (pt == lastPt)
    {
      a_columnBegin = lasti;
      found = true;
    }
    lastPt = pt;
    lasti = i;
    i = (i + 1) % facePointsSize;
  }

  if (!found)
  {
    // should have found first pair of points with matching x/y values
    return false;
  }

  // find end of column of points
  found = false;
  while (!found && i != a_columnBegin)
  {
    pt = GetPointXy0(a_facePoints[i]);
    if (pt != lastPt)
    {
      a_columnEnd = lasti;
      found = true;
    }
    else
    {
      lastPt = pt;
      lasti = i;
      i = (i + 1) % facePointsSize;
    }
  }

  if (!found)
  {
    // should have found end of points with matching x/y values
    return false;
  }

  return true;
} // XmUGrid::Impl::GetNextFaceColumn

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Impl::GetFacePointSegments(const VecInt& a_facePts,
                                         int a_columnBegin,
                                         int a_columnEnd,
                                         VecPt3d& a_segments) const
{
  int i = a_columnBegin;
  while (i != a_columnEnd)
  {
    Pt3d p = GetPointXy0(a_facePts[i]);
    a_segments.push_back(p);
    if (i != a_columnBegin)
      a_segments.push_back(p);
    i = (i + 1) % a_facePts.size();
  }
  a_segments.push_back(GetPointXy0(a_facePts[i]));
} // XmUGrid::Impl::GetFacePointSegments
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::FaceOrientationWithFail(int a_cellIdx, int a_faceIdx) const
{
  // if first face point has same x and y as any other face point then
  // the face is a side face
  if (IsSideFace(a_cellIdx, a_faceIdx))
    return XMU_ORIENTATION_SIDE;

  // else if face connected to cell lower idx cell then top face
  // else if face connected to cell higher idx cell then bottom face
  XmUGridFaceOrientation orientation = ConnectedTopOrBottom(a_cellIdx, a_faceIdx);
  if (orientation != XMU_ORIENTATION_UNKNOWN)
    return orientation;

  // else use sign of area to determine up or down
  orientation = GetOrientationFromArea(a_cellIdx, a_faceIdx);
  if (orientation != XMU_ORIENTATION_UNKNOWN)
    return orientation;

  return XMU_ORIENTATION_UNKNOWN;
} // XmUGrid::Impl::FaceOrientationWithFail
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::VerticalOrientationFromOpposing(int a_cellIdx,
                                                                      int a_faceIdx) const
{
  XmUGridFaceOrientation orientation = XMU_ORIENTATION_UNKNOWN;

  // Assume only 1 top and bottom.  Find other face that's top or bottom and
  // give answer as opposite.  Otherwise first top other is bottom.
  bool firstUnknown = false;
  bool foundTop = false;
  bool foundBot = false;
  for (int face = 0; face < GetCell3dFaceCount(a_cellIdx); ++face)
  {
    int orientation = FaceOrientationWithFail(a_cellIdx, face);
    if (orientation == XMU_ORIENTATION_UNKNOWN)
    {
      if (a_faceIdx == face)
        firstUnknown = true;
    }
    else if (orientation == XMU_ORIENTATION_TOP)
      foundTop = true;
    else if (orientation == XMU_ORIENTATION_BOTTOM)
      foundBot = true;
  }

  if (foundTop)
    orientation = XMU_ORIENTATION_BOTTOM;
  else if (foundBot)
    orientation = XMU_ORIENTATION_TOP;
  else if (firstUnknown)
    orientation = XMU_ORIENTATION_TOP;
  else
    orientation = XMU_ORIENTATION_BOTTOM;
  return orientation;
} // XmUGrid::Impl::VerticalOrientationFromOpposing
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::Impl::IsSideFace(int a_cellIdx, int a_faceIdx) const
{
  VecInt facePts;
  GetCell3dFacePoints(a_cellIdx, a_faceIdx, facePts);

  // if any face point has same x and y as any other face point then
  // the face is a side face
  if (!facePts.empty())
  {
    Pt3d ptLast = GetPointXy0(facePts[0]);
    for (size_t facePtIdx = 1; facePtIdx < facePts.size(); ++facePtIdx)
    {
      Pt3d ptCurr = GetPointXy0(facePts[facePtIdx]);
      if (ptLast == ptCurr)
        return true;

      ptLast = ptCurr;
    }
  }

  return false;
} // XmUGrid::Impl::IsSideFace
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::ConnectedTopOrBottom(int a_cellIdx, int a_faceIdx) const
{
  int adjacentCellIdx = GetCell3dFaceAdjacentCell(a_cellIdx, a_faceIdx);
  if (adjacentCellIdx != XM_NONE && adjacentCellIdx < a_cellIdx)
    return XMU_ORIENTATION_TOP;
  else if (adjacentCellIdx != XM_NONE && adjacentCellIdx > a_cellIdx)
    return XMU_ORIENTATION_BOTTOM;
  return XMU_ORIENTATION_UNKNOWN;
} // ConnectedTopOrBottom
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::Impl::GetOrientationFromArea(int a_cellIdx, int a_faceIdx) const
{
  VecInt facePtIdxs;
  GetCell3dFacePoints(a_cellIdx, a_faceIdx, facePtIdxs);
  VecPt3d facePts;
  for (int facePtIdx : facePtIdxs)
    facePts.push_back(GetPointXy0(facePtIdx));
  double area = gmPolygonArea(&facePts[0], facePts.size());
  if (area > 0.0)
    return XMU_ORIENTATION_TOP;
  else if (area < 0.0)
    return XMU_ORIENTATION_BOTTOM;
  return XMU_ORIENTATION_UNKNOWN;
} // XmUGrid::Impl::GetOrientationFromArea

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> XmUGrid::New()
{
  std::shared_ptr<XmUGrid> ugrid(new XmUGrid());
  return ugrid;
} // XmUGrid::New
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> XmUGrid::New(const VecPt3d& a_locations, const VecInt& a_cellstream)
{
  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New();
  ugrid->SetLocations(a_locations);
  ugrid->SetCellstream(a_cellstream);
  return ugrid;
} // XmUGrid::New
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGrid::XmUGrid()
: m_impl(new Impl)
{
} // XmUGrid::XmUGrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGrid::XmUGrid(const XmUGrid& a_xmUGrid)
: m_impl(new Impl)
{
  *m_impl = *a_xmUGrid.m_impl;
} // XmUGrid::XmUGrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGrid::XmUGrid(XmUGrid&& a_xmUGrid)
: XmUGrid()
{
  Swap(a_xmUGrid);
} // XmUGrid::XmUGrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGrid::~XmUGrid()
{
} // XmUGrid::~XmUGrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGrid& XmUGrid::operator=(XmUGrid a_xmUGrid)
{
  Swap(a_xmUGrid);
  return *this;
} // XmUGrid::operator=
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::Swap(XmUGrid& a_xmUGrid)
{
  std::swap(m_impl, a_xmUGrid.m_impl);
} // XmUGrid::Swap
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::IsValidCellstream(const VecInt& a_cellstream)
{
  if (a_cellstream.empty())
    return true;

  int numItems = (int)a_cellstream.size();
  int currIdx = 0;
  while (currIdx < numItems)
  {
    // get cell type
    int cellType = a_cellstream[currIdx++];
    if (currIdx >= numItems)
      return false;

    // get the number of points
    int numPoints = a_cellstream[currIdx++];
    if (currIdx >= numItems)
      return false;

    switch (cellType)
    {
    case XMU_EMPTY_CELL:
      if (numPoints != 0)
        return false;
      break;
    case XMU_VERTEX:
      if (numPoints != 1)
        return false;
      currIdx += numPoints;
      break;
    case XMU_LINE:
      if (numPoints != 2)
        return false;
      currIdx += numPoints;
      break;
    case XMU_TRIANGLE:
      if (numPoints != 3)
        return false;
      currIdx += numPoints;
      break;
    case XMU_PIXEL:
    case XMU_QUAD:
    case XMU_TETRA:
      if (numPoints != 4)
        return false;
      currIdx += numPoints;
      break;
    case XMU_PYRAMID:
      if (numPoints != 5)
        return false;
      currIdx += numPoints;
      break;
    case XMU_WEDGE:
    case XMU_PENTAGONAL_PRISM:
      if (numPoints != 6)
        return false;
      currIdx += numPoints;
      break;
    case XMU_HEXAGONAL_PRISM:
      if (numPoints != 7)
        return false;
      currIdx += numPoints;
      break;
    case XMU_VOXEL:
    case XMU_HEXAHEDRON:
      if (numPoints != 8)
        return false;
      currIdx += numPoints;
      break;
    case XMU_POLY_VERTEX:
      if (numPoints < 1)
        return false;
      currIdx += numPoints;
      break;
    case XMU_POLY_LINE:
      if (numPoints < 2)
        return false;
      currIdx += numPoints;
      break;
    case XMU_TRIANGLE_STRIP:
    case XMU_POLYGON:
      if (numPoints < 3)
        return false;
      currIdx += numPoints;
      break;
    case XMU_POLYHEDRON:
    {
      int numFaces = numPoints;
      if (numFaces < 3)
        return false;
      for (int faceIdx = 0; faceIdx < numFaces; ++faceIdx)
      {
        if (currIdx >= numItems)
          return false;
        int numPoints = a_cellstream[currIdx++];
        if (numPoints < 3)
          return false;
        currIdx += numPoints;
      }
    }
    break;
    default:
      currIdx += numPoints;
      break;
    }

    if (currIdx > numItems)
      return false;
  }

  return true;
} // XmUGrid::IsValidCellstream

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::GetModified() const
{
  return m_impl->GetModified();
} // XmUGrid::GetModified

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::SetUnmodified()
{
  m_impl->SetUnmodified();
} // XmUGrid::SetUnmodified

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::SetUseCache(bool a_useCache)
{
  m_impl->SetUseCache(a_useCache);
} // XmUGrid::SetUseCache

// Points

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetPointCount() const
{
  return m_impl->GetPointCount();
} // XmUGrid::GetPointCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecPt3d& XmUGrid::GetLocations() const
{
  return m_impl->GetLocations();
} // XmUGrid::GetLocations

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::SetLocations(const VecPt3d& a_locations)
{
  m_impl->SetLocations(a_locations);
} // XmUGrid::SetLocations

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
Pt3d XmUGrid::GetPointLocation(int a_pointIdx) const
{
  return m_impl->GetPointLocation(a_pointIdx);
} // XmUGrid::GetPointLocation

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::SetPointLocation(int a_pointIdx, const Pt3d& a_location)
{
  return m_impl->SetPointLocation(a_pointIdx, a_location);
} // XmUGrid::SetPointLocation

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
Pt3d XmUGrid::GetPointXy0(int a_pointIdx) const
{
  return m_impl->GetPointXy0(a_pointIdx);
} // XmUGrid::GetPointXy0

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecPt3d XmUGrid::GetPointsLocations(const VecInt& a_points) const
{
  return m_impl->GetPointsLocations(a_points);
} // XmUGrid::GetPointsLocations

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetExtents(Pt3d& a_min, Pt3d& a_max) const
{
  m_impl->GetExtents(a_min, a_max);
} // XmUGrid::GetExtents

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetPointAdjacentCellCount(int a_pointIdx) const
{
  return m_impl->GetPointAdjacentCellCount(a_pointIdx);
} // XmUGrid::GetPointAdjacentCellCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetPointAdjacentCells(int a_pointIdx) const
{
  return m_impl->GetPointAdjacentCells(a_pointIdx);
} // XmUGrid::GetPointAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetPointAdjacentCells(int a_pointIdx, VecInt& a_adjacentCells) const
{
  m_impl->GetPointAdjacentCells(a_pointIdx, a_adjacentCells);
} // XmUGrid::GetPointAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetPointsAdjacentCells(const VecInt& a_points) const
{
  return m_impl->GetPointsAdjacentCells(a_points);
} // XmUGrid::GetPointsAdjacentCells

//------------------------------------------------------------------------------
void XmUGrid::GetPointsAdjacentCells(
  const VecInt& a_pointIdxs,
  //------------------------------------------------------------------------------
  VecInt& a_adjacentCellIdxs) const
{
  m_impl->GetPointsAdjacentCells(a_pointIdxs, a_adjacentCellIdxs);
} // XmUGrid::GetPointsAdjacentCells

//------------------------------------------------------------------------------
void XmUGrid::GetPointsAdjacentCells(
  int a_pointIdx1,
  int a_pointIdx2,
  //------------------------------------------------------------------------------
  VecInt& a_adjacentCellIdxs) const
{
  m_impl->GetPointsAdjacentCells(a_pointIdx1, a_pointIdx2, a_adjacentCellIdxs);
} // XmUGrid::GetPointsAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::IsValidPointChange(int a_changedPtIdx, const Pt3d& a_newPosition) const
{
  return m_impl->IsValidPointChange(a_changedPtIdx, a_newPosition);
} // XmUGrid::IsValidPointChange

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCellCount() const
{
  return m_impl->GetCellCount();
} // XmUGrid::GetCellCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCellPointCount(int a_cellIdx) const
{
  return m_impl->GetCellPointCount(a_cellIdx);
} // XmUGrid::GetCellPointCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetCellPoints(int a_cellIdx) const
{
  return m_impl->GetCellPoints(a_cellIdx);
} // XmUGrid::GetCellPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::GetCellPoints(int a_cellIdx, VecInt& a_cellPoints) const
{
  return m_impl->GetCellPoints(a_cellIdx, a_cellPoints);
} // XmUGrid::GetCellPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetCellLocations(int a_cellIdx, VecPt3d& a_cellLocations) const
{
  m_impl->GetCellLocations(a_cellIdx, a_cellLocations);
} // XmUGrid::GetCellLocations

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridCellType XmUGrid::GetCellType(int a_cellIdx) const
{
  return m_impl->GetCellType(a_cellIdx);
} // XmUGrid::GetCellType

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::vector<int> XmUGrid::GetDimensionCounts() const
{
  return m_impl->GetDimensionCounts();
} // XmUGrid::GetDimensionCounts

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCellDimension(int a_cellIdx) const
{
  return m_impl->GetCellDimension(a_cellIdx);
} // XmUGrid::GetCellDimension

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetCellExtents(int a_cellIdx, Pt3d& a_min, Pt3d& a_max) const
{
  m_impl->GetCellExtents(a_cellIdx, a_min, a_max);
} // XmUGrid::GetCellExtents

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecInt& XmUGrid::GetCellstream() const
{
  return m_impl->GetCellstream();
} // XmUGrid::GetCellstream

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::SetCellstream(const VecInt& a_cellstream)
{
  return m_impl->SetCellstream(a_cellstream);
} // XmUGrid::SetCellstream

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::GetCellCellstream(int a_cellIdx, VecInt& a_cellstream) const
{
  return m_impl->GetCellCellstream(a_cellIdx, a_cellstream);
} // XmUGrid::GetCellCellstream

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetCellAdjacentCells(int a_cellIdx) const
{
  return m_impl->GetCellAdjacentCells(a_cellIdx);
} // XmUGrid::GetCellAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetCellAdjacentCells(int a_cellIdx, VecInt& a_cellNeighbors) const
{
  m_impl->GetCellAdjacentCells(a_cellIdx, a_cellNeighbors);
} // XmUGrid::GetCellAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::GetCellPlanViewPolygon(int a_cellIdx, VecPt3d& a_polygon) const
{
  return m_impl->GetCellPlanViewPolygon(a_cellIdx, a_polygon);
} // XmUGrid::GetCellPlanViewPolygon

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::GetCellCentroid(int a_cellIdx, Pt3d& a_centroid) const
{
  return m_impl->GetCellCentroid(a_cellIdx, a_centroid);
} // XmUGrid::GetCellCentroid

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCellEdgeCount(int a_cellIdx) const
{
  return m_impl->GetCellEdgeCount(a_cellIdx);
} // XmUGrid::GetCellEdgeCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmEdge XmUGrid::GetCellEdge(int a_cellIdx, int a_edgeIdx) const
{
  return m_impl->GetCellEdge(a_cellIdx, a_edgeIdx);
} // XmUGrid::GetCellEdge

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetCellEdgeAdjacentCells(int a_cellIdx, int a_edgeIdx) const
{
  return m_impl->GetCellEdgeAdjacentCells(a_cellIdx, a_edgeIdx);
} // XmUGrid::GetCellEdgeAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetCellEdgeAdjacentCells(int a_cellIdx,
                                       int a_edgeIdx,
                                       VecInt& a_adjacentCellIdxs) const
{
  m_impl->GetCellEdgeAdjacentCells(a_cellIdx, a_edgeIdx, a_adjacentCellIdxs);
} // XmUGrid::GetCellEdgeAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCell2dEdgeAdjacentCell(int a_cellIdx, int a_edgeIdx) const
{
  return m_impl->GetCell2dEdgeAdjacentCell(a_cellIdx, a_edgeIdx);
} // XmUGrid::GetCell2dEdgeAdjacentCell

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetEdgeAdjacentCells(const XmEdge& a_edge, VecInt& a_adjacentCellIdxs) const
{
  m_impl->GetEdgeAdjacentCells(a_edge, a_adjacentCellIdxs);
} // XmUGrid::GetEdgeAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetEdgeAdjacentCells(const XmEdge& a_edge) const
{
  return m_impl->GetEdgeAdjacentCells(a_edge);
} // XmUGrid::GetEdgeAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::vector<XmEdge> XmUGrid::GetCellEdges(int a_cellIdx) const
{
  return m_impl->GetCellEdges(a_cellIdx);
} // XmUGrid::GetCellEdges

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetCellEdges(int a_cellIdx, std::vector<XmEdge>& a_edges) const
{
  m_impl->GetCellEdges(a_cellIdx, a_edges);
} // XmUGrid::GetCellEdges

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetPointAdjacentPoints(int a_pointIdx, VecInt& a_edgePoints) const
{
  m_impl->GetPointAdjacentPoints(a_pointIdx, a_edgePoints);
} // XmUGrid::GetPointAdjacentPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetPointAdjacentLocations(int a_pointIdx, VecPt3d& a_edgePoints) const
{
  m_impl->GetPointAdjacentLocations(a_pointIdx, a_edgePoints);
} // XmUGrid::GetPointAdjacentLocations

// Faces

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCell3dFaceCount(int a_cellIdx) const
{
  return m_impl->GetCell3dFaceCount(a_cellIdx);
} // XmUGrid::GetCell3dFaceCount
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCell3dFacePointCount(int a_cellIdx, int a_faceIdx) const
{
  return m_impl->GetCell3dFacePointCount(a_cellIdx, a_faceIdx);
} // XmUGrid::GetCell3dFacePointCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt XmUGrid::GetCell3dFacePoints(int a_cellIdx, int a_faceIdx) const
{
  return m_impl->GetCell3dFacePoints(a_cellIdx, a_faceIdx);
} // XmUGrid::GetCell3dFacePoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGrid::GetCell3dFacePoints(int a_cellIdx, int a_faceIdx, VecInt& a_facePtIdxs) const
{
  m_impl->GetCell3dFacePoints(a_cellIdx, a_faceIdx, a_facePtIdxs);
} // XmUGrid::GetCell3dFacePoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt2d XmUGrid::GetCell3dFacesPoints(int a_cellIdx) const
{
  return m_impl->GetCell3dFacesPoints(a_cellIdx);
} // XmUGrid::GetCell3dFacesPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGrid::GetCell3dFaceAdjacentCell(int a_cellIdx, int a_faceIdx) const
{
  return m_impl->GetCell3dFaceAdjacentCell(a_cellIdx, a_faceIdx);
} // XmUGrid::GetCell3dFaceAdjacentCell

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool XmUGrid::GetCell3dFaceAdjacentCell(int a_cellIdx,
                                        int a_faceIdx,
                                        int& a_neighborCell,
                                        int& a_neighborFace) const
{
  return m_impl->GetCell3dFaceAdjacentCell(a_cellIdx, a_faceIdx, a_neighborCell, a_neighborFace);
} // XmUGrid::GetCell3dFaceAdjacentCell

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridFaceOrientation XmUGrid::GetCell3dFaceOrientation(int a_cellIdx, int a_faceIdx) const
{
  return m_impl->GetCell3dFaceOrientation(a_cellIdx, a_faceIdx);
} // XmUGrid::GetCell3dFaceOrientation

//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> TEST_XmUGrid1Left90Tri()
{
  VecPt3d points = {{0, 0, 0}, {20, 0, 0}, {0, 20, 0}};

  // Cell type (5), number of points (3), point numbers, counterclockwise
  std::vector<int> cells = {5, 3, 0, 1, 2};

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);

  return ugrid;
} // TEST_UGrid1Left90Tri
//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> TEST_XmUGridSimpleQuad()
{
  VecPt3d points = {{0, 10, 0}, {10, 10, 0}, {20, 10, 0},  {0, 0, 0},   {10, 0, 0},
                    {20, 0, 0}, {0, -10, 0}, {10, -10, 0}, {20, -10, 0}};

  // Cell type (9), number of points (4), point numbers, counterclockwise
  VecInt cellstream = {XMU_QUAD, 4, 0, 3, 4, 1, XMU_QUAD, 4, 1, 4, 5, 2,
                       XMU_QUAD, 4, 3, 6, 7, 4, XMU_QUAD, 4, 4, 7, 8, 5};

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cellstream);
  return ugrid;
} // TEST_XmUGridSimpleQuad

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> TEST_XmUGrid2dLinear()
{
  // clang-format off
  VecPt3d points = { { 0, 0, 0 }, { 10, 0, 0 }, { 20, 0, 0 }, { 30, 0, 0 }, { 40, 0, 0 },
                    {0, 10, 0},  {10, 10, 0}, {20, 10, 0}, {40, 10, 0}, {0, 20, 0},
                    {10, 20, 0}, {20, 20, 0}, {30, 20, 0}, {40, 20, 0}};

  // Cell type (5), number of points (3), point numbers, counterclockwise
  std::vector<int> cells = {(int)XMU_QUAD, 4, 0, 1, 6, 5,            // 0
                            (int)XMU_PIXEL, 4, 1, 2, 6, 7,           // 1
                            (int)XMU_TRIANGLE, 3, 2, 3, 7,           // 2
                            (int)XMU_POLYGON, 6, 3, 4, 8, 13, 12, 7, // 3
                            (int)XMU_POLY_LINE, 3, 7, 11, 10,        // 4
                            (int)XMU_LINE, 2, 5, 9};                 // 5
  // clang-format on

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);

  return ugrid;
} // TEST_XmUGrid2dLinear
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> TEST_XmUGrid3dLinear()
{
  // clang-format off
  VecPt3d points = {{0, 0, 0},   {10, 0, 0},   {20, 0, 0},   {30, 0, 0},   {40, 0, 0},
                    {0, 10, 0},  {10, 10, 0},  {20, 10, 0},  {30, 10, 0},  {40, 10, 0},
                    {0, 20, 0},  {10, 20, 0},  {20, 20, 0},  {30, 20, 0},  {40, 20, 0},
                    {0, 0, 10},  {10, 0, 10},  {20, 0, 10},  {30, 0, 10},  {40, 0, 10},
                    {0, 10, 10}, {10, 10, 10}, {20, 10, 10}, {30, 10, 10}, {40, 10, 10},
                    {0, 20, 10}, {10, 20, 10}, {20, 20, 10}, {30, 20, 10}, {40, 20, 10}};

  // Cell type (5), number of points (3), point numbers, counterclockwise
  std::vector<int> cells = {(int)XMU_TETRA, 4, 0, 1, 5, 15,
                            (int)XMU_VOXEL, 8, 1, 2, 6, 7, 16, 17, 21, 22,
                            (int)XMU_HEXAHEDRON, 8, 2, 3, 8, 7, 17, 18, 23, 22,
                            (int)XMU_POLYHEDRON, 6, // A polyhedron with 6 faces
                            4, 9, 8, 13, 14, // Bottom face with 4 points : 9, 8, 13, 14
                            4, 8, 9, 24, 23, // Front face with 4 points : 8, 9, 24, 23
                            4, 9, 14, 29, 24, // Right face with 4 points : 9, 14, 29, 28
                            4, 14, 13, 28, 29, // Back face with 4 points : 14, 13, 28, 29
                            4, 8, 13, 28, 23, // Left face with 4 points : 13, 8, 23, 28
                            4, 23, 24, 29, 28, // Top face with 4 points : 23, 24, 29, 28
                            (int)XMU_WEDGE, 6, 3, 4, 18, 8, 9, 23,
                            (int)XMU_PYRAMID, 5, 5, 6, 11, 10, 20};
  // clang-format on

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);

  return ugrid;
} // TEST_XmUGrid3dLinear
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<XmUGrid> TEST_XmUGridHexagonalPolyhedron()
{
  VecPt3d points = {
    {0, 0, 10}, {10, 0, 10}, {10, 10, 10}, {0, 10, 10},
    {0, 0, 0},  {10, 0, 0},  {10, 10, 0},  {0, 10, 0},
  };

  // clang-format off

  // Cell type (5), number of points (3), point numbers, counterclockwise
  std::vector<int> cells = {(int)XMU_POLYHEDRON, 6,
                            4, 0, 1, 2, 3,
                            4, 4, 5, 7, 2,
                            4, 5, 6, 2, 1,
                            4, 6, 7, 3, 2,
                            4, 7, 4, 0, 3,
                            4, 4, 7, 6, 5};
  // clang-format on

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);

  return ugrid;
} // TEST_XmUGridHexagonalPolyhedron
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuildQuadUGrid(int a_rows, int a_cols)
{
  Pt3d origin(0.0, 0.0, 0.0);
  return TEST_XmUBuildQuadUGrid(a_rows, a_cols, origin);
} // TEST_XmUBuildQuadUGrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuildQuadUGrid(int a_rows, int a_cols, const xms::Pt3d& a_origin)
{
  VecPt3d points;
  points.reserve(a_rows * a_cols);
  for (int i = 0; i < a_rows; ++i)
  {
    for (int j = 0; j < a_cols; ++j)
    {
      points.push_back(Pt3d(j + a_origin.x, a_rows - i + a_origin.y));
    }
  }

  VecInt cells;
  cells.reserve((a_rows - 1) * (a_cols - 1) * 6);
  for (int i = 0; i < a_rows - 1; ++i)
  {
    for (int j = 0; j < a_cols - 1; ++j)
    {
      cells.push_back(XMU_QUAD);
      cells.push_back(4);
      cells.push_back(j + a_cols * i);
      cells.push_back(j + a_cols * (i + 1));
      cells.push_back(j + 1 + a_cols * (i + 1));
      cells.push_back(j + 1 + a_cols * i);
    }
  }

  return XmUGrid::New(points, cells);
} // TEST_XmUBuildQuadUGrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuildHexahedronUgrid(int a_rows, int a_cols, int a_lays)
{
  Pt3d origin(0.0, 0.0, 0.0);
  return TEST_XmUBuildHexahedronUgrid(a_rows, a_cols, a_lays, origin);
} // TEST_XmUBuildHexahedronUgrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuildHexahedronUgrid(int a_rows,
                                                int a_cols,
                                                int a_lays,
                                                const xms::Pt3d& a_origin)
{
  VecPt3d points;
  points.reserve(a_rows * a_cols);
  for (int k = 0; k < a_lays; ++k)
  {
    for (int i = 0; i < a_rows; ++i)
    {
      for (int j = 0; j < a_cols; ++j)
      {
        points.push_back(
          Pt3d(j + a_origin.x, a_rows - i - 1 + a_origin.y, a_lays - k - 1 + a_origin.z));
      }
    }
  }

  VecInt cells;
  cells.reserve((a_rows - 1) * (a_cols - 1) * (a_lays - 1) * 10);
  int numInLayer = a_rows * a_cols;
  for (int k = 0; k < a_lays - 1; ++k)
  {
    int layOffset = numInLayer * k;
    for (int i = 0; i < a_rows - 1; ++i)
    {
      for (int j = 0; j < a_cols - 1; ++j)
      {
        cells.push_back(XMU_HEXAHEDRON);
        cells.push_back(8);

        cells.push_back(j + a_cols * i + layOffset);
        cells.push_back(j + 1 + a_cols * i + layOffset);
        cells.push_back(j + 1 + a_cols * (i + 1) + layOffset);
        cells.push_back(j + a_cols * (i + 1) + layOffset);

        cells.push_back(j + a_cols * i + layOffset + numInLayer);
        cells.push_back(j + 1 + a_cols * i + layOffset + numInLayer);
        cells.push_back(j + 1 + a_cols * (i + 1) + layOffset + numInLayer);
        cells.push_back(j + a_cols * (i + 1) + layOffset + numInLayer);
      }
    }
  }
  return XmUGrid::New(points, cells);
} // TEST_XmUBuildHexahedronUgrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuildPolyhedronUgrid(int a_rows, int a_cols, int a_lays)
{
  Pt3d origin(0.0, 0.0, 0.0);
  return TEST_XmUBuildPolyhedronUgrid(a_rows, a_cols, a_lays, origin);
} // TEST_XmUBuildPolyhedronUgrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuildPolyhedronUgrid(int a_rows,
                                                int a_cols,
                                                int a_lays,
                                                const xms::Pt3d& a_origin)
{
  VecPt3d points;
  points.reserve(a_rows * a_cols);
  for (int k = 0; k < a_lays; ++k)
  {
    for (int i = 0; i < a_rows; ++i)
    {
      for (int j = 0; j < a_cols; ++j)
      {
        points.push_back(Pt3d(j + a_origin.x, a_rows - i + a_origin.y, a_lays - k + a_origin.z));
      }
    }
  }

  VecInt cells;
  cells.reserve((a_rows - 1) * (a_cols - 1) * (a_lays - 1) * 10);
  int numInLayer = a_rows * a_cols;
  for (int k = 0; k < a_lays - 1; ++k)
  {
    int layOffset = numInLayer * k;
    for (int i = 0; i < a_rows - 1; ++i)
    {
      for (int j = 0; j < a_cols - 1; ++j)
      {
        int pt0 = j + a_cols * i + layOffset;
        int pt1 = j + 1 + a_cols * i + layOffset;
        int pt2 = j + a_cols * (i + 1) + layOffset;
        int pt3 = j + 1 + a_cols * (i + 1) + layOffset;
        int pt4 = pt0 + numInLayer;
        int pt5 = pt1 + numInLayer;
        int pt6 = pt2 + numInLayer;
        int pt7 = pt3 + numInLayer;

        cells.push_back(XMU_POLYHEDRON);
        cells.push_back(6);
        // front
        cells.push_back(4);
        cells.push_back(pt2);
        cells.push_back(pt6);
        cells.push_back(pt7);
        cells.push_back(pt3);
        // right
        cells.push_back(4);
        cells.push_back(pt1);
        cells.push_back(pt5);
        cells.push_back(pt7);
        cells.push_back(pt3);
        // back
        cells.push_back(4);
        cells.push_back(pt0);
        cells.push_back(pt1);
        cells.push_back(pt5);
        cells.push_back(pt4);
        // left
        cells.push_back(4);
        cells.push_back(pt0);
        cells.push_back(pt4);
        cells.push_back(pt6);
        cells.push_back(pt2);
        // top
        cells.push_back(4);
        cells.push_back(pt0);
        cells.push_back(pt2);
        cells.push_back(pt3);
        cells.push_back(pt1);
        // bottom
        cells.push_back(4);
        cells.push_back(pt4);
        cells.push_back(pt5);
        cells.push_back(pt7);
        cells.push_back(pt6);
      }
    }
  }
  return XmUGrid::New(points, cells);
} // TEST_XmUBuildPolyhedronUgrid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::shared_ptr<xms::XmUGrid> TEST_XmUBuild3DChevronUgrid()
{
  VecPt3d ugridPoints = {{0.0, 0.0, 0.0},   {20.0, 10.0, 0.0}, {40.0, 0.0, 0.0},
                         {20.0, 50.0, 0.0}, {0.0, 0.0, 10.0},  {20.0, 10.0, 10.0},
                         {40.0, 0.0, 10.0}, {20.0, 50.0, 10.0}};
  // clang-format off
  VecInt ugridCell = { (int)XMU_POLYHEDRON, 6, 
                      4, 0, 3, 2, 1,
                      4, 4, 5, 6, 7,
                      4, 0, 1, 5, 4,
                      4, 1, 2, 6, 5,
                      4, 2, 3, 7, 6,
                      4, 3, 0, 4, 7};
  // clang-format on

  return XmUGrid::New(ugridPoints, ugridCell);
} // TEST_XmUBuild3DChevronUgrid

} // namespace xms

#ifdef CXX_TEST
//------------------------------------------------------------------------------
// Unit Tests
//------------------------------------------------------------------------------
#include <xmsgrid/ugrid/XmUGrid.t.h>

#include <xmscore/testing/TestTools.h>

using namespace xms;


//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testOperators()
{
  std::shared_ptr<XmUGrid> xmUGrid = TEST_XmUGrid1Left90Tri();

  // test copy constructor
  XmUGrid grid1(*xmUGrid);
  TS_ASSERT_EQUALS(xmUGrid->GetLocations(), grid1.GetLocations());

  // test move constructor
  XmUGrid grid2(std::move(grid1));
  TS_ASSERT(grid1.GetLocations().empty());
  TS_ASSERT_EQUALS(xmUGrid->GetLocations(), grid2.GetLocations());

  // test assignment operator
  grid1 = grid2;
  TS_ASSERT_EQUALS(xmUGrid->GetLocations(), grid1.GetLocations());
  TS_ASSERT_EQUALS(xmUGrid->GetLocations(), grid2.GetLocations());
} // XmUGridUnitTests::testOperators

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testUGridStreams()
{
  // test empty UGrid
  std::shared_ptr<XmUGrid> emptyUGrid = XmUGrid::New();
  VecPt3d points = emptyUGrid->GetLocations();
  TS_ASSERT_EQUALS(0, points.size());

  VecInt cellstream = emptyUGrid->GetCellstream();
  TS_ASSERT(XmUGrid::IsValidCellstream(cellstream));
  TS_ASSERT_EQUALS(0, cellstream.size());

  // test adding points and cell stream

  //     0----1----2
  //     |    |    |
  //     3----4----5
  //     |    |    |
  //     6----7----8

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New();

  points = {{0, 10, 0}, {10, 10, 0}, {20, 10, 0},  {0, 0, 0},   {10, 0, 0},
            {20, 0, 0}, {0, -10, 0}, {10, -10, 0}, {20, -10, 0}};

  // Cell type (9), number of points (4), point numbers, counterclockwise
  cellstream = {9, 4, 0, 3, 4, 1, 9, 4, 1, 4, 5, 2, 9, 4, 3, 6, 7, 4, 9, 4, 4, 7, 8, 5};

  TS_ASSERT(!ugrid->GetModified());
  ugrid->SetLocations(points);
  VecPt3d pointsOut = ugrid->GetLocations();
  TS_ASSERT_EQUALS(points, pointsOut);
  TS_ASSERT(ugrid->GetModified());

  ugrid->SetUnmodified();
  TS_ASSERT(!ugrid->GetModified());
  TS_ASSERT(ugrid->SetCellstream(cellstream));
  VecInt cellstreamOut = ugrid->GetCellstream();
  TS_ASSERT_EQUALS(cellstream, cellstreamOut);
  TS_ASSERT(ugrid->GetModified());

  // Test invalid cell streams
  cellstream = {-1};
  TS_ASSERT(!XmUGrid::IsValidCellstream(cellstream));
  cellstream = {9, 4, 0, 3, 4};
  TS_ASSERT(!XmUGrid::IsValidCellstream(cellstream));
  cellstream = {9, 3, 0, 3, 4, 1};
  TS_ASSERT(!XmUGrid::IsValidCellstream(cellstream));

  // Test adding cellstream then points (should fail)
  std::shared_ptr<XmUGrid> ugridBadOrder = XmUGrid::New();
  TS_ASSERT(!ugridBadOrder->SetCellstream(cellstream));

} // XmUGridUnitTests::testUGridStreams

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetSetPoint()
{
  //     0----1----2
  //     |    |    |
  //     3----4----5
  //     |    |    |
  //     6----7----8

  VecPt3d points = {{0, 10, 0}, {10, 10, 10}, {20, 10, 0},  {0, 0, 0},   {10, 0, 0},
                    {20, 0, 0}, {0, -10, 0},  {10, -10, 0}, {20, -10, 0}};

  // Cell type (9), number of points (4), point numbers, counterclockwise
  VecInt cellstream = {9, 4, 0, 3, 4, 1, 9, 4, 1, 4, 5, 2, 9, 4, 3, 6, 7, 4, 9, 4, 4, 7, 8, 5};

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cellstream);

  TS_ASSERT_EQUALS(Pt3d(), ugrid->GetPointLocation(-1));
  for (int i = 0; i < (int)points.size(); i++)
  {
    TS_ASSERT_EQUALS(points[i], ugrid->GetPointLocation(i));
  }
  TS_ASSERT_EQUALS(Pt3d(10, 10, 0), ugrid->GetPointXy0(1));
  TS_ASSERT_EQUALS(Pt3d(), ugrid->GetPointLocation((int)points.size()));
  TS_ASSERT(ugrid->GetModified());

  ugrid->SetUnmodified();
  TS_ASSERT(!ugrid->GetModified());
  TS_ASSERT(!ugrid->SetPointLocation(-1, Pt3d()));
  TS_ASSERT(!ugrid->SetPointLocation((int)points.size(), Pt3d()));
  TS_ASSERT(!ugrid->GetModified());

  TS_ASSERT(ugrid->SetPointLocation(0, Pt3d(-10, 10, 0)));
  TS_ASSERT_EQUALS(Pt3d(-10, 10, 0), ugrid->GetPointLocation(0));
  TS_ASSERT(ugrid->GetModified());
} // XmUGridUnitTests::testGetSetPoint
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellCellstream()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  VecPt3d points = {{0, 10, 0}, {10, 10, 0}, {20, 10, 0},  {0, 0, 0},   {10, 0, 0},
                    {20, 0, 0}, {0, -10, 0}, {10, -10, 0}, {20, -10, 0}};

  VecInt cell0 = {9, 4, 0, 3, 4, 1}, cell1 = {9, 4, 1, 4, 5, 2}, cell2 = {9, 4, 3, 6, 7, 4},
         cell3 = {9, 4, 4, 7, 8, 5};
  // Cell type (9), number of points (4), point numbers, counterclockwise
  VecInt cellstream;
  cellstream.insert(cellstream.end(), cell0.begin(), cell0.end());
  cellstream.insert(cellstream.end(), cell1.begin(), cell1.end());
  cellstream.insert(cellstream.end(), cell2.begin(), cell2.end());
  cellstream.insert(cellstream.end(), cell3.begin(), cell3.end());

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cellstream);
  VecInt cellResult;

  TS_ASSERT_EQUALS(false, ugrid->GetCellCellstream(-1, cellResult));

  TS_ASSERT_EQUALS(true, ugrid->GetCellCellstream(0, cellResult));
  TS_ASSERT_EQUALS(cell0, cellResult);
  TS_ASSERT_EQUALS(true, ugrid->GetCellCellstream(1, cellResult));
  TS_ASSERT_EQUALS(cell1, cellResult);
  TS_ASSERT_EQUALS(true, ugrid->GetCellCellstream(2, cellResult));
  TS_ASSERT_EQUALS(cell2, cellResult);
  TS_ASSERT_EQUALS(true, ugrid->GetCellCellstream(3, cellResult));
  TS_ASSERT_EQUALS(cell3, cellResult);

  TS_ASSERT_EQUALS(false, ugrid->GetCellCellstream((int)points.size(), cellResult));

} // XmUGridUnitTests::testGetCellCellstream
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellType()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);
  TS_ASSERT_EQUALS(14, ugrid2d->GetPointCount());
  TS_ASSERT_EQUALS(6, ugrid2d->GetCellCount());
  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid2d->GetCellType(-1));
  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid2d->GetCellType(6));
  TS_ASSERT_EQUALS(XMU_QUAD, ugrid2d->GetCellType(0));
  TS_ASSERT_EQUALS(XMU_PIXEL, ugrid2d->GetCellType(1));
  TS_ASSERT_EQUALS(XMU_TRIANGLE, ugrid2d->GetCellType(2));
  TS_ASSERT_EQUALS(XMU_POLYGON, ugrid2d->GetCellType(3));
  TS_ASSERT_EQUALS(XMU_POLY_LINE, ugrid2d->GetCellType(4));
  TS_ASSERT_EQUALS(XMU_LINE, ugrid2d->GetCellType(5));

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  TS_ASSERT_EQUALS(30, ugrid3d->GetPointCount());
  TS_ASSERT_EQUALS(6, ugrid3d->GetCellCount());
  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid3d->GetCellType(-1));
  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid3d->GetCellType(6));
  TS_ASSERT_EQUALS(XMU_TETRA, ugrid3d->GetCellType(0));
  TS_ASSERT_EQUALS(XMU_VOXEL, ugrid3d->GetCellType(1));
  TS_ASSERT_EQUALS(XMU_HEXAHEDRON, ugrid3d->GetCellType(2));
  TS_ASSERT_EQUALS(XMU_POLYHEDRON, ugrid3d->GetCellType(3));
  TS_ASSERT_EQUALS(XMU_WEDGE, ugrid3d->GetCellType(4));
  TS_ASSERT_EQUALS(XMU_PYRAMID, ugrid3d->GetCellType(5));
} // XmUGridUnitTests::testGetCellType
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellDimension()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid2d->GetCellDimension(-1));
  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid2d->GetCellDimension(6));
  TS_ASSERT_EQUALS(2, ugrid2d->GetCellDimension(0));
  TS_ASSERT_EQUALS(2, ugrid2d->GetCellDimension(1));
  TS_ASSERT_EQUALS(2, ugrid2d->GetCellDimension(2));
  TS_ASSERT_EQUALS(2, ugrid2d->GetCellDimension(3));
  TS_ASSERT_EQUALS(1, ugrid2d->GetCellDimension(4));
  TS_ASSERT_EQUALS(1, ugrid2d->GetCellDimension(5));

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid3d->GetCellDimension(-1));
  TS_ASSERT_EQUALS(XMU_INVALID_CELL_TYPE, ugrid3d->GetCellDimension(6));
  TS_ASSERT_EQUALS(3, ugrid3d->GetCellDimension(0));
  TS_ASSERT_EQUALS(3, ugrid3d->GetCellDimension(1));
  TS_ASSERT_EQUALS(3, ugrid3d->GetCellDimension(2));
  TS_ASSERT_EQUALS(3, ugrid3d->GetCellDimension(3));
  TS_ASSERT_EQUALS(3, ugrid3d->GetCellDimension(4));
  TS_ASSERT_EQUALS(3, ugrid3d->GetCellDimension(5));

  // Test Dimension Counts
  std::vector<int> twoDResults(4, 0), threeDResults(4, 0);
  twoDResults[1] = 2;
  twoDResults[2] = 4;
  threeDResults[3] = 6;
  TS_ASSERT_EQUALS(twoDResults, ugrid2d->GetDimensionCounts());
  TS_ASSERT_EQUALS(threeDResults, ugrid3d->GetDimensionCounts());

  // Test GetDimensionCounts with an empty grid
  std::shared_ptr<XmUGrid> emptyUgrid = XmUGrid::New();
  std::vector<int> emptyResults = emptyUgrid->GetDimensionCounts();
  TS_ASSERT_EQUALS(VecInt(4, 0), emptyResults);
} // XmUGridUnitTests::testGetCellDimension
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetExtents()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  Pt3d min, max;
  Pt3d expectedMin = {0.0, 0.0, 0.0}, expectMax = {40.0, 20.0, 0.0};
  ugrid2d->GetExtents(min, max);
  TS_ASSERT_EQUALS(expectedMin, min);
  TS_ASSERT_EQUALS(expectMax, max);

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  expectedMin = {0.0, 0.0, 0.0};
  expectMax = {40.0, 20.0, 10.0};
  ugrid3d->GetExtents(min, max);
  TS_ASSERT_EQUALS(expectedMin, min);
  TS_ASSERT_EQUALS(expectMax, max);

  Pt3d origin(-25.0, -25.0, -15.0);
  std::shared_ptr<XmUGrid> ugridBuild = TEST_XmUBuildHexahedronUgrid(51, 51, 31, origin);
  TS_REQUIRE_NOT_NULL(ugridBuild);

  expectedMin = {-25.0, -25.0, -15.0};
  expectMax = {25.0, 25.0, 15.0};
  ugridBuild->GetExtents(min, max);
  TS_ASSERT_EQUALS(expectedMin, min);
  TS_ASSERT_EQUALS(expectMax, max);
} // XmUGridUnitTests::testGetExtents
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellEdgeCount()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  TS_ASSERT_EQUALS(-1, ugrid2d->GetCellEdgeCount(-1));
  TS_ASSERT_EQUALS(-1, ugrid2d->GetCellEdgeCount(6));
  TS_ASSERT_EQUALS(4, ugrid2d->GetCellEdgeCount(0));
  TS_ASSERT_EQUALS(4, ugrid2d->GetCellEdgeCount(1));
  TS_ASSERT_EQUALS(3, ugrid2d->GetCellEdgeCount(2));
  TS_ASSERT_EQUALS(6, ugrid2d->GetCellEdgeCount(3));
  TS_ASSERT_EQUALS(2, ugrid2d->GetCellEdgeCount(4));
  TS_ASSERT_EQUALS(1, ugrid2d->GetCellEdgeCount(5));

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  TS_ASSERT_EQUALS(-1, ugrid3d->GetCellEdgeCount(-1));
  TS_ASSERT_EQUALS(-1, ugrid3d->GetCellEdgeCount(6));
  TS_ASSERT_EQUALS(6, ugrid3d->GetCellEdgeCount(0));
  TS_ASSERT_EQUALS(12, ugrid3d->GetCellEdgeCount(1));
  TS_ASSERT_EQUALS(12, ugrid3d->GetCellEdgeCount(2));
  TS_ASSERT_EQUALS(12, ugrid3d->GetCellEdgeCount(3));
  TS_ASSERT_EQUALS(9, ugrid3d->GetCellEdgeCount(4));
  TS_ASSERT_EQUALS(8, ugrid3d->GetCellEdgeCount(5));

} // XmUGridUnitTests::testGetCellEdgeCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCell3dFaceCount()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  // test 1D and 2D cells
  TS_ASSERT_EQUALS(-1, ugrid2d->GetCell3dFaceCount(-1));
  for (int cellIdx = 0; cellIdx < ugrid2d->GetCellCount(); ++cellIdx)
  {
    TS_ASSERT_EQUALS(0, ugrid2d->GetCell3dFaceCount(cellIdx));
  }
  TS_ASSERT_EQUALS(-1, ugrid2d->GetCell3dFaceCount(ugrid2d->GetCellCount()));

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  // test 3D cells
  TS_ASSERT_EQUALS(-1, ugrid3d->GetCell3dFaceCount(-1));
  TS_ASSERT_EQUALS(-1, ugrid3d->GetCell3dFaceCount(6));
  TS_ASSERT_EQUALS(4, ugrid3d->GetCell3dFaceCount(0));
  TS_ASSERT_EQUALS(6, ugrid3d->GetCell3dFaceCount(1));
  TS_ASSERT_EQUALS(6, ugrid3d->GetCell3dFaceCount(2));
  TS_ASSERT_EQUALS(6, ugrid3d->GetCell3dFaceCount(3));
  TS_ASSERT_EQUALS(5, ugrid3d->GetCell3dFaceCount(4));
  TS_ASSERT_EQUALS(5, ugrid3d->GetCell3dFaceCount(5));
} // XmUGridUnitTests::testGetCell3dFaceCount

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCell3dFacePointCount()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  // test 1D and 2D cells - should return 0 for valid cell index
  TS_ASSERT_EQUALS(-1, ugrid2d->GetCell3dFacePointCount(-1, 0));
  TS_ASSERT_EQUALS(0, ugrid2d->GetCell3dFacePointCount(1, 0));
  TS_ASSERT_EQUALS(-1, ugrid2d->GetCell3dFacePointCount(ugrid2d->GetCellCount(), 0));

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  // test 3D cells
  TS_ASSERT_EQUALS(-1, ugrid3d->GetCell3dFacePointCount(-1, 0));
  TS_ASSERT_EQUALS(-1, ugrid3d->GetCell3dFacePointCount(6, 0));
  TS_ASSERT_EQUALS(4, ugrid3d->GetCell3dFacePointCount(1, 3));
  TS_ASSERT_EQUALS(4, ugrid3d->GetCell3dFacePointCount(3, 2));
} // XmUGridUnitTests::testGetCell3dFacePointCount
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetPointAdjacentCellsSimple()
{
  //     0----1----2
  //     |    |    |
  //     3----4----5

  VecPt3d points = {{0, 10, 0}, {10, 10, 0}, {20, 10, 0}, {0, 0, 0}, {10, 0, 0}, {20, 0, 0}};

  VecInt cellstream = {XMU_QUAD, 4, 0, 3, 4, 1, XMU_QUAD, 4, 1, 4, 5, 2};

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cellstream);

  VecInt cellEmpty = {};
  VecInt cellZero = {0};
  VecInt cellZeroAndOne = {0, 1};
  VecInt cellOne = {1};

  TS_ASSERT_EQUALS(cellEmpty, ugrid->GetPointAdjacentCells(-1));
  TS_ASSERT_EQUALS(cellZero, ugrid->GetPointAdjacentCells(0));
  TS_ASSERT_EQUALS(cellZeroAndOne, ugrid->GetPointAdjacentCells(1));
  TS_ASSERT_EQUALS(cellOne, ugrid->GetPointAdjacentCells(2));
  TS_ASSERT_EQUALS(cellZero, ugrid->GetPointAdjacentCells(3));
  TS_ASSERT_EQUALS(cellZeroAndOne, ugrid->GetPointAdjacentCells(4));
  TS_ASSERT_EQUALS(cellOne, ugrid->GetPointAdjacentCells(5));
  TS_ASSERT_EQUALS(cellEmpty, ugrid->GetPointAdjacentCells(6));

  TS_ASSERT_EQUALS(0, ugrid->GetPointAdjacentCellCount(-1));
  TS_ASSERT_EQUALS(1, ugrid->GetPointAdjacentCellCount(0));
  TS_ASSERT_EQUALS(2, ugrid->GetPointAdjacentCellCount(1));
  TS_ASSERT_EQUALS(1, ugrid->GetPointAdjacentCellCount(2));
  TS_ASSERT_EQUALS(1, ugrid->GetPointAdjacentCellCount(3));
  TS_ASSERT_EQUALS(2, ugrid->GetPointAdjacentCellCount(4));
  TS_ASSERT_EQUALS(1, ugrid->GetPointAdjacentCellCount(5));
  TS_ASSERT_EQUALS(0, ugrid->GetPointAdjacentCellCount(6));

} // XmUGridUnitTests::testGetPointAdjacentCellsSimple

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetPointAdjacentCells()
{
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  VecInt2d cellsFor2DPoints = {{0},          // point 0
                               {0, 1},       // point 1
                               {1, 2},       // point 2
                               {2, 3},       // point 3
                               {3},          // point 4
                               {0, 5},       // point 5
                               {0, 1},       // point 6
                               {1, 2, 3, 4}, // point 7
                               {3},          // point 8
                               {5},          // point 9
                               {4},          // point 10
                               {4},          // point 11
                               {3},          // point 12
                               {3}};         // point 13

  // test 1D and 2D cells
  for (int i(0); i < cellsFor2DPoints.size(); i++)
  {
    TS_ASSERT_EQUALS(cellsFor2DPoints[i], ugrid2d->GetPointAdjacentCells(i));
  }

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  VecInt2d cellsFor3DPoints = {{0},       // point 0
                               {0, 1},    // point 1
                               {1, 2},    // point 2
                               {2, 4},    // point 3
                               {4},       // point 4
                               {0, 5},    // point 5
                               {1, 5},    // point 6
                               {1, 2},    // point 7
                               {2, 3, 4}, // point 8
                               {3, 4},    // point 9
                               {5},       // point 10
                               {5},       // point 11
                               {},        // point 12
                               {3},       // point 13
                               {3},       // point 14
                               {0},       // point 15
                               {1},       // point 16
                               {1, 2},    // point 17
                               {2, 4},    // point 18
                               {},        // point 19
                               {5},       // point 20
                               {1},       // point 21
                               {1, 2},    // point 22
                               {2, 3, 4}, // point 23
                               {3},       // point 24
                               {},        // point 25
                               {},        // point 26
                               {},        // point 27
                               {3},       // point 28
                               {3},       // point 29
                               {}};       // point 30

  // test 3D cells
  for (int i(0); i < cellsFor3DPoints.size(); i++)
  {
    if (cellsFor3DPoints[i] != ugrid3d->GetPointAdjacentCells(i))
    {
      TS_FAIL(i);
    }
    TS_ASSERT_EQUALS(cellsFor3DPoints[i], ugrid3d->GetPointAdjacentCells(i));
  }

} // XmUGridUnitTests::testGetPointAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellPoints()
{
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  VecInt2d expectedGetCellPoints = {{0, 3, 4, 1}, {1, 4, 5, 2}, {3, 6, 7, 4}, {4, 7, 8, 5}};
  VecInt2d cellPoints;
  VecInt cellPoint1D;
  TS_ASSERT(!ugrid->GetCellPoints(-1, cellPoint1D));
  TS_ASSERT(!ugrid->GetCellPoints(5, cellPoint1D));

  for (int i(0); i < ugrid->GetCellCount(); i++)
  {
    TS_ASSERT(ugrid->GetCellPoints(i, cellPoint1D));
    cellPoints.push_back(cellPoint1D);
    TS_ASSERT_EQUALS(expectedGetCellPoints[i], cellPoints[i]);
  }

  expectedGetCellPoints = {{0, 1, 6, 5},         {1, 2, 7, 6}, {2, 3, 7},
                           {3, 4, 8, 13, 12, 7}, {7, 11, 10},  {5, 9}};
  cellPoints.clear();
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  for (int i(0); i < ugrid2d->GetCellCount(); i++)
  {
    TS_ASSERT(ugrid2d->GetCellPoints(i, cellPoint1D));
    cellPoints.push_back(cellPoint1D);
    TS_ASSERT_EQUALS(expectedGetCellPoints[i], cellPoints[i]);
  }

  expectedGetCellPoints = {{0, 1, 5, 15},
                           {1, 2, 6, 7, 16, 17, 21, 22},
                           {2, 3, 8, 7, 17, 18, 23, 22},
                           {9, 8, 13, 14, 24, 23, 29, 28},
                           {3, 4, 18, 8, 9, 23},
                           {5, 6, 11, 10, 20}};
  VecInt expectedPointCounts = {4, 8, 8, 8, 6, 5};
  cellPoints.clear();
  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  for (int i(0); i < ugrid3d->GetCellCount(); i++)
  {
    TS_ASSERT(ugrid3d->GetCellPoints(i, cellPoint1D));
    cellPoints.push_back(cellPoint1D);
    TS_ASSERT_EQUALS(expectedGetCellPoints[i], cellPoints[i]);
    TS_ASSERT_EQUALS(expectedPointCounts[i], ugrid3d->GetCellPointCount(i));
  }

} // XmUGridUnitTests::testGetCellPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellEdge()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  // 2D shape (simple)
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  {
    XmEdge edge0(0, 3), edge1(3, 4), edge2(4, 1), edge3(1, 0);     // cell 0
    XmEdge edge4(1, 4), edge5(4, 5), edge6(5, 2), edge7(2, 1);     // cell 1
    XmEdge edge8(3, 6), edge9(6, 7), edge10(7, 4), edge11(4, 3);   // cell 2
    XmEdge edge12(4, 7), edge13(7, 8), edge14(8, 5), edge15(5, 4); // cell 3
    XmEdge emptyEdge;
    std::vector<std::vector<XmEdge>> expectedCellsCellEdges;
    std::vector<XmEdge> emptyVec = {emptyEdge};
    std::vector<XmEdge> cellEdges;

    cellEdges.push_back(ugrid->GetCellEdge(0, -1));
    TS_ASSERT_EQUALS(emptyVec, cellEdges);

    cellEdges.clear();
    cellEdges.push_back(ugrid->GetCellEdge(0, ugrid->GetCellEdgeCount(0)));
    TS_ASSERT_EQUALS(emptyVec, cellEdges);

    cellEdges.clear();
    expectedCellsCellEdges = {{edge0, edge1, edge2, edge3},
                              {edge4, edge5, edge6, edge7},
                              {edge8, edge9, edge10, edge11},
                              {edge12, edge13, edge14, edge15}};
    for (int i(0); i < ugrid->GetCellCount(); i++)
    {
      for (int j(0); j < ugrid->GetCellEdgeCount(i); j++)
      {
        cellEdges.push_back(ugrid->GetCellEdge(i, j));
      }
      TS_ASSERT_EQUALS(expectedCellsCellEdges[i], cellEdges);
      cellEdges.clear();
    }
  }
  {
    XmEdge edge0(0, 1), edge1(1, 6), edge2(6, 5), edge3(5, 0); // Quad
    XmEdge edge4(1, 2), edge5(2, 7), edge6(7, 6), edge7(6, 1); // Pixel
    XmEdge edge8(2, 3), edge9(3, 7), edge10(7, 2);             // Triangle
    XmEdge edge11(3, 4), edge12(4, 8), edge13(8, 13), edge14(13, 12), edge15(12, 7),
      edge16(7, 3);                       // Polygon
    XmEdge edge17(7, 11), edge18(11, 10); // PolyLine
    XmEdge edge19(5, 9);                  // Line

    // 2D Shapes
    std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
    TS_REQUIRE_NOT_NULL(ugrid2d);

    std::vector<std::vector<XmEdge>> expectedCellsCellEdges;
    expectedCellsCellEdges = {{edge0, edge1, edge2, edge3},
                              {edge4, edge5, edge6, edge7},
                              {edge8, edge9, edge10},
                              {edge11, edge12, edge13, edge14, edge15, edge16},
                              {edge17, edge18},
                              {edge19}};
    std::vector<XmEdge> cellEdges;
    for (int i(0); i < ugrid2d->GetCellCount(); i++)
    {
      for (int j(0); j < ugrid2d->GetCellEdgeCount(i); j++)
      {
        cellEdges.push_back(ugrid2d->GetCellEdge(i, j));
      }
      TS_ASSERT_EQUALS(expectedCellsCellEdges[i], cellEdges);
      cellEdges.clear();
    }
  }

  {
    // 3D Shapes
    std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
    TS_REQUIRE_NOT_NULL(ugrid3d);

    // clang-format off
    std::vector<std::vector<XmEdge>> expectedCellsCellEdges =
    {
      {{0, 1}, {1, 5}, {5, 0}, {0, 15}, {1, 15}, {5, 15}}, // XMU_TETRA
      {{1, 2}, {2, 7}, {6, 7}, {1, 6}, {16, 17}, {17, 22}, {21, 22}, {16, 21}, {1, 16}, {2, 17}, {6, 21}, {7, 22}}, // XMU_VOXEL
      {{2, 3}, {3, 8}, {7, 8}, {2, 7}, {17, 18}, {18, 23}, {22, 23}, {17, 22}, {2, 17}, {3, 18}, {7, 22}, {8, 23}}, // XMU_HEXAHEDRON
      {{8, 9}, {8, 13}, {8, 23}, {9, 14}, {9, 24}, {13, 14}, {13, 28}, {14, 29}, {23, 24}, {23, 28}, {24, 29}, {28, 29}}, // XMU_POLYHEDRON
      {{3, 4}, {4, 18}, {18, 3}, {8, 9}, {9, 23}, {23, 8}, {3, 8}, {4, 9}, {18, 23}}, // XMU_WEDGE
      {{5, 6}, {6, 11}, {11, 10}, {10, 5}, {5, 20}, {6, 20}, {11, 20}, {10, 20}} // XMU_PYRAMID
    };
    // clang-format on
    std::vector<XmEdge> cellEdges;
    for (int i(0); i < ugrid3d->GetCellCount(); i++)
    {
      for (int j(0); j < ugrid3d->GetCellEdgeCount(i); j++)
      {
        cellEdges.push_back(ugrid3d->GetCellEdge(i, j));
      }
      TS_ASSERT_EQUALS(expectedCellsCellEdges[i], cellEdges);
      cellEdges.clear();
    }
  }
} // XmUGridUnitTests::testGetCellEdge

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetPointsAdjacentCells()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  VecInt expectedCells;
  VecInt points = {-1, 0};
  VecInt retrievedCells;
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);

  points = {0, -1};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);

  points = {0, 8};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  points.clear();
  retrievedCells.clear();

  expectedCells = {0, 2};
  points = {3, 4};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  points.clear();
  retrievedCells.clear();

  expectedCells = {0};
  points = {0, 3};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  points.clear();
  retrievedCells.clear();

  points = {3, 4, 5};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  points.clear();
  retrievedCells.clear();

  points = {ugrid->GetPointCount(), 0};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);

  points = {0, ugrid->GetPointCount()};
  retrievedCells = ugrid->GetPointsAdjacentCells(points);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);

} // XmUGridUnitTests::testGetPointsAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellAdjacentCells()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  VecInt expectedCells;

  VecInt retrievedCells;
  retrievedCells = ugrid->GetCellAdjacentCells(-1);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  retrievedCells = ugrid->GetCellAdjacentCells(0);
  expectedCells = {2, 1, 3};
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  retrievedCells = ugrid->GetCellAdjacentCells(4);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  retrievedCells = ugrid->GetCellAdjacentCells(5);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells.clear();
  retrievedCells.clear();

  // 2D Shapes
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  expectedCells = {1, 5};
  retrievedCells = ugrid2d->GetCellAdjacentCells(0);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {0, 2, 3, 4};
  retrievedCells = ugrid2d->GetCellAdjacentCells(1);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {1, 3, 4};
  retrievedCells = ugrid2d->GetCellAdjacentCells(2);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {2, 1, 4};
  retrievedCells = ugrid2d->GetCellAdjacentCells(3);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {1, 2, 3};
  retrievedCells = ugrid2d->GetCellAdjacentCells(4);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {0};
  retrievedCells = ugrid2d->GetCellAdjacentCells(5);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);

  // 3D Shapes
  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  expectedCells = {1, 5};
  retrievedCells = ugrid3d->GetCellAdjacentCells(0);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {0, 2, 5};
  retrievedCells = ugrid3d->GetCellAdjacentCells(1);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {1, 4, 3};
  retrievedCells = ugrid3d->GetCellAdjacentCells(2);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {4, 2};
  retrievedCells = ugrid3d->GetCellAdjacentCells(3);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {2, 3};
  retrievedCells = ugrid3d->GetCellAdjacentCells(4);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);
  expectedCells = {0, 1};
  retrievedCells = ugrid3d->GetCellAdjacentCells(5);
  TS_ASSERT_EQUALS(expectedCells, retrievedCells);

} // XmUGridUnitTests::testGetCellAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellEdgeAdjacentCells()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  VecInt expectedFail;
  VecInt2d expectedCells = {{}, {2}, {1}, {}};
  VecInt expected2dCells = {-1, 2, 1, -1};

  int adjacentCell;
  VecInt adjacentCells = ugrid->GetCellEdgeAdjacentCells(0, -1);
  TS_ASSERT_EQUALS(expectedFail, adjacentCells);
  adjacentCell = ugrid->GetCell2dEdgeAdjacentCell(0, -1);
  TS_ASSERT_EQUALS(-1, adjacentCell);

  for (int i(0); i < ugrid->GetCellEdgeCount(0); i++)
  {
    adjacentCells = ugrid->GetCellEdgeAdjacentCells(0, i);
    TS_ASSERT_EQUALS(expectedCells[i], adjacentCells);

    adjacentCell = ugrid->GetCell2dEdgeAdjacentCell(0, i);
    TS_ASSERT_EQUALS(expected2dCells[i], adjacentCell);
  }
  adjacentCells = ugrid->GetCellEdgeAdjacentCells(0, ugrid->GetCellEdgeCount(0));
  TS_ASSERT_EQUALS(expectedFail, adjacentCells);
  adjacentCell = ugrid->GetCell2dEdgeAdjacentCell(0, ugrid->GetCellEdgeCount(0));
  TS_ASSERT_EQUALS(-1, adjacentCell);

  // Test a cell in the middle of an hexahedron grid
  std::shared_ptr<XmUGrid> hexUgrid = TEST_XmUBuildHexahedronUgrid(5, 5, 5);
  TS_REQUIRE_NOT_NULL(hexUgrid);

  expectedCells = {{2, 6, 18},   {6, 7, 23},   {6, 10, 26},  {5, 6, 21},
                   {18, 34, 38}, {23, 38, 39}, {26, 38, 42}, {21, 37, 38},
                   {17, 18, 21}, {18, 19, 23}, {21, 26, 25}, {23, 26, 27}};

  for (int i(0); i < hexUgrid->GetCellEdgeCount(22); i++)
  {
    adjacentCells = hexUgrid->GetCellEdgeAdjacentCells(22, i);
    TS_ASSERT_EQUALS(expectedCells[i], adjacentCells);
  }
} // XmUGridUnitTests::testGetCellEdgeAdjacentCells
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testEdgeAdjacentCells()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  VecInt expectedFail;
  VecInt2d expectedCells = {{0, 1}, {0, 2}, {0}};
  std::vector<XmEdge> edges = {{1, 4}, {3, 4}, {0, 3}};

  VecInt adjacentCells = ugrid->GetEdgeAdjacentCells(XmEdge(-1, -1));
  TS_ASSERT_EQUALS(expectedFail, adjacentCells);
  for (int i(0); i < edges.size(); i++)
  {
    adjacentCells = ugrid->GetEdgeAdjacentCells(edges[i]);
    TS_ASSERT_EQUALS(expectedCells[i], adjacentCells);
  }

  XmEdge badEdge(0, ugrid->GetPointCount());
  adjacentCells = ugrid->GetEdgeAdjacentCells(badEdge);
  TS_ASSERT_EQUALS(expectedFail, adjacentCells);

} // XmUGridUnitTests::testEdgeAdjacentCells

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCell3dFacePoints()
{
  // 2D Tests, include the bounds check
  VecInt emptyCellFaces = {};
  VecInt2d expectedCellFaces = {// XMU_QUAD
                                {},
                                // XMU_PIXEL
                                {},
                                // XMU_TRIANGLE
                                {},
                                // XMU_POLYGON
                                {},
                                // XMU_POLY_LINE
                                {},
                                // XMU_LINE
                                {}};
  VecInt cellFaces;

  // 2D Shapes
  std::shared_ptr<XmUGrid> ugrid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(ugrid2d);

  cellFaces = ugrid2d->GetCell3dFacePoints(-1, 0);
  TS_ASSERT_EQUALS(emptyCellFaces, cellFaces);
  cellFaces = ugrid2d->GetCell3dFacePoints(0, -1);
  TS_ASSERT_EQUALS(emptyCellFaces, cellFaces);

  for (int i(0); i < ugrid2d->GetCellCount(); i++)
  {
    cellFaces = ugrid2d->GetCell3dFacePoints(i, 0);
    TS_ASSERT_EQUALS(expectedCellFaces[i], cellFaces);
  }

  cellFaces = ugrid2d->GetCell3dFacePoints(ugrid2d->GetCellCount(), 0);
  TS_ASSERT_EQUALS(emptyCellFaces, cellFaces);
  cellFaces = ugrid2d->GetCell3dFacePoints(0, 1);
  TS_ASSERT_EQUALS(emptyCellFaces, cellFaces);

  // 3D Shapes
  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  expectedCellFaces = {// Tetra
                       {0, 1, 15},
                       {1, 5, 15},
                       {5, 0, 15},
                       {0, 5, 1},
                       // Voxel
                       {1, 16, 21, 6},
                       {2, 7, 22, 17},
                       {1, 2, 17, 16},
                       {6, 21, 22, 7},
                       {1, 6, 7, 2},
                       {16, 17, 22, 21},
                       // Hexahedron
                       {2, 17, 22, 7},
                       {3, 8, 23, 18},
                       {2, 3, 18, 17},
                       {7, 22, 23, 8},
                       {2, 7, 8, 3},
                       {17, 18, 23, 22},
                       // Polyhedron
                       {9, 8, 13, 14},
                       {8, 9, 24, 23},
                       {9, 14, 29, 24},
                       {14, 13, 28, 29},
                       {8, 13, 28, 23},
                       {23, 24, 29, 28},
                       // Wedge
                       {3, 4, 18},
                       {8, 23, 9},
                       {3, 8, 9, 4},
                       {4, 9, 23, 18},
                       {18, 23, 8, 3},
                       // Pyramid
                       {5, 10, 11, 6},
                       {5, 6, 20},
                       {6, 11, 20},
                       {11, 10, 20},
                       {10, 5, 20}};
  int currId(0);
  for (int i(0); i < ugrid3d->GetCellCount(); i++)
  {
    for (int j(0); j < ugrid3d->GetCell3dFaceCount(i); j++, currId++)
    {
      cellFaces = ugrid3d->GetCell3dFacePoints(i, j);
      TS_ASSERT_EQUALS(expectedCellFaces[currId], cellFaces);
    }
  }
} // XmUGridUnitTests::testGetCell3dFacePoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCell3dFaceAdjacentCell()
{
  // 2 hexahedrons
  int rows = 3;
  int cols = 2;
  int lays = 2;

  std::shared_ptr<xms::XmUGrid> grid = TEST_XmUBuildHexahedronUgrid(rows, cols, lays);
  TS_REQUIRE_NOT_NULL(grid);

  // clang-format off
  VecInt expectedNeighbor = {-1, -1, -1, 1, -1, -1, -1, -1, 0, -1, -1, -1};
  VecInt expectedFace =     {-1, -1, -1, 2, -1, -1, -1, -1, 3, -1, -1, -1};
  // clang-format on

  int expectedIdx = 0;
  int neighborCellIdx;

  neighborCellIdx = grid->GetCell3dFaceAdjacentCell(-1, 0);
  TS_ASSERT_EQUALS(-1, neighborCellIdx);
  neighborCellIdx = grid->GetCell3dFaceAdjacentCell(0, -1);
  TS_ASSERT_EQUALS(-1, neighborCellIdx);

  for (int cellIdx = 0; cellIdx < grid->GetCellCount(); cellIdx++)
  {
    for (int faceIdx = 0; faceIdx < grid->GetCell3dFaceCount(cellIdx); faceIdx++, expectedIdx++)
    {
      // Check first function
      neighborCellIdx = grid->GetCell3dFaceAdjacentCell(cellIdx, faceIdx);
      TS_ASSERT_EQUALS(expectedNeighbor[expectedIdx], neighborCellIdx);

      // Check overload function
      int neighborFaceIdx = -1;
      grid->GetCell3dFaceAdjacentCell(cellIdx, faceIdx, neighborCellIdx, neighborFaceIdx);
      TS_ASSERT_EQUALS(expectedNeighbor[expectedIdx], neighborCellIdx);
      TS_ASSERT_EQUALS(expectedFace[expectedIdx], neighborFaceIdx);

      // Check that faces are the same
      VecInt expectedFacePt = grid->GetCell3dFacePoints(cellIdx, faceIdx);
      std::sort(expectedFacePt.begin(), expectedFacePt.end());
      VecInt neighborFacePt = grid->GetCell3dFacePoints(neighborCellIdx, neighborFaceIdx);
      std::sort(neighborFacePt.begin(), neighborFacePt.end());
      if (!expectedFacePt.empty() && !neighborFacePt.empty())
      {
        TS_ASSERT_EQUALS(expectedFacePt, neighborFacePt);
      }
    }
  }
  neighborCellIdx = grid->GetCell3dFaceAdjacentCell(grid->GetCellCount(), 0);
  TS_ASSERT_EQUALS(-1, neighborCellIdx);
  neighborCellIdx = grid->GetCell3dFaceAdjacentCell(0, grid->GetCellCount());
  TS_ASSERT_EQUALS(-1, neighborCellIdx);

  rows = 3;
  cols = 3;
  lays = 3;
  grid = TEST_XmUBuildHexahedronUgrid(rows, cols, lays);
  TS_REQUIRE_NOT_NULL(grid);

  expectedNeighbor.clear();
  // clang-format off
  expectedNeighbor = {-1, 1, -1, 2, -1, 4, 0, -1, -1, 3, -1, 5, -1, 3, 0, -1, -1, 6, 2, -1, 1, -1, -1, 7, -1, 5, -1, 6, 0, -1, 4, -1, -1, 7, 1, -1, -1, 7, 4, -1, 2, -1, 6, -1, 5, -1, 3, -1};
  expectedFace =     {-1, 0, -1, 2, -1, 4, 1, -1, -1, 2, -1, 4, -1, 0, 3, -1, -1, 4, 1, -1, 3, -1, -1, 4, -1, 0, -1, 2, 5, -1, 1, -1, -1, 2, 5, -1, -1, 0, 3, -1, 5, -1, 1, -1, 3, -1, 5, -1};
  // clang-format on

  expectedIdx = 0;
  for (int cellIdx = 0; cellIdx < grid->GetCellCount(); ++cellIdx)
  {
    for (int faceIdx = 0; faceIdx < grid->GetCell3dFaceCount(cellIdx); faceIdx++, expectedIdx++)
    {
      // Check first function
      neighborCellIdx = grid->GetCell3dFaceAdjacentCell(cellIdx, faceIdx);
      TS_ASSERT_EQUALS(expectedNeighbor[expectedIdx], neighborCellIdx);

      // Check overload function
      int neighborFaceIdx = -1;
      grid->GetCell3dFaceAdjacentCell(cellIdx, faceIdx, neighborCellIdx, neighborFaceIdx);
      TS_ASSERT_EQUALS(expectedNeighbor[expectedIdx], neighborCellIdx);
      TS_ASSERT_EQUALS(expectedFace[expectedIdx], neighborFaceIdx);

      // Check that faces are the same
      VecInt expectedFacePt = grid->GetCell3dFacePoints(cellIdx, faceIdx);
      std::sort(expectedFacePt.begin(), expectedFacePt.end());
      VecInt neighborFacePt = grid->GetCell3dFacePoints(neighborCellIdx, neighborFaceIdx);
      std::sort(neighborFacePt.begin(), neighborFacePt.end());
      if (!expectedFacePt.empty() && !neighborFacePt.empty())
      {
        TS_ASSERT_EQUALS(expectedFacePt, neighborFacePt);
      }
    }
  }
} // XmUGridUnitTests::testGetCell3dFaceAdjacentCell

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetPointAdjacentPoints()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  std::shared_ptr<xms::XmUGrid> grid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(grid);

  VecInt attachedIdxs;
  grid->GetPointAdjacentPoints(0, attachedIdxs);
  VecInt expectedIdxs = {1, 3};
  TS_ASSERT_EQUALS(expectedIdxs, attachedIdxs);
  grid->GetPointAdjacentPoints(3, attachedIdxs);
  expectedIdxs = {0, 4, 6};
  TS_ASSERT_EQUALS(expectedIdxs, attachedIdxs);
  grid->GetPointAdjacentPoints(4, attachedIdxs);
  expectedIdxs = {1, 3, 5, 7};
  TS_ASSERT_EQUALS(expectedIdxs, attachedIdxs);
  VecPt3d attachedPts;
  grid->GetPointAdjacentLocations(0, attachedPts);
  VecPt3d expectedPts = {{10, 10, 0}, {0, 0, 0}};
  TS_ASSERT_EQUALS(expectedPts, attachedPts);
  grid->GetPointAdjacentLocations(3, attachedPts);
  expectedPts = {{0, 10, 0}, {10, 0, 0}, {0, -10, 0}};
  TS_ASSERT_EQUALS(expectedPts, attachedPts);
  grid->GetPointAdjacentLocations(4, attachedPts);
  expectedPts = {{10, 10, 0}, {0, 0, 0}, {20, 0, 0}, {10, -10, 0}};
  TS_ASSERT_EQUALS(expectedPts, attachedPts);
} // XmUGridUnitTests::testGetPointAdjacentPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellPlanViewPolygon()
{
  // 2d
  std::shared_ptr<xms::XmUGrid> grid2d = TEST_XmUGrid2dLinear();
  TS_REQUIRE_NOT_NULL(grid2d);

  VecPt3d2d expectedPolygons{
    {{0, 0, 0}, {10, 0, 0}, {10, 10, 0}, {0, 10, 0}},
    {{10, 0, 0}, {20, 0, 0}, {20, 10, 0}, {10, 10, 0}},
    {{20, 0, 0}, {30, 0, 0}, {20, 10, 0}},
    {{30, 0, 0}, {40, 0, 0}, {40, 10, 0}, {40, 20, 0}, {30, 20, 0}, {20, 10, 0}},
    {},
    {}};
  VecPt3d viewPolygon;
  VecPt3d empty;
  // Bounds testing
  TS_ASSERT(!grid2d->GetCellPlanViewPolygon(-1, viewPolygon));
  TS_ASSERT_EQUALS(empty, viewPolygon);
  TS_ASSERT(!grid2d->GetCellPlanViewPolygon(grid2d->GetCellCount(), viewPolygon));
  TS_ASSERT_EQUALS(empty, viewPolygon);
  // Regular testing
  for (int i = 0; i < grid2d->GetCellCount(); ++i)
  {
    // Should return true for the first 4 cells, and false for the last 2
    TS_ASSERT(grid2d->GetCellPlanViewPolygon(i, viewPolygon) != i >= 4);
    TS_ASSERT_EQUALS(expectedPolygons[i], viewPolygon);
  }

  // 3d
  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  expectedPolygons = {{{0, 0, 0}, {10, 0, 0}, {0, 10, 0}},
                      {{20, 0, 0}, {20, 10, 0}, {10, 10, 0}, {10, 0, 0}},
                      {{30, 0, 0}, {30, 10, 0}, {20, 10, 0}, {20, 0, 0}},
                      {{40, 10, 0}, {40, 20, 0}, {30, 20, 0}, {30, 10, 0}},
                      {{30, 0, 0}, {40, 0, 0}, {40, 10, 0}, {30, 10, 0}},
                      {{0, 10, 0}, {10, 10, 0}, {10, 20, 0}, {0, 20, 0}}};

  for (int i = 0; i < ugrid3d->GetCellCount(); ++i)
  {
    TS_ASSERT(ugrid3d->GetCellPlanViewPolygon(i, viewPolygon));
    TS_ASSERT_EQUALS(expectedPolygons[i], viewPolygon);
  }

  std::shared_ptr<xms::XmUGrid> chevronUgrid = TEST_XmUBuild3DChevronUgrid();
  TS_REQUIRE_NOT_NULL(chevronUgrid);

  expectedPolygons = {{{0.0, 0.0, 0.0}, {20.0, 10.0, 0.0}, {40.0, 0.0, 0.0}, {20.0, 50.0, 0.0}}};
  TS_ASSERT(chevronUgrid->GetCellPlanViewPolygon(0, viewPolygon));
  TS_ASSERT_EQUALS(expectedPolygons[0], viewPolygon);

} // XmUGridUnitTests::testGetCellPlanViewPolygon

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCellPlanViewPolygonMultiSideFace()
{
  // 3d - For vertically prismatic cells returns the unique plan view points for side faces
  //      in sorted order. For cells that aren't vertically prismatic returns an empty vector.
  VecPt3d points = {
    {0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {0, 0, 1}, {1, 0, 1},
    {2, 0, 1}, {0, 0, 2}, {2, 0, 2},                       // front face
    {0, 1, 2}, {1, 1, 2}, {2, 1, 2},                       // points midway back
    {0, 2, 0}, {2, 2, 0}, {0, 2, 2}, {1, 2, 2}, {2, 2, 2}, // back points
    {4, 2, 0}, {3, 1, 2}                                   // tetrahedron
  };
  // clang-format off
  VecInt cells = {
    XMU_POLYHEDRON, 10,     // hexahedron split into smaller faces
    4, 0, 1, 4, 3,          // hex front - lower left face
    4, 1, 2, 5, 4,          // hex front - lower right face
    5, 3, 4, 5, 7, 6,       // hex front - top face
    5, 6, 7, 10, 9, 8,      // hex top - front face
    4, 8, 9, 14, 13,        // hex top - left back face
    4, 9, 10, 15, 14,       // hex top - right back face
    6, 0, 3, 6, 8, 13, 11,  // hex left
    6, 2, 12, 15, 10, 7, 5, // hex right
    5, 0, 11, 12, 2, 1,     // hex bottom
    5, 11, 13, 14, 15, 12,  // hex back
    XMU_TETRA, 4, 2, 16, 12, 17 // tetrahedron with 
  };
  // clang-format on
  std::shared_ptr<XmUGrid> ugrid3d = XmUGrid::New(points, cells);
  TS_REQUIRE_NOT_NULL(ugrid3d);

  // vertically prismatic
  VecPt3d cellPoints;
  ugrid3d->GetCellPlanViewPolygon(0, cellPoints);
  VecPt3d expectedGetCellPoints = {{1, 0, 0}, {2, 0, 0}, {2, 1, 0}, {2, 2, 0},
                                   {1, 2, 0}, {0, 2, 0}, {0, 1, 0}, {0, 0, 0}};
  TS_ASSERT_EQUALS(expectedGetCellPoints, cellPoints);

  // not vertically prismatic
  ugrid3d->GetCellPlanViewPolygon(1, cellPoints);
  expectedGetCellPoints = {{2, 0, 0}, {4, 2, 0}, {2, 2, 0}};
  TS_ASSERT_EQUALS(expectedGetCellPoints, cellPoints);
} // XmUGridUnitTests::testGetCellPlanViewPolygonMultiSideFace
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testIsCellValidWithPointChange()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8

  // VecPt3d points = { { 0, 10, 0 }, { 10, 10, 0 }, { 20, 10, 0 }, { 0, 0, 0 }, { 10, 0, 0 },
  //{ 20, 0, 0 }, { 0, -10, 0 }, { 10, -10, 0 }, { 20, -10, 0 } };

  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  Pt3d validPt = {5.0, 5.0, 0.0};
  Pt3d invalid = {500.0, 500.0, 0.0};

  TS_ASSERT(ugrid->IsValidPointChange(4, validPt));
  TS_ASSERT(!ugrid->IsValidPointChange(4, invalid));

  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUBuildHexahedronUgrid(4, 4, 4);

  validPt = ugrid3d->GetPointLocation(21);
  validPt.x += 0.5;
  validPt.y += 0.5;
  validPt.z += 0.5;
  TS_ASSERT(ugrid3d->IsValidPointChange(21, validPt));
  TS_ASSERT(!ugrid3d->IsValidPointChange(21, invalid));

  TS_ASSERT(ugrid->SetPointLocation(4, validPt));
  TS_ASSERT(ugrid->SetPointLocation(4, invalid));
  TS_ASSERT(ugrid3d->SetPointLocation(21, validPt));
  TS_ASSERT(ugrid3d->SetPointLocation(21, invalid));
} // XmUGridUnitTests::testIsCellValidWithPointChange

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testPointFunctions()
{
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  // Test GetPointCount
  TS_ASSERT_EQUALS(ugrid->GetPointCount(), 9);
  VecPt3d points = {{0, 10, 0}, {10, 10, 0}, {20, 10, 0},  {0, 0, 0},   {10, 0, 0},
                    {20, 0, 0}, {0, -10, 0}, {10, -10, 0}, {20, -10, 0}};

  // Test Locations
  TS_ASSERT_EQUALS(ugrid->GetLocations(), points);
  for (int i = 0; i < points.size(); ++i)
  {
    TS_ASSERT_EQUALS(points[i], ugrid->GetPointLocation(i));
  }

  // Test SetPointLocation & GetPointLocation
  Pt3d invalid = {100, 100, 100};
  Pt3d valid = {5, 5, 5};
  TS_ASSERT(!ugrid->IsValidPointChange(4, invalid));
  TS_ASSERT(ugrid->SetPointLocation(4, valid));
  TS_ASSERT_EQUALS(valid, ugrid->GetPointLocation(4));

  VecInt pointIndices = {0, 3, 6};
  VecPt3d expectedPoints = {{0, 10, 0}, {0, 0, 0}, {0, -10, 0}};
  VecPt3d pointsPt3d;

  // Test GetPointsLocations
  pointsPt3d = ugrid->GetPointsLocations(pointIndices);
  TS_ASSERT_EQUALS(expectedPoints, pointsPt3d);

  // Test GetExtents
  Pt3d min, max, expectedMin = {0, -10, 0}, expectedMax = {20, 10, 5};
  ugrid->GetExtents(min, max);
  TS_ASSERT_EQUALS(expectedMin, min);
  TS_ASSERT_EQUALS(expectedMax, max);

  // Test GetPointAdjacentCells
  VecInt expectedCells = {0, 2};
  VecInt cellsAssociated = ugrid->GetPointAdjacentCells(3);
  TS_ASSERT_EQUALS(expectedCells, cellsAssociated);

  // Test GetPointsAdjacentCells
  VecInt expectedCommonCells = {0, 1};
  VecInt pointIndex = {1};
  VecInt commonCells = ugrid->GetPointsAdjacentCells(pointIndex);
  TS_ASSERT_EQUALS(expectedCells, cellsAssociated);

} // XmUGridUnitTests::testPointFunctions

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testCellFunctions()
{
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  // Test GetCellPoints
  VecInt pointsOfCell = ugrid->GetCellPoints(0);
  VecInt expectedPoints = {0, 3, 4, 1};
  TS_ASSERT_EQUALS(expectedPoints, pointsOfCell);
  // Test overload
  ugrid->GetCellPoints(0, pointsOfCell);
  TS_ASSERT_EQUALS(expectedPoints, pointsOfCell);
  VecPt3d locations;
  ugrid->GetCellLocations(0, locations);
  VecPt3d expectedLocations = {{0, 10, 0}, {0, 0, 0}, {10, 0, 0}, {10, 10, 0}};
  TS_ASSERT_EQUALS(expectedLocations, locations);

  // Test GetCellExtents
  Pt3d min, max;
  ugrid->GetCellExtents(0, min, max);
  TS_ASSERT_EQUALS(Pt3d(0, 0, 0), min);
  TS_ASSERT_EQUALS(Pt3d(10, 10, 0), max);

  // Test GetCellType
  TS_ASSERT_EQUALS(XMU_QUAD, ugrid->GetCellType(0));

  // Test GetDimensionCounts
  VecInt expectedDimensions = {0, 0, 4, 0};
  TS_ASSERT_EQUALS(expectedDimensions, ugrid->GetDimensionCounts());

  // Test GetCellDimension
  TS_ASSERT_EQUALS(2, ugrid->GetCellDimension(0));

  // Test GetCellCentroid
  Pt3d centroid;
  TS_ASSERT(ugrid->GetCellCentroid(0, centroid));
  TS_ASSERT_EQUALS(Pt3d(5, 5, 0), centroid);
} // XmUGridUnitTests::testCellFunctions
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testCellstreamFunctions()
{
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  // Test GetCellstream
  VecInt cellstream = ugrid->GetCellstream();
  VecInt expectedCellstream = {XMU_QUAD, 4, 0, 3, 4, 1, XMU_QUAD, 4, 1, 4, 5, 2,
                               XMU_QUAD, 4, 3, 6, 7, 4, XMU_QUAD, 4, 4, 7, 8, 5};
  TS_ASSERT_EQUALS(expectedCellstream, cellstream);

  // Test GetCellCellstream
  ugrid->GetCellCellstream(0, cellstream);
  expectedCellstream = {XMU_QUAD, 4, 0, 3, 4, 1};
  TS_ASSERT_EQUALS(expectedCellstream, cellstream);

} // XmUGridUnitTests::testCellstreamFunctions
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testCellEdgeAdjacentCellFunctions()
{
  //     0-----1-----2
  //     |  0  |  1  |
  //     3-----4-----5
  //     |  2  |  3  |
  //     6-----7-----8
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  // Get Adjacent Cells from cell and edge
  VecInt2d expectedCells = {{}, {2}, {1}, {}};
  VecInt expected2dCells = {-1, 2, 1, -1};
  int adjacentCell;
  VecInt adjacentCells;
  for (int i(0); i < ugrid->GetCellEdgeCount(0); i++)
  {
    adjacentCells = ugrid->GetCellEdgeAdjacentCells(0, i);
    TS_ASSERT_EQUALS(expectedCells[i], adjacentCells);
    // For 2D cells only
    adjacentCell = ugrid->GetCell2dEdgeAdjacentCell(0, i);
    TS_ASSERT_EQUALS(expected2dCells[i], adjacentCell);
  }

  // Get Adjacent cells from given Edge
  VecInt2d expectedCellsFromEdge = {{0, 1}, {0, 2}, {0}};
  std::vector<XmEdge> edges = {{1, 4}, {3, 4}, {0, 3}};
  adjacentCells = ugrid->GetEdgeAdjacentCells(XmEdge(-1, -1));
  for (int i(0); i < edges.size(); i++)
  {
    adjacentCells = ugrid->GetEdgeAdjacentCells(edges[i]);
    TS_ASSERT_EQUALS(expectedCellsFromEdge[i], adjacentCells);
  }
} // XmUGridUnitTests::testCellAdjacentCellFunctions
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testCellEdges()
{
  std::shared_ptr<XmUGrid> ugrid = TEST_XmUGridSimpleQuad();
  TS_REQUIRE_NOT_NULL(ugrid);

  // Test GetCellEdges
  std::vector<XmEdge> edges = ugrid->GetCellEdges(0);
  std::vector<XmEdge> expectededges = {{0, 3}, {3, 4}, {4, 1}, {1, 0}};
  TS_ASSERT_EQUALS(expectededges, edges);

} // XmUGridUnitTests::testCellEdges
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testCell3dFaceFunctions()
{
  // 3D Shapes
  std::shared_ptr<XmUGrid> ugrid3d = TEST_XmUGrid3dLinear();
  TS_REQUIRE_NOT_NULL(ugrid3d);

  VecInt cellFace;
  VecInt2d expectedCellFaces = {// Tetra
                                {0, 1, 15},
                                {1, 5, 15},
                                {5, 0, 15},
                                {0, 5, 1},
                                // Voxel
                                {1, 16, 21, 6},
                                {2, 7, 22, 17},
                                {1, 2, 17, 16},
                                {6, 21, 22, 7},
                                {1, 6, 7, 2},
                                {16, 17, 22, 21},
                                // Hexahedron
                                {2, 17, 22, 7},
                                {3, 8, 23, 18},
                                {2, 3, 18, 17},
                                {7, 22, 23, 8},
                                {2, 7, 8, 3},
                                {17, 18, 23, 22},
                                // Polyhedron
                                {9, 8, 13, 14},
                                {8, 9, 24, 23},
                                {9, 14, 29, 24},
                                {14, 13, 28, 29},
                                {8, 13, 28, 23},
                                {23, 24, 29, 28},
                                // Wedge
                                {3, 4, 18},
                                {8, 23, 9},
                                {3, 8, 9, 4},
                                {4, 9, 23, 18},
                                {18, 23, 8, 3},
                                // Pyramid
                                {5, 10, 11, 6},
                                {5, 6, 20},
                                {6, 11, 20},
                                {11, 10, 20},
                                {10, 5, 20}};
  int expectedIdx = 0;
  for (int i(0); i < ugrid3d->GetCellCount(); i++)
  {
    for (int j(0); j < ugrid3d->GetCell3dFaceCount(i); j++, expectedIdx++)
    {
      cellFace = ugrid3d->GetCell3dFacePoints(i, j);
      TS_ASSERT_EQUALS(expectedCellFaces[expectedIdx], cellFace);
    }
  }

  VecInt2d cellFaces;
  expectedIdx = 0;
  for (int i(0); i < ugrid3d->GetCellCount(); i++)
  {
    cellFaces = ugrid3d->GetCell3dFacesPoints(i);
    TS_ASSERT_EQUALS(cellFaces.size(), ugrid3d->GetCell3dFaceCount(i));
    for (int j(0); j < ugrid3d->GetCell3dFaceCount(i); j++, expectedIdx++)
    {
      TS_ASSERT_EQUALS(expectedCellFaces[expectedIdx], cellFaces[j]);
    }
  }

  VecInt faceOrientations;
  for (int cellIdx = 0; cellIdx < ugrid3d->GetCellCount(); ++cellIdx)
  {
    for (int faceIdx = 0; faceIdx < ugrid3d->GetCell3dFaceCount(cellIdx); ++faceIdx)
    {
      faceOrientations.push_back(ugrid3d->GetCell3dFaceOrientation(cellIdx, faceIdx));
    }
  }

  VecInt expectedFaceOrientations = {
    // Tetra
    XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_TOP, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_BOTTOM,
    // Voxel
    XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
    XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_TOP,
    // Hexahedron
    XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
    XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_TOP,
    // Polyhedron
    XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
    XMU_ORIENTATION_SIDE, XMU_ORIENTATION_TOP,
    // Wedge
    XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_TOP,
    XMU_ORIENTATION_SIDE,
    // Pyramid
    XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_TOP, XMU_ORIENTATION_TOP,
    XMU_ORIENTATION_SIDE};

  TS_ASSERT_EQUALS(expectedFaceOrientations, faceOrientations);
} // XmUGridUnitTests::testCell3dFaceFunctions
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCell3dFaceOrientationHexahedrons()
{
  std::shared_ptr<XmUGrid> xmUGrid = TEST_XmUBuildHexahedronUgrid(2, 2, 2, Pt3d());
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 0));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 1));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 2));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 3));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_TOP, xmUGrid->GetCell3dFaceOrientation(0, 4));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_BOTTOM, xmUGrid->GetCell3dFaceOrientation(0, 5));
  xmUGrid = TEST_XmUBuildPolyhedronUgrid(2, 2, 2, Pt3d());
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 0));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 1));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 2));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_SIDE, xmUGrid->GetCell3dFaceOrientation(0, 3));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_TOP, xmUGrid->GetCell3dFaceOrientation(0, 4));
  TS_ASSERT_EQUALS(XMU_ORIENTATION_BOTTOM, xmUGrid->GetCell3dFaceOrientation(0, 5));
} // XmUGridUnitTests::testGetCell3dFaceOrientationHexahedrons
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testGetCell3dFaceOrientationConcaveCell()
{
  // clang-format off
  // vertically prismatic concave cell
  VecPt3d nodes = {
    Pt3d(1620022.8468, 6134363.759, 0),    Pt3d(1620009.9411, 6134414.9476, 0),
    Pt3d(1619994.9996, 6134289.4991, 0),   Pt3d(1619866.1047, 6134542.8755, 0),
    Pt3d(1619745.5374, 6134167.0467, 0),   Pt3d(1619829.9996, 6134192.9991, 0),
    Pt3d(1619773.3077, 6134545.2322, 0),   Pt3d(1619710.0815, 6134182.8542, 0),
    Pt3d(1619693.1618, 6134208.2547, 0),   Pt3d(1619645.5529, 6134438.0278, 0),
    Pt3d(1619774.9993, 6134371.9982, 0),   Pt3d(1620022.8468, 6134363.759, -10),
    Pt3d(1620009.9411, 6134414.9476, -10), Pt3d(1619994.9996, 6134289.4991, -10),
    Pt3d(1619866.1047, 6134542.8755, -10), Pt3d(1619745.5374, 6134167.0467, -10),
    Pt3d(1619829.9996, 6134192.9991, -10), Pt3d(1619773.3077, 6134545.2322, -10),
    Pt3d(1619710.0815, 6134182.8542, -10), Pt3d(1619693.1618, 6134208.2547, -10),
    Pt3d(1619645.5529, 6134438.0278, -10), Pt3d(1619774.9993, 6134371.9982, -10) };
  VecInt elements = {
    XMU_POLYHEDRON,
    13, // number of faces
    11, 5, 10, 2, 0, 1, 3, 6, 9, 8, 7, 4, // top
    11, 16, 15, 18, 19, 20, 17, 14, 12, 11, 13, 21, // bottom
    4, 5, 16, 21, 10, // sides
    4, 10, 21, 13, 2,
    4, 2, 13, 11, 0,
    4, 0, 11, 12, 1,
    4, 1, 12, 14, 3,
    4, 3, 14, 17, 6,
    4, 6, 17, 20, 9,
    4, 9, 20, 19, 8,
    4, 8, 19, 18, 7,
    4, 7, 18, 15, 4,
    4, 4, 1, 16, 5
  };
  // clang-format on

  std::shared_ptr<XmUGrid> xmUGrid = XmUGrid::New(nodes, elements);
  VecInt actual;
  for (int faceIdx = 0; faceIdx < xmUGrid->GetCell3dFaceCount(0); ++faceIdx)
  {
    actual.push_back(xmUGrid->GetCell3dFaceOrientation(0, faceIdx));
  }
  VecInt expected = {XMU_ORIENTATION_TOP,  XMU_ORIENTATION_BOTTOM, XMU_ORIENTATION_SIDE,
                     XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,   XMU_ORIENTATION_SIDE,
                     XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,   XMU_ORIENTATION_SIDE,
                     XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,   XMU_ORIENTATION_SIDE,
                     XMU_ORIENTATION_SIDE};
  TS_ASSERT_EQUALS(expected, actual);
} // XmUGridUnitTests::testFaceOrientationConcaveCell
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridUnitTests::testCell3dFunctionCaching()
{
  std::shared_ptr<XmUGrid> xmUGrid = TEST_XmUBuildHexahedronUgrid(2, 3, 2);

  VecInt2d expectedNeighbors = {{-1, 1, -1, -1, -1, -1}, {0, -1, -1, -1, -1, -1}};
  VecInt2d expectedOrientations = {
    {XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
     XMU_ORIENTATION_TOP, XMU_ORIENTATION_BOTTOM},
    {XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
     XMU_ORIENTATION_TOP, XMU_ORIENTATION_BOTTOM}};

  // test four times, to fill the cache, read from the cache, turn off the cache,
  // and turn on the cache
  bool useCacheValues[] = {true, true, false, true};
  for (bool useCacheValue : useCacheValues)
  {
    xmUGrid->SetUseCache(useCacheValue);
    for (int cellIdx = 0; cellIdx < xmUGrid->GetCellCount(); ++cellIdx)
    {
      TS_ASSERT_EQUALS(6, xmUGrid->GetCell3dFaceCount(cellIdx));
      for (int faceIdx = 0; faceIdx < xmUGrid->GetCell3dFaceCount(cellIdx); ++faceIdx)
      {
        int neighbor = xmUGrid->GetCell3dFaceAdjacentCell(cellIdx, faceIdx);
        TS_ASSERT_EQUALS(expectedNeighbors[cellIdx][faceIdx], neighbor);
        int orientation = (int)xmUGrid->GetCell3dFaceOrientation(cellIdx, faceIdx);
        TS_ASSERT_EQUALS(expectedOrientations[cellIdx][faceIdx], orientation);
      }
    }
  }

  // change to different points and cells and test again
  expectedNeighbors = {{-1, 1, -1, -1, -1, -1}, {0, 2, -1, -1, -1, -1}, {1, -1, -1, -1, -1, -1}};
  expectedOrientations = {{XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
                           XMU_ORIENTATION_SIDE, XMU_ORIENTATION_TOP, XMU_ORIENTATION_BOTTOM},
                          {XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
                           XMU_ORIENTATION_SIDE, XMU_ORIENTATION_TOP, XMU_ORIENTATION_BOTTOM},
                          {XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE, XMU_ORIENTATION_SIDE,
                           XMU_ORIENTATION_SIDE, XMU_ORIENTATION_TOP, XMU_ORIENTATION_BOTTOM}};
  std::shared_ptr<XmUGrid> newXmUGrid = TEST_XmUBuildHexahedronUgrid(2, 4, 2);
  xmUGrid->SetLocations(newXmUGrid->GetLocations());
  xmUGrid->SetCellstream(newXmUGrid->GetCellstream());
  for (bool useCacheValue : useCacheValues)
  {
    xmUGrid->SetUseCache(useCacheValue);
    for (int cellIdx = 0; cellIdx < xmUGrid->GetCellCount(); ++cellIdx)
    {
      TS_ASSERT_EQUALS(6, xmUGrid->GetCell3dFaceCount(cellIdx));
      for (int faceIdx = 0; faceIdx < xmUGrid->GetCell3dFaceCount(cellIdx); ++faceIdx)
      {
        int neighbor = xmUGrid->GetCell3dFaceAdjacentCell(cellIdx, faceIdx);
        TS_ASSERT_EQUALS(expectedNeighbors[cellIdx][faceIdx], neighbor);
        int orientation = (int)xmUGrid->GetCell3dFaceOrientation(cellIdx, faceIdx);
        TS_ASSERT_EQUALS(expectedOrientations[cellIdx][faceIdx], orientation);
      }
    }
  }
} // XmUGridUnitTests::testCell3dFunctionCaching
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//#define SPEEDTEST
#ifdef SPEEDTEST
#include <boost/timer/timer.hpp>
#endif
void XmUGridUnitTests::testLargeUGridLinkSpeed()
{
#ifdef SPEEDTEST
  int rows = 1000;
  int cols = 500;
  int lays = 4;

  // std::shared_ptr<xms::XmUGrid> grid = TEST_XmUBuildQuadUGrid(rows, cols);
  // std::shared_ptr<xms::XmUGrid> grid = TEST_XmUBuildHexahedronUgrid(rows, cols, lays);
  std::shared_ptr<xms::XmUGrid> grid = TEST_XmUBuildPolyhedronUgrid(rows, cols, lays);

  {
    boost::timer::cpu_timer timer;
    std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(grid->GetLocations(), grid->GetCellstream());
    std::cerr << timer.format() << '\n';
  }

#endif
} // XmUGridUnitTests::testLargeUGridLinkSpeed

#endif