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

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

#include "vtkBitArray.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataSetAttributes.h"
#include "vtkExtentTranslator.h"
#include "vtkHyperTreeGrid.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkStreamingDemandDrivenPipeline.h"

#include "vtkDoubleArray.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"

#include "vtkIncrementalPointLocator.h"
#include "vtkMergePoints.h"

#include "vtkHyperTreeGridNonOrientedGeometryCursor.h"
#include "vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight.h"

#include <set>
#include <vector>

static const unsigned int VonNeumannCursors3D[] = { 0, 1, 2, 4, 5, 6 };
static const unsigned int VonNeumannOrientations3D[] = { 2, 1, 0, 0, 1, 2 };
static const unsigned int VonNeumannOffsets3D[] = { 0, 0, 0, 1, 1, 1 };

static const unsigned int EdgeIndices[3][2][4] = { { { 3, 11, 7, 8 }, { 1, 10, 5, 9 } },
  { { 0, 9, 4, 8 }, { 2, 10, 6, 11 } }, { { 0, 1, 2, 3 }, { 4, 5, 6, 7 } } };

const unsigned char FULL_WORK_FACES = (unsigned char)255;

vtkStandardNewMacro(vtkHyperTreeGridGeometry);

//-----------------------------------------------------------------------------
vtkHyperTreeGridGeometry::vtkHyperTreeGridGeometry()
{
  // Create storage for corners of leaf cells
  this->Points = vtkPoints::New();

  // Create storage for untructured leaf cells
  this->Cells = vtkCellArray::New();

  // Default dimension is 0
  this->Dimension = 0;

  // Default orientation is 0
  this->Orientation = 0;

  // Default orientation is 0
  this->BranchFactor = 0;

  // Default Locator is 0
  this->Merging = false;
  this->Locator = nullptr;

  // Default interface values
  this->HasInterface = false;
  this->Normals = nullptr;
  this->Intercepts = nullptr;
  this->FaceIDs = vtkIdList::New();
  this->FacePoints = vtkPoints::New();
  this->FacePoints->SetNumberOfPoints(4);
  this->FacesA = vtkIdTypeArray::New();
  this->FacesA->SetNumberOfComponents(2);
  this->FacesB = vtkIdTypeArray::New();
  this->FacesB->SetNumberOfComponents(2);
  this->FaceScalarsA = vtkDoubleArray::New();
  this->FaceScalarsA->SetNumberOfTuples(4);
  this->FaceScalarsB = vtkDoubleArray::New();
  this->FaceScalarsB->SetNumberOfTuples(4);
}

//-----------------------------------------------------------------------------
vtkHyperTreeGridGeometry::~vtkHyperTreeGridGeometry()
{
  if (this->Points)
  {
    this->Points->Delete();
    this->Points = nullptr;
  }

  if (this->Cells)
  {
    this->Cells->Delete();
    this->Cells = nullptr;
  }

  if (this->FacePoints)
  {
    this->FacePoints->Delete();
    this->FacePoints = nullptr;
  }

  if (this->FaceIDs)
  {
    this->FaceIDs->Delete();
    this->FaceIDs = nullptr;
  }

  if (this->FacesA)
  {
    this->FacesA->Delete();
    this->FacesA = nullptr;
  }

  if (this->FacesB)
  {
    this->FacesB->Delete();
    this->FacesB = nullptr;
  }

  if (this->FaceScalarsA)
  {
    this->FaceScalarsA->Delete();
    this->FaceScalarsA = nullptr;
  }

  if (this->FaceScalarsB)
  {
    this->FaceScalarsB->Delete();
    this->FaceScalarsB = nullptr;
  }
}

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

  if (this->Points)
  {
    os << indent << "Points:\n";
    this->Points->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "Points: ( none )\n";
  }

  if (this->Cells)
  {
    os << indent << "Cells:\n";
    this->Cells->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "Cells: ( none )\n";
  }

  os << indent << "Dimension: " << this->Dimension << endl;
  os << indent << "Orientation: " << this->Orientation << endl;
  os << indent << "Merging: " << this->Merging << endl;
  os << indent << "HasInterface: " << this->HasInterface << endl;
  if (this->Normals)
  {
    os << indent << ":\n";
    this->Normals->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "Normals: ( none )\n";
  }
  if (this->Intercepts)
  {
    os << indent << ":\n";
    this->Intercepts->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "Intercepts: ( none )\n";
  }
  if (this->FacePoints)
  {
    os << indent << ":\n";
    this->FacePoints->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "FacePoints: ( none )\n";
  }
  if (this->FaceIDs)
  {
    os << indent << ":\n";
    this->FaceIDs->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "FaceIDs: ( none )\n";
  }
  os << indent << "EdgesA:";
  for (unsigned int i = 0; i < 12; ++i)
  {
    os << " " << EdgesA[i];
  }
  os << endl;
  os << indent << "EdgesB:";
  for (unsigned int i = 0; i < 12; ++i)
  {
    os << " " << EdgesB[i];
  }
  os << endl;
  if (this->FacesA)
  {
    os << indent << ":\n";
    this->FacesA->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "FacesA: ( none )\n";
  }
  if (this->FacesB)
  {
    os << indent << ":\n";
    this->FacesB->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "FacesB: ( none )\n";
  }
  if (this->FaceScalarsA)
  {
    os << indent << ":\n";
    this->FaceScalarsA->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "FaceScalarsA: ( none )\n";
  }
  if (this->FaceScalarsB)
  {
    os << indent << ":\n";
    this->FaceScalarsB->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "FaceScalarsB: ( none )\n";
  }
}

//----------------------------------------------------------------------------
int vtkHyperTreeGridGeometry::FillOutputPortInformation(int, vtkInformation* info)
{
  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
  return 1;
}

//-----------------------------------------------------------------------------
int vtkHyperTreeGridGeometry::ProcessTrees(vtkHyperTreeGrid* input, vtkDataObject* outputDO)
{
  // Downcast output data object to polygonal data set
  vtkPolyData* output = vtkPolyData::SafeDownCast(outputDO);
  if (!output)
  {
    vtkErrorMacro("Incorrect type of output: " << outputDO->GetClassName());
    return 0;
  }

  // Retrieve useful grid parameters for speed of access
  this->Dimension = input->GetDimension();
  this->Orientation = input->GetOrientation();

  this->BranchFactor = static_cast<int>(input->GetBranchFactor());

  // Initialize output cell data
  this->InData = input->GetPointData();
  this->OutData = output->GetCellData();
  this->OutData->CopyAllocate(this->InData);

  // Retrieve material mask
  this->Mask = input->HasMask() ? input->GetMask() : nullptr;

  // Retrieve pure material mask
  this->PureMask = input->GetPureMask();

  // Retrieve interface data when relevant
  this->HasInterface = input->GetHasInterface();
  if (this->HasInterface)
  {
    this->Normals =
      vtkDoubleArray::SafeDownCast(this->InData->GetArray(input->GetInterfaceNormalsName()));
    this->Intercepts =
      vtkDoubleArray::SafeDownCast(this->InData->GetArray(input->GetInterfaceInterceptsName()));
  } // this->HasInterface

  // Create storage for corners of leaf cells
  if (this->Points)
  {
    this->Points->Delete();
  }
  this->Points = vtkPoints::New();

  // Create storage for untructured leaf cells
  if (this->Cells)
  {
    this->Cells->Delete();
  }
  this->Cells = vtkCellArray::New();

  // JB Initialize a Locator
  if (this->Merging)
  {
    if (this->Locator)
    {
      this->Locator->Delete();
    }
    this->Locator = vtkMergePoints::New();
    this->Locator->InitPointInsertion(this->Points, input->GetBounds());
  }

  // Iterate over all hyper trees
  vtkIdType index;
  vtkHyperTreeGrid::vtkHyperTreeGridIterator it;
  input->InitializeTreeIterator(it);
  if (this->Dimension == 3)
  {
    vtkNew<vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight> cursor;
    while (it.GetNextTree(index))
    {
      // Initialize new cursor at root of current tree
      // In 3 dimensions, von Neumann neighborhood information is needed
      input->InitializeNonOrientedVonNeumannSuperCursorLight(cursor, index);
      // Build geometry recursively
      this->RecursivelyProcessTree3D(cursor, FULL_WORK_FACES);
    } // it
  }
  else
  {
    vtkNew<vtkHyperTreeGridNonOrientedGeometryCursor> cursor;
    while (it.GetNextTree(index))
    {
      // Initialize new cursor at root of current tree
      // Otherwise, geometric properties of the cells suffice
      input->InitializeNonOrientedGeometryCursor(cursor, index);
      // Build geometry recursively
      this->RecursivelyProcessTreeNot3D(cursor);
    } // it
  }   // else

  // Set output geometry and topology
  output->SetPoints(this->Points);
  if (this->Dimension == 1)
  {
    output->SetLines(this->Cells);
  }
  else
  {
    output->SetPolys(this->Cells);
  }

  if (this->Points)
  {
    this->Points->Delete();
    this->Points = nullptr;
  }
  if (this->Cells)
  {
    this->Cells->Delete();
    this->Cells = nullptr;
  }

  if (this->Locator)
  {
    this->Locator->Delete();
    this->Locator = nullptr;
  }

  return 1;
}

//----------------------------------------------------------------------------
void vtkHyperTreeGridGeometry::RecursivelyProcessTreeNot3D(
  vtkHyperTreeGridNonOrientedGeometryCursor* cursor)
{
  if (this->Mask ? this->Mask->GetValue(cursor->GetGlobalNodeIndex()) : false)
  {
    return;
  }

  // Create geometry output if cursor is at leaf
  if (cursor->IsLeaf())
  {
    // Cursor is at leaf, process it depending on its dimension
    switch (this->Dimension)
    {
      case 1:
        this->ProcessLeaf1D(cursor);
        break;
      case 2:
        this->ProcessLeaf2D(cursor);
        break;
      default:
        break;
    } // switch ( this->Dimension )
    return;
  } // if ( cursor->IsLeaf() )

  // FR Il existe une filles qui n'est pas dans le mat, on la chercher partout
  unsigned int numChildren = cursor->GetNumberOfChildren();
  for (unsigned int ichild = 0; ichild < numChildren; ++ichild)
  {
    cursor->ToChild(ichild);
    // Recurse
    this->RecursivelyProcessTreeNot3D(cursor);
    cursor->ToParent();
  } // ichild
}

//----------------------------------------------------------------------------
// JB Meme code que vtkAdaptativeDataSetSurfaceFiltre ??
void vtkHyperTreeGridGeometry::ProcessLeaf1D(vtkHyperTreeGridNonOrientedGeometryCursor* cursor)
{
  // Cell at cursor center is a leaf, retrieve its global index
  vtkIdType inId = cursor->GetGlobalNodeIndex();
  if (inId < 0)
  {
    return;
  }

  // In 1D the geometry is composed of edges, create storage for endpoint IDs
  vtkIdType ids[2];

  // First endpoint is at origin of cursor
  const double* origin = cursor->GetOrigin();

  // Second endpoint is at origin of cursor plus its length
  double pt[3];
  memcpy(pt, origin, 3 * sizeof(double));
  pt[this->Orientation] += cursor->GetSize()[this->Orientation];

  if (this->Locator)
  {
    this->Locator->InsertUniquePoint(origin, ids[0]);
    this->Locator->InsertUniquePoint(pt, ids[1]);
  }
  else
  {
    ids[0] = this->Points->InsertNextPoint(origin);
    ids[1] = this->Points->InsertNextPoint(pt);
  }

  // Insert edge into 1D geometry
  vtkIdType outId = this->Cells->InsertNextCell(2, ids);

  // Copy edge data from that of the cell from which it comes
  this->OutData->CopyData(this->InData, inId, outId);
}

//----------------------------------------------------------------------------
// JB Meme code que vtkAdaptativeDataSetSurfaceFiltre ??
void vtkHyperTreeGridGeometry::ProcessLeaf2D(vtkHyperTreeGridNonOrientedGeometryCursor* cursor)

{
  // Cell at cursor center is a leaf, retrieve its global index
  vtkIdType inId = cursor->GetGlobalNodeIndex();
  assert(inId >= 0);

  // Reset interface variables if needed
  if (this->HasInterface)
  {
    size_t int12sz = 12 * sizeof(vtkIdType);
    memset(this->EdgesA, -1, int12sz);
    memset(this->EdgesB, -1, int12sz);
    this->FacesA->Reset();
    this->FacesB->Reset();
  } // if ( this->HasInterface )

  // Insert face into 2D geometry depending on orientation
  // this->AddFace( inId, cursor->GetOrigin(), cursor->GetSize(), 0, this->Orientation );
  this->AddFace2(inId, inId, cursor->GetOrigin(), cursor->GetSize(), 0, this->Orientation);
}

//----------------------------------------------------------------------------
void vtkHyperTreeGridGeometry::RecursivelyProcessTree3D(
  vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight* cursor, unsigned char crtWorkFaces)
{
  // FR Traitement specifique pour la maille fille centrale en raffinement 3
  // Create geometry output if cursor is at leaf
  if (cursor->IsLeaf() || cursor->IsMasked())
  {
    // Cursor is at leaf, process it depending on its dimension
    this->ProcessLeaf3D(cursor); // JBVTK9 ProcessLeaf3D2 prends en compte les interfaces... :)?
    return;
  } // if ( cursor->IsLeaf() )

  // FR Parce que le curseur est un super curseur
  bool pureMask = false;
  if (this->Mask != nullptr)
  {
    // JB Question : que fait le PureMaterialMask quand un masque est mis sur un coarse et pas
    // toutes les filles
    pureMask = this->PureMask->GetValue(cursor->GetGlobalNodeIndex()) != 0;
  }
  if (!pureMask)
  {
    // FR La maille courante est donc pure (pureMask == false car GetPureMask()->GetValue(..)
    // retourne false)
    std::set<int> childList;

    const unsigned int numChildren = cursor->GetNumberOfChildren();
    std::vector<unsigned char> workFaces(numChildren, 0);

    // FR Toutes les filles sont dans le materiau, on traite que les filles du bord
    for (unsigned int f = 0; f < 3; ++f) // dimension
    {
      for (unsigned int o = 0; o < 2; ++o) // gauche, centre, droite
      {
        int neighborIdx = (2 * o - 1) * (f + 1);
        if ((crtWorkFaces & (1 << (this->Dimension + neighborIdx))))
        {
          bool isValidN = cursor->HasTree(this->Dimension + neighborIdx);
          vtkIdType idN = 0;
          if (isValidN)
          {
            idN = cursor->GetGlobalNodeIndex(this->Dimension + neighborIdx);
          }
          if (!isValidN || (this->Mask && this->PureMask->GetValue(idN)))
          {
            // FR La maille voisine n'existe pas ou n'est pas pure (PureMask->GetValue(id) retourne
            // false) FR Fille du bord
            int iMin = (f == 0 && o == 1) ? this->BranchFactor - 1 : 0;
            int iMax = (f == 0 && o == 0) ? 1 : this->BranchFactor;
            int jMin = (f == 1 && o == 1) ? this->BranchFactor - 1 : 0;
            int jMax = (f == 1 && o == 0) ? 1 : this->BranchFactor;
            int kMin = (f == 2 && o == 1) ? this->BranchFactor - 1 : 0;
            int kMax = (f == 2 && o == 0) ? 1 : this->BranchFactor;
            for (int i = iMin; i < iMax; ++i)
            {
              for (int j = jMin; j < jMax; ++j)
              {
                for (int k = kMin; k < kMax; ++k)
                {
                  unsigned int ichild = i + this->BranchFactor * (j + this->BranchFactor * k);

                  // FR Les mailles de coin peuvent etre sollicitees plusieurs fois suivant chacune
                  // des faces
                  childList.insert(ichild);
                  workFaces[ichild] |= (1 << (this->Dimension + neighborIdx));
                } // k
              }   // j
            }     // i
          }       // if ...
        }
      } // o
    }   // f
    for (std::set<int>::iterator it = childList.begin(); it != childList.end(); ++it)
    {
      cursor->ToChild(*it);
      this->RecursivelyProcessTree3D(cursor, workFaces[*it]);
      cursor->ToParent();
    } // ichild
    return;
  }
  // FR Il existe une filles qui n'est pas dans le mat, on la chercher partout
  unsigned int numChildren = cursor->GetNumberOfChildren();
  for (unsigned int ichild = 0; ichild < numChildren; ++ichild)
  {
    cursor->ToChild(ichild);
    this->RecursivelyProcessTree3D(cursor, FULL_WORK_FACES);
    cursor->ToParent();
  } // ichild
}

//----------------------------------------------------------------------------
// JB Meme code que vtkAdaptativeDataSetSurfaceFiltre ??
void vtkHyperTreeGridGeometry::ProcessLeaf3D(
  vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight* superCursor)
{
  // Cell at cursor center is a leaf, retrieve its global index, and mask
  vtkIdType inId = superCursor->GetGlobalNodeIndex();
  if (inId < 0)
  {
    return;
  }
  unsigned level = superCursor->GetLevel();

  // Reset interface variables if needed
  if (this->HasInterface)
  {
    size_t int12sz = 12 * sizeof(vtkIdType);
    memset(this->EdgesA, -1, int12sz);
    memset(this->EdgesB, -1, int12sz);
    this->FacesA->Reset();
    this->FacesB->Reset();

    // Retrieve intercept type
    this->Intercepts->GetComponent(inId, 2);
  } // if ( this->HasInterface )

  // Iterate over all cursors of Von Neumann neighborhood around center
  unsigned int nc = superCursor->GetNumberOfCursors() - 1;
  for (unsigned int c = 0; c < nc; ++c)
  {
    // Retrieve cursor to neighbor across face
    // Retrieve tree, leaf flag, and mask of neighbor cursor
    bool leafN;
    vtkIdType idN;
    unsigned int levelN;
    vtkHyperTree* treeN = superCursor->GetInformation(VonNeumannCursors3D[c], levelN, leafN, idN);
    int maskedN = superCursor->IsMasked(VonNeumannCursors3D[c]);

    // In 3D masked and unmasked cells are handled differently:
    // . If cell is unmasked, and face neighbor is a masked leaf, or no such neighbor
    //   exists, then generate face.
    // . If cell is masked, and face neighbor exists and is an unmasked leaf, then
    //   generate face, breaking ties at same level. This ensures that faces between
    //   unmasked and masked cells will be generated once and only once.
    if ((!superCursor->IsMasked() && (!treeN || maskedN)) ||
      (superCursor->IsMasked() && treeN && leafN && levelN < level && !maskedN))
    {
      // Generate face with corresponding normal and offset
      this->AddFace(superCursor->IsMasked() ? idN : inId, superCursor->GetOrigin(),
        superCursor->GetSize(), VonNeumannOffsets3D[c], VonNeumannOrientations3D[c]);
    }
    /*
        // In 3D masked and unmasked cells are handled differently
        if ( superCursor->IsMasked() )
        {
          // If cell is masked, and face neighbor exists and is an unmasked leaf,
          // generate face, breaking ties at same level. This ensures that faces
          // between unmasked and masked cells will be generated once and only once.
          if ( treeN && leafN && levelN < level && ! maskedN )
          {
            // Generate face with corresponding normal and offset
            this->AddFace2( inId, ?, superCursor->GetOrigin(), superCursor->GetSize(),
                           VonNeumannOffsets3D[c], VonNeumannOrientations3D[c] );
          }
        } // if ( superCursor->IsMasked() )
        else
        {
          // If cell is unmasked, and face neighbor is a masked leaf, or no
          // such neighbor exists, generate face except in some interface cases.
          bool addFace = ! treeN || ( leafN && maskedN );
          bool create = true;
          if ( ! addFace && type < 2. )
          {
            // Mixed cells must be handled but not created
            create  = false;
            addFace = true;
          }
          if ( ( ! addFace || ! create ) && this->HasInterface && leafN )
          {
            // Face must be created if neighbor is a mixed cell
            create  = this->Intercepts->GetComponent( idN, 2 ) < 2.;
            addFace |= create;
          }
          if ( addFace )
          {
            // Generate or handle face with corresponding normal and offset
            this->AddFace2( inId, ?, superCursor->GetOrigin(), superCursor->GetSize(),
                           VonNeumannOffsets3D[c], VonNeumannOrientations3D[c],
                           create );
          } // if ( addFace )
        } // else
    */
  } // c

  // Handle interfaces separately
  if (this->HasInterface)
  {
    // Create face A when its edges are present
    vtkIdType nA = this->FacesA->GetNumberOfTuples();
    if (nA > 0)
    {
      this->FaceIDs->Reset();
      vtkIdType i0 = 0;
      vtkIdType edge0[2];
      this->FacesA->GetTypedTuple(i0, edge0);
      this->FaceIDs->InsertNextId(this->EdgesA[edge0[1]]);
      while (edge0[0] != edge0[1])
      {
        // Iterate over edges of face A
        for (vtkIdType i = 0; i < nA; ++i)
        {
          // Seek next edge then break out from loop
          vtkIdType edge[2];
          this->FacesA->GetTypedTuple(i, edge);
          if (i0 != i)
          {
            if (edge[0] == edge0[1])
            {
              edge0[1] = edge[1];
              i0 = i;
              break;
            }
            if (edge[1] == edge0[1])
            {
              edge0[1] = edge[0];
              i0 = i;
              break;
            }
          } // if ( i != i0 )
        }   // nA
        this->FaceIDs->InsertNextId(this->EdgesA[edge0[1]]);
      } // while ( edge0[0] != edge0[1] )

      // Create new face
      vtkIdType outId = this->Cells->InsertNextCell(this->FaceIDs);

      // Copy face data from that of the cell from which it comes
      this->OutData->CopyData(this->InData, inId, outId);
    } // if ( nA > 0 )

    // Create face B when its vertices are present
    vtkIdType nB = this->FacesB->GetNumberOfTuples();
    if (nB > 0)
    {
      this->FaceIDs->Reset();
      int i0 = 0;
      vtkIdType edge0[2];
      this->FacesB->GetTypedTuple(i0, edge0);
      this->FaceIDs->InsertNextId(this->EdgesB[edge0[1]]);
      while (edge0[0] != edge0[1])
      {
        // Iterate over faces B
        for (vtkIdType i = 0; i < nB; ++i)
        {
          // Seek next edge then break out from loop
          vtkIdType edge[2];
          this->FacesB->GetTypedTuple(i, edge);
          if (i0 != i)
          {
            if (edge[0] == edge0[1])
            {
              edge0[1] = edge[1];
              i0 = i;
              break;
            }
            if (edge[1] == edge0[1])
            {
              edge0[1] = edge[0];
              i0 = i;
              break;
            }
          } // if ( i0 != i )
        }   // nB
        this->FaceIDs->InsertNextId(this->EdgesB[edge0[1]]);
      } // while ( edge0[0] != edge0[1] )

      // Create new face
      vtkIdType outId = this->Cells->InsertNextCell(this->FaceIDs);

      // Copy face data from that of the cell from which it comes
      this->OutData->CopyData(this->InData, inId, outId);
    } // if ( nB > 0 )
  }   // if ( this->HasInterface )
}

//----------------------------------------------------------------------------
void vtkHyperTreeGridGeometry::AddFace(vtkIdType useId, const double* origin, const double* size,
  unsigned int offset, unsigned int orientation)
{
  // cerr << "AddFace" << std::endl;

  // Storage for point coordinates
  double pt[] = { 0., 0., 0. };

  // Storage for face vertex IDs
  vtkIdType ids[4];

  // First cell vertex is always at origin of cursor
  memcpy(pt, origin, 3 * sizeof(double));

  if (this->Locator)
  {
    if (offset)
    {
      // Offset point coordinate as needed
      pt[orientation] += size[orientation];
    }
    this->Locator->InsertUniquePoint(pt, ids[0]);
    // Create other face vertices depending on orientation
    unsigned int axis1 = orientation ? 0 : 1;
    unsigned int axis2 = orientation == 2 ? 1 : 2;
    pt[axis1] += size[axis1];
    this->Locator->InsertUniquePoint(pt, ids[1]);
    pt[axis2] += size[axis2];
    this->Locator->InsertUniquePoint(pt, ids[2]);
    pt[axis1] = origin[axis1];
    this->Locator->InsertUniquePoint(pt, ids[3]);
  }
  else
  {
    if (offset)
    {
      // Offset point coordinate as needed
      pt[orientation] += size[orientation];
    }
    ids[0] = this->Points->InsertNextPoint(pt);
#ifdef TRACE
    cerr << "Point #" << ids[0] << " : ";
    for (unsigned int ipt = 0; ipt < 3; ++ipt)
    {
      cerr << pt[ipt] << " ";
    }
    cerr << std::endl;
#endif
    // Create other face vertices depending on orientation
    unsigned int axis1 = orientation ? 0 : 1;
    unsigned int axis2 = orientation == 2 ? 1 : 2;
    pt[axis1] += size[axis1];
    ids[1] = this->Points->InsertNextPoint(pt);
#ifdef TRACE
    cerr << "Point #" << ids[1] << " : ";
    for (unsigned int ipt = 0; ipt < 3; ++ipt)
    {
      cerr << pt[ipt] << " ";
    }
    cerr << std::endl;
#endif
    pt[axis2] += size[axis2];
    ids[2] = this->Points->InsertNextPoint(pt);
#ifdef TRACE
    cerr << "Point #" << ids[2] << " : ";
    for (unsigned int ipt = 0; ipt < 3; ++ipt)
    {
      cerr << pt[ipt] << " ";
    }
    cerr << std::endl;
#endif
    pt[axis1] = origin[axis1];
    ids[3] = this->Points->InsertNextPoint(pt);
#ifdef TRACE
    cerr << "Point #" << ids[3] << " : ";
    for (unsigned int ipt = 0; ipt < 3; ++ipt)
    {
      cerr << pt[ipt] << " ";
    }
    cerr << std::endl;
    cerr << "Face #" << this->Cells->GetNumberOfCells() << " : ";
    for (unsigned int ipt = 0; ipt < 3; ++ipt)
    {
      cerr << ids[ipt] << " ";
    }
    cerr << std::endl;
#endif
  }

  // Insert next face
  vtkIdType outId = this->Cells->InsertNextCell(4, ids);

  // Copy face data from that of the cell from which it comes
  this->OutData->CopyData(this->InData, useId, outId);
}
//----------------------------------------------------------------------------
void vtkHyperTreeGridGeometry::AddFace2(vtkIdType inId, vtkIdType useId, const double* origin,
  const double* size, unsigned int offset, unsigned int orientation, bool create)
{
  // First cell vertex is always at origin of cursor
  double pt[3];
  memcpy(pt, origin, 3 * sizeof(double));
  if (offset)
  {
    // Offset point coordinate as needed
    pt[orientation] += size[orientation];
  }

  // Storage for face vertex IDs
  vtkIdType ids[4];
  unsigned int nPts = 4;

  // Keep track of face axes
  unsigned int axis1 = orientation ? 0 : 1;
  unsigned int axis2 = orientation == 2 ? 1 : 2;

  if (this->HasInterface)
  {
    // Retrieve intercept tuple and type
    double* inter = this->Intercepts->GetTuple(inId);
    double type = inter[2];

    // Distinguish cases depending on intercept type
    if (type < 2)
    {
      // Create interface intersection points
      this->FacePoints->SetPoint(0, pt);
      pt[axis1] += size[axis1];
      this->FacePoints->SetPoint(1, pt);
      pt[axis2] += size[axis2];
      this->FacePoints->SetPoint(2, pt);
      pt[axis1] = origin[axis1];
      this->FacePoints->SetPoint(3, pt);

      // Create interface intersection faces
      double coordsA[3];
      double* normal = this->Normals->GetTuple(inId);
      for (vtkIdType pId = 0; pId < 4; ++pId)
      {
        // Retrieve vertex coordinates
        this->FacePoints->GetPoint(pId, coordsA);

        // Set face scalars
        if (type != 1.)
        {
          double val =
            inter[0] + normal[0] * coordsA[0] + normal[1] * coordsA[1] + normal[2] * coordsA[2];
          this->FaceScalarsA->SetTuple1(pId, val);
        } // if ( type != 1. )
        if (type != -1.)
        {
          double val =
            inter[1] + normal[0] * coordsA[0] + normal[1] * coordsA[1] + normal[2] * coordsA[2];
          this->FaceScalarsB->SetTuple1(pId, val);
        } // if ( type != -1. )
      }   // p

      // Storage for points
      double coordsB[3];
      double coordsC[3];
      nPts = 0;

      // Distinguish between relevant types
      if (type == 1)
      {
        // Take negative values in B
        vtkIdType pair[2];
        unsigned int indPair = 0;

        // Loop over face vertices
        for (int p = 0; p < 4; ++p)
        {
          // Retrieve vertex coordinates
          this->FacePoints->GetPoint(p, coordsA);

          // Retrieve vertex scalars
          double A = this->FaceScalarsB->GetTuple1(p);
          double B = this->FaceScalarsB->GetTuple1((p + 1) % 4);

          // Add point when necessary
          if (create && A <= 0.)
          {
            ids[nPts] = this->Points->InsertNextPoint(coordsA);
            ++nPts;
          }
          if (A * B < 0)
          {
            unsigned int i = EdgeIndices[orientation][offset][p];
            if (this->EdgesB[i] == -1)
            {
              // Compute barycenter of A and B
              this->FacePoints->GetPoint((p + 1) % 4, coordsB);
              for (int j = 0; j < 3; ++j)
              {
                coordsC[j] = (B * coordsA[j] - A * coordsB[j]) / (B - A);
              }
              this->EdgesB[i] = this->Points->InsertNextPoint(coordsC);
            } // if ( this->EdgesB[i] == -1 )

            // Update points
            ids[nPts] = this->EdgesB[i];
            ++nPts;
            if (indPair)
            {
              pair[1] = i;
            }
            else
            {
              pair[0] = i;
            };
            ++indPair;
          } // if ( A * B < 0 )
        }   // p

        // Insert pair only if it makes sense
        if (indPair == 2)
        {
          this->FacesB->InsertNextTypedTuple(pair);
        }
      } // if (type = 1 )

      else if (!type)
      {
        // Take positive values in A
        vtkIdType pairA[2];
        unsigned int indPairA = 0;

        // Take negative values in B
        vtkIdType pairB[2];
        unsigned int indPairB = 0;

        // Loop over face vertices
        for (int p = 0; p < 4; ++p)
        {
          // Retrieve vertex coordinates
          this->FacePoints->GetPoint(p, coordsA);

          // Retrieve vertex scalars
          double A1 = this->FaceScalarsA->GetTuple1(p);
          double B1 = this->FaceScalarsA->GetTuple1((p + 1) % 4);
          double A2 = this->FaceScalarsB->GetTuple1(p);
          double B2 = this->FaceScalarsB->GetTuple1((p + 1) % 4);

          // Add point when necessary
          if (create && A1 >= 0. && A2 <= 0.)
          {
            ids[nPts] = this->Points->InsertNextPoint(coordsA);
            ++nPts;
          }
          if (A1 < 0. && A1 * B1 < 0.)
          {
            unsigned int i = EdgeIndices[orientation][offset][p];
            if (this->EdgesA[i] == -1)
            {
              // Compute barycenter of A and B
              this->FacePoints->GetPoint((p + 1) % 4, coordsB);
              for (int j = 0; j < 3; ++j)
              {
                coordsC[j] = (B1 * coordsA[j] - A1 * coordsB[j]) / (B1 - A1);
              }
              this->EdgesA[i] = this->Points->InsertNextPoint(coordsC);
            } // if ( this->EdgesB[i] == -1 )

            // Update points
            ids[nPts] = this->EdgesA[i];
            ++nPts;
            if (indPairA)
            {
              pairA[1] = i;
            }
            else
            {
              pairA[0] = i;
            };
            ++indPairA;
          } // if ( A1 < 0. && A1 * B1 < 0. )
          if (A2 * B2 < 0.)
          {
            unsigned int i = EdgeIndices[orientation][offset][p];
            if (this->EdgesA[i] == -1)
            {
              // Compute barycenter of A and B
              this->FacePoints->GetPoint((p + 1) % 4, coordsB);
              for (int j = 0; j < 3; ++j)
              {
                coordsC[j] = (B2 * coordsA[j] - A2 * coordsB[j]) / (B2 - A2);
              }
              this->EdgesB[i] = this->Points->InsertNextPoint(coordsC);
            } // if ( this->EdgesA[i] == -1 )

            // Update points
            ids[nPts] = this->EdgesB[i];
            ++nPts;
            if (indPairB)
            {
              pairB[1] = i;
            }
            else
            {
              pairB[0] = i;
            };
            ++indPairB;
          } // if ( A2 * B2 < 0. )
          if (A1 > 0. && A1 * B1 < 0.)
          {
            unsigned int i = EdgeIndices[orientation][offset][p];
            if (this->EdgesA[i] == -1)
            {
              // Compute barycenter of A and B
              this->FacePoints->GetPoint((p + 1) % 4, coordsB);
              for (int j = 0; j < 3; ++j)
              {
                coordsC[j] = (B1 * coordsA[j] - A1 * coordsB[j]) / (B1 - A1);
              }
              this->EdgesA[i] = this->Points->InsertNextPoint(coordsC);
            } // if ( this->EdgesA[i] == -1 )

            // Update points
            ids[nPts] = this->EdgesA[i];
            ++nPts;
            if (indPairA)
            {
              pairA[1] = i;
            }
            else
            {
              pairA[0] = i;
            };
            ++indPairA;
          } // if ( A1 > 0. && A1 * B1 < 0. )
        }   // p

        // Insert pairs only if it makes sense
        if (indPairA == 2)
        {
          this->FacesA->InsertNextTypedTuple(pairA);
        }
        if (indPairB == 2)
        {
          this->FacesB->InsertNextTypedTuple(pairB);
        }
      } // else if( ! type )
      else if (type == -1.)
      {
        // Take positive values in A
        vtkIdType pair[] = { -1, -1 };
        unsigned int indPair = 0;

        // Loop over face vertices
        for (int p = 0; p < 4; ++p)
        {
          // Retrieve vertex coordinates
          this->FacePoints->GetPoint(p, coordsA);

          // Retrieve vertex scalars
          double A = this->FaceScalarsA->GetTuple1(p);
          double B = this->FaceScalarsA->GetTuple1((p + 1) % 4);

          // Add point when necessary
          if (create && A >= 0.)
          {
            ids[nPts] = this->Points->InsertNextPoint(coordsA);
            ++nPts;
          }
          if (A * B < 0.)
          {
            unsigned int i = EdgeIndices[orientation][offset][p];
            if (this->EdgesA[i] == -1)
            {
              // Compute barycenter of A and B
              this->FacePoints->GetPoint((p + 1) % 4, coordsB);
              for (int j = 0; j < 3; ++j)
              {
                coordsC[j] = (B * coordsA[j] - A * coordsB[j]) / (B - A);
              }
              this->EdgesA[i] = this->Points->InsertNextPoint(coordsC);
            } // if ( this->EdgesB[i] == -1 )

            // Update points
            ids[nPts] = this->EdgesA[i];
            ++nPts;
            if (indPair)
            {
              pair[1] = i;
            }
            else
            {
              pair[0] = i;
            };
            ++indPair;
          } // if ( A * B < 0. )
        }   // p

        // Insert pair only if it makes sense
        if (indPair == 2)
        {
          this->FacesA->InsertNextTypedTuple(pair);
        }
      } // else if ( type == -1. )
    }   // if ( type < 2 )
    else
    {
      // Create quadrangle vertices depending on orientation
      ids[0] = this->Points->InsertNextPoint(pt);
      pt[axis1] += size[axis1];
      ids[1] = this->Points->InsertNextPoint(pt);
      pt[axis2] += size[axis2];
      ids[2] = this->Points->InsertNextPoint(pt);
      pt[axis1] = origin[axis1];
      ids[3] = this->Points->InsertNextPoint(pt);
    } // else
  }   // if ( this->HasInterface )
  else
  {
    // Create quadrangle vertices depending on orientation
    ids[0] = this->Points->InsertNextPoint(pt);
    pt[axis1] += size[axis1];
    ids[1] = this->Points->InsertNextPoint(pt);
    pt[axis2] += size[axis2];
    ids[2] = this->Points->InsertNextPoint(pt);
    pt[axis1] = origin[axis1];
    ids[3] = this->Points->InsertNextPoint(pt);
  } // else

  // Insert next face if needed
  if (create)
  {
    // Create cell and corresponding ID
    vtkIdType outId = this->Cells->InsertNextCell(nPts, ids);

    // Copy face data from that of the cell from which it comes
    this->OutData->CopyData(this->InData, useId, outId);
  } // if ( create )
}
