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

  Program:   Visualization Toolkit
  Module:    vtkQuad.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 "vtkQuad.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkLine.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkTriangle.h"

vtkStandardNewMacro(vtkQuad);

static const double VTK_DIVERGED = 1.e6;

//------------------------------------------------------------------------------
struct IntersectionStruct
{
  bool Intersected = false;
  int SubId = -1;
  double X[3] = { 0.0, 0.0, 0.0 };
  double PCoords[3] = { 0.0, 0.0, 0.0 };
  double T = -1.0;

  operator bool() { return this->Intersected; }

  void CopyValues(double& t, double* x, double* pcoords, int& subId) const
  {
    t = this->T;
    subId = this->SubId;
    for (int i = 0; i < 3; ++i)
    {
      x[i] = this->X[i];
      pcoords[i] = this->PCoords[i];
    }
  }

  static IntersectionStruct CellIntersectWithLine(
    vtkCell* cell, const double* p1, const double* p2, double tol)
  {
    IntersectionStruct res;
    res.Intersected = cell->IntersectWithLine(p1, p2, tol, res.T, res.X, res.PCoords, res.SubId);
    return res;
  }
};

//------------------------------------------------------------------------------
// Construct the quad with four points.
vtkQuad::vtkQuad()
{
  this->Points->SetNumberOfPoints(4);
  this->PointIds->SetNumberOfIds(4);
  for (int i = 0; i < 4; i++)
  {
    this->Points->SetPoint(i, 0.0, 0.0, 0.0);
    this->PointIds->SetId(i, 0);
  }
  this->Line = vtkLine::New();
  this->Triangle = vtkTriangle::New();
}

//------------------------------------------------------------------------------
vtkQuad::~vtkQuad()
{
  this->Line->Delete();
  this->Triangle->Delete();
}

//------------------------------------------------------------------------------
static const int VTK_QUAD_MAX_ITERATION = 20;
static const double VTK_QUAD_CONVERGED = 1.e-04;

inline static void ComputeNormal(
  vtkQuad* self, double pt1[3], double pt2[3], double pt3[3], double n[3])
{
  vtkTriangle::ComputeNormal(pt1, pt2, pt3, n);

  // If first three points are co-linear, then use fourth point
  //
  double pt4[3];
  if (n[0] == 0.0 && n[1] == 0.0 && n[2] == 0.0)
  {
    self->Points->GetPoint(3, pt4);
    vtkTriangle::ComputeNormal(pt2, pt3, pt4, n);
  }
}

//------------------------------------------------------------------------------
int vtkQuad::EvaluatePosition(const double x[3], double closestPoint[3], int& subId,
  double pcoords[3], double& dist2, double weights[])
{
  int i, j;
  double pt1[3], pt2[3], pt3[3], pt[3], n[3];
  double det;
  double maxComponent;
  int idx = 0, indices[2];
  int iteration, converged;
  double params[2];
  double fcol[2], rcol[2], scol[2], cp[3];
  double derivs[8];

  subId = 0;
  pcoords[0] = pcoords[1] = params[0] = params[1] = 0.5;
  pcoords[2] = 0.0;

  // Get normal for quadrilateral
  //
  this->Points->GetPoint(0, pt1);
  this->Points->GetPoint(1, pt2);
  this->Points->GetPoint(2, pt3);
  ComputeNormal(this, pt1, pt2, pt3, n);

  // Project point to plane
  //
  vtkPlane::ProjectPoint(x, pt1, n, cp);

  // Construct matrices.  Since we have over determined system, need to find
  // which 2 out of 3 equations to use to develop equations. (Any 2 should
  // work since we've projected point to plane.)
  //
  for (maxComponent = 0.0, i = 0; i < 3; i++)
  {
    if (fabs(n[i]) > maxComponent)
    {
      maxComponent = fabs(n[i]);
      idx = i;
    }
  }
  for (j = 0, i = 0; i < 3; i++)
  {
    if (i != idx)
    {
      indices[j++] = i;
    }
  }

  // Use Newton's method to solve for parametric coordinates
  //
  for (iteration = converged = 0; !converged && (iteration < VTK_QUAD_MAX_ITERATION); iteration++)
  {
    //  calculate element interpolation functions and derivatives
    //
    vtkQuad::InterpolationFunctions(pcoords, weights);
    vtkQuad::InterpolationDerivs(pcoords, derivs);

    //  calculate newton functions
    //
    for (i = 0; i < 2; i++)
    {
      fcol[i] = rcol[i] = scol[i] = 0.0;
    }
    for (i = 0; i < 4; i++)
    {
      this->Points->GetPoint(i, pt);
      for (j = 0; j < 2; j++)
      {
        fcol[j] += pt[indices[j]] * weights[i];
        rcol[j] += pt[indices[j]] * derivs[i];
        scol[j] += pt[indices[j]] * derivs[i + 4];
      }
    }

    for (j = 0; j < 2; j++)
    {
      fcol[j] -= cp[indices[j]];
    }

    //  compute determinants and generate improvements
    //
    if ((det = vtkMath::Determinant2x2(rcol, scol)) == 0.0)
    {
      return -1;
    }

    pcoords[0] = params[0] - vtkMath::Determinant2x2(fcol, scol) / det;
    pcoords[1] = params[1] - vtkMath::Determinant2x2(rcol, fcol) / det;

    //  check for convergence
    //
    if (((fabs(pcoords[0] - params[0])) < VTK_QUAD_CONVERGED) &&
      ((fabs(pcoords[1] - params[1])) < VTK_QUAD_CONVERGED))
    {
      converged = 1;
    }
    // Test for bad divergence (S.Hirschberg 11.12.2001)
    else if ((fabs(pcoords[0]) > VTK_DIVERGED) || (fabs(pcoords[1]) > VTK_DIVERGED))
    {
      return -1;
    }

    //  if not converged, repeat
    //
    else
    {
      params[0] = pcoords[0];
      params[1] = pcoords[1];
    }
  }

  //  if not converged, set the parametric coordinates to arbitrary values
  //  outside of element
  //
  if (!converged)
  {
    return -1;
  }

  vtkQuad::InterpolationFunctions(pcoords, weights);

  if (pcoords[0] >= -0.001 && pcoords[0] <= 1.001 && pcoords[1] >= -0.001 && pcoords[1] <= 1.001)
  {
    if (closestPoint)
    {
      dist2 = vtkMath::Distance2BetweenPoints(cp, x); // projection distance
      closestPoint[0] = cp[0];
      closestPoint[1] = cp[1];
      closestPoint[2] = cp[2];
    }
    return 1;
  }
  else
  {
    double t;
    double pt4[3];

    if (closestPoint)
    {
      this->Points->GetPoint(3, pt4);

      if (pcoords[0] < 0.0 && pcoords[1] < 0.0)
      {
        dist2 = vtkMath::Distance2BetweenPoints(x, pt1);
        for (i = 0; i < 3; i++)
        {
          closestPoint[i] = pt1[i];
        }
      }
      else if (pcoords[0] > 1.0 && pcoords[1] < 0.0)
      {
        dist2 = vtkMath::Distance2BetweenPoints(x, pt2);
        for (i = 0; i < 3; i++)
        {
          closestPoint[i] = pt2[i];
        }
      }
      else if (pcoords[0] > 1.0 && pcoords[1] > 1.0)
      {
        dist2 = vtkMath::Distance2BetweenPoints(x, pt3);
        for (i = 0; i < 3; i++)
        {
          closestPoint[i] = pt3[i];
        }
      }
      else if (pcoords[0] < 0.0 && pcoords[1] > 1.0)
      {
        dist2 = vtkMath::Distance2BetweenPoints(x, pt4);
        for (i = 0; i < 3; i++)
        {
          closestPoint[i] = pt4[i];
        }
      }
      else if (pcoords[0] < 0.0)
      {
        dist2 = vtkLine::DistanceToLine(x, pt1, pt4, t, closestPoint);
      }
      else if (pcoords[0] > 1.0)
      {
        dist2 = vtkLine::DistanceToLine(x, pt2, pt3, t, closestPoint);
      }
      else if (pcoords[1] < 0.0)
      {
        dist2 = vtkLine::DistanceToLine(x, pt1, pt2, t, closestPoint);
      }
      else if (pcoords[1] > 1.0)
      {
        dist2 = vtkLine::DistanceToLine(x, pt3, pt4, t, closestPoint);
      }
    }
    return 0;
  }
}

//------------------------------------------------------------------------------
void vtkQuad::EvaluateLocation(
  int& vtkNotUsed(subId), const double pcoords[3], double x[3], double* weights)
{
  int i, j;
  double pt[3];

  vtkQuad::InterpolationFunctions(pcoords, weights);

  x[0] = x[1] = x[2] = 0.0;
  for (i = 0; i < 4; i++)
  {
    this->Points->GetPoint(i, pt);
    for (j = 0; j < 3; j++)
    {
      x[j] += pt[j] * weights[i];
    }
  }
}

//------------------------------------------------------------------------------
// Compute iso-parametric interpolation functions
//
void vtkQuad::InterpolationFunctions(const double pcoords[3], double sf[4])
{
  double rm, sm;

  rm = 1. - pcoords[0];
  sm = 1. - pcoords[1];

  sf[0] = rm * sm;
  sf[1] = pcoords[0] * sm;
  sf[2] = pcoords[0] * pcoords[1];
  sf[3] = rm * pcoords[1];
}

//------------------------------------------------------------------------------
void vtkQuad::InterpolationDerivs(const double pcoords[3], double derivs[8])
{
  double rm, sm;

  rm = 1. - pcoords[0];
  sm = 1. - pcoords[1];

  derivs[0] = -sm;
  derivs[1] = sm;
  derivs[2] = pcoords[1];
  derivs[3] = -pcoords[1];
  derivs[4] = -rm;
  derivs[5] = -pcoords[0];
  derivs[6] = pcoords[0];
  derivs[7] = rm;
}

//------------------------------------------------------------------------------
int vtkQuad::CellBoundary(int vtkNotUsed(subId), const double pcoords[3], vtkIdList* pts)
{
  double t1 = pcoords[0] - pcoords[1];
  double t2 = 1.0 - pcoords[0] - pcoords[1];

  pts->SetNumberOfIds(2);

  // compare against two lines in parametric space that divide element
  // into four pieces.
  if (t1 >= 0.0 && t2 >= 0.0)
  {
    pts->SetId(0, this->PointIds->GetId(0));
    pts->SetId(1, this->PointIds->GetId(1));
  }

  else if (t1 >= 0.0 && t2 < 0.0)
  {
    pts->SetId(0, this->PointIds->GetId(1));
    pts->SetId(1, this->PointIds->GetId(2));
  }

  else if (t1 < 0.0 && t2 < 0.0)
  {
    pts->SetId(0, this->PointIds->GetId(2));
    pts->SetId(1, this->PointIds->GetId(3));
  }

  else //( t1 < 0.0 && t2 >= 0.0 )
  {
    pts->SetId(0, this->PointIds->GetId(3));
    pts->SetId(1, this->PointIds->GetId(0));
  }

  if (pcoords[0] < 0.0 || pcoords[0] > 1.0 || pcoords[1] < 0.0 || pcoords[1] > 1.0)
  {
    return 0;
  }
  else
  {
    return 1;
  }
}

//------------------------------------------------------------------------------
// Marching (convex) quadrilaterals
//
namespace
{ // required so we don't violate ODR
constexpr vtkIdType edges[4][2] = {
  { 0, 1 },
  { 1, 2 },
  { 3, 2 },
  { 0, 3 },
};

typedef int EDGE_LIST;
struct LINE_CASES_t
{
  EDGE_LIST edges[5];
};
using LINE_CASES = struct LINE_CASES_t;

LINE_CASES lineCases[] = {
  { { -1, -1, -1, -1, -1 } },
  { { 0, 3, -1, -1, -1 } },
  { { 1, 0, -1, -1, -1 } },
  { { 1, 3, -1, -1, -1 } },
  { { 2, 1, -1, -1, -1 } },
  { { 0, 3, 2, 1, -1 } },
  { { 2, 0, -1, -1, -1 } },
  { { 2, 3, -1, -1, -1 } },
  { { 3, 2, -1, -1, -1 } },
  { { 0, 2, -1, -1, -1 } },
  { { 1, 0, 3, 2, -1 } },
  { { 1, 2, -1, -1, -1 } },
  { { 3, 1, -1, -1, -1 } },
  { { 0, 1, -1, -1, -1 } },
  { { 3, 0, -1, -1, -1 } },
  { { -1, -1, -1, -1, -1 } },
};
}

//------------------------------------------------------------------------------
const vtkIdType* vtkQuad::GetEdgeArray(vtkIdType edgeId)
{
  return edges[edgeId];
}

//------------------------------------------------------------------------------
void vtkQuad::Contour(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
  vtkCellArray* verts, vtkCellArray* lines, vtkCellArray* vtkNotUsed(polys), vtkPointData* inPd,
  vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId, vtkCellData* outCd)
{
  static const int CASE_MASK[4] = { 1, 2, 4, 8 };
  LINE_CASES* lineCase;
  EDGE_LIST* edge;
  int i, j, index;
  const vtkIdType* vert;
  int newCellId;
  vtkIdType pts[2];
  int e1, e2;
  double t, x1[3], x2[3], x[3], deltaScalar;
  vtkIdType offset = verts->GetNumberOfCells();

  // Build the case table
  for (i = 0, index = 0; i < 4; i++)
  {
    if (cellScalars->GetComponent(i, 0) >= value)
    {
      index |= CASE_MASK[i];
    }
  }

  lineCase = lineCases + index;
  edge = lineCase->edges;

  for (; edge[0] > -1; edge += 2)
  {
    for (i = 0; i < 2; i++) // insert line
    {
      vert = edges[edge[i]];
      // calculate a preferred interpolation direction
      deltaScalar = (cellScalars->GetComponent(vert[1], 0) - cellScalars->GetComponent(vert[0], 0));
      if (deltaScalar > 0)
      {
        e1 = vert[0];
        e2 = vert[1];
      }
      else
      {
        e1 = vert[1];
        e2 = vert[0];
        deltaScalar = -deltaScalar;
      }

      // linear interpolation
      if (deltaScalar == 0.0)
      {
        t = 0.0;
      }
      else
      {
        t = (value - cellScalars->GetComponent(e1, 0)) / deltaScalar;
      }

      this->Points->GetPoint(e1, x1);
      this->Points->GetPoint(e2, x2);

      for (j = 0; j < 3; j++)
      {
        x[j] = x1[j] + t * (x2[j] - x1[j]);
      }
      if (locator->InsertUniquePoint(x, pts[i]))
      {
        if (outPd)
        {
          vtkIdType p1 = this->PointIds->GetId(e1);
          vtkIdType p2 = this->PointIds->GetId(e2);
          outPd->InterpolateEdge(inPd, pts[i], p1, p2, t);
        }
      }
    }
    // check for degenerate line
    if (pts[0] != pts[1])
    {
      newCellId = offset + lines->InsertNextCell(2, pts);
      if (outCd)
      {
        outCd->CopyData(inCd, cellId, newCellId);
      }
    }
  }
}

//------------------------------------------------------------------------------
vtkCell* vtkQuad::GetEdge(int edgeId)
{
  int edgeIdPlus1 = edgeId + 1;

  if (edgeIdPlus1 > 3)
  {
    edgeIdPlus1 = 0;
  }

  // load point id's
  this->Line->PointIds->SetId(0, this->PointIds->GetId(edgeId));
  this->Line->PointIds->SetId(1, this->PointIds->GetId(edgeIdPlus1));

  // load coordinates
  this->Line->Points->SetPoint(0, this->Points->GetPoint(edgeId));
  this->Line->Points->SetPoint(1, this->Points->GetPoint(edgeIdPlus1));

  return this->Line;
}

//------------------------------------------------------------------------------
// Intersect plane; see whether point is in quadrilateral. This code
// splits the quad into two triangles and intersects them (because the
// quad may be non-planar).
//
int vtkQuad::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
  double x[3], double pcoords[3], int& subId)
{
  int diagonalCase;
  double d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0), this->Points->GetPoint(2));
  double d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1), this->Points->GetPoint(3));
  subId = 0;

  // Figure out how to uniquely tessellate the quad. Watch out for
  // equivalent triangulations (i.e., the triangulation is equivalent
  // no matter where the diagonal). In this case use the point ids as
  // a tie breaker to ensure unique triangulation across the quad.
  //
  if (d1 == d2) // rare case; discriminate based on point id
  {
    int i, id, maxId = 0, maxIdx = 0;
    for (i = 0; i < 4; i++) // find the maximum id
    {
      if ((id = this->PointIds->GetId(i)) > maxId)
      {
        maxId = id;
        maxIdx = i;
      }
    }
    if (maxIdx == 0 || maxIdx == 2)
      diagonalCase = 0;
    else
      diagonalCase = 1;
  }
  else if (d1 < d2)
  {
    diagonalCase = 0;
  }
  else // d2 < d1
  {
    diagonalCase = 1;
  }

  // Note: in the following code the parametric coords must be adjusted to
  // reflect the use of the triangle parametric coordinate system.
  IntersectionStruct res;
  switch (diagonalCase)
  {
    case 0:
    {
      this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
      this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
      this->Triangle->Points->SetPoint(2, this->Points->GetPoint(2));
      IntersectionStruct firstIntersect =
        IntersectionStruct::CellIntersectWithLine(this->Triangle, p1, p2, tol);

      this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
      this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
      this->Triangle->Points->SetPoint(2, this->Points->GetPoint(0));
      IntersectionStruct secondIntersect =
        IntersectionStruct::CellIntersectWithLine(this->Triangle, p1, p2, tol);

      bool useFirstIntersection = (firstIntersect && secondIntersect)
        ? (firstIntersect.T <= secondIntersect.T)
        : firstIntersect;
      bool useSecondIntersection = (firstIntersect && secondIntersect)
        ? (secondIntersect.T < firstIntersect.T)
        : secondIntersect;

      if (useFirstIntersection)
      {
        res = firstIntersect;
        res.PCoords[0] += res.PCoords[1];
      }
      else if (useSecondIntersection)
      {
        res = secondIntersect;
        res.PCoords[0] = 1.0 - (res.PCoords[0] + res.PCoords[1]);
        res.PCoords[1] = 1.0 - res.PCoords[1];
      }
    }
    break;

    case 1:
    {
      this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
      this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
      this->Triangle->Points->SetPoint(2, this->Points->GetPoint(3));
      IntersectionStruct firstIntersect =
        IntersectionStruct::CellIntersectWithLine(this->Triangle, p1, p2, tol);

      this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
      this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
      this->Triangle->Points->SetPoint(2, this->Points->GetPoint(1));
      IntersectionStruct secondIntersect =
        IntersectionStruct::CellIntersectWithLine(this->Triangle, p1, p2, tol);

      bool useFirstIntersection = (firstIntersect && secondIntersect)
        ? (firstIntersect.T <= secondIntersect.T)
        : firstIntersect;
      bool useSecondIntersection = (firstIntersect && secondIntersect)
        ? (secondIntersect.T < firstIntersect.T)
        : secondIntersect;

      if (useFirstIntersection)
      {
        res = firstIntersect;
      }
      else if (useSecondIntersection)
      {
        res = secondIntersect;
        res.PCoords[0] = 1.0 - res.PCoords[0];
        res.PCoords[1] = 1.0 - res.PCoords[1];
      }
    }
    break;
  }

  if (res)
  {
    res.CopyValues(t, x, pcoords, subId);
  }

  return res.Intersected;
}

//------------------------------------------------------------------------------
int vtkQuad::Triangulate(int vtkNotUsed(index), vtkIdList* ptIds, vtkPoints* pts)
{
  double d1, d2;

  pts->Reset();
  ptIds->Reset();

  // use minimum diagonal (Delaunay triangles) - assumed convex
  d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0), this->Points->GetPoint(2));
  d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1), this->Points->GetPoint(3));

  if (d1 <= d2)
  {
    ptIds->InsertId(0, this->PointIds->GetId(0));
    pts->InsertPoint(0, this->Points->GetPoint(0));
    ptIds->InsertId(1, this->PointIds->GetId(1));
    pts->InsertPoint(1, this->Points->GetPoint(1));
    ptIds->InsertId(2, this->PointIds->GetId(2));
    pts->InsertPoint(2, this->Points->GetPoint(2));

    ptIds->InsertId(3, this->PointIds->GetId(0));
    pts->InsertPoint(3, this->Points->GetPoint(0));
    ptIds->InsertId(4, this->PointIds->GetId(2));
    pts->InsertPoint(4, this->Points->GetPoint(2));
    ptIds->InsertId(5, this->PointIds->GetId(3));
    pts->InsertPoint(5, this->Points->GetPoint(3));
  }
  else
  {
    ptIds->InsertId(0, this->PointIds->GetId(0));
    pts->InsertPoint(0, this->Points->GetPoint(0));
    ptIds->InsertId(1, this->PointIds->GetId(1));
    pts->InsertPoint(1, this->Points->GetPoint(1));
    ptIds->InsertId(2, this->PointIds->GetId(3));
    pts->InsertPoint(2, this->Points->GetPoint(3));

    ptIds->InsertId(3, this->PointIds->GetId(1));
    pts->InsertPoint(3, this->Points->GetPoint(1));
    ptIds->InsertId(4, this->PointIds->GetId(2));
    pts->InsertPoint(4, this->Points->GetPoint(2));
    ptIds->InsertId(5, this->PointIds->GetId(3));
    pts->InsertPoint(5, this->Points->GetPoint(3));
  }

  return 1;
}

//------------------------------------------------------------------------------
void vtkQuad::Derivatives(
  int vtkNotUsed(subId), const double pcoords[3], const double* values, int dim, double* derivs)
{
  double v0[2], v1[2], v2[2], v3[2], v10[3], v20[3], lenX;
  double x0[3], x1[3], x2[3], x3[3], n[3], vec20[3], vec30[3];
  double *J[2], J0[2], J1[2];
  double *JI[2], JI0[2], JI1[2];
  double funcDerivs[8], sum[2], dBydx, dBydy;
  int i, j;

  // Project points of quad into 2D system
  this->Points->GetPoint(0, x0);
  this->Points->GetPoint(1, x1);
  this->Points->GetPoint(2, x2);
  ComputeNormal(this, x0, x1, x2, n);
  this->Points->GetPoint(3, x3);

  for (i = 0; i < 3; i++)
  {
    v10[i] = x1[i] - x0[i];
    vec20[i] = x2[i] - x0[i];
    vec30[i] = x3[i] - x0[i];
  }

  vtkMath::Cross(n, v10, v20); // creates local y' axis

  if ((lenX = vtkMath::Normalize(v10)) <= 0.0 || vtkMath::Normalize(v20) <= 0.0) // degenerate
  {
    for (j = 0; j < dim; j++)
    {
      for (i = 0; i < 3; i++)
      {
        derivs[j * dim + i] = 0.0;
      }
    }
    return;
  }

  v0[0] = v0[1] = 0.0; // convert points to 2D (i.e., local system)
  v1[0] = lenX;
  v1[1] = 0.0;
  v2[0] = vtkMath::Dot(vec20, v10);
  v2[1] = vtkMath::Dot(vec20, v20);
  v3[0] = vtkMath::Dot(vec30, v10);
  v3[1] = vtkMath::Dot(vec30, v20);

  this->InterpolationDerivs(pcoords, funcDerivs);

  // Compute Jacobian and inverse Jacobian
  J[0] = J0;
  J[1] = J1;
  JI[0] = JI0;
  JI[1] = JI1;

  J[0][0] =
    v0[0] * funcDerivs[0] + v1[0] * funcDerivs[1] + v2[0] * funcDerivs[2] + v3[0] * funcDerivs[3];
  J[0][1] =
    v0[1] * funcDerivs[0] + v1[1] * funcDerivs[1] + v2[1] * funcDerivs[2] + v3[1] * funcDerivs[3];
  J[1][0] =
    v0[0] * funcDerivs[4] + v1[0] * funcDerivs[5] + v2[0] * funcDerivs[6] + v3[0] * funcDerivs[7];
  J[1][1] =
    v0[1] * funcDerivs[4] + v1[1] * funcDerivs[5] + v2[1] * funcDerivs[6] + v3[1] * funcDerivs[7];

  // Compute inverse Jacobian, return if Jacobian is singular
  if (!vtkMath::InvertMatrix(J, JI, 2))
  {
    for (j = 0; j < dim; j++)
    {
      for (i = 0; i < 3; i++)
      {
        derivs[j * dim + i] = 0.0;
      }
    }
    return;
  }

  // Loop over "dim" derivative values. For each set of values,
  // compute derivatives
  // in local system and then transform into modelling system.
  // First compute derivatives in local x'-y' coordinate system
  for (j = 0; j < dim; j++)
  {
    sum[0] = sum[1] = 0.0;
    for (i = 0; i < 4; i++) // loop over interp. function derivatives
    {
      sum[0] += funcDerivs[i] * values[dim * i + j];
      sum[1] += funcDerivs[4 + i] * values[dim * i + j];
    }
    dBydx = sum[0] * JI[0][0] + sum[1] * JI[0][1];
    dBydy = sum[0] * JI[1][0] + sum[1] * JI[1][1];

    // Transform into global system (dot product with global axes)
    derivs[3 * j] = dBydx * v10[0] + dBydy * v20[0];
    derivs[3 * j + 1] = dBydx * v10[1] + dBydy * v20[1];
    derivs[3 * j + 2] = dBydx * v10[2] + dBydy * v20[2];
  }
}

//------------------------------------------------------------------------------
// support quad clipping
typedef int QUAD_EDGE_LIST;
struct QUAD_CASES_t
{
  QUAD_EDGE_LIST edges[14];
};
using QUAD_CASES = struct QUAD_CASES_t;

static QUAD_CASES quadCases[] = {
  { { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },    // 0
  { { 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 1
  { { 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 2
  { { 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 3
  { { 3, 102, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 4
  { { 3, 100, 0, 3, 3, 102, 2, 1, 4, 0, 1, 2, 3, -1 } },             // 5
  { { 4, 101, 102, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 6
  { { 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 102, 2, -1, -1 } },      // 7
  { { 3, 103, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 8
  { { 4, 100, 0, 2, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 9
  { { 3, 101, 1, 0, 3, 103, 3, 2, 4, 0, 1, 2, 3, -1 } },             // 10
  { { 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 103, -1, -1 } },      // 11
  { { 4, 102, 103, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 12
  { { 3, 100, 0, 103, 3, 0, 1, 103, 3, 1, 102, 103, -1, -1 } },      // 13
  { { 3, 0, 101, 102, 3, 0, 102, 3, 3, 102, 103, 3, -1, -1 } },      // 14
  { { 4, 100, 101, 102, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 15
};

static QUAD_CASES quadCasesComplement[] = {
  { { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },    // 0
  { { 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 1
  { { 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 2
  { { 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 3
  { { 3, 102, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 4
  { { 3, 100, 0, 3, 3, 102, 2, 1, -1, -1, -1, -1, -1, -1 } },        // 5
  { { 4, 101, 102, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 6
  { { 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 102, 2, -1, -1 } },      // 7
  { { 3, 103, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },      // 8
  { { 4, 100, 0, 2, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 9
  { { 3, 101, 1, 0, 3, 103, 3, 2, -1, -1, -1, -1, -1, -1 } },        // 10
  { { 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 103, -1, -1 } },      // 11
  { { 4, 102, 103, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 } },     // 12
  { { 3, 100, 0, 103, 3, 0, 1, 103, 3, 1, 102, 103, -1, -1 } },      // 13
  { { 3, 0, 101, 102, 3, 0, 102, 3, 3, 102, 103, 3, -1, -1 } },      // 14
  { { 4, 100, 101, 102, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1 } }, // 15
};

//------------------------------------------------------------------------------
// Clip this quad using scalar value provided. Like contouring, except
// that it cuts the quad to produce other quads and/or triangles.
void vtkQuad::Clip(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
  vtkCellArray* polys, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId,
  vtkCellData* outCd, int insideOut)
{
  static const int CASE_MASK[4] = { 1, 2, 4, 8 };
  QUAD_CASES* quadCase;
  QUAD_EDGE_LIST* edge;
  int i, j, index;
  const vtkIdType* vert;
  int e1, e2;
  int newCellId;
  vtkIdType pts[4];
  int vertexId;
  double t, x1[3], x2[3], x[3], deltaScalar;
  double scalar0, scalar1, e1Scalar;

  // Build the index into the case table
  if (insideOut)
  {
    for (i = 0, index = 0; i < 4; i++)
    {
      if (cellScalars->GetComponent(i, 0) <= value)
      {
        index |= CASE_MASK[i];
      }
    }
    // Select case based on the index and get the list of edges for this case
    quadCase = quadCases + index;
  }
  else
  {
    for (i = 0, index = 0; i < 4; i++)
    {
      if (cellScalars->GetComponent(i, 0) > value)
      {
        index |= CASE_MASK[i];
      }
    }
    // Select case based on the index and get the list of edges for this case
    quadCase = quadCasesComplement + index;
  }

  edge = quadCase->edges;

  // generate each quad
  for (; edge[0] > -1; edge += edge[0] + 1)
  {
    for (i = 0; i < edge[0]; i++) // insert quad or triangle
    {
      // vertex exists, and need not be interpolated
      if (edge[i + 1] >= 100)
      {
        vertexId = edge[i + 1] - 100;
        this->Points->GetPoint(vertexId, x);
        if (locator->InsertUniquePoint(x, pts[i]))
        {
          outPd->CopyData(inPd, this->PointIds->GetId(vertexId), pts[i]);
        }
      }

      else // new vertex, interpolate
      {
        vert = edges[edge[i + 1]];

        // calculate a preferred interpolation direction
        scalar0 = cellScalars->GetComponent(vert[0], 0);
        scalar1 = cellScalars->GetComponent(vert[1], 0);
        deltaScalar = scalar1 - scalar0;

        if (deltaScalar > 0)
        {
          e1 = vert[0];
          e2 = vert[1];
          e1Scalar = scalar0;
        }
        else
        {
          e1 = vert[1];
          e2 = vert[0];
          e1Scalar = scalar1;
          deltaScalar = -deltaScalar;
        }

        // linear interpolation
        if (deltaScalar == 0.0)
        {
          t = 0.0;
        }
        else
        {
          t = (value - e1Scalar) / deltaScalar;
        }

        this->Points->GetPoint(e1, x1);
        this->Points->GetPoint(e2, x2);

        for (j = 0; j < 3; j++)
        {
          x[j] = x1[j] + t * (x2[j] - x1[j]);
        }

        if (locator->InsertUniquePoint(x, pts[i]))
        {
          vtkIdType p1 = this->PointIds->GetId(e1);
          vtkIdType p2 = this->PointIds->GetId(e2);
          outPd->InterpolateEdge(inPd, pts[i], p1, p2, t);
        }
      }
    }
    // check for degenerate output
    if (edge[0] == 3) // i.e., a triangle
    {
      if (pts[0] == pts[1] || pts[0] == pts[2] || pts[1] == pts[2])
      {
        continue;
      }
    }
    else // a quad
    {
      if ((pts[0] == pts[3] && pts[1] == pts[2]) || (pts[0] == pts[1] && pts[3] == pts[2]))
      {
        continue;
      }
    }

    newCellId = polys->InsertNextCell(edge[0], pts);
    outCd->CopyData(inCd, cellId, newCellId);
  }
}

//------------------------------------------------------------------------------
static double vtkQuadCellPCoords[12] = {
  0.0, 0.0, 0.0, //
  1.0, 0.0, 0.0, //
  1.0, 1.0, 0.0, //
  0.0, 1.0, 0.0  //
};
double* vtkQuad::GetParametricCoords()
{
  return vtkQuadCellPCoords;
}

//------------------------------------------------------------------------------
void vtkQuad::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "Line:\n";
  this->Line->PrintSelf(os, indent.GetNextIndent());
  os << indent << "Triangle:\n";
  this->Triangle->PrintSelf(os, indent.GetNextIndent());
}