/*=========================================================================

  Program:   Visualization Toolkit
  Module:    UnitTestCells.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#include "vtkCellType.h"
#include "vtkSmartPointer.h"

#include "vtkEmptyCell.h"
#include "vtkHexagonalPrism.h"
#include "vtkHexahedron.h"
#include "vtkLine.h"
#include "vtkPentagonalPrism.h"
#include "vtkPixel.h"
#include "vtkPolyLine.h"
#include "vtkPolyVertex.h"
#include "vtkPolygon.h"
#include "vtkPolyhedron.h"
#include "vtkPyramid.h"
#include "vtkQuad.h"
#include "vtkTetra.h"
#include "vtkTriangle.h"
#include "vtkTriangleStrip.h"
#include "vtkVertex.h"
#include "vtkVoxel.h"
#include "vtkWedge.h"

#include "vtkQuadraticEdge.h"
#include "vtkQuadraticHexahedron.h"
#include "vtkQuadraticLinearQuad.h"
#include "vtkQuadraticLinearWedge.h"
#include "vtkQuadraticPolygon.h"
#include "vtkQuadraticPyramid.h"
#include "vtkQuadraticQuad.h"
#include "vtkQuadraticTetra.h"
#include "vtkQuadraticTriangle.h"
#include "vtkQuadraticWedge.h"

#include "vtkBiQuadraticQuad.h"
#include "vtkBiQuadraticQuadraticHexahedron.h"
#include "vtkBiQuadraticQuadraticWedge.h"
#include "vtkBiQuadraticTriangle.h"
#include "vtkTriQuadraticHexahedron.h"
#include "vtkTriQuadraticPyramid.h"

#include "vtkCubicLine.h"

#include "vtkCellArray.h"
#include "vtkMath.h"
#include "vtkMathUtilities.h"
#include "vtkPoints.h"
#include <map>
#include <sstream>
#include <string>
#include <vector>

static vtkSmartPointer<vtkEmptyCell> MakeEmptyCell();
static vtkSmartPointer<vtkVertex> MakeVertex();
static vtkSmartPointer<vtkPolyVertex> MakePolyVertex();
static vtkSmartPointer<vtkLine> MakeLine();
static vtkSmartPointer<vtkPolyLine> MakePolyLine();
static vtkSmartPointer<vtkTriangle> MakeTriangle();
static vtkSmartPointer<vtkTriangleStrip> MakeTriangleStrip();
static vtkSmartPointer<vtkPolygon> MakePolygon();
static vtkSmartPointer<vtkQuad> MakeQuad();
static vtkSmartPointer<vtkPixel> MakePixel();
static vtkSmartPointer<vtkVoxel> MakeVoxel();
static vtkSmartPointer<vtkHexahedron> MakeHexahedron();
static vtkSmartPointer<vtkPyramid> MakePyramid();
static vtkSmartPointer<vtkTetra> MakeTetra();
static vtkSmartPointer<vtkWedge> MakeWedge();
static vtkSmartPointer<vtkPentagonalPrism> MakePentagonalPrism();
static vtkSmartPointer<vtkHexagonalPrism> MakeHexagonalPrism();
static vtkSmartPointer<vtkPolyhedron> MakeCube();
static vtkSmartPointer<vtkPolyhedron> MakeDodecahedron();

static vtkSmartPointer<vtkQuadraticEdge> MakeQuadraticEdge();
static vtkSmartPointer<vtkQuadraticHexahedron> MakeQuadraticHexahedron();
static vtkSmartPointer<vtkQuadraticPolygon> MakeQuadraticPolygon();
static vtkSmartPointer<vtkQuadraticLinearQuad> MakeQuadraticLinearQuad();
static vtkSmartPointer<vtkQuadraticLinearWedge> MakeQuadraticLinearWedge();
static vtkSmartPointer<vtkQuadraticPyramid> MakeQuadraticPyramid();
static vtkSmartPointer<vtkQuadraticQuad> MakeQuadraticQuad();
static vtkSmartPointer<vtkQuadraticTetra> MakeQuadraticTetra();
static vtkSmartPointer<vtkQuadraticTriangle> MakeQuadraticTriangle();
static vtkSmartPointer<vtkQuadraticWedge> MakeQuadraticWedge();

static vtkSmartPointer<vtkBiQuadraticQuad> MakeBiQuadraticQuad();
static vtkSmartPointer<vtkBiQuadraticQuadraticHexahedron> MakeBiQuadraticQuadraticHexahedron();
static vtkSmartPointer<vtkBiQuadraticQuadraticWedge> MakeBiQuadraticQuadraticWedge();
static vtkSmartPointer<vtkBiQuadraticTriangle> MakeBiQuadraticTriangle();
static vtkSmartPointer<vtkTriQuadraticHexahedron> MakeTriQuadraticHexahedron();
static vtkSmartPointer<vtkTriQuadraticPyramid> MakeTriQuadraticPyramid();
static vtkSmartPointer<vtkCubicLine> MakeCubicLine();

template <typename T>
int TestOneCell(VTKCellType cellType, vtkSmartPointer<T> cell, int linear = 1);
//------------------------------------------------------------------------------
int UnitTestCells(int, char*[])
{
  std::map<std::string, int> results;

  results["EmptyCell"] = TestOneCell<vtkEmptyCell>(VTK_EMPTY_CELL, MakeEmptyCell());
  results["Vertex"] = TestOneCell<vtkVertex>(VTK_VERTEX, MakeVertex());
  results["PolyVertex"] = TestOneCell<vtkPolyVertex>(VTK_POLY_VERTEX, MakePolyVertex());
  results["Line"] = TestOneCell<vtkLine>(VTK_LINE, MakeLine());
  results["PolyLine"] = TestOneCell<vtkPolyLine>(VTK_POLY_LINE, MakePolyLine());
  results["Triangle"] = TestOneCell<vtkTriangle>(VTK_TRIANGLE, MakeTriangle());
  results["TriangleStrip"] = TestOneCell<vtkTriangleStrip>(VTK_TRIANGLE_STRIP, MakeTriangleStrip());
  results["Polygon"] = TestOneCell<vtkPolygon>(VTK_POLYGON, MakePolygon());
  results["Pixel"] = TestOneCell<vtkPixel>(VTK_PIXEL, MakePixel());
  results["Quad"] = TestOneCell<vtkQuad>(VTK_QUAD, MakeQuad());
  results["Tetra"] = TestOneCell<vtkTetra>(VTK_TETRA, MakeTetra());
  results["Voxel"] = TestOneCell<vtkVoxel>(VTK_VOXEL, MakeVoxel());
  results["Hexahedron"] = TestOneCell<vtkHexahedron>(VTK_HEXAHEDRON, MakeHexahedron());
  results["Wedge"] = TestOneCell<vtkWedge>(VTK_WEDGE, MakeWedge());
  results["Pyramid"] = TestOneCell<vtkPyramid>(VTK_PYRAMID, MakePyramid());
  results["PentagonalPrism"] =
    TestOneCell<vtkPentagonalPrism>(VTK_PENTAGONAL_PRISM, MakePentagonalPrism());
  results["HexagonalPrism"] =
    TestOneCell<vtkHexagonalPrism>(VTK_HEXAGONAL_PRISM, MakeHexagonalPrism());
  results["Polyhedron(Cube)"] = TestOneCell<vtkPolyhedron>(VTK_POLYHEDRON, MakeCube());
  results["Polyhedron(Dodecahedron)"] =
    TestOneCell<vtkPolyhedron>(VTK_POLYHEDRON, MakeDodecahedron());

  results["QuadraticEdge"] =
    TestOneCell<vtkQuadraticEdge>(VTK_QUADRATIC_EDGE, MakeQuadraticEdge(), 0);
  results["QuadraticHexahedron"] =
    TestOneCell<vtkQuadraticHexahedron>(VTK_QUADRATIC_HEXAHEDRON, MakeQuadraticHexahedron(), 0);
  results["QuadraticPolygon"] =
    TestOneCell<vtkQuadraticPolygon>(VTK_QUADRATIC_POLYGON, MakeQuadraticPolygon(), 0);
  results["QuadraticLinearQuad"] =
    TestOneCell<vtkQuadraticLinearQuad>(VTK_QUADRATIC_LINEAR_QUAD, MakeQuadraticLinearQuad(), 0);
  results["QuadraticLinearWedge"] =
    TestOneCell<vtkQuadraticLinearWedge>(VTK_QUADRATIC_LINEAR_WEDGE, MakeQuadraticLinearWedge(), 0);
  results["QuadraticPyramid"] =
    TestOneCell<vtkQuadraticPyramid>(VTK_QUADRATIC_PYRAMID, MakeQuadraticPyramid(), 0);
  results["QuadraticQuad"] =
    TestOneCell<vtkQuadraticQuad>(VTK_QUADRATIC_QUAD, MakeQuadraticQuad(), 0);
  results["QuadraticTetra"] =
    TestOneCell<vtkQuadraticTetra>(VTK_QUADRATIC_TETRA, MakeQuadraticTetra(), 0);
  results["QuadraticTrangle"] =
    TestOneCell<vtkQuadraticTriangle>(VTK_QUADRATIC_TRIANGLE, MakeQuadraticTriangle(), 0);
  results["QuadraticWedge"] =
    TestOneCell<vtkQuadraticWedge>(VTK_QUADRATIC_WEDGE, MakeQuadraticWedge(), 0);

  results["BiQuadraticQuad"] =
    TestOneCell<vtkBiQuadraticQuad>(VTK_BIQUADRATIC_QUAD, MakeBiQuadraticQuad(), 0);
  results["BiQuadraticQuadraticHexahedron"] = TestOneCell<vtkBiQuadraticQuadraticHexahedron>(
    VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON, MakeBiQuadraticQuadraticHexahedron(), 0);
  results["BiQuadraticQuadraticWedge"] = TestOneCell<vtkBiQuadraticQuadraticWedge>(
    VTK_BIQUADRATIC_QUADRATIC_WEDGE, MakeBiQuadraticQuadraticWedge(), 0);
  results["BiQuadraticTrangle"] =
    TestOneCell<vtkBiQuadraticTriangle>(VTK_BIQUADRATIC_TRIANGLE, MakeBiQuadraticTriangle(), 0);
  results["CubicLine"] = TestOneCell<vtkCubicLine>(VTK_CUBIC_LINE, MakeCubicLine(), 0);

  results["TriQuadraticHexahedron"] = TestOneCell<vtkTriQuadraticHexahedron>(
    VTK_TRIQUADRATIC_HEXAHEDRON, MakeTriQuadraticHexahedron(), 0);
  results["TriQuadraticPyramid"] =
    TestOneCell<vtkTriQuadraticPyramid>(VTK_TRIQUADRATIC_PYRAMID, MakeTriQuadraticPyramid(), 0);

  int status = 0;
  std::cout << "----- Unit Test Summary -----" << std::endl;
  std::map<std::string, int>::iterator it;
  for (it = results.begin(); it != results.end(); ++it)
  {
    std::cout << std::setw(25) << it->first << " " << (it->second ? " FAILED" : " OK") << std::endl;
    if (it->second != 0)
    {
      ++status;
    }
  }
  if (status)
  {
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

vtkSmartPointer<vtkEmptyCell> MakeEmptyCell()
{
  auto anEmptyCell = vtkSmartPointer<vtkEmptyCell>::New();

  return anEmptyCell;
}

vtkSmartPointer<vtkVertex> MakeVertex()
{
  auto aVertex = vtkSmartPointer<vtkVertex>::New();

  aVertex->GetPointIds()->SetId(0, 0);
  aVertex->GetPoints()->SetPoint(0, 10.0, 20.0, 30.0);

  return aVertex;
}

vtkSmartPointer<vtkPolyVertex> MakePolyVertex()
{
  auto aPolyVertex = vtkSmartPointer<vtkPolyVertex>::New();

  aPolyVertex->GetPointIds()->SetNumberOfIds(2);
  aPolyVertex->GetPointIds()->SetId(0, 0);
  aPolyVertex->GetPointIds()->SetId(1, 1);

  aPolyVertex->GetPoints()->SetNumberOfPoints(2);
  aPolyVertex->GetPoints()->SetPoint(0, 10.0, 20.0, 30.0);
  aPolyVertex->GetPoints()->SetPoint(1, 30.0, 20.0, 10.0);

  return aPolyVertex;
}

vtkSmartPointer<vtkLine> MakeLine()
{
  auto aLine = vtkSmartPointer<vtkLine>::New();

  aLine->GetPointIds()->SetId(0, 0);
  aLine->GetPointIds()->SetId(1, 1);
  aLine->GetPoints()->SetPoint(0, 10.0, 20.0, 30.0);
  aLine->GetPoints()->SetPoint(1, 30.0, 20.0, 10.0);

  return aLine;
}

vtkSmartPointer<vtkPolyLine> MakePolyLine()
{
  auto aPolyLine = vtkSmartPointer<vtkPolyLine>::New();

  aPolyLine->GetPointIds()->SetNumberOfIds(3);
  aPolyLine->GetPointIds()->SetId(0, 0);
  aPolyLine->GetPointIds()->SetId(1, 1);
  aPolyLine->GetPointIds()->SetId(2, 2);

  aPolyLine->GetPoints()->SetNumberOfPoints(3);
  aPolyLine->GetPoints()->SetPoint(0, 10.0, 20.0, 30.0);
  aPolyLine->GetPoints()->SetPoint(1, 10.0, 30.0, 30.0);
  aPolyLine->GetPoints()->SetPoint(2, 10.0, 30.0, 40.0);

  return aPolyLine;
}

vtkSmartPointer<vtkTriangle> MakeTriangle()
{
  auto aTriangle = vtkSmartPointer<vtkTriangle>::New();

  for (int i = 0; i < aTriangle->GetNumberOfPoints(); ++i)
  {
    aTriangle->GetPointIds()->SetId(i, i);
  }

  aTriangle->GetPoints()->SetPoint(0, -10.0, -10.0, 0.0);
  aTriangle->GetPoints()->SetPoint(1, 10.0, -10.0, 0.0);
  aTriangle->GetPoints()->SetPoint(2, 10.0, 10.0, 0.0);

  return aTriangle;
}

vtkSmartPointer<vtkTriangleStrip> MakeTriangleStrip()
{
  auto aTriangleStrip = vtkSmartPointer<vtkTriangleStrip>::New();

  aTriangleStrip->GetPointIds()->SetNumberOfIds(4);
  aTriangleStrip->GetPointIds()->SetId(0, 0);
  aTriangleStrip->GetPointIds()->SetId(1, 1);
  aTriangleStrip->GetPointIds()->SetId(2, 2);
  aTriangleStrip->GetPointIds()->SetId(3, 3);

  aTriangleStrip->GetPoints()->SetNumberOfPoints(4);
  aTriangleStrip->GetPoints()->SetPoint(0, 10.0, 10.0, 10.0);
  aTriangleStrip->GetPoints()->SetPoint(1, 12.0, 10.0, 10.0);
  aTriangleStrip->GetPoints()->SetPoint(2, 11.0, 12.0, 10.0);
  aTriangleStrip->GetPoints()->SetPoint(3, 13.0, 10.0, 10.0);

  return aTriangleStrip;
}

vtkSmartPointer<vtkPolygon> MakePolygon()
{
  auto aPolygon = vtkSmartPointer<vtkPolygon>::New();

  aPolygon->GetPointIds()->SetNumberOfIds(4);
  aPolygon->GetPointIds()->SetId(0, 0);
  aPolygon->GetPointIds()->SetId(1, 1);
  aPolygon->GetPointIds()->SetId(2, 2);
  aPolygon->GetPointIds()->SetId(3, 3);

  aPolygon->GetPoints()->SetNumberOfPoints(4);
  aPolygon->GetPoints()->SetPoint(0, 0.0, 0.0, 0.0);
  aPolygon->GetPoints()->SetPoint(1, 10.0, 0.0, 0.0);
  aPolygon->GetPoints()->SetPoint(2, 10.0, 10.0, 0.0);
  aPolygon->GetPoints()->SetPoint(3, 0.0, 10.0, 0.0);

  return aPolygon;
}

vtkSmartPointer<vtkQuad> MakeQuad()
{
  auto aQuad = vtkSmartPointer<vtkQuad>::New();

  for (int i = 0; i < aQuad->GetNumberOfPoints(); ++i)
  {
    aQuad->GetPointIds()->SetId(i, i);
  }

  aQuad->GetPoints()->SetPoint(0, -10.0, -10.0, 0.0);
  aQuad->GetPoints()->SetPoint(1, 10.0, -10.0, 0.0);
  aQuad->GetPoints()->SetPoint(2, 10.0, 10.0, 0.0);
  aQuad->GetPoints()->SetPoint(3, -10.0, 10.0, 0.0);

  return aQuad;
}

vtkSmartPointer<vtkPixel> MakePixel()
{
  auto aPixel = vtkSmartPointer<vtkPixel>::New();

  for (int i = 0; i < aPixel->GetNumberOfPoints(); ++i)
  {
    aPixel->GetPointIds()->SetId(i, i);
  }

  aPixel->GetPoints()->SetPoint(0, 10.0, 10.0, 10.0);
  aPixel->GetPoints()->SetPoint(1, 12.0, 10.0, 10.0);
  aPixel->GetPoints()->SetPoint(2, 10.0, 12.0, 10.0);
  aPixel->GetPoints()->SetPoint(3, 12.0, 12.0, 10.0);

  return aPixel;
}

vtkSmartPointer<vtkVoxel> MakeVoxel()
{
  auto aVoxel = vtkSmartPointer<vtkVoxel>::New();

  for (int i = 0; i < aVoxel->GetNumberOfPoints(); ++i)
  {
    aVoxel->GetPointIds()->SetId(i, i);
  }

  aVoxel->GetPoints()->SetPoint(0, 10, 10, 10);
  aVoxel->GetPoints()->SetPoint(1, 12, 10, 10);
  aVoxel->GetPoints()->SetPoint(2, 10, 12, 10);
  aVoxel->GetPoints()->SetPoint(3, 12, 12, 10);
  aVoxel->GetPoints()->SetPoint(4, 10, 10, 12);
  aVoxel->GetPoints()->SetPoint(5, 12, 10, 12);
  aVoxel->GetPoints()->SetPoint(6, 10, 12, 12);
  aVoxel->GetPoints()->SetPoint(7, 12, 12, 12);

  return aVoxel;
}

vtkSmartPointer<vtkHexahedron> MakeHexahedron()
{
  auto aHexahedron = vtkSmartPointer<vtkHexahedron>::New();

  for (int i = 0; i < aHexahedron->GetNumberOfPoints(); ++i)
  {
    aHexahedron->GetPointIds()->SetId(i, i);
  }

  aHexahedron->GetPoints()->SetPoint(0, 10, 10, 10);
  aHexahedron->GetPoints()->SetPoint(1, 12, 10, 10);
  aHexahedron->GetPoints()->SetPoint(2, 12, 12, 10);
  aHexahedron->GetPoints()->SetPoint(3, 10, 12, 10);
  aHexahedron->GetPoints()->SetPoint(4, 10, 10, 12);
  aHexahedron->GetPoints()->SetPoint(5, 12, 10, 12);
  aHexahedron->GetPoints()->SetPoint(6, 12, 12, 12);
  aHexahedron->GetPoints()->SetPoint(7, 10, 12, 12);

  return aHexahedron;
}

vtkSmartPointer<vtkPyramid> MakePyramid()
{
  auto aPyramid = vtkSmartPointer<vtkPyramid>::New();

  for (int i = 0; i < aPyramid->GetNumberOfPoints(); ++i)
  {
    aPyramid->GetPointIds()->SetId(i, i);
  }

  aPyramid->GetPoints()->SetPoint(0, 0.0, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(1, 1.0, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(2, 1.0, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(3, 0.0, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(4, 0.5, 0.5, 1.0);

  return aPyramid;
}

vtkSmartPointer<vtkQuadraticPyramid> MakeQuadraticPyramid()
{
  auto aPyramid = vtkSmartPointer<vtkQuadraticPyramid>::New();

  for (int i = 0; i < aPyramid->GetNumberOfPoints(); ++i)
  {
    aPyramid->GetPointIds()->SetId(i, i);
  }

  aPyramid->GetPoints()->SetPoint(0, 0.0, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(1, 1.0, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(2, 1.0, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(3, 0.0, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(4, 0.5, 0.5, 1.0);

  aPyramid->GetPoints()->SetPoint(5, 0.5, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(6, 1.0, 0.5, 0.0);
  aPyramid->GetPoints()->SetPoint(7, 0.5, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(8, 0.0, 0.5, 0.0);
  aPyramid->GetPoints()->SetPoint(9, 0.25, 0.25, 0.5);
  aPyramid->GetPoints()->SetPoint(10, 0.75, 0.25, 0.5);
  aPyramid->GetPoints()->SetPoint(11, 0.75, 0.75, 0.5);
  aPyramid->GetPoints()->SetPoint(12, 0.25, 0.75, 0.5);

  return aPyramid;
}

vtkSmartPointer<vtkTriQuadraticPyramid> MakeTriQuadraticPyramid()
{
  auto aPyramid = vtkSmartPointer<vtkTriQuadraticPyramid>::New();

  for (int i = 0; i < aPyramid->GetNumberOfPoints(); ++i)
  {
    aPyramid->GetPointIds()->SetId(i, i);
  }

  aPyramid->GetPoints()->SetPoint(0, 0.0, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(1, 1.0, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(2, 1.0, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(3, 0.0, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(4, 0.5, 0.5, 1.0);

  aPyramid->GetPoints()->SetPoint(5, 0.5, 0.0, 0.0);
  aPyramid->GetPoints()->SetPoint(6, 1.0, 0.5, 0.0);
  aPyramid->GetPoints()->SetPoint(7, 0.5, 1.0, 0.0);
  aPyramid->GetPoints()->SetPoint(8, 0.0, 0.5, 0.0);
  aPyramid->GetPoints()->SetPoint(9, 0.25, 0.25, 0.5);
  aPyramid->GetPoints()->SetPoint(10, 0.75, 0.25, 0.5);
  aPyramid->GetPoints()->SetPoint(11, 0.75, 0.75, 0.5);
  aPyramid->GetPoints()->SetPoint(12, 0.25, 0.75, 0.5);

  aPyramid->GetPoints()->SetPoint(13, 0.5, 0.5, 0);
  aPyramid->GetPoints()->SetPoint(14, 0.5, 1.0 / 6.0, 1.0 / 3.0);
  aPyramid->GetPoints()->SetPoint(15, 5.0 / 6.0, 0.5, 1.0 / 3.0);
  aPyramid->GetPoints()->SetPoint(16, 0.5, 5.0 / 6.0, 1.0 / 3.0);
  aPyramid->GetPoints()->SetPoint(17, 1.0 / 6.0, 0.5, 1.0 / 3.0);

  aPyramid->GetPoints()->SetPoint(18, 0.5, 0.5, 0.2);

  return aPyramid;
}

vtkSmartPointer<vtkQuadraticEdge> MakeQuadraticEdge()
{
  auto anEdge = vtkSmartPointer<vtkQuadraticEdge>::New();

  for (int i = 0; i < anEdge->GetNumberOfPoints(); ++i)
  {
    anEdge->GetPointIds()->SetId(i, i);
  }

  anEdge->GetPoints()->SetPoint(0, 0, 0, 0);
  anEdge->GetPoints()->SetPoint(1, 1, 0, 0);
  anEdge->GetPoints()->SetPoint(2, .5, 0, 0);

  return anEdge;
}

vtkSmartPointer<vtkQuadraticHexahedron> MakeQuadraticHexahedron()
{
  auto aHexahedron = vtkSmartPointer<vtkQuadraticHexahedron>::New();

  double* pcoords = aHexahedron->GetParametricCoords();
  for (int i = 0; i < aHexahedron->GetNumberOfPoints(); ++i)
  {
    aHexahedron->GetPointIds()->SetId(i, i);
    aHexahedron->GetPoints()->SetPoint(i, *(pcoords + 3 * i) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 1) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 2) + vtkMath::Random(-.1, .1));
  }

  return aHexahedron;
}

vtkSmartPointer<vtkBiQuadraticQuadraticHexahedron> MakeBiQuadraticQuadraticHexahedron()
{
  auto aHexahedron = vtkSmartPointer<vtkBiQuadraticQuadraticHexahedron>::New();

  double* pcoords = aHexahedron->GetParametricCoords();
  for (int i = 0; i < aHexahedron->GetNumberOfPoints(); ++i)
  {
    aHexahedron->GetPointIds()->SetId(i, i);
    aHexahedron->GetPoints()->SetPoint(i, *(pcoords + 3 * i) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 1) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 2) + vtkMath::Random(-.1, .1));
  }

  return aHexahedron;
}

vtkSmartPointer<vtkTriQuadraticHexahedron> MakeTriQuadraticHexahedron()
{
  auto aHexahedron = vtkSmartPointer<vtkTriQuadraticHexahedron>::New();

  double* pcoords = aHexahedron->GetParametricCoords();
  for (int i = 0; i < aHexahedron->GetNumberOfPoints(); ++i)
  {
    aHexahedron->GetPointIds()->SetId(i, i);
    aHexahedron->GetPoints()->SetPoint(i, *(pcoords + 3 * i) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 1) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 2) + vtkMath::Random(-.1, .1));
  }

  return aHexahedron;
}

vtkSmartPointer<vtkQuadraticPolygon> MakeQuadraticPolygon()
{
  auto aPolygon = vtkSmartPointer<vtkQuadraticPolygon>::New();

  aPolygon->GetPointIds()->SetNumberOfIds(8);
  aPolygon->GetPointIds()->SetId(0, 0);
  aPolygon->GetPointIds()->SetId(1, 1);
  aPolygon->GetPointIds()->SetId(2, 2);
  aPolygon->GetPointIds()->SetId(3, 3);
  aPolygon->GetPointIds()->SetId(4, 4);
  aPolygon->GetPointIds()->SetId(5, 5);
  aPolygon->GetPointIds()->SetId(6, 6);
  aPolygon->GetPointIds()->SetId(7, 7);

  aPolygon->GetPoints()->SetNumberOfPoints(8);
  aPolygon->GetPoints()->SetPoint(0, 0.0, 0.0, 0.0);
  aPolygon->GetPoints()->SetPoint(1, 2.0, 0.0, 0.0);
  aPolygon->GetPoints()->SetPoint(2, 2.0, 2.0, 0.0);
  aPolygon->GetPoints()->SetPoint(3, 0.0, 2.0, 0.0);
  aPolygon->GetPoints()->SetPoint(4, 1.0, 0.0, 0.0);
  aPolygon->GetPoints()->SetPoint(5, 2.0, 1.0, 0.0);
  aPolygon->GetPoints()->SetPoint(6, 1.0, 2.0, 0.0);
  aPolygon->GetPoints()->SetPoint(7, 0.0, 1.0, 0.0);
  aPolygon->GetPoints()->SetPoint(5, 3.0, 1.0, 0.0);

  return aPolygon;
}

vtkSmartPointer<vtkQuadraticLinearQuad> MakeQuadraticLinearQuad()
{
  auto aLinearQuad = vtkSmartPointer<vtkQuadraticLinearQuad>::New();

  double* pcoords = aLinearQuad->GetParametricCoords();
  for (int i = 0; i < aLinearQuad->GetNumberOfPoints(); ++i)
  {
    aLinearQuad->GetPointIds()->SetId(i, i);
    aLinearQuad->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aLinearQuad;
}

vtkSmartPointer<vtkQuadraticLinearWedge> MakeQuadraticLinearWedge()
{
  auto aLinearWedge = vtkSmartPointer<vtkQuadraticLinearWedge>::New();

  double* pcoords = aLinearWedge->GetParametricCoords();
  for (int i = 0; i < 12; ++i)
  {
    aLinearWedge->GetPointIds()->SetId(i, i);
    aLinearWedge->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aLinearWedge;
}

vtkSmartPointer<vtkQuadraticQuad> MakeQuadraticQuad()
{
  auto aQuad = vtkSmartPointer<vtkQuadraticQuad>::New();

  double* pcoords = aQuad->GetParametricCoords();
  for (int i = 0; i < 8; ++i)
  {
    aQuad->GetPointIds()->SetId(i, i);
    aQuad->GetPoints()->SetPoint(i, *(pcoords + 3 * i) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 1) + vtkMath::Random(-.1, .1), *(pcoords + 3 * i + 2));
  }

  return aQuad;
}

vtkSmartPointer<vtkQuadraticTetra> MakeQuadraticTetra()
{
  auto aTetra = vtkSmartPointer<vtkQuadraticTetra>::New();

  double* pcoords = aTetra->GetParametricCoords();
  for (int i = 0; i < 10; ++i)
  {
    aTetra->GetPointIds()->SetId(i, i);
    aTetra->GetPoints()->SetPoint(i, *(pcoords + 3 * i) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 1) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 2) + vtkMath::Random(-.1, .1));
  }

  return aTetra;
}

vtkSmartPointer<vtkQuadraticTriangle> MakeQuadraticTriangle()
{
  auto aTriangle = vtkSmartPointer<vtkQuadraticTriangle>::New();

  double* pcoords = aTriangle->GetParametricCoords();
  for (int i = 0; i < aTriangle->GetNumberOfPoints(); ++i)
  {
    aTriangle->GetPointIds()->SetId(i, i);
    aTriangle->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aTriangle;
}

vtkSmartPointer<vtkBiQuadraticTriangle> MakeBiQuadraticTriangle()
{
  auto aTriangle = vtkSmartPointer<vtkBiQuadraticTriangle>::New();

  double* pcoords = aTriangle->GetParametricCoords();
  for (int i = 0; i < aTriangle->GetNumberOfPoints(); ++i)
  {
    aTriangle->GetPointIds()->SetId(i, i);
    aTriangle->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aTriangle;
}

vtkSmartPointer<vtkBiQuadraticQuad> MakeBiQuadraticQuad()
{
  auto aQuad = vtkSmartPointer<vtkBiQuadraticQuad>::New();

  double* pcoords = aQuad->GetParametricCoords();
  for (int i = 0; i < aQuad->GetNumberOfPoints(); ++i)
  {
    aQuad->GetPointIds()->SetId(i, i);
    aQuad->GetPoints()->SetPoint(i, *(pcoords + 3 * i) + vtkMath::Random(-.1, .1),
      *(pcoords + 3 * i + 1) + vtkMath::Random(-.1, .1), *(pcoords + 3 * i + 2));
  }

  return aQuad;
}

vtkSmartPointer<vtkCubicLine> MakeCubicLine()
{
  auto aLine = vtkSmartPointer<vtkCubicLine>::New();

  double* pcoords = aLine->GetParametricCoords();
  for (int i = 0; i < aLine->GetNumberOfPoints(); ++i)
  {
    aLine->GetPointIds()->SetId(i, i);
    aLine->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aLine;
}

vtkSmartPointer<vtkQuadraticWedge> MakeQuadraticWedge()
{
  auto aWedge = vtkSmartPointer<vtkQuadraticWedge>::New();

  double* pcoords = aWedge->GetParametricCoords();
  for (int i = 0; i < aWedge->GetNumberOfPoints(); ++i)
  {
    aWedge->GetPointIds()->SetId(i, i);
    aWedge->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aWedge;
}

vtkSmartPointer<vtkBiQuadraticQuadraticWedge> MakeBiQuadraticQuadraticWedge()
{
  auto aWedge = vtkSmartPointer<vtkBiQuadraticQuadraticWedge>::New();

  double* pcoords = aWedge->GetParametricCoords();
  for (int i = 0; i < aWedge->GetNumberOfPoints(); ++i)
  {
    aWedge->GetPointIds()->SetId(i, i);
    aWedge->GetPoints()->SetPoint(
      i, *(pcoords + 3 * i), *(pcoords + 3 * i + 1), *(pcoords + 3 * i + 2));
  }

  return aWedge;
}

vtkSmartPointer<vtkTetra> MakeTetra()
{
  auto aTetra = vtkSmartPointer<vtkTetra>::New();

  for (int i = 0; i < aTetra->GetNumberOfPoints(); ++i)
  {
    aTetra->GetPointIds()->SetId(i, i);
  }

  aTetra->GetPoints()->SetPoint(0, 10.0, 10.0, 10.0);
  aTetra->GetPoints()->SetPoint(1, 12.0, 10.0, 10.0);
  aTetra->GetPoints()->SetPoint(2, 11.0, 12.0, 10.0);
  aTetra->GetPoints()->SetPoint(3, 11.0, 11.0, 12.0);

  return aTetra;
}

vtkSmartPointer<vtkWedge> MakeWedge()
{
  auto aWedge = vtkSmartPointer<vtkWedge>::New();

  for (int i = 0; i < aWedge->GetNumberOfPoints(); ++i)
  {
    aWedge->GetPointIds()->SetId(i, i);
  }

  aWedge->GetPoints()->SetPoint(0, 0.0, 1.0, 0.0);
  aWedge->GetPoints()->SetPoint(1, 0.0, 0.0, 0.0);
  aWedge->GetPoints()->SetPoint(2, 0.0, 0.5, 0.5);
  aWedge->GetPoints()->SetPoint(3, 1.0, 1.0, 0.0);
  aWedge->GetPoints()->SetPoint(4, 1.0, 0.0, 0.0);
  aWedge->GetPoints()->SetPoint(5, 1.0, 0.5, 0.5);

  return aWedge;
}

vtkSmartPointer<vtkPolyhedron> MakeCube()
{
  auto aCube = vtkSmartPointer<vtkPolyhedron>::New();

  // create polyhedron (cube)
  auto points = vtkSmartPointer<vtkPoints>::New();

  aCube->GetPointIds()->SetNumberOfIds(8);
  aCube->GetPointIds()->SetId(0, 0);
  aCube->GetPointIds()->SetId(1, 1);
  aCube->GetPointIds()->SetId(2, 2);
  aCube->GetPointIds()->SetId(3, 3);
  aCube->GetPointIds()->SetId(4, 4);
  aCube->GetPointIds()->SetId(5, 5);
  aCube->GetPointIds()->SetId(6, 6);
  aCube->GetPointIds()->SetId(7, 7);

  aCube->GetPoints()->SetNumberOfPoints(8);
  aCube->GetPoints()->SetPoint(0, -1.0, -1.0, -1.0);
  aCube->GetPoints()->SetPoint(1, 1.0, -1.0, -1.0);
  aCube->GetPoints()->SetPoint(2, 1.0, 1.0, -1.0);
  aCube->GetPoints()->SetPoint(3, -1.0, 1.0, -1.0);
  aCube->GetPoints()->SetPoint(4, -1.0, -1.0, 1.0);
  aCube->GetPoints()->SetPoint(5, 1.0, -1.0, 1.0);
  aCube->GetPoints()->SetPoint(6, 1.0, 1.0, 1.0);
  aCube->GetPoints()->SetPoint(7, -1.0, 1.0, 1.0);

  vtkIdType faces[31] = {
    6,             // number of faces
    4, 0, 3, 2, 1, //
    4, 0, 4, 7, 3, //
    4, 4, 5, 6, 7, //
    4, 5, 1, 2, 6, //
    4, 0, 1, 5, 4, //
    4, 2, 3, 7, 6  //
  };

  aCube->SetFaces(faces);
  aCube->Initialize();

  return aCube;
}

vtkSmartPointer<vtkPolyhedron> MakeDodecahedron()
{
  auto aDodecahedron = vtkSmartPointer<vtkPolyhedron>::New();

  // create polyhedron (dodecahedron)
  auto points = vtkSmartPointer<vtkPoints>::New();

  for (int i = 0; i < 20; ++i)
  {
    aDodecahedron->GetPointIds()->InsertNextId(i);
  }

  aDodecahedron->GetPoints()->InsertNextPoint(1.21412, 0, 1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(0.375185, 1.1547, 1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(-0.982247, 0.713644, 1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(-0.982247, -0.713644, 1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(0.375185, -1.1547, 1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(1.96449, 0, 0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(0.607062, 1.86835, 0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(-1.58931, 1.1547, 0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(-1.58931, -1.1547, 0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(0.607062, -1.86835, 0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(1.58931, 1.1547, -0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(-0.607062, 1.86835, -0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(-1.96449, 0, -0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(-0.607062, -1.86835, -0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(1.58931, -1.1547, -0.375185);
  aDodecahedron->GetPoints()->InsertNextPoint(0.982247, 0.713644, -1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(-0.375185, 1.1547, -1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(-1.21412, 0, -1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(-0.375185, -1.1547, -1.58931);
  aDodecahedron->GetPoints()->InsertNextPoint(0.982247, -0.713644, -1.58931);

  vtkIdType faces[73] = {
    12,                   // number of faces
    5, 0, 1, 2, 3, 4,     // number of ids on face, ids
    5, 0, 5, 10, 6, 1,    //
    5, 1, 6, 11, 7, 2,    //
    5, 2, 7, 12, 8, 3,    //
    5, 3, 8, 13, 9, 4,    //
    5, 4, 9, 14, 5, 0,    //
    5, 15, 10, 5, 14, 19, //
    5, 16, 11, 6, 10, 15, //
    5, 17, 12, 7, 11, 16, //
    5, 18, 13, 8, 12, 17, //
    5, 19, 14, 9, 13, 18, //
    5, 19, 18, 17, 16, 15 //
  };

  aDodecahedron->SetFaces(faces);
  aDodecahedron->Initialize();

  return aDodecahedron;
}

vtkSmartPointer<vtkPentagonalPrism> MakePentagonalPrism()
{
  auto aPentagonalPrism = vtkSmartPointer<vtkPentagonalPrism>::New();

  for (int i = 0; i < aPentagonalPrism->GetNumberOfPoints(); ++i)
  {
    aPentagonalPrism->GetPointIds()->SetId(i, i);
  }

  aPentagonalPrism->GetPoints()->SetPoint(0, 11, 10, 10);
  aPentagonalPrism->GetPoints()->SetPoint(1, 13, 10, 10);
  aPentagonalPrism->GetPoints()->SetPoint(2, 14, 12, 10);
  aPentagonalPrism->GetPoints()->SetPoint(3, 12, 14, 10);
  aPentagonalPrism->GetPoints()->SetPoint(4, 10, 12, 10);
  aPentagonalPrism->GetPoints()->SetPoint(5, 11, 10, 14);
  aPentagonalPrism->GetPoints()->SetPoint(6, 13, 10, 14);
  aPentagonalPrism->GetPoints()->SetPoint(7, 14, 12, 14);
  aPentagonalPrism->GetPoints()->SetPoint(8, 12, 14, 14);
  aPentagonalPrism->GetPoints()->SetPoint(9, 10, 12, 14);

  return aPentagonalPrism;
}

vtkSmartPointer<vtkHexagonalPrism> MakeHexagonalPrism()
{
  auto aHexagonalPrism = vtkSmartPointer<vtkHexagonalPrism>::New();

  for (int i = 0; i < aHexagonalPrism->GetNumberOfPoints(); ++i)
  {
    aHexagonalPrism->GetPointIds()->SetId(i, i);
  }

  aHexagonalPrism->GetPoints()->SetPoint(0, 11, 10, 10);
  aHexagonalPrism->GetPoints()->SetPoint(1, 13, 10, 10);
  aHexagonalPrism->GetPoints()->SetPoint(2, 14, 12, 10);
  aHexagonalPrism->GetPoints()->SetPoint(3, 13, 14, 10);
  aHexagonalPrism->GetPoints()->SetPoint(4, 11, 14, 10);
  aHexagonalPrism->GetPoints()->SetPoint(5, 10, 12, 10);
  aHexagonalPrism->GetPoints()->SetPoint(6, 11, 10, 14);
  aHexagonalPrism->GetPoints()->SetPoint(7, 13, 10, 14);
  aHexagonalPrism->GetPoints()->SetPoint(8, 14, 12, 14);
  aHexagonalPrism->GetPoints()->SetPoint(9, 13, 14, 14);
  aHexagonalPrism->GetPoints()->SetPoint(10, 11, 14, 14);
  aHexagonalPrism->GetPoints()->SetPoint(11, 10, 12, 14);

  return aHexagonalPrism;
}

template <typename T>
int TestOneCell(VTKCellType cellType, vtkSmartPointer<T> aCell, int linear)
{
  int status = 0;
  std::cout << "Testing " << aCell->GetClassName() << std::endl;

  std::cout << "  Testing Print of an uninitialized cell...";
  std::ostringstream cellPrint;
  aCell->Print(cellPrint);
  std::cout << "PASSED" << std::endl;

  std::cout << "  Testing GetCellType...";
  if (cellType != aCell->GetCellType())
  {
    std::cout << "Expected " << cellType << " but got " << aCell->GetCellType() << " FAILED"
              << std::endl;
    ++status;
  }
  else
  {
    std::cout << "PASSED" << std::endl;
  }

  std::cout << "  Testing GetCellDimension...";
  std::cout << aCell->GetCellDimension();
  std::cout << "...PASSED" << std::endl;

  std::cout << "  Testing IsLinear...";
  if (aCell->IsLinear() != 1 && linear)
  {
    ++status;
    std::cout << "...FAILED" << std::endl;
  }
  else
  {
    std::cout << "...PASSED" << std::endl;
  }

  std::cout << "  Testing IsPrimaryCell...";
  std::cout << aCell->IsPrimaryCell();
  std::cout << "...PASSED" << std::endl;

  std::cout << "  Testing IsExplicitCell...";
  std::cout << aCell->IsExplicitCell();
  std::cout << "...PASSED" << std::endl;

  std::cout << "  Testing RequiresInitialization...";
  std::cout << aCell->RequiresInitialization();
  std::cout << "...PASSED" << std::endl;

  std::cout << "  Testing RequiresExplicitFaceRepresentation...";
  std::cout << aCell->RequiresExplicitFaceRepresentation();
  std::cout << "...PASSED" << std::endl;

  if (aCell->RequiresInitialization())
  {
    aCell->Initialize();
  }
  std::cout << "  Testing GetNumberOfPoints...";
  std::cout << aCell->GetNumberOfPoints();
  std::cout << "...PASSED" << std::endl;

  std::cout << "  Testing GetNumberOfEdges...";
  std::cout << aCell->GetNumberOfEdges();
  std::cout << "...PASSED" << std::endl;

  std::cout << "  Testing GetNumberOfFaces...";
  std::cout << aCell->GetNumberOfFaces();
  std::cout << "...PASSED" << std::endl;

  if (std::string(aCell->GetClassName()) != "vtkEmptyCell" &&
    std::string(aCell->GetClassName()) != "vtkVertex" &&
    std::string(aCell->GetClassName()) != "vtkPolyhedron")
  {
    std::cout << "  Testing GetParametricCoords...";
    double* parametricCoords = aCell->GetParametricCoords();
    if (aCell->IsPrimaryCell() && parametricCoords == nullptr)
    {
      ++status;
      std::cout << "...FAILED" << std::endl;
    }
    else if (parametricCoords)
    {
      std::vector<double> pweights(aCell->GetNumberOfPoints());
      // The pcoords should correspond to the cell points
      for (int p = 0; p < aCell->GetNumberOfPoints(); ++p)
      {
        double vertex[3];
        aCell->GetPoints()->GetPoint(p, vertex);
        int subId = 0;
        double x[3];
        aCell->EvaluateLocation(subId, parametricCoords + 3 * p, x, pweights.data());
        if (!vtkMathUtilities::FuzzyCompare(x[0], vertex[0], 1.e-3) ||
          !vtkMathUtilities::FuzzyCompare(x[1], vertex[1], 1.e-3) ||
          !vtkMathUtilities::FuzzyCompare(x[2], vertex[2], 1.e-3))
        {
          std::cout << "EvaluateLocation failed...";
          std::cout << "pcoords[" << p << "]: " << parametricCoords[3 * p] << " "
                    << parametricCoords[3 * p + 1] << " " << parametricCoords[3 * p + 2]
                    << std::endl;
          std::cout << "x[" << p << "]: " << x[0] << " " << x[1] << " " << x[2] << std::endl;
          std::cout << "...FAILED" << std::endl;
          ++status;
        }
      }
      std::cout << "...PASSED" << std::endl;
    }
  }
  std::cout << "  Testing GetBounds...";
  double bounds[6];
  aCell->GetBounds(bounds);
  std::cout << bounds[0] << "," << bounds[1] << " " << bounds[2] << "," << bounds[3] << " "
            << bounds[4] << "," << bounds[5];
  std::cout << "...PASSED" << std::endl;

  if (aCell->GetNumberOfPoints() > 0)
  {
    std::cout << "  Testing GetParametricCenter...";
    double pcenter[3], center[3];
    pcenter[0] = pcenter[1] = pcenter[2] = -12345.0;
    aCell->GetParametricCenter(pcenter);
    std::cout << pcenter[0] << ", " << pcenter[1] << ", " << pcenter[2];
    std::vector<double> cweights(aCell->GetNumberOfPoints());
    int pSubId = 0;
    aCell->EvaluateLocation(pSubId, pcenter, center, cweights.data());
    std::cout << " -> " << center[0] << ", " << center[1] << ", " << center[2];
    if (center[0] < bounds[0] || center[0] > bounds[1] || center[1] < bounds[2] ||
      center[1] > bounds[3] || center[2] < bounds[4] || center[2] > bounds[5])
    {
      std::cout << " The computed center is not within the bounds of the cell" << std::endl;
      std::cout << "bounds: " << bounds[0] << "," << bounds[1] << " " << bounds[2] << ","
                << bounds[3] << " " << bounds[4] << "," << bounds[5] << std::endl;
      std::cout << "parametric center " << pcenter[0] << ", " << pcenter[1] << ", " << pcenter[2]
                << " "
                << "center: " << center[0] << ", " << center[1] << ", " << center[2] << std::endl;
      std::cout << "...FAILED" << std::endl;
    }
    else
    {
      std::cout << "...PASSED" << std::endl;
    }
  }

  std::cout << "  Testing GetParametricDistance...";
  double pcenter[3];
  aCell->GetParametricCenter(pcenter);
  double pd = aCell->GetParametricDistance(pcenter);
  if (pd == 0.0)
  {
    std::cout << "...PASSED" << std::endl;
  }
  else
  {
    ++status;
    std::cout << "...FAILED" << std::endl;
  }

  std::cout << "  Testing CellBoundaries...";
  vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
  int cellStatus = aCell->CellBoundary(0, pcenter, cellIds);
  if (aCell->GetCellDimension() > 0 && cellStatus != 1)
  {
    ++status;
    std::cout << "FAILED" << std::endl;
  }
  else
  {
    for (int c = 0; c < cellIds->GetNumberOfIds(); ++c)
    {
      std::cout << " " << cellIds->GetId(c) << ", ";
    }
    std::cout << "PASSED" << std::endl;
  }

  if (aCell->GetNumberOfPoints() > 0 && strcmp(aCell->GetClassName(), "vtkQuadraticEdge") != 0)
  {
    std::cout << "  Testing Derivatives...";
    // Create scalars and set first scalar to 1.0
    std::vector<double> scalars(aCell->GetNumberOfPoints());
    scalars[0] = 1.0;
    for (int s = 1; s < aCell->GetNumberOfPoints(); ++s)
    {
      scalars[s] = 0.0;
    }
    std::vector<double> derivs(3, -12345.0);
    aCell->Derivatives(0, pcenter, &(*scalars.begin()), 1, &(*derivs.begin()));
    if (derivs[0] == -12345. && derivs[1] == -12345. && derivs[2] == -12345.)
    {
      std::cout << " not computed";
    }
    else
    {
      std::cout << " " << derivs[0] << " " << derivs[1] << " " << derivs[2] << " ";
    }
    std::cout << "...PASSED" << std::endl;
  }

  std::cout << "  Testing EvaluateLocation vertex matches pcoord...";
  int status5 = 0;
  double* locations = aCell->GetParametricCoords();
  if (locations)
  {
    std::vector<double> lweights(aCell->GetNumberOfPoints());
    for (int l = 0; l < aCell->GetNumberOfPoints(); ++l)
    {
      double point[3];
      double vertex[3];
      aCell->GetPoints()->GetPoint(l, vertex);
      int subId = 0;
      aCell->EvaluateLocation(subId, locations + 3 * l, point, lweights.data());
      for (int v = 0; v < 3; ++v)
      {
        if (!vtkMathUtilities::FuzzyCompare(point[v], vertex[v], 1.e-3))
        {
          std::cout << " " << point[0] << ", " << point[1] << ", " << point[2]
                    << " != " << vertex[0] << ", " << vertex[1] << ", " << vertex[2] << " ";
          std::cout << "eps ratio is: "
                    << (point[v] - vertex[v]) / std::numeric_limits<double>::epsilon() << std::endl;

          ++status5;
          break;
        }
      }
    }
  }
  if (status5)
  {
    std::cout << "...FAILED" << std::endl;
    ++status;
  }
  else
  {
    std::cout << "...PASSED" << std::endl;
  }

  std::cout << "  Testing EvaluatePosition pcoord matches vertex...";
  // Each vertex should correspond to a pcoord.
  int subId = 0;
  int status6 = 0;
  std::vector<double> weights(aCell->GetNumberOfPoints());
  double* vlocations = aCell->GetParametricCoords();
  if (vlocations)
  {
    for (int i = 0; i < aCell->GetNumberOfPoints(); ++i)
    {
      int status61 = 0;
      double closestPoint[3];
      double point[3];
      double pcoords[3];
      double dist2;
      aCell->GetPoints()->GetPoint(i, point);
      aCell->EvaluatePosition(point, closestPoint, subId, pcoords, dist2, weights.data());
      for (int v = 0; v < 3; ++v)
      {
        if (!vtkMathUtilities::FuzzyCompare(*(vlocations + 3 * i + v), pcoords[v], 1.e-3))
        {
          ++status61;
        }
      }
      if (status61)
      {
        std::cout << std::endl
                  << "Point id: " << i << ": " << *(vlocations + 3 * i + 0) << ", "
                  << *(vlocations + 3 * i + 1) << ", " << *(vlocations + 3 * i + 2)
                  << " != " << pcoords[0] << ", " << pcoords[1] << ", " << pcoords[2] << " ";
        ++status6;
      }
    }
  }
  if (status6)
  {
    ++status;
    std::cout << "...FAILED" << std::endl;
  }
  else
  {
    std::cout << "...PASSED" << std::endl;
  }

  std::cout << "  Testing EvaluatePosition in/out test...";

  int status2 = 0;
  std::vector<std::vector<double>> testPoints;
  std::vector<int> inOuts;
  std::vector<std::string> typePoint;

  // First test cell points
  for (int i = 0; i < aCell->GetNumberOfPoints(); ++i)
  {
    std::vector<double> point(3);
    aCell->GetPoints()->GetPoint(i, &(*point.begin()));
    testPoints.push_back(point);
    inOuts.push_back(1);
    typePoint.emplace_back("cell point");
  }
  // Then test center of cell
  if (aCell->GetNumberOfPoints() > 0)
  {
    std::vector<double> tCenter(3);
    aCell->EvaluateLocation(subId, pcenter, &(*tCenter.begin()), weights.data());
    testPoints.push_back(tCenter);
    inOuts.push_back(1);
    typePoint.emplace_back("cell center");
    // Test a point above the cell
    if (aCell->GetCellDimension() == 2)
    {
      std::vector<double> above(3);
      above[0] = tCenter[0];
      above[1] = tCenter[1];
      above[2] = tCenter[2] + aCell->GetLength2();
      testPoints.push_back(above);
      inOuts.push_back(0);
      typePoint.emplace_back("point above cell");
    }
  }

  // Test points at the center of each edge
  for (int e = 0; e < aCell->GetNumberOfEdges(); ++e)
  {
    std::vector<double> eCenter(3);
    vtkCell* c = aCell->GetEdge(e);
    c->GetParametricCenter(pcenter);
    c->EvaluateLocation(subId, pcenter, &(*eCenter.begin()), weights.data());
    testPoints.push_back(eCenter);
    typePoint.emplace_back("edge center");
    inOuts.push_back(1);
  }

  // Test points at the center of each face
  for (int f = 0; f < aCell->GetNumberOfFaces(); ++f)
  {
    std::vector<double> fCenter(3);
    vtkCell* c = aCell->GetFace(f);
    c->GetParametricCenter(pcenter);
    c->EvaluateLocation(subId, pcenter, &(*fCenter.begin()), weights.data());
    testPoints.push_back(fCenter);
    inOuts.push_back(1);
    typePoint.emplace_back("face center");
  }

  // Test a point outside the cell
  if (aCell->GetNumberOfPoints() > 0)
  {
    std::vector<double> outside(3, -12345.0);
    testPoints.push_back(outside);
    inOuts.push_back(0);
    typePoint.emplace_back("outside point");
  }
  for (size_t p = 0; p < testPoints.size(); ++p)
  {
    double closestPoint[3], pcoords[3], dist2;
    int inOut = aCell->EvaluatePosition(
      &(*testPoints[p].begin()), closestPoint, subId, pcoords, dist2, weights.data());
    if ((inOut == 0 || inOut == -1) && inOuts[p] == 0)
    {
      continue;
    }
    else if (inOut == 1 && dist2 == 0.0 && inOuts[p] == 1)
    {
      continue;
    }
    else if (inOut == 1 && dist2 != 0.0 && inOuts[p] == 0)
    {
      continue;
    }
    // inOut failed
    std::cout << typePoint[p] << " failed inOut: " << inOut << " "
              << "point: " << testPoints[p][0] << ", " << testPoints[p][1] << ", "
              << testPoints[p][2] << "-> "
              << "pcoords: " << pcoords[0] << ", " << pcoords[1] << ", " << pcoords[2] << ": "
              << "closestPoint: " << closestPoint[0] << ", " << closestPoint[1] << ", "
              << closestPoint[2] << " "
              << "dist2: " << dist2;
    std::cout << " weights: ";
    for (int w = 0; w < aCell->GetNumberOfPoints(); ++w)
    {
      std::cout << weights[w] << " ";
    }
    std::cout << std::endl;
    status2 += 1;
  }
  if (status2)
  {
    ++status;
    std::cout << "FAILED" << std::endl;
  }
  else
  {
    std::cout << "PASSED" << std::endl;
  }

  if (aCell->GetNumberOfPoints() > 0 && aCell->GetCellDimension() > 0)
  {
    std::cout << "  Testing IntersectWithLine...";
    double tol = 1.e-5;
    double t;
    double startPoint[3];
    double endPoint[3];
    double intersection[3], pintersection[3];
    aCell->GetParametricCenter(pcenter);
    aCell->EvaluateLocation(subId, pcenter, startPoint, weights.data());
    endPoint[0] = startPoint[0];
    endPoint[1] = startPoint[1];
    endPoint[2] = startPoint[2] + aCell->GetLength2();
    startPoint[2] = startPoint[2] - aCell->GetLength2();
    int status3 = 0;
    int result =
      aCell->IntersectWithLine(startPoint, endPoint, tol, t, intersection, pintersection, subId);
    if (result == 0)
    {
      ++status3;
    }
    else
    {
      std::cout << " t: " << t << " ";
    }
    startPoint[2] = endPoint[2] + aCell->GetLength2();
    result =
      aCell->IntersectWithLine(startPoint, endPoint, tol, t, intersection, pintersection, subId);
    if (result == 1)
    {
      ++status3;
    }

    if (status3 != 0)
    {
      ++status;
      std::cout << "...FAILED" << std::endl;
    }
    else
    {
      std::cout << "...PASSED" << std::endl;
    }
  }

  // Triangulate
  std::cout << "  Testing Triangulate...";
  int index = 0;
  vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
  ptIds->SetNumberOfIds(100);
  vtkSmartPointer<vtkPoints> triPoints = vtkSmartPointer<vtkPoints>::New();
  aCell->Triangulate(index, ptIds, triPoints);
  int pts = ptIds->GetNumberOfIds();
  if (aCell->GetCellDimension() == 0)
  {
    std::cout << "Generated " << pts << " Points";
  }
  else if (aCell->GetCellDimension() == 1)
  {
    std::cout << "Generated " << pts / 2 << " Lines";
  }
  else if (aCell->GetCellDimension() == 2)
  {
    std::cout << "Generated " << pts / 3 << " Triangles";
  }
  else if (aCell->GetCellDimension() == 3)
  {
    std::cout << "Generated " << pts / 4 << " Tetra";
  }
  std::cout << "...PASSED" << std::endl;

  if (status)
  {
    std::cout << aCell->GetClassName() << " FAILED" << std::endl;
  }
  else
  {
    std::cout << aCell->GetClassName() << " PASSED" << std::endl;
  }
  return status;
}