// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-FileCopyrightText: Copyright 2003 Sandia Corporation
// SPDX-License-Identifier: LicenseRef-BSD-3-Clause-Sandia-USGov

/* Do not edit this file! It was generated by hand. Mostly.
 * Edit vtkStreamingTessellatorGenerator.py instead.
 */

#include "vtkObjectFactory.h"

#include "vtkEdgeSubdivisionCriterion.h"
#include "vtkStreamingTessellator.h"


#undef UGLY_ASPECT_RATIO_HACK
#undef DBG_MIDPTS

#include <algorithm>
#include <stack>

#ifdef PARAVIEW_DEBUG_TESSELLATOR
#  define VTK_TESSELLATOR_INCR_CASE_COUNT(cs) this->CaseCounts[cs]++
#  define VTK_TESSELLATOR_INCR_SUBCASE_COUNT(cs,sc) this->SubcaseCounts[cs][sc]++
#else // PARAVIEW_DEBUG_TESSELLATOR
#  define VTK_TESSELLATOR_INCR_CASE_COUNT(cs)
#  define VTK_TESSELLATOR_INCR_SUBCASE_COUNT(cs,sc)
#endif // PARAVIEW_DEBUG_TESSELLATOR

VTK_ABI_NAMESPACE_BEGIN
static void DefaultFacet3Callback(
  const double* a, const double* b, const double* c, const double* d,
  vtkEdgeSubdivisionCriterion*, void* pd, const void* )
{
  (void)a;
  (void)b;
  (void)c;
  (void)d;
  (void)pd;
}

static void DefaultFacet2Callback(
  const double* a, const double* b, const double* c,
  vtkEdgeSubdivisionCriterion*, void* pd, const void* )
{
  (void)a;
  (void)b;
  (void)c;
  (void)pd;
}

static void DefaultFacet1Callback(
  const double* a, const double* b,
  vtkEdgeSubdivisionCriterion*, void* pd, const void* )
{
  (void)a;
  (void)b;
  (void)pd;
}

static void DefaultFacet0Callback(
  const double* a,
  vtkEdgeSubdivisionCriterion*, void* pd, const void* )
{
  (void)a;
  (void)pd;
}

vtkStandardNewMacro(vtkStreamingTessellator);

void vtkStreamingTessellator::PrintSelf( ostream& os, vtkIndent indent )
{
  this->Superclass::PrintSelf( os, indent );
  os << indent << "PointDimension:       "
     << this->PointDimension[1] << " " << this->PointDimension[2] << " " << this->PointDimension[3] << endl;
  os << indent << "EmbeddingDimension:   "
     << this->EmbeddingDimension[1] << " " << this->EmbeddingDimension[2] << " " << this->EmbeddingDimension[3] << endl;
  os << indent << "PrivateData:          " << this->PrivateData << endl;
  os << indent << "ConstPrivateData:     " << this->ConstPrivateData << endl;
  os << indent << "SubdivisionAlgorithm: " << this->Algorithm << endl;
  os << indent << "VertexCallback:       " << this->Callback0 << endl;
  os << indent << "EdgeCallback:         " << this->Callback1 << endl;
  os << indent << "TriangleCallback:     " << this->Callback2 << endl;
  os << indent << "TetrahedronCallback:  " << this->Callback3 << endl;
#ifdef PARAVIEW_DEBUG_TESSELLATOR
  os << indent << "CaseCounts:           " << this->CaseCounts << endl;
  os << indent << "SubcaseCounts:        " << this->SubcaseCounts << endl;
#endif // PARAVIEW_DEBUG_TESSELLATOR
}

vtkStreamingTessellator::vtkStreamingTessellator()
{
  this->PrivateData = nullptr;
  this->ConstPrivateData = nullptr;
  this->Algorithm = nullptr;
  this->Callback0 = DefaultFacet0Callback;
  this->Callback1 = DefaultFacet1Callback;
  this->Callback2 = DefaultFacet2Callback;
  this->Callback3 = DefaultFacet3Callback;
  this->MaximumNumberOfSubdivisions = 3;
  for ( int i=0; i<4; ++i )
  {
    this->EmbeddingDimension[i] = i;
    this->PointDimension[i] = i+3; // By default, FieldSize = 0
  }

}

vtkStreamingTessellator::~vtkStreamingTessellator()
{
  if ( this->Algorithm )
    this->Algorithm->UnRegister( this );
}

void vtkStreamingTessellator::SetEmbeddingDimension( int k, int d )
{
  if ( d > 8 )
  {
    vtkErrorMacro( "Embedding dimension may not be > 8. (You asked for " << d << "). Truncating to 8." );
    d = 8;
  }

  if ( k == 0 || k < -1 || k >= 4 )
  {
    vtkWarningMacro( "Invalid argument k=" << k );
    return;
  }

  if ( k < 0 )
  {
    for ( k=0; k<4; k++ )
      if ( this->EmbeddingDimension[k] != d )
      {
        this->PointDimension[k] += d - this->EmbeddingDimension[k] ;
        this->EmbeddingDimension[k] = d;
        this->Modified();
      }
    return;
  }
  if ( this->EmbeddingDimension[k] != d )
  {
    this->PointDimension[k] += d - this->EmbeddingDimension[k] ;
    this->EmbeddingDimension[k] = d;
    this->Modified();
  }
}

void vtkStreamingTessellator::SetFieldSize( int k, int s )
{
  if ( s > vtkStreamingTessellator::MaxFieldSize )
  {
    vtkErrorMacro( "Embedding dimension may not be > " << MaxFieldSize << ". (You asked for " << s << "). Truncating to " << MaxFieldSize );
    s = vtkStreamingTessellator::MaxFieldSize;
  }

  if ( k == 0 || k < -1 || k >= 4 )
  {
    vtkWarningMacro( "Invalid argument k=" << k );
    return;
  }

  if ( k < 0 )
  {
    // Use field size for all facet types (point, line, triangle, tet, ...)
    for ( k=0; k<4; k++ )
      if ( this->PointDimension[k] != s + this->EmbeddingDimension[k] + 3 )
      {
        this->PointDimension[k] = s + this->EmbeddingDimension[k] + 3;
        this->Modified();
      }
    return;
  }

  if ( this->PointDimension[k] != s + this->EmbeddingDimension[k] + 3 )
  {
    this->PointDimension[k] = s + this->EmbeddingDimension[k] + 3;
    this->Modified();
  }
}

void vtkStreamingTessellator::SetMaximumNumberOfSubdivisions( int num_subdiv_in )
{
  if ( this->MaximumNumberOfSubdivisions == num_subdiv_in )
    return;

  if ( num_subdiv_in < 0 )
  {
    vtkErrorMacro( "MaximumNumberOfSubdivisions must be 0 or greater (you requested " << num_subdiv_in << ")" );
    return;
  }

  this->MaximumNumberOfSubdivisions = num_subdiv_in;
  this->Modified();
}

void vtkStreamingTessellator::SetTetrahedronCallback( TetrahedronProcessorFunction f )
{
  this->Callback3 = f ? f : DefaultFacet3Callback;
}

vtkStreamingTessellator::TetrahedronProcessorFunction vtkStreamingTessellator::GetTetrahedronCallback() const
{
  return this->Callback3;
}

void vtkStreamingTessellator::SetTriangleCallback( TriangleProcessorFunction f )
{
  this->Callback2 = f ? f : DefaultFacet2Callback;
}

vtkStreamingTessellator::TriangleProcessorFunction vtkStreamingTessellator::GetTriangleCallback() const
{
  return this->Callback2;
}

void vtkStreamingTessellator::SetEdgeCallback( EdgeProcessorFunction f )
{
  this->Callback1 = f ? f : DefaultFacet1Callback;
}

vtkStreamingTessellator::EdgeProcessorFunction vtkStreamingTessellator::GetEdgeCallback() const
{
  return this->Callback1;
}

void vtkStreamingTessellator::SetVertexCallback( VertexProcessorFunction f )
{
  this->Callback0 = f ? f : DefaultFacet0Callback;
}

vtkStreamingTessellator::VertexProcessorFunction vtkStreamingTessellator::GetVertexCallback() const
{
  return this->Callback0;
}

void vtkStreamingTessellator::SetPrivateData( void* Private )
{
  this->PrivateData = Private;
}

void* vtkStreamingTessellator::GetPrivateData() const
{
  return this->PrivateData;
}

void vtkStreamingTessellator::SetConstPrivateData( const void* ConstPrivate )
{
  this->ConstPrivateData = ConstPrivate;
}

const void* vtkStreamingTessellator::GetConstPrivateData() const
{
  return this->ConstPrivateData;
}

void vtkStreamingTessellator::SetSubdivisionAlgorithm( vtkEdgeSubdivisionCriterion* a )
{
  if ( a != this->Algorithm )
  {
    if ( this->Algorithm )
      this->Algorithm->UnRegister( this );

    this->Algorithm = a;
    this->Modified();

    if ( this->Algorithm )
      this->Algorithm->Register( this );
  }
}

vtkEdgeSubdivisionCriterion* vtkStreamingTessellator::GetSubdivisionAlgorithm()
{
  return this->Algorithm;
}

const vtkEdgeSubdivisionCriterion* vtkStreamingTessellator::GetSubdivisionAlgorithm() const
{
  return this->Algorithm;
}


// Returns true if || a0a1 || < || b0b1 ||
// We use this to test which triangulation has the best
// aspect ratio when there are 2 to choose from.
static bool compareHopfCrossStringDist( const double* a0, const double* a1, const double* b0, const double* b1 )
{
  double SqMagA = 0.;
  double SqMagB = 0.;
  for (int i=0; i<3; i++)
  {
    double tmp;
    tmp = a0[i] - a1[i];
    SqMagA += tmp * tmp;
    tmp = b0[i] - b1[i];
    SqMagB += tmp * tmp;
  }
  return SqMagA < SqMagB;
}


int vtkStreamingTessellator::BestTets( int* vtkNotUsed(connOffsets), double** vtkNotUsed(verts), int vtkNotUsed(permOffset), int vtkNotUsed(sgn) ) const
{
  // Re-run vtkStreamingTessellatorGenerator.py with QualityThang=1
  // to get this implemented (along with on-the-fly quality improvement)
  return 1;
}

void vtkStreamingTessellator::AdaptivelySample3FacetLinear(
  double* v0, double* v1, double* v2, double* v3) const
{
  Callback3(v0, v1, v2, v3, this->Algorithm, this->PrivateData, this->ConstPrivateData);
}

void vtkStreamingTessellator::AdaptivelySample2FacetLinear(double* v0, double* v1, double* v2) const
{
  Callback2(v0, v1, v2, this->Algorithm, this->PrivateData, this->ConstPrivateData);
}

void vtkStreamingTessellator::AdaptivelySample1FacetLinear(double* v0, double* v1) const
{
  Callback1(v0, v1, this->Algorithm, this->PrivateData, this->ConstPrivateData);
}

void vtkStreamingTessellator::AdaptivelySample0Facet( double* v0 ) const
{
  Callback0( v0, this->Algorithm, this->PrivateData, this->ConstPrivateData );
}

void vtkStreamingTessellator::AdaptivelySample1Facet( double* v0, double* v1, int maxDepth ) const
{
  int edgeCode = 0;

  double midpt0[11+vtkStreamingTessellator::MaxFieldSize];
  // make valgrind happy
  std::fill(midpt0,midpt0+this->PointDimension[1],0.);

  if ( maxDepth-- > 0 )
  {
      for ( int i=0; i<this->PointDimension[1]; i++ )
        midpt0[i] = (v0[i] + v1[i])/2.;

      if ( this->Algorithm->EvaluateLocationAndFields(midpt0, 3+this->EmbeddingDimension[1] ) )
        edgeCode += 1;
  }

  switch (edgeCode) {
    // No edges to subdivide
    case 0:
      Callback1( v0, v1, this->Algorithm, this->PrivateData, this->ConstPrivateData );
      break ;

      // One edge to subdivide
    case 1:
      this->AdaptivelySample1Facet( v0, midpt0, maxDepth );
      this->AdaptivelySample1Facet( midpt0, v1, maxDepth );
      break;
  }
}

void vtkStreamingTessellator::AdaptivelySample2Facet( double* v0, double* v1, double* v2, int maxDepth, int move ) const
{
  int edgeCode = 0;

  double midpt0[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt1[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt2[11+vtkStreamingTessellator::MaxFieldSize];

  // Make valgrind happy
  std::fill(midpt0,midpt0+this->PointDimension[2],0.);
  std::fill(midpt1,midpt1+this->PointDimension[2],0.);
  std::fill(midpt2,midpt2+this->PointDimension[2],0.);

  if ( maxDepth-- > 0 )
  {

    for ( int i=0; i<this->PointDimension[2]; i++ )
    {
      midpt0[i] = (v0[i] + v1[i])/2.;
      midpt1[i] = (v1[i] + v2[i])/2.;
      midpt2[i] = (v2[i] + v0[i])/2.;
    }

    if ( (move & 1) && Algorithm->EvaluateLocationAndFields( midpt0, 3+this->EmbeddingDimension[2] ) )
      edgeCode += 1;
    if ( (move & 2) && Algorithm->EvaluateLocationAndFields( midpt1, 3+this->EmbeddingDimension[2] ) )
      edgeCode += 2;
    if ( (move & 4) && Algorithm->EvaluateLocationAndFields( midpt2, 3+this->EmbeddingDimension[2] ) )
      edgeCode += 4;
#ifdef UGLY_ASPECT_RATIO_HACK
    double dist0=0.;
    double dist1=0.;
    double dist2=0.;
    double tmp;
    for ( int j=0; j<3; ++j )
    {
      tmp = v0[j] - v1[j];
      dist0 += tmp*tmp;
      tmp = v1[j] - v2[j];
      dist1 += tmp*tmp;
      tmp = v2[j] - v0[j];
      dist2 += tmp*tmp;
    }

    if ( edgeCode & 1 ) dist0 /= 2.;
    if ( edgeCode & 2 ) dist1 /= 2.;
    if ( edgeCode & 4 ) dist2 /= 2.;

#define MAR2 2.25
    if ( (!(edgeCode & 1)) && (move&1) && ((dist0/dist1 > MAR2) || (dist0/dist2 > MAR2)) )
    {
      edgeCode += 1;
      move &= 6;
    }
    if ( (!(edgeCode & 2)) && (move&2) && ((dist1/dist0 > MAR2) || (dist1/dist2 > MAR2)) )
    {
      edgeCode += 2;
      move &= 5;
    }
    if ( (!(edgeCode & 4)) && (move&4) && ((dist2/dist1 > MAR2) || (dist2/dist0 > MAR2)) )
    {
      edgeCode += 4;
      move &= 3;
    }
#endif // UGLY_ASPECT_RATIO_HACK
  }

#ifdef DBG_MIDPTS
  if ( maxDepth == 0 )
  {
    fprintf( stderr, "midpoint of v%d (%g %g %g/%g %g %g)-v%d (%g %g %g/%g %g %g) = (%g %g %g/%g %g %g)\n",
      0, v0[0], v0[1], v0[2], v0[3], v0[4], v0[5],
      1, v1[0], v1[1], v1[2], v1[3], v1[4], v1[5],
         midpt0[0], midpt0[1], midpt0[2], midpt0[3], midpt0[4], midpt0[5]
    );

    fprintf( stderr, "midpoint of v%d (%g %g %g/%g %g %g)-v%d (%g %g %g/%g %g %g) = (%g %g %g/%g %g %g)\n",
      1, v1[0], v1[1], v1[2], v1[3], v1[4], v1[5],
      2, v2[0], v2[1], v2[2], v2[3], v2[4], v2[5],
         midpt1[0], midpt1[1], midpt1[2], midpt1[3], midpt1[4], midpt1[5]
    );

    fprintf( stderr, "midpoint of v%d (%g %g %g/%g %g %g)-v%d (%g %g %g/%g %g %g) = (%g %g %g/%g %g %g)\n\n",
      2, v2[0], v2[1], v2[2], v2[3], v2[4], v2[5],
      0, v0[0], v0[1], v0[2], v0[3], v0[4], v0[5],
         midpt2[0], midpt2[1], midpt2[2], midpt2[3], midpt2[4], midpt2[5]
    );
  }
#endif // DBG_MIDPTS

  switch (edgeCode)
  {
    // No edges to subdivide
  case 0:
    Callback2( v0, v1, v2, this->Algorithm, this->PrivateData, this->ConstPrivateData );
    break ;

    // One edge to subdivide
  case 1:
    this->AdaptivelySample2Facet( v0, midpt0, v2, maxDepth, move | 2 );
    this->AdaptivelySample2Facet( midpt0, v1, v2, maxDepth, move | 4 );
    break;
  case 2:
    this->AdaptivelySample2Facet( v0, v1, midpt1, maxDepth, move | 4 );
    this->AdaptivelySample2Facet( v0, midpt1, v2, maxDepth, move | 1 );
    break;
  case 4:
    this->AdaptivelySample2Facet( v0, v1, midpt2, maxDepth, move | 2 );
    this->AdaptivelySample2Facet( midpt2, v1, v2, maxDepth, move | 1 );
    break;

    // Two edges to subdivide
  case 3:
    this->AdaptivelySample2Facet( midpt0, v1, midpt1, maxDepth, move | 4 );
    if ( compareHopfCrossStringDist( v2, midpt0, v0, midpt1 ) )
    {
      this->AdaptivelySample2Facet( midpt0, midpt1,   v2  , maxDepth, move | 5 );
      this->AdaptivelySample2Facet(   v0,   midpt0,   v2  , maxDepth, move | 2 );
    }
    else
    {
      this->AdaptivelySample2Facet(   v0  , midpt0, midpt1, maxDepth, move | 6 );
      this->AdaptivelySample2Facet(   v0,   midpt1,   v2  , maxDepth, move | 1 );
    }
    break;
  case 5:
    this->AdaptivelySample2Facet( v0, midpt0, midpt2, maxDepth, move | 2 );
    if ( compareHopfCrossStringDist( v2, midpt0, v1, midpt2 ) )
    {
      this->AdaptivelySample2Facet( midpt0,   v1,     v2  , maxDepth, move | 4 );
      this->AdaptivelySample2Facet( midpt2, midpt0,   v2  , maxDepth, move | 3 );
    }
    else
    {
      this->AdaptivelySample2Facet( midpt0,   v1,   midpt2, maxDepth, move | 6 );
      this->AdaptivelySample2Facet( midpt2,   v1,     v2,   maxDepth, move | 1 );
    }
    break;
  case 6:
    this->AdaptivelySample2Facet( midpt2, midpt1, v2, maxDepth, move | 1 );
    if ( compareHopfCrossStringDist( v0, midpt1, v1, midpt2 ) )
    {
      this->AdaptivelySample2Facet(   v0,   midpt1, midpt2, maxDepth, move | 3 );
      this->AdaptivelySample2Facet(   v0,     v1,   midpt1, maxDepth, move | 4 );
    }
    else
    {
      this->AdaptivelySample2Facet(   v0,     v1,   midpt2, maxDepth, move | 2 );
      this->AdaptivelySample2Facet( midpt2,   v1,   midpt1, maxDepth, move | 5 );
    }
    break;

    // Three edges to subdivide
  case 7:
    this->AdaptivelySample2Facet( midpt0, midpt1, midpt2, maxDepth, 7 );
    this->AdaptivelySample2Facet(   v0  , midpt0, midpt2, maxDepth, move | 2 );
    this->AdaptivelySample2Facet( midpt0,   v1  , midpt1, maxDepth, move | 4 );
    this->AdaptivelySample2Facet( midpt2, midpt1,   v2  , maxDepth, move | 1 );
    break;
  }
}

void vtkStreamingTessellator::AdaptivelySample3Facet( double* v0, double* v1, double* v2, double* v3, int maxDepth ) const
{
  int edgeCode = 0;

  double midpt0[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt1[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt2[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt3[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt4[11+vtkStreamingTessellator::MaxFieldSize];
  double midpt5[11+vtkStreamingTessellator::MaxFieldSize];
  double facept0[11+vtkStreamingTessellator::MaxFieldSize];
  double facept1[11+vtkStreamingTessellator::MaxFieldSize];
  double facept2[11+vtkStreamingTessellator::MaxFieldSize];
  double facept3[11+vtkStreamingTessellator::MaxFieldSize];

  // Make valgrind happy
  std::fill(midpt0,midpt0+this->PointDimension[3],0.);
  std::fill(midpt1,midpt1+this->PointDimension[3],0.);
  std::fill(midpt2,midpt2+this->PointDimension[3],0.);
  std::fill(midpt3,midpt3+this->PointDimension[3],0.);
  std::fill(midpt4,midpt4+this->PointDimension[3],0.);
  std::fill(midpt5,midpt5+this->PointDimension[3],0.);

  double edgeLength2[6];
  if ( maxDepth-- > 0 )
  {
    for ( int i=0; i<this->PointDimension[3]; i++ )
    {
      midpt0[i] = (v0[i] + v1[i])/2.;
      midpt1[i] = (v1[i] + v2[i])/2.;
      midpt2[i] = (v2[i] + v0[i])/2.;
      midpt3[i] = (v0[i] + v3[i])/2.;
      midpt4[i] = (v1[i] + v3[i])/2.;
      midpt5[i] = (v2[i] + v3[i])/2.;
    }

    if ( Algorithm->EvaluateLocationAndFields(midpt0, 3+this->EmbeddingDimension[3] ) )
      edgeCode |=  1;
    if ( Algorithm->EvaluateLocationAndFields(midpt1, 3+this->EmbeddingDimension[3] ) )
      edgeCode |=  2;
    if ( Algorithm->EvaluateLocationAndFields(midpt2, 3+this->EmbeddingDimension[3] ) )
      edgeCode |=  4;

    if ( Algorithm->EvaluateLocationAndFields(midpt3, 3+this->EmbeddingDimension[3] ) )
      edgeCode |=  8;
    if ( Algorithm->EvaluateLocationAndFields(midpt4, 3+this->EmbeddingDimension[3] ) )
      edgeCode |= 16;
    if ( Algorithm->EvaluateLocationAndFields(midpt5, 3+this->EmbeddingDimension[3] ) )
      edgeCode |= 32;

    edgeLength2[0] = edgeLength2[1] = edgeLength2[2] = edgeLength2[3]
      = edgeLength2[4] = edgeLength2[5] = 0;
    for ( int c=0; c<3; ++c )
    {
      double tmp;
      tmp = v1[c] - v0[c];
      edgeLength2[0] += tmp*tmp;
      tmp = v2[c] - v1[c];
      edgeLength2[1] += tmp*tmp;
      tmp = v2[c] - v0[c];
      edgeLength2[2] += tmp*tmp;
      tmp = v3[c] - v0[c];
      edgeLength2[3] += tmp*tmp;
      tmp = v3[c] - v1[c];
      edgeLength2[4] += tmp*tmp;
      tmp = v3[c] - v2[c];
      edgeLength2[5] += tmp*tmp;
    }
  }

  if ( edgeCode == 0 )
  {
    // No edges to subdivide
    Callback3( v0, v1, v2, v3, this->Algorithm, this->PrivateData, this->ConstPrivateData );
  }
  else
  {
    // Do the subdivision
    double* vertices[10] =
    {
      v0, v1, v2, v3,     midpt0, midpt1, midpt2, midpt3, midpt4, midpt5
    };

    // Generate tetrahedra that are compatible except when edge
    // lengths are equal on indeterminately subdivided faces.
    double* permuted[14];
    double permlen[6]; // permuted edge lengths
    int C = vtkStreamingTessellator::EdgeCodesToCaseCodesPlusPermutation[ edgeCode ][0];
    int P = vtkStreamingTessellator::EdgeCodesToCaseCodesPlusPermutation[ edgeCode ][1];
    int i;

    // 1. Permute the tetrahedron into our canonical configuration
    for ( i=0; i<4; ++i )
    {
      permuted[i] = vertices[ vtkStreamingTessellator::PermutationsFromIndex[P][i] ];
    }
    for ( i=4; i<10; ++i )
    {
      // permute the edge lengths, too
      permuted[i] = vertices[ vtkStreamingTessellator::PermutationsFromIndex[P][i] ];
      permlen[i-4]  = edgeLength2[ vtkStreamingTessellator::PermutationsFromIndex[P][i] - 4 ];
    }
    // Add our local (heap) storage for face points to the list.
    permuted[10] = facept0;
    permuted[11] = facept1;
    permuted[12] = facept2;
    permuted[13] = facept3;

    int comparisonBits;
    std::stack<vtkIdType*> outputTets;
    std::stack<vtkIdType*> outputPerm;
    std::stack<int>        outputSign;

    // cout << "Case " << C << "  Permutation " << P << endl;
    // 2. Generate tetrahedra based on the configuration.
    //    Note that case 0 is handled above (edgeCode == 0).
    switch (C)
    {

    case 1: // Ruprecht-Müller Case 1
      VTK_TESSELLATOR_INCR_CASE_COUNT(0);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 0 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(0,0);
      break;
    case 2: // Ruprecht-Müller Case 2a
      comparisonBits =
        (permlen[0] <= permlen[1] ? 1 : 0) | (permlen[0] >= permlen[1] ? 2 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[10][i] = (permuted[0][i] + permuted[2][i])*0.375 + permuted[1][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[10], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(1);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 9 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(1,0);
      switch (comparisonBits)
      {
        case 2: // 0>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 14 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(1,1);
          break;
        case 1: // 1>0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 23 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(1,2);
          break;
        case 3: // 0=1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 32 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(1,3);
          break;
      }
      break;
    case 3: // Ruprecht-Müller Case 2b
      VTK_TESSELLATOR_INCR_CASE_COUNT(2);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 49 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(2,0);
      break;
    case 4: // Ruprecht-Müller Case 3a
      comparisonBits =
        (permlen[0] <= permlen[2] ? 1 : 0) | (permlen[0] >= permlen[2] ? 2 : 0) |
        (permlen[0] <= permlen[3] ? 4 : 0) | (permlen[0] >= permlen[3] ? 8 : 0) |
        (permlen[2] <= permlen[3] ? 16 : 0) | (permlen[2] >= permlen[3] ? 32 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[10][i] = (permuted[1][i] + permuted[2][i])*0.375 + permuted[0][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[10], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 12) == 12 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[11][i] = (permuted[1][i] + permuted[3][i])*0.375 + permuted[0][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[11], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 48) == 48 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[13][i] = (permuted[2][i] + permuted[3][i])*0.375 + permuted[0][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[13], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(3);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 66 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,0);
      switch (comparisonBits)
      {
        case 42: // 0>2>3<0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 71 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,1);
          break;
        case 37: // 2>3>0<2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 71 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[11] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,2);
          break;
        case 22: // 3>0>2<3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 71 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[3] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,3);
          break;
        case 21: // 3>2>0<3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 71 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[22] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,4);
          break;
        case 41: // 2>0>3<2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 71 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[12] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,5);
          break;
        case 26: // 0>3>2<0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 71 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,6);
          break;
        case 43: // 0=2>3<0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 84 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,7);
          break;
        case 53: // 2=3>0<2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 84 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[11] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,8);
          break;
        case 30: // 0=3>2<0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 84 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[3] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,9);
          break;
        case 23: // 3>0=2<3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 105 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,10);
          break;
        case 58: // 0>2=3<0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 105 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[11] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,11);
          break;
        case 45: // 2>0=3<2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 105 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[3] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,12);
          break;
        case 63: // 0=2=3=0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 126 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(3,13);
          break;
      }
      break;
    case 5: // Ruprecht-Müller Case 3b
      VTK_TESSELLATOR_INCR_CASE_COUNT(4);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 171 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(4,0);
      break;
    case 6: // Ruprecht-Müller Case 3c
      comparisonBits =
        (permlen[0] <= permlen[3] ? 1 : 0) | (permlen[0] >= permlen[3] ? 2 : 0) |
        (permlen[0] <= permlen[1] ? 4 : 0) | (permlen[0] >= permlen[1] ? 8 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[11][i] = (permuted[1][i] + permuted[3][i])*0.375 + permuted[0][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[11], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 12) == 12 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[10][i] = (permuted[0][i] + permuted[2][i])*0.375 + permuted[1][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[10], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(5);
      switch (comparisonBits)
      {
        case 10: // 0>1,0>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 188 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,0);
          break;
        case 5: // 1>0,3>0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 209 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,1);
          break;
        case 9: // 0>1,3>0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 230 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,2);
          break;
        case 6: // 1>0,0>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 251 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,3);
          break;
        case 14: // 0=1,0>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 272 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,4);
          break;
        case 11: // 0=3,0>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 272 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,5);
          break;
        case 13: // 3>0,0=1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 301 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,6);
          break;
        case 7: // 1>0,0=3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 301 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,7);
          break;
        case 15: // 0=1,0=3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 330 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(5,8);
          break;
      }
      break;
    case 7: // Ruprecht-Müller Case 3d
      comparisonBits =
        (permlen[0] <= permlen[2] ? 1 : 0) | (permlen[0] >= permlen[2] ? 2 : 0) |
        (permlen[0] <= permlen[4] ? 4 : 0) | (permlen[0] >= permlen[4] ? 8 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[10][i] = (permuted[1][i] + permuted[2][i])*0.375 + permuted[0][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[10], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 12) == 12 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[11][i] = (permuted[0][i] + permuted[3][i])*0.375 + permuted[1][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[11], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(6);
      switch (comparisonBits)
      {
        case 10: // 0>4,0>2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 371 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,0);
          break;
        case 5: // 4>0,2>0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 392 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,1);
          break;
        case 9: // 0>4,2>0
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 413 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,2);
          break;
        case 6: // 4>0,0>2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 434 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,3);
          break;
        case 14: // 0=4,0>2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 455 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,4);
          break;
        case 11: // 0=2,0>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 455 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,5);
          break;
        case 13: // 2>0,0=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 484 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,6);
          break;
        case 7: // 4>0,0=2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 484 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,7);
          break;
        case 15: // 0=4,0=2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 513 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(6,8);
          break;
      }
      break;
    case 8: // Ruprecht-Müller Case 4a
      comparisonBits =
        (permlen[4] <= permlen[5] ? 1 : 0) | (permlen[4] >= permlen[5] ? 2 : 0) |
        (permlen[3] <= permlen[4] ? 4 : 0) | (permlen[3] >= permlen[4] ? 8 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[12][i] = (permuted[1][i] + permuted[2][i])*0.375 + permuted[3][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[12], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 12) == 12 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[11][i] = (permuted[0][i] + permuted[1][i])*0.375 + permuted[3][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[11], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(7);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 554 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,0);
      switch (comparisonBits)
      {
        case 5: // 5>4>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 563 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,1);
          break;
        case 10: // 3>4>5
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 563 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[13] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,2);
          break;
        case 6: // 3<4>5
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 580 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,3);
          break;
        case 9: // 3>4<5
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 597 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,4);
          break;
        case 14: // 3=4>5
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 614 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,5);
          break;
        case 7: // 4=5,4>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 614 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[13] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,6);
          break;
        case 13: // 5>4,3=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 639 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,7);
          break;
        case 11: // 3>4=5
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 639 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[13] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,8);
          break;
        case 15: // 3=4=5
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 664 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(7,9);
          break;
      }
      break;
    case 9: // Ruprecht-Müller Case 4b
      comparisonBits =
        (permlen[2] <= permlen[3] ? 1 : 0) | (permlen[2] >= permlen[3] ? 2 : 0) |
        (permlen[1] <= permlen[2] ? 4 : 0) | (permlen[1] >= permlen[2] ? 8 : 0) |
        (permlen[1] <= permlen[4] ? 16 : 0) | (permlen[1] >= permlen[4] ? 32 : 0) |
        (permlen[3] <= permlen[4] ? 64 : 0) | (permlen[3] >= permlen[4] ? 128 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[13][i] = (permuted[2][i] + permuted[3][i])*0.375 + permuted[0][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[13], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 12) == 12 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[10][i] = (permuted[1][i] + permuted[0][i])*0.375 + permuted[2][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[10], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 48) == 48 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[12][i] = (permuted[2][i] + permuted[3][i])*0.375 + permuted[1][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[12], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 192) == 192 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[11][i] = (permuted[0][i] + permuted[1][i])*0.375 + permuted[3][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[11], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(8);
      switch (comparisonBits)
      {
        case 85: // 2>1,3>2,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,0);
          break;
        case 153: // 1>2,3>2,3>4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,1);
          break;
        case 170: // 1>2,2>3,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,2);
          break;
        case 102: // 2>1,2>3,4>3,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,3);
          break;
        case 90: // 1>2,2>3,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[9] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,4);
          break;
        case 150: // 2>1,2>3,3>4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[7] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,5);
          break;
        case 165: // 2>1,3>2,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[19] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,6);
          break;
        case 105: // 1>2,3>2,4>3,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 697 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[23] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,7);
          break;
        case 149: // 2>1,3>2,3>4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 722 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,8);
          break;
        case 89: // 1>2,3>2,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 722 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,9);
          break;
        case 106: // 1>2,2>3,4>3,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 722 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,10);
          break;
        case 166: // 2>1,2>3,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 722 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,11);
          break;
        case 86: // 2>1,2>3,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 747 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,12);
          break;
        case 169: // 1>2,3>2,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 747 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,13);
          break;
        case 157: // 1=2,3>2,3>4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,14);
          break;
        case 93: // 1=2,3>2,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,15);
          break;
        case 230: // 2>1,2>3,3=4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,16);
          break;
        case 234: // 1>2,2>3,3=4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,17);
          break;
        case 107: // 1>2,2=3,4>3,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[9] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,18);
          break;
        case 181: // 2>1,3>2,3>4,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[7] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,19);
          break;
        case 182: // 2>1,2>3,3>4,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[19] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,20);
          break;
        case 91: // 1>2,2=3,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 772 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[23] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,21);
          break;
        case 110: // 1=2,2>3,4>3,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,22);
          break;
        case 174: // 1=2,2>3,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,23);
          break;
        case 217: // 1>2,3>2,3=4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,24);
          break;
        case 213: // 2>1,3>2,3=4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,25);
          break;
        case 151: // 2>1,2=3,3>4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[9] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,26);
          break;
        case 122: // 1>2,2>3,4>3,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[7] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,27);
          break;
        case 121: // 1>2,3>2,4>3,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[19] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,28);
          break;
        case 167: // 2>1,2=3,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 809 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[23] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,29);
          break;
        case 94: // 1=2,2>3,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,30);
          break;
        case 173: // 1=2,3>2,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,31);
          break;
        case 233: // 1>2,3>2,3=4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,32);
          break;
        case 214: // 2>1,2>3,3=4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,33);
          break;
        case 87: // 2>1,2=3,4>3,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[9] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,34);
          break;
        case 118: // 2>1,2>3,4>3,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[7] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,35);
          break;
        case 185: // 1>2,3>2,3>4,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[19] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,36);
          break;
        case 171: // 1>2,2=3,3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 846 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[23] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,37);
          break;
        case 221: // 1=2,3>2,3=4,4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 883 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,38);
          break;
        case 238: // 1=2>3=4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 883 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,39);
          break;
        case 123: // 1>2,2=3,4>3,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 883 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[9] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,40);
          break;
        case 183: // 2>1,2=3,3>4,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 883 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[7] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,41);
          break;
        case 95: // 4>1=2=3,4>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 928 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,42);
          break;
        case 189: // 1=2,3>2,3>4,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 928 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,43);
          break;
        case 235: // 1>2=3=4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 928 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,44);
          break;
        case 246: // 2>1,2>3,3=4,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 928 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,45);
          break;
        case 175: // 1=2=3>4,1>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 977 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,46);
          break;
        case 126: // 1=2>3,4>3,1=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 977 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,47);
          break;
        case 215: // 2>1,2=3=4>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 977 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,48);
          break;
        case 249: // 1>2,3>2,3=4=1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 977 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,49);
          break;
        case 255: // 1=2=3=4=1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1026 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(8,50);
          break;
      }
      break;
    case 10: // Ruprecht-Müller Case 5
      comparisonBits =
        (permlen[1] <= permlen[2] ? 1 : 0) | (permlen[1] >= permlen[2] ? 2 : 0) |
        (permlen[3] <= permlen[4] ? 4 : 0) | (permlen[3] >= permlen[4] ? 8 : 0) |
        0;
      if ( (comparisonBits & 3) == 3 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[10][i] = (permuted[1][i] + permuted[0][i])*0.375 + permuted[2][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[10], 3+this->EmbeddingDimension[3] );
      }
      if ( (comparisonBits & 12) == 12 )
      {
        // Compute face point
        for ( i=0; i<this->PointDimension[3]; i++ )
        {
          permuted[11][i] = (permuted[0][i] + permuted[1][i])*0.375 + permuted[3][i]/4.;
        }
        Algorithm->EvaluateLocationAndFields(permuted[11], 3+this->EmbeddingDimension[3] );
      }
      VTK_TESSELLATOR_INCR_CASE_COUNT(9);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1091 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,0);
      switch (comparisonBits)
      {
        case 10: // 1>2,3>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1100 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,1);
          break;
        case 5: // 2>1,4>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1100 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,2);
          break;
        case 6: // 1>2,4>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1121 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,3);
          break;
        case 9: // 2>1,3>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1121 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,4);
          break;
        case 11: // 1=2,3>4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1142 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,5);
          break;
        case 7: // 1=2,4>3
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1142 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[14] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,6);
          break;
        case 14: // 3=4,1>2
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1142 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[5] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,7);
          break;
        case 13: // 3=4,2>1
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1142 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[15] );
          outputSign.push( -1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,8);
          break;
        case 15: // 1=2,3=4
          outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1175 );
          outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
          outputSign.push( 1 );
          VTK_TESSELLATOR_INCR_SUBCASE_COUNT(9,9);
          break;
      }
      break;
    case 11: // Ruprecht-Müller Case 6
      VTK_TESSELLATOR_INCR_CASE_COUNT(10);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1216 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(10,0);
      outputTets.push( vtkStreamingTessellator::TetrahedralDecompositions + 1233 );
      outputPerm.push( vtkStreamingTessellator::PermutationsFromIndex[0] );
      outputSign.push( 1 );
      VTK_TESSELLATOR_INCR_SUBCASE_COUNT(10,1);
      break;
    }

    vtkIdType* tets;
    vtkIdType  ntets;
    vtkIdType* perm;
    int        sgn;
#ifdef PARAVIEW_DEBUG_TESSELLATOR
    if ( outputTets.empty() )
    {
    cout << "Argh! Case " << C << " Perm " << P << " has no output!" << endl;
    }
#endif // PARAVIEW_DEBUG_TESSELLATOR
    while ( ! outputTets.empty() )
    {
      tets = outputTets.top();
      ntets = *tets;
      tets++;
      perm = outputPerm.top();
      sgn = outputSign.top();

      outputTets.pop();
      outputPerm.pop();
      outputSign.pop();

      int t;
      if ( sgn > 0 )
      {
        for ( t = 0; t < ntets; ++t )
        {
          this->AdaptivelySample3Facet(
            permuted[ perm[ tets[0] ] ], permuted[ perm[ tets[1] ] ],
            permuted[ perm[ tets[2] ] ], permuted[ perm[ tets[3] ] ],
            maxDepth );
          tets += 4;
        }
      }
      else
      {
        // we have an inverted tet... reverse the first 2 vertices
        // so the orientation is positive.
        for ( t = 0; t < ntets; ++t )
        {
          this->AdaptivelySample3Facet(
            permuted[ perm[ tets[1] ] ], permuted[ perm[ tets[0] ] ],
            permuted[ perm[ tets[2] ] ], permuted[ perm[ tets[3] ] ],
            maxDepth );
          tets += 4;
        }
      }
    }
  }
}


/*
 * The array below is indexed by the edge code for a tetrahedron.
 * Looking up a row with a tet's edge code will return C and P.
 * C is a configuration number and P is a permutation index.
 *
 * C is based on the case number from Ruprecht and
 * Müller's (1998) paper on adaptive tetrahedra. (The case
 * numbers are shown to the left of the row in the column
 * labeled case. The only difference is that we introduce
 * a case 3d which is part of case 3c in the paper.)
 *
 * P is an index into the vtkStreamingTessellator::PermutationsFromIndex array below,
 * and is used to transform the current tetrahedron into
 * the canonical configuration associated with C.
 *
 * The 6-digit binary number to the left (which is shown in
 * the horribly UNconventional LSB->MSB order) is the edge
 * code for the row. The 6 digits correspond to the 6 edges
 * of the tetrahedron; a '0' implies no subdivision while
 * a '1' implies subdivision should occur. The ordering of
 * the bits is
 *
 * Edge 0-1, Edge 1-2, Edge 2-0, Edge 0-3, Edge 1-3, Edge 2-3,
 *
 * where the numbers are vertices of the tetrahedron 0-1-2-3.
 * Note that Tet 0-1-2-3 must be positive (i.e., the plane
 * specified by Triangle 0-1-2 must have a normal pointing
 * towards vertex 3, and Triangle 0-1-2's normal must be
 * calculated using the cross-product (Edge 0-1) x (Edge 0-2)).
 *
 * ===========
 * References:
 * (Ruprect and Müller, 1998) A Scheme for Edge-based Adaptive
 *   Tetrahedron Subdivision, Mathematical Visualization (eds.
 *   Hege and Polthier), pp. 61--70. Springer-Verlag. 1998.
 */
int vtkStreamingTessellator::EdgeCodesToCaseCodesPlusPermutation[64][2] =
{
  /*      code case      C    P */
  /* 000000  0  0  */ {  0,   0 },
  /* 100000  1  1  */ {  1,   0 },
  /* 010000  2  1  */ {  1,   1 },
  /* 110000  3  2a */ {  2,   0 },
  /* 001000  4  1  */ {  1,   2 },
  /* 101000  5  2a */ {  2,   2 },
  /* 011000  6  2a */ {  2,   1 },
  /* 111000  7  3b */ {  5,  11 },
  /* 000100  8  1  */ {  1,  10 },
  /* 100100  9  2a */ {  2,   5 },
  /* 010100 10  2b */ {  3,   1 },
  /* 110100 11  3c */ {  6,   0 },
  /* 001100 12  2a */ {  2,  10 },
  /* 101100 13  3a */ {  4,   0 },
  /* 011100 14  3d */ {  7,   2 },
  /* 111100 15  4a */ {  8,   6 },
  /* 000010 16  1  */ {  1,   6 },
  /* 100010 17  2a */ {  2,   4 },
  /* 010010 18  2a */ {  2,   8 },
  /* 110010 19  3a */ {  4,   1 },
  /* 001010 20  2b */ {  3,   2 },
  /* 101010 21  3d */ {  7,   0 },
  /* 011010 22  3c */ {  6,   1 },
  /* 111010 23  4a */ {  8,   9 },
  /* 000110 24  2a */ {  2,   3 },
  /* 100110 25  3b */ {  5,   0 },
  /* 010110 26  3d */ {  7,   4 },
  /* 110110 27  4a */ {  8,  11 },
  /* 001110 28  3c */ {  6,  10 },
  /* 101110 29  4a */ {  8,   7 },
  /* 011110 30  4b */ {  9,   0 },
  /* 111110 31  5  */ { 10,   7 },
  /* 000001 32  1  */ {  1,   7 },
  /* 100001 33  2b */ {  3,   0 },
  /* 010001 34  2a */ {  2,   7 },
  /* 110001 35  3d */ {  7,   1 },
  /* 001001 36  2a */ {  2,  11 },
  /* 101001 37  3c */ {  6,   2 },
  /* 011001 38  3a */ {  4,   2 },
  /* 111001 39  4a */ {  8,   3 },
  /* 000101 40  2a */ {  2,   9 },
  /* 100101 41  3d */ {  7,  10 },
  /* 010101 42  3c */ {  6,   7 },
  /* 110101 43  4b */ {  9,   2 },
  /* 001101 44  3b */ {  5,   7 },
  /* 101101 45  4a */ {  8,   8 },
  /* 011101 46  4a */ {  8,   4 },
  /* 111101 47  5  */ { 10,   6 },
  /* 000011 48  2a */ {  2,   6 },
  /* 100011 49  3c */ {  6,   4 },
  /* 010011 50  3b */ {  5,   1 },
  /* 110011 51  4a */ {  8,  10 },
  /* 001011 52  3d */ {  7,   7 },
  /* 101011 53  4b */ {  9,   1 },
  /* 011011 54  4a */ {  8,   5 },
  /* 111011 55  5  */ { 10,  10 },
  /* 000111 56  3a */ {  4,  10 },
  /* 100111 57  4a */ {  8,   1 },
  /* 010111 58  4a */ {  8,   2 },
  /* 110111 59  5  */ { 10,   2 },
  /* 001111 60  4a */ {  8,   0 },
  /* 101111 61  5  */ { 10,   1 },
  /* 011111 62  5  */ { 10,   0 },
  /* 111111 63  6  */ { 11,   0 },
};


/* Does this mean anything? If so, then you are either
 * superstitious or much more clever than I (or both?).
 */
/* permutation index, P:  0  1  2  3  4  5  6  7  8  9 10 11 */
/* number of references: 12  9  9  3  4  2  5  6  2  3  7  2 */


/*
 * The array below is a list of all the _positive_
 * permutations of Tetrahedron 0-1-2-3. Given a
 * permutation index, it returns a row of 14 values:
 * these are the vertex numbers of the permuted
 * tetrahedron. The first 4 values are the permuted
 * corner indices, the next 6 values are the
 * permuted edge midpoint indices, and the final
 * entries reference mid-face points inserted
 * to maintain a compatible tetrahedralization.
 *
 * There are 24 entries, 6 for each of the 4 faces of
 * the tetrahedron.
 */
vtkIdType vtkStreamingTessellator::PermutationsFromIndex[24][14] =
{
  /* corners      midpoints          face points   */
  /* POSITIVE ARRANGEMENTS                         */
  { 0, 1, 2, 3,   4, 5, 6, 7, 8, 9,  10, 11, 12, 13 }, /* Face 0-1-2 */
  { 1, 2, 0, 3,   5, 6, 4, 8, 9, 7,  10, 12, 13, 11 },
  { 2, 0, 1, 3,   6, 4, 5, 9, 7, 8,  10, 13, 11, 12 },

  { 0, 3, 1, 2,   7, 8, 4, 6, 9, 5,  11, 13, 12, 10 }, /* Face 0-3-1 */
  { 3, 1, 0, 2,   8, 4, 7, 9, 5, 6,  11, 12, 10, 13 },
  { 1, 0, 3, 2,   4, 7, 8, 5, 6, 9,  11, 10, 13, 12 },

  { 1, 3, 2, 0,   8, 9, 5, 4, 7, 6,  12, 11, 13, 10 }, /* Face 1-3-2 */
  { 3, 2, 1, 0,   9, 5, 8, 7, 6, 4,  12, 13, 10, 11 },
  { 2, 1, 3, 0,   5, 8, 9, 6, 4, 7,  12, 10, 11, 13 },

  { 2, 3, 0, 1,   9, 7, 6, 5, 8, 4,  13, 12, 11, 10 }, /* Face 2-3-0 */
  { 3, 0, 2, 1,   7, 6, 9, 8, 4, 5,  13, 11, 10, 12 },
  { 0, 2, 3, 1,   6, 9, 7, 4, 5, 8,  13, 10, 12, 11 },

  /* NEGATIVE ARRANGEMENTS                         */
  { 0, 2, 1, 3,   6, 5, 4, 7, 9, 8,  10, 13, 12, 11 }, /* Face 0-1-2 */
  { 2, 1, 0, 3,   5, 4, 6, 9, 8, 7,  10, 12, 11, 13 },
  { 1, 0, 2, 3,   4, 6, 5, 8, 7, 9,  10, 11, 13, 12 },

  { 0, 1, 3, 2,   4, 8, 7, 6, 5, 9,  11, 10, 12, 13 }, /* Face 0-3-1 */
  { 1, 3, 0, 2,   8, 7, 4, 5, 9, 6,  11, 12, 13, 10 },
  { 3, 0, 1, 2,   7, 4, 8, 9, 6, 5,  11, 13, 10, 12 },

  { 1, 2, 3, 0,   5, 9, 8, 4, 6, 7,  12, 10, 13, 11 }, /* Face 1-3-2 */
  { 2, 3, 1, 0,   9, 8, 5, 6, 7, 4,  12, 13, 11, 10 },
  { 3, 1, 2, 0,   8, 5, 9, 7, 4, 6,  12, 11, 10, 13 },

  { 2, 0, 3, 1,   6, 7, 9, 5, 4, 8,  13, 10, 11, 12 }, /* Face 2-3-0 */
  { 0, 3, 2, 1,   7, 9, 6, 4, 8, 5,  13, 11, 12, 10 },
  { 3, 2, 0, 1,   9, 6, 7, 8, 5, 4,  13, 12, 10, 11 }
};

/*
 * Below is a list of output tetrahedra. The array is
 * generated by vtkStreamingCaseGenerator.py
 * which also generates the code that references it.
 * Each set of tetrahedra begins with a single integer
 * that is the number of tetrahedra for that particular
 * case. It is followed by 5 integers for each output
 * tetrahedron; the first four numbers on each row are
 * indices of the output tetrahedron. The final number
 * is a bit vector specifying which edges of the
 * tetrahedron are internal to the parent tetrahedron
 * being decomposed.
 *
 * Multiple lists of output tetrahedra may be
 * combined to create the tessellation of a single
 * input tetrahedron.
 */

vtkIdType vtkStreamingTessellator::TetrahedralDecompositions[] =
{
// case 1_0
   2,
   0,  4,  2,  3,
   4,  1,  2,  3,

// case 2a_0
   1,
   3,  4,  5,  1,

// case 2a, 0>1
   2,
   0,  4,  2,  3,
   4,  5,  2,  3,

// case 2a, 1>0
   2,
   0,  4,  5,  3,
   0,  5,  2,  3,

// case 2a, 0=1
   4,
  10,  3,  0,  4,
  10,  3,  4,  5,
  10,  3,  5,  2,
  10,  3,  2,  0,

// case 2b_0
   4,
   0,  4,  9,  3,
   4,  1,  9,  3,
   0,  4,  2,  9,
   4,  1,  2,  9,

// case 3a_0
   1,
   4,  7,  6,  0,

// case 3a, 0>2>3<0
   3,
   1,  3,  2,  4,
   4,  6,  3,  2,
   4,  6,  7,  3,

// case 3a, 0=2>3<0
   5,
   4,  6,  7,  3,
  10,  1,  2,  3,
  10,  2,  6,  3,
  10,  6,  4,  3,
  10,  4,  1,  3,

// case 3a, 3>0=2<3
   5,
   1,  3,  2,  7,
  10,  1,  2,  7,
  10,  2,  6,  7,
  10,  6,  4,  7,
  10,  4,  1,  7,

// case 3a, 0=2=3=0
  11,
   2,  6, 10, 13,
   3,  7, 13, 11,
   4,  1, 10, 11,
  11,  6, 10,  4,
  11,  6, 13, 10,
  11,  6,  7, 13,
  11,  6,  4,  7,
   2, 10, 11, 13,
   1, 10, 11,  2,
   2, 11,  3, 13,
   3,  2,  1, 11,

// case 3b_0
   4,
   0,  7,  4,  2,
   4,  7,  8,  2,
   4,  8,  1,  2,
   7,  3,  8,  2,

// case 3c, 0>1,0>3
   5,
   4,  2,  7,  5,
   4,  2,  0,  7,
   4,  3,  1,  5,
   4,  3,  5,  7,
   3,  5,  7,  2,

// case 3c, 1>0,3>0
   5,
   0,  5,  2,  7,
   0,  5,  7,  4,
   7,  1,  4,  5,
   7,  1,  5,  3,
   3,  5,  7,  2,

// case 3c, 0>1,3>0
   5,
   4,  2,  7,  5,
   4,  2,  0,  7,
   7,  1,  4,  5,
   7,  1,  5,  3,
   3,  5,  7,  2,

// case 3c, 1>0,0>3
   5,
   0,  5,  2,  7,
   0,  5,  7,  4,
   4,  3,  1,  5,
   4,  3,  5,  7,
   3,  5,  7,  2,

// case 3c, 0=1,0>3
   7,
   4,  1,  5,  3,
  10,  0,  4,  7,
  10,  2,  0,  7,
  10,  7,  4,  3,
  10,  2,  7,  3,
  10,  5,  2,  3,
  10,  4,  5,  3,

// case 3c, 3>0,0=1
   7,
   7,  1,  5,  3,
   7,  5,  2,  3,
  10,  0,  4,  7,
  10,  2,  0,  7,
  10,  5,  2,  7,
  10,  4,  5,  7,
   1,  5,  4,  7,

// case 3c, 0=1,0=3
  10,
   4,  1,  5, 11,
  11,  1,  5,  3,
  10,  0,  4,  7,
  10,  2,  0,  7,
  10,  5,  2,  3,
  10,  2,  7,  3,
  10,  7,  4, 11,
  10,  7, 11,  3,
  10,  4,  5, 11,
  10, 11,  5,  3,

// case 3d, 0>4,0>2
   5,
   4,  3,  6,  0,
   4,  3,  8,  6,
   4,  2,  8,  1,
   4,  2,  6,  8,
   2,  3,  6,  8,

// case 3d, 4>0,2>0
   5,
   8,  0,  6,  4,
   8,  0,  3,  6,
   6,  1,  8,  4,
   6,  1,  2,  8,
   2,  3,  6,  8,

// case 3d, 0>4,2>0
   5,
   4,  3,  6,  0,
   4,  3,  8,  6,
   6,  1,  8,  4,
   6,  1,  2,  8,
   2,  3,  6,  8,

// case 3d, 4>0,0>2
   5,
   8,  0,  6,  4,
   8,  0,  3,  6,
   4,  2,  8,  1,
   4,  2,  6,  8,
   2,  3,  6,  8,

// case 3d, 0=4,0>2
   7,
   4,  1,  2,  8,
  11,  4,  0,  6,
  11,  0,  3,  6,
  11,  2,  4,  6,
  11,  3,  2,  6,
  11,  3,  8,  2,
  11,  8,  4,  2,

// case 3d, 2>0,0=4
   7,
   6,  2,  8,  1,
   6,  8,  2,  3,
  11,  4,  0,  6,
  11,  0,  3,  6,
   8, 11,  3,  6,
   8,  4, 11,  6,
   1,  6,  4,  8,

// case 3d, 0=4,0=2
  10,
   4,  1, 10,  8,
  10,  2,  8,  1,
  11,  4,  0,  6,
  11,  0,  3,  6,
  11,  3,  8,  2,
  11,  3,  2,  6,
  11, 10,  4,  6,
  11, 10,  6,  2,
   8,  4, 11, 10,
  11, 10,  2,  8,

// case 4a_0
   2,
   7,  8,  9,  3,
   7,  9,  8,  6,

// case 4a, 5>4>3
   4,
   8,  0,  6,  1,
   8,  0,  7,  6,
   9,  1,  6,  2,
   9,  1,  8,  6,

// case 4a, 3<4>5
   4,
   8,  0,  6,  1,
   8,  0,  7,  6,
   8,  2,  6,  9,
   8,  2,  1,  6,

// case 4a, 3>4<5
   4,
   6,  9,  8,  1,
   6,  9,  1,  2,
   6,  7,  0,  1,
   6,  7,  1,  8,

// case 4a, 3=4>5
   6,
   6,  7,  0, 11,
   6,  0,  1, 11,
   6,  7, 11,  8,
   6, 11,  1,  8,
   1,  2,  6,  8,
   2,  6,  8,  9,

// case 4a, 5>4,3=4
   6,
   6,  7,  0, 11,
   6,  0,  1, 11,
   6,  7, 11,  8,
   6, 11,  1,  8,
   1,  2,  6,  9,
   1,  6,  8,  9,

// case 4a, 3=4=5
   8,
   6,  7,  0, 11,
   6,  0,  1, 11,
   6,  7, 11,  8,
   6, 11,  1,  8,
   6,  1,  2, 12,
   6,  2,  9, 12,
   6,  9,  8, 12,
   6,  8,  1, 12,

// case 4b, 2>1,3>2,4>3,4>1
   6,
   6,  8,  1,  5,
   6,  8,  0,  1,
   6,  8,  7,  0,
   6,  8,  2,  7,
   7,  8,  2,  3,
   6,  8,  5,  2,

// case 4b, 2>1,3>2,3>4,4>1
   6,
   6,  8,  1,  5,
   6,  8,  7,  1,
   6,  7,  0,  1,
   8,  7,  3,  2,
   6,  8,  5,  2,
   6,  8,  2,  7,

// case 4b, 2>1,2>3,4>3,4>1
   6,
   6,  8,  5,  2,
   6,  8,  2,  3,
   6,  8,  3,  7,
   6,  8,  7,  0,
   6,  8,  0,  1,
   6,  8,  1,  5,

// case 4b, 1=2,3>2,3>4,4>1
   9,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1,  7,
  10,  7,  1,  8,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6,  2,  7,  8,
   6,  5,  2,  8,
   7,  8,  2,  3,

// case 4b, 1=2,2>3,4>3,1>4
   9,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1,  8,
  10,  7,  0,  8,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6,  3,  7,  8,
   6,  5,  3,  8,
   6,  5,  2,  3,

// case 4b, 1=2,2>3,4>3,4>1
   9,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1,  8,
  10,  7,  0,  8,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6,  3,  7,  8,
   6,  5,  2,  8,
   6,  2,  3,  8,

// case 4b, 1=2,3>2,3=4,4>1
  11,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1, 11,
  10, 11,  1,  8,
  10,  0, 11,  7,
  10,  7, 11,  8,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6,  2,  7,  8,
   6,  5,  2,  8,
   7,  8,  2,  3,

// case 4b, 4>1=2=3,4>3
  12,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1,  8,
  10,  7,  0,  8,
  13,  6,  2,  5,
  13,  3,  7,  8,
  13,  2,  3,  8,
  13,  2,  8,  5,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6, 13,  7,  8,
   6,  5, 13,  8,

// case 4b, 1=2=3>4,1>4
  12,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1,  7,
  10,  7,  1,  8,
  13,  6,  2,  5,
  13,  3,  7,  8,
  13,  2,  3,  5,
  13,  3,  8,  5,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6, 13,  7,  8,
   6,  5, 13,  8,

// case 4b, 1=2=3=4=1
  16,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1, 11,
  10, 11,  1,  8,
  10,  0, 11,  7,
  10,  7, 11,  8,
  13,  6,  2,  5,
  13,  3,  7,  8,
  13,  2,  3, 12,
  13,  2, 12,  5,
  13, 12,  3,  8,
  13, 12,  5,  8,
   6,  7, 10,  8,
   6, 10,  5,  8,
   6,  5, 13,  8,
   6, 13,  7,  8,

// case 5_0
   2,
   7,  8,  9,  3,
   6,  5,  2,  9,

// case 5, 1>2,3>4
   5,
   5,  7,  1,  8,
   5,  7,  0,  1,
   5,  7,  6,  0,
   5,  7,  9,  6,
   5,  7,  8,  9,

// case 5, 1>2,4>3
   5,
   0,  5,  6,  7,
   0,  5,  7,  8,
   0,  5,  8,  1,
   5,  7,  9,  6,
   5,  7,  8,  9,

// case 5, 1=2,3>4
   8,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1,  7,
  10,  7,  1,  8,
  10,  8,  5,  9,
  10,  6,  7,  9,
  10,  7,  8,  9,
  10,  5,  6,  9,

// case 5, 1=2,3=4
  10,
  10,  6,  0,  7,
  10,  1,  5,  8,
  10,  0,  1, 11,
  10, 11,  1,  8,
  10,  0, 11,  7,
  10,  7, 11,  8,
  10,  8,  5,  9,
  10,  6,  7,  9,
  10,  7,  8,  9,
  10,  5,  6,  9,

// case 6_0
   4,
   7,  8,  9,  3,
   6,  5,  2,  9,
   4,  1,  5,  8,
   0,  4,  6,  7,

// case 6_1
   4,
   6,  4,  5,  8,
   6,  5,  9,  8,
   6,  9,  7,  8,
   6,  7,  4,  8,

};

VTK_ABI_NAMESPACE_END
