//=========================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt 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 "vtkDiscreteModelFace.h"

#include "ModelVertexClassification.h"
#include "vtkBitArray.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkConnectivityFilter.h"
#include "vtkDiscreteModel.h"
#include "vtkDiscreteModelEdge.h"
#include "vtkDiscreteModelEntityGroup.h"
#include "vtkDiscreteModelVertex.h"
#include "vtkFeatureEdges.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkModelFaceUse.h"
#include "vtkModelItemIterator.h"
#include "vtkModelShellUse.h"
#include "vtkNew.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPointSet.h"
#include "vtkPolyData.h"
#include "vtkPolyDataNormals.h"
#include "vtkSerializer.h"
#include "vtkSmartPointer.h"
#include "vtkSplitEventData.h"

vtkDiscreteModelFace* vtkDiscreteModelFace::New()
{
  vtkObject* ret = vtkObjectFactory::CreateInstance("vtkDiscreteModelFace");
  if (ret)
  {
    return static_cast<vtkDiscreteModelFace*>(ret);
  }
  return new vtkDiscreteModelFace;
}

vtkDiscreteModelFace::vtkDiscreteModelFace()
{
}

vtkDiscreteModelFace::~vtkDiscreteModelFace()
{
}

vtkModelEntity* vtkDiscreteModelFace::GetThisModelEntity()
{
  return this;
}

bool vtkDiscreteModelFace::Split(
  double splitAngle, std::map<vtkIdType, FaceEdgeSplitInfo>& FaceSplitInfo)
{
  vtkObject* geometry = this->GetGeometry();
  vtkPolyData* poly = vtkPolyData::SafeDownCast(geometry);

  if (poly == 0)
  {
    if (geometry)
    {
      // we are on the server...
      return 0;
    }
    // BoundaryRep has not been set -> return error
    return 0;
  }

  // on server, go ahead and perform the split
  vtkNew<vtkPolyDataNormals> pdNormals;
  pdNormals->SetFeatureAngle(splitAngle);
  pdNormals->SetInputData(0, poly);

  vtkNew<vtkConnectivityFilter> pdConnectivity;
  pdConnectivity->SetInputConnection(pdNormals->GetOutputPort());
  pdConnectivity->SetColorRegions(1);
  pdConnectivity->Update();

  vtkIdTypeArray* newFaceTags =
    vtkIdTypeArray::SafeDownCast(vtkDataSet::SafeDownCast(pdConnectivity->GetOutputDataObject(0))
                                   ->GetCellData()
                                   ->GetArray("RegionId"));

  // filter the output
  // actually, newModelFaces also contains the old/source model face

  // the ids in newModelFaces are with respect to the master poly data
  std::map<vtkIdType, vtkSmartPointer<vtkIdList> > newModelFaces;

  vtkIdTypeArray* masterGeometryCellIndex = this->GetReverseClassificationArray();
  for (vtkIdType i = 0; i < newFaceTags->GetNumberOfTuples(); i++)
  {
    vtkIdType tag = newFaceTags->GetValue(i);
    if (newModelFaces.find(tag) == newModelFaces.end())
    {
      newModelFaces[tag] = vtkSmartPointer<vtkIdList>::New();
    }
    newModelFaces[tag]->InsertNextId(masterGeometryCellIndex->GetValue(i));
  }

  if (newModelFaces.size() > 1)
  {
    vtkSplitEventData* splitEventData = vtkSplitEventData::New();
    splitEventData->SetSourceEntity(this);
    vtkIdList* newFaceIds = vtkIdList::New();
    newFaceIds->SetNumberOfIds(newModelFaces.size() - 1); // the first one is the original
    std::map<vtkIdType, vtkSmartPointer<vtkIdList> >::iterator mit = newModelFaces.begin();
    mit++; // skip the first model face as it is the old/source model face

    bool saveExistingLoop = false;
    for (vtkIdType i = 0; mit != newModelFaces.end(); mit++, i++)
    {
      FaceEdgeSplitInfo splitInfo;
      // if this is the last new face, save the loopinfo for existing face
      saveExistingLoop = (i == (newFaceIds->GetNumberOfIds() - 1));
      vtkDiscreteModelFace* face =
        this->BuildFromExistingModelFace(mit->second, splitInfo, saveExistingLoop);
      vtkIdType newId = face->GetUniquePersistentId();
      FaceSplitInfo.insert(std::make_pair(newId, splitInfo));
      newFaceIds->SetId(i, newId);
    }
    splitEventData->SetCreatedModelEntityIds(newFaceIds);
    newFaceIds->Delete();
    this->GetModel()->InvokeModelGeometricEntityEvent(ModelGeometricEntitySplit, splitEventData);
    splitEventData->Delete();
  }

  return 1;
}

vtkDiscreteModelFace* vtkDiscreteModelFace::BuildFromExistingModelFace(
  vtkIdList* masterCellIds, FaceEdgeSplitInfo& splitInfo, bool saveLoopForExistingFace)
{
  vtkDiscreteModel* thisModel = vtkDiscreteModel::SafeDownCast(this->GetModel());
  bool blockEvent = thisModel->GetBlockModelGeometricEntityEvent();
  thisModel->SetBlockModelGeometricEntityEvent(true);

  vtkDiscreteModelFace* newModelFace =
    vtkDiscreteModelFace::SafeDownCast(thisModel->BuildModelFace(0, 0, 0));
  thisModel->SetBlockModelGeometricEntityEvent(blockEvent);

  newModelFace->AddCellsToGeometry(masterCellIds);

  // Handle splits needed for existing edges
  if (this->GetNumberOfModelEdges())
  {
    this->SplitEdges(newModelFace, splitInfo);
  }

  // put in the adjacencies for the model face use
  for (int i = 0; i < 2; i++)
  {
    vtkModelFaceUse* sourceFaceUse = this->GetModelFaceUse(i);
    vtkModelShellUse* sourceShellUse = sourceFaceUse->GetModelShellUse();
    vtkModelFaceUse* targetFaceUse = newModelFace->GetModelFaceUse(i);
    if (sourceShellUse)
    {
      sourceShellUse->AddModelFaceUse(targetFaceUse);
    }
  }

  // new model face will belong to same groups as old model face
  newModelFace->CopyModelEntityGroups(this);
  thisModel->InvokeModelGeometricEntityEvent(ModelGeometricEntityCreated, newModelFace);
  for (int i = 0; i < 2; i++)
  {
    if (vtkModelRegion* region = newModelFace->GetModelRegion(i))
    {
      thisModel->InvokeModelGeometricEntityEvent(ModelGeometricEntityBoundaryModified, region);
    }
  }

  // update loops and edges info if needed
  if (this->GetNumberOfModelEdges())
  {
    // destroy the existing loop uses
    this->DestroyLoopUses();
    newModelFace->BuildEdges(true, splitInfo);
    this->BuildEdges(true, splitInfo, saveLoopForExistingFace);
  }

  return newModelFace;
}

void vtkDiscreteModelFace::GetAllPointIds(vtkIdList* ptsList)
{
  vtkBitArray* Points = vtkBitArray::New();
  this->GatherAllPointIdsMask(Points);
  for (vtkIdType i = 0; i < Points->GetNumberOfTuples(); i++)
  {
    if (Points->GetValue(i) != 0)
    {
      ptsList->InsertNextId(i);
    }
  }
  Points->Delete();
}

void vtkDiscreteModelFace::GetInteriorPointIds(vtkIdList* ptsList)
{
  vtkBitArray* allModelFacePoints = vtkBitArray::New();
  this->GatherAllPointIdsMask(allModelFacePoints);
  vtkBitArray* boundaryModelFacePoints = vtkBitArray::New();
  this->GatherBoundaryPointIdsMask(boundaryModelFacePoints);
  for (vtkIdType i = 0; i < allModelFacePoints->GetNumberOfTuples(); i++)
  {
    if (allModelFacePoints->GetValue(i) != 0 && boundaryModelFacePoints->GetValue(i) == 0)
    {
      ptsList->InsertNextId(i);
    }
  }
  allModelFacePoints->Delete();
  boundaryModelFacePoints->Delete();
}

void vtkDiscreteModelFace::GetBoundaryPointIds(vtkIdList* ptsList)
{
  vtkBitArray* points = vtkBitArray::New();
  this->GatherBoundaryPointIdsMask(points);
  for (vtkIdType i = 0; i < points->GetNumberOfTuples(); i++)
  {
    if (points->GetValue(i) != 0)
    {
      ptsList->InsertNextId(i);
    }
  }
  points->Delete();
}

void vtkDiscreteModelFace::GatherAllPointIdsMask(vtkBitArray* points)
{
  vtkPointSet* grid = vtkPointSet::SafeDownCast(this->GetGeometry());
  if (!grid)
  {
    vtkWarningMacro("Cannot access model face's grid.");
    return;
  }
  // Since not all points in Grid's points are attached to cells, we
  // iterate over the cells and get their points to set them in Points.
  points->SetNumberOfComponents(1);
  points->SetNumberOfTuples(grid->GetNumberOfPoints());
  for (vtkIdType i = 0; i < grid->GetNumberOfPoints(); i++)
  {
    points->SetValue(i, 0);
  }
  vtkIdType numberOfCells = grid->GetNumberOfCells();
  vtkIdList* cellPoints = vtkIdList::New();
  for (vtkIdType i = 0; i < numberOfCells; i++)
  {
    grid->GetCellPoints(i, cellPoints);
    for (vtkIdType j = 0; j < cellPoints->GetNumberOfIds(); j++)
    {
      vtkIdType id = cellPoints->GetId(j);
      points->SetValue(id, 1);
    }
  }
  cellPoints->Delete();
}

void vtkDiscreteModelFace::GatherBoundaryPointIdsMask(vtkBitArray* points)
{
  vtkPointSet* grid = vtkPointSet::SafeDownCast(this->GetGeometry());
  // add in an array of the original point ids to grid
  vtkIdTypeArray* originalPointIds = vtkIdTypeArray::New();
  originalPointIds->SetNumberOfComponents(1);
  originalPointIds->SetNumberOfTuples(grid->GetNumberOfPoints());
  for (vtkIdType i = 0; i < grid->GetNumberOfPoints(); i++)
  {
    originalPointIds->SetValue(i, i);
  }
  const char arrayName[] = "vtkModelNodalGroupPointIdArray";
  originalPointIds->SetName(arrayName);
  grid->GetPointData()->AddArray(originalPointIds);
  originalPointIds->Delete();

  vtkNew<vtkFeatureEdges> featureEdges;
  featureEdges->BoundaryEdgesOn();
  featureEdges->NonManifoldEdgesOff();
  featureEdges->ManifoldEdgesOff();
  featureEdges->FeatureEdgesOff();
  featureEdges->SetInputData(grid);
  featureEdges->Update();
  vtkPointSet* gridEdges = vtkPointSet::SafeDownCast(featureEdges->GetOutput(0));

  // Since not all points in GridEdges's points are attached to cells, we
  // iterate over the cells and get their points to set them in Points.
  originalPointIds = vtkIdTypeArray::SafeDownCast(gridEdges->GetPointData()->GetArray(arrayName));

  points->SetNumberOfComponents(1);
  points->SetNumberOfTuples(grid->GetNumberOfPoints());
  for (vtkIdType i = 0; i < grid->GetNumberOfPoints(); i++)
  {
    points->SetValue(i, 0);
  }
  vtkIdType numberOfCells = gridEdges->GetNumberOfCells();
  vtkNew<vtkIdList> cellPoints;
  for (vtkIdType i = 0; i < numberOfCells; i++)
  {
    gridEdges->GetCellPoints(i, cellPoints.GetPointer());
    for (vtkIdType j = 0; j < cellPoints->GetNumberOfIds(); j++)
    {
      vtkIdType id = cellPoints->GetId(j);
      vtkIdType masterGridId = originalPointIds->GetValue(id);
      points->SetValue(masterGridId, 1);
    }
  }
  grid->GetPointData()->RemoveArray(arrayName);
}
void vtkDiscreteModelFace::ExtractEdges(vtkPolyData* result)
{
  vtkPolyData* grid = vtkPolyData::SafeDownCast(this->GetGeometry());
  result->Reset();
  // add in an array for storing the cell ID that "owns" the edge
  vtkNew<vtkIdTypeArray> cellIds;
  cellIds->SetNumberOfComponents(1);
  cellIds->Allocate(grid->GetNumberOfCells() * 0.1);
  result->SetPoints(grid->GetPoints());
  const char arrayName[] = "OriginalFacetCellID";
  cellIds->SetName(arrayName);
  result->GetCellData()->AddArray(cellIds.GetPointer());

  // Setup the resulting edges
  vtkNew<vtkCellArray> newLines;
  newLines->Allocate(grid->GetNumberOfCells() * 0.1);
  result->SetLines(newLines.GetPointer());

  // Get the original cell IDs w/r to the mesh
  vtkIdTypeArray* masterGeometryCellIndex = this->GetReverseClassificationArray();

  // Iterate over all of the cells of the face and find any boundary edge
  vtkIdType cellId, npts, *pts, i, linePts[2];
  vtkNew<vtkIdList> neighbors;
  vtkCellArray* facePolys = grid->GetPolys();
  neighbors->Allocate(VTK_CELL_SIZE);
  grid->BuildLinks();
  for (cellId = 0, facePolys->InitTraversal(); facePolys->GetNextCell(npts, pts); cellId++)
  {
    if (npts < 3)
    {
      continue; // Degenerate Cell!
    }
    // Set the "first line point to be the last one in the poly
    linePts[0] = pts[npts - 1];
    for (i = 0; i < npts; i++)
    {
      grid->GetCellEdgeNeighbors(cellId, linePts[0], pts[i], neighbors.GetPointer());
      if (!(neighbors->GetNumberOfIds()))
      {
        linePts[1] = pts[i];
        newLines->InsertNextCell(2, linePts);
        cellIds->InsertNextValue(masterGeometryCellIndex->GetValue(cellId));
      }
      linePts[0] = pts[i]; // Set this to be the new start point
    }
  }
}

void vtkDiscreteModelFace::BuildEdges(
  bool showEdge, FaceEdgeSplitInfo& splitInfo, bool saveLoopInfo)
{
  vtkModelEdge* gedge;
  vtkNew<vtkPolyData> edges;
  std::vector<LoopInfo> loops;
  std::vector<bool> visited;
  NewModelEdgeInfo newEdgesInfo;
  int i, nlines, dummyInt, nLoops, j, nEdges;
  size_t currentLoop;
  vtkIdType gid;
  vtkIdType idLoops, idEdges, idOrient, fid = this->GetUniquePersistentId();

  vtkCellArray* lines;

  this->ExtractEdges(edges.GetPointer());
  vtkIdTypeArray* facetIds =
    dynamic_cast<vtkIdTypeArray*>(edges->GetCellData()->GetArray("OriginalFacetCellID", dummyInt));
  if (!facetIds)
  {
    std::cout << "Could not find Facet Mapping!!\n";
    return;
  }

  edges->BuildLinks(); // We need to "walk" the loops

  lines = edges->GetLines();
  nlines = lines->GetNumberOfCells();
  if (!nlines)
  {
    std::cout << "\tFace " << this->GetUniquePersistentId() << "has no edges!\n";
    return;
  }
  visited.assign(nlines, false);
  loops.clear();
  newEdgesInfo.Reset();
  for (i = 0; i < nlines; i++)
  {
    if (visited[i])
    {
      continue;
    }
    //Create a new loop
    currentLoop = loops.size();
    loops.resize(currentLoop + 1);

    this->WalkLoop(i, edges.GetPointer(), visited, facetIds, newEdgesInfo, loops[currentLoop]);
  }

  // OK so we have now processed all of the loops we now need to create
  // any new model edges
  std::map<int, vtkDiscreteModelEdge*> newModelEdges;
  this->CreateModelEdges(newEdgesInfo, newModelEdges, showEdge, splitInfo);

  // Now we are ready to add the model edges to the model face
  nLoops = static_cast<int>(loops.size());
  if (!nLoops)
  {
    // This face has no loops - like a sphere
    return;
  }

  if (saveLoopInfo)
  {
    // The map of <faceId, nloops, nEdges, (gedges[n]...), (orientations[n]...)>
    splitInfo.FaceEdgeLoopIDs->InsertNextTypedTuple(&fid);
    idLoops = nLoops;
    splitInfo.FaceEdgeLoopIDs->InsertNextTypedTuple(&idLoops);
  }
  std::vector<int> orientations;
  std::vector<vtkModelEdge*> gedges;
  for (i = 0; i < nLoops; i++)
  {
    // std::cout << "\tLoop " << i << ":";
    nEdges = static_cast<int>(loops[i].loop.size());
    orientations.resize(nEdges);
    gedges.resize(nEdges);
    if (saveLoopInfo)
    {
      idEdges = nEdges;
      splitInfo.FaceEdgeLoopIDs->InsertNextTypedTuple(&idEdges);
    }

    for (j = 0; j < nEdges; j++)
    {
      // Is this a new model edge or a prexisting one
      gid = loops[i].loop[j].first;
      if (gid < 0)
      {
        gedge = newModelEdges[gid];
      }
      else
      {
        gedge =
          dynamic_cast<vtkModelEdge*>(this->GetModel()->GetModelEntity(vtkModelEdgeType, gid));
      }
      if (!gedge)
      {
        std::cout << "Could not find model edge: " << gid << "\n";
        continue;
      }
      if (loops[i].loop[j].second)
      {
        orientations[j] = 1;
        // std::cout << gedge->GetUniquePersistentId() << " ";
      }
      else
      {
        orientations[j] = 0;
        // std::cout << "~" << gedge->GetUniquePersistentId() << " ";
      }
      gedges[j] = gedge;

      if (saveLoopInfo)
      {
        gid = gedge->GetUniquePersistentId();
        splitInfo.FaceEdgeLoopIDs->InsertNextTypedTuple(&gid);
        idOrient = orientations[j];
        splitInfo.FaceEdgeLoopIDs->InsertNextTypedTuple(&idOrient);
      }
    }
    // std::cout << "\n";
    this->AddLoop(nEdges, &gedges.front(), &orientations.front());
  }
}

void vtkDiscreteModelFace::WalkLoop(vtkIdType startingEdge, vtkPolyData* edges,
  std::vector<bool>& visited, vtkIdTypeArray* facetIds, NewModelEdgeInfo& newEdgesInfo,
  LoopInfo& loopInfo)
{
  vtkIdType currentEdge = startingEdge;
  vtkNew<vtkIdList> cellIds, pointIds;
  vtkIdType currentPoint, nextPoint;
  vtkIdType gedge;
  std::string currentFaceInfo;
  vtkDiscreteModel* thisModel = vtkDiscreteModel::SafeDownCast(this->GetModel());
  const DiscreteMesh& mesh = thisModel->GetMesh();
  vtkDiscreteModel::ClassificationType& classificationInfo = thisModel->GetMeshClassification();
  vtkDiscreteModelEdge* dedge;
  // Setup the newEdgesInfo to tag the first mesh that needs to be
  // created and to force a new model edge to be "defined" when that
  // happens
  newEdgesInfo.ClearCurrentFaceInfo();
  newEdgesInfo.ClearTaggedFaceInfo();
  // Walk the loop
  currentPoint = -1; // Indicates the point has not been set

  while (!visited[currentEdge])
  {
    visited[currentEdge] = true;
    // Get the points of the edge
    edges->GetCellPoints(currentEdge, pointIds.GetPointer());
    if (pointIds->GetNumberOfIds() != 2)
    {
      std::cout << "Found incorrect number of points: " << pointIds->GetNumberOfIds() << "\n";
      break;
    }
    if (currentPoint == -1)
    {
      currentPoint = pointIds->GetId(0);
      nextPoint = pointIds->GetId(1);
    }
    else if (currentPoint != pointIds->GetId(0))
    {
      if (currentPoint == pointIds->GetId(1))
      {
        std::cout << "Found flipped Edge!\n";
        nextPoint = pointIds->GetId(0);
      }
      else
      {
        std::cout << "Found Disconnected Edge!\n";
      }
    }
    else
    {
      nextPoint = pointIds->GetId(1);
    }
    // Does the edge exist?
    bool createdEdge, orientation;
    vtkIdType medge = mesh.AddEdgeIfNotExisting(currentPoint, nextPoint, orientation, createdEdge);
    // std::cout << "Edge: " << currentEdge << " Verts: "
    //           << currentPoint << ", " << nextPoint << " MEdge: "
    //           << medge << " Created: " << createdEdge << "\n";
    if (createdEdge)
    {
      //Need to get the face information bounding the edge
      currentFaceInfo =
        this->EncodeModelFaces(facetIds->GetValue(currentEdge), currentPoint, nextPoint);
      gedge = newEdgesInfo.InsertMeshEdge(medge, currentFaceInfo);
      loopInfo.InsertModelEdge(gedge, true);
    }
    else
    {
      // The model edge exists so the next time we find an edge
      // that doesn't exist it must be on a new model edge
      newEdgesInfo.ClearCurrentFaceInfo();
      // Need to get the model edge classified on the mesh edge
      dedge = dynamic_cast<vtkDiscreteModelEdge*>(classificationInfo.GetEntity(medge));

      gedge = dedge->GetUniquePersistentId();
      loopInfo.InsertModelEdge(gedge, orientation);
    }
    // Lets get the next edge in the loop
    currentPoint = nextPoint;
    edges->GetPointCells(currentPoint, cellIds.GetPointer());
    // Should be 2 lines comming into the point
    if (cellIds->GetNumberOfIds() < 2)
    {
      std::cout << "Found unconnected edge!\n";
    }
    else if (cellIds->GetNumberOfIds() > 2)
    {
      std::cout << "Found non-manifold edge!\n";
      // Is there only one edge left to visit?
      vtkIdType numEdgesToVisit = 0, lastEdgeToBeVisited = 0, j;
      vtkIdType testMEdge, testEdge;
      for (j = 0; j < cellIds->GetNumberOfIds(); j++)
      {
        if (!visited[cellIds->GetId(j)])
        {
          ++numEdgesToVisit;
          lastEdgeToBeVisited = cellIds->GetId(j);
        }
      }
      if (!numEdgesToVisit)
      {
        continue; // We have reached the end of the loop
      }
      if (numEdgesToVisit == 1)
      {
        // There is only one way to go - take it!
        currentEdge = lastEdgeToBeVisited;
        continue;
      }
      // OK we need to determine how to go - if we created a new model
      // edge then we should look for an edge with the same face classification
      // else we should see if we have an edge on the same model edge
      if (createdEdge)
      {
        // For each edge see if we can find the same bounding faces
        std::string nextFaceInfo;
        for (j = 0; j < cellIds->GetNumberOfIds(); j++)
        {
          testEdge = cellIds->GetId(j);
          // Skip if we have already visited the edge or if the edge is on
          // an existing model edge
          if (visited[testEdge])
          {
            continue;
          }
          edges->GetCellPoints(testEdge, pointIds.GetPointer());
          if (mesh.EdgeExists(pointIds->GetId(0), pointIds->GetId(1), testMEdge))
          {
            continue;
          }
          nextFaceInfo = this->EncodeModelFaces(
            facetIds->GetValue(testEdge), pointIds->GetId(0), pointIds->GetId(1));
          if (currentFaceInfo == nextFaceInfo)
          {
            // Found next edge!
            currentEdge = testEdge;
            break;
          }
        }
      }
      else
      {
        // This is the reverse of the previous method - here we want to follow
        // only existing model edges
        for (j = 0; j < cellIds->GetNumberOfIds(); j++)
        {
          testEdge = cellIds->GetId(j);
          // Skip if we have already visited the edge or if the edge is not on
          // an existing model edge
          if (visited[testEdge])
          {
            continue;
          }
          edges->GetCellPoints(testEdge, pointIds.GetPointer());
          if (!mesh.EdgeExists(pointIds->GetId(0), pointIds->GetId(1), testMEdge))
          {
            continue;
          }
          // Need to get the model edge classified on the mesh edge
          dedge = dynamic_cast<vtkDiscreteModelEdge*>(classificationInfo.GetEntity(testMEdge));
          if (dedge->GetUniquePersistentId() == gedge)
          {
            // Found next edge!
            currentEdge = testEdge;
            break;
          }
        }
      }
    }
    else
    {
      if (cellIds->GetId(0) == currentEdge)
      {
        currentEdge = cellIds->GetId(1);
      }
      else if (cellIds->GetId(1) != currentEdge)
      {
        std::cout << "Invalid Cell Structure!!\n";
      }
      else
      {
        currentEdge = cellIds->GetId(0);
      }
    }
  }
  // Do we need to fix the first edge in the loop?  This would occur if
  // we didn't  not start at a "model vertex".  The condition we need to
  // check for is if the first and last edges are either both "new" or both
  // existing
  size_t numLoops = loopInfo.loop.size();
  if (numLoops == 1)
  {
    // There is only a single model edge so there is no problem
    return;
  }
  // Does the loop start and end with the same model edge?  If so remove one
  if (loopInfo.loop[0].first == loopInfo.loop[numLoops - 1].first)
  {
    loopInfo.loop.erase(loopInfo.loop.begin());
  }
  else if ((loopInfo.loop[0].first < 0) && (loopInfo.loop[numLoops - 1].first < 0))
  {
    // OK so we know that the first and last edges of the loop need to be
    // created.  If they are bounding the same model faces they need to be
    // combined.
    if (newEdgesInfo.CheckLastModelEdge(gedge))
    {
      loopInfo.RemoveModelEdge(gedge);
    }
  }
}

std::string vtkDiscreteModelFace::EncodeModelFaces(vtkIdType facetId, vtkIdType v0, vtkIdType v1)
{
  vtkDiscreteModel* thisModel = vtkDiscreteModel::SafeDownCast(this->GetModel());

  vtkNew<vtkIdList> cellIds;
  vtkDiscreteModel::ClassificationType& classificationInfo = thisModel->GetMeshClassification();
  int j, nids;
  vtkDiscreteModelFace* dface;
  vtkIdType gface;
  thisModel->GetMesh().GetCellEdgeNeighbors(facetId, v0, v1, cellIds.GetPointer());
  // Convert the facet Ids to model Ids
  nids = cellIds->GetNumberOfIds();
  if (nids == 0)
  {
    cellIds->InsertNextId(this->GetUniquePersistentId());
  }
  else
  {
    for (j = 0; j < nids; j++)
    {
      dface = dynamic_cast<vtkDiscreteModelFace*>(classificationInfo.GetEntity(cellIds->GetId(j)));
      gface = dface->GetUniquePersistentId();
      cellIds->SetId(j, gface);
    }
  }

  return NewModelEdgeInfo::to_key(cellIds.GetPointer());
}

void vtkDiscreteModelFace::CreateModelEdges(NewModelEdgeInfo& newEdgesInfo,
  std::map<int, vtkDiscreteModelEdge*>& newEdges, bool bShow, FaceEdgeSplitInfo& splitInfo)
{
  // Lets walk the new edge information and create the appropriate model
  // topology
  vtkDiscreteModel* thisModel = vtkDiscreteModel::SafeDownCast(this->GetModel());

  newEdges.clear();
  vtkIdType currentModelEdgeId = 0, numEdges;
  vtkNew<vtkIdList> edgeCells, edgePnts;
  const DiscreteMesh& mesh = thisModel->GetMesh();
  vtkDiscreteModel::ClassificationType& classificationInfo = thisModel->GetMeshClassification();
  ModelVertexClassification vertexInfo(thisModel);
  vtkDiscreteModelVertex *v0 = NULL, *v1;
  vtkDiscreteModelEdge* gedge;
  std::vector<std::pair<vtkIdType, vtkIdType> >& info = newEdgesInfo.info;
  std::size_t i, numMeshEdges = info.size();

  //First we need to extend the classification information to include the
  // mesh edges we created
  numEdges = classificationInfo.size(vtkDiscreteModel::ClassificationType::EDGE_DATA);
  classificationInfo.resize(
    numEdges + numMeshEdges, vtkDiscreteModel::ClassificationType::EDGE_DATA);

  for (i = 0; i < numMeshEdges; i++)
  {
    //Are we on the same model edge?
    if (info[i].second == currentModelEdgeId)
    {
      // Just add the cell to the list
      edgeCells->InsertNextId(info[i].first);
      continue;
    }

    // Have we found any mesh edges to create a model edge yet?
    // Note that the first time into the loop this will not be true
    numEdges = edgeCells->GetNumberOfIds();
    if (numEdges)
    {
      // We need to create v1 for this model edge by looking at the
      // last mesh edge added
      mesh.GetCellPointIds(edgeCells->GetId(numEdges - 1), edgePnts.GetPointer());

      vtkIdType modelId;
      v1 = vertexInfo.AddModelVertex(edgePnts->GetId(1), true, modelId);
      if (v0 && v1)
      {
        gedge = dynamic_cast<vtkDiscreteModelEdge*>(thisModel->BuildModelEdge(v0, v1));
        // Don't show edges by default
        gedge->SetVisibility(bShow);
        v0->SetVisibility(bShow);
        v1->SetVisibility(bShow);

        // The map of <NewEdgeId, VertId1, VertId2>
        vtkIdType newEdgeVV[3] = { gedge->GetUniquePersistentId(), v0->GetUniquePersistentId(),
          v1->GetUniquePersistentId() };
        splitInfo.CreatedModelEdgeVertIDs->InsertNextTypedTuple(newEdgeVV);

        gedge->AddCellsToGeometry(edgeCells.GetPointer());
        newEdges[currentModelEdgeId] = gedge;
        // std::cout << "Create Model Edge: " << gedge->GetUniquePersistentId()
        //           << " (ID = " << currentModelEdgeId
        //           << ") From Vertices ("
        //           << v0->GetUniquePersistentId() << ","
        //           << v1->GetUniquePersistentId() << ")\n";
      }
      else
      {
        std::cout << "Could not build model edge!!\n";
      }
      edgeCells->Reset();
    }

    currentModelEdgeId = info[i].second;
    // Get the first vertex of the edge
    mesh.GetCellPointIds(info[i].first, edgePnts.GetPointer());

    vtkIdType modelId;
    v0 = vertexInfo.AddModelVertex(edgePnts->GetId(0), true, modelId);
    // Add the mesh edge to the list
    edgeCells->InsertNextId(info[i].first);
  }

  // Finally we need to process the last model edge
  // Have we found any mesh edges to create a model edge yet?
  // Note that the first time into the loop this will not be true
  numEdges = edgeCells->GetNumberOfIds();
  if (numEdges)
  {
    // We need to create v1 for this model edge by looking at the
    // last mesh edge added
    mesh.GetCellPointIds(edgeCells->GetId(numEdges - 1), edgePnts.GetPointer());

    vtkIdType modelId;
    v1 = vertexInfo.AddModelVertex(edgePnts->GetId(1), true, modelId);
    if (v0 && v1)
    {
      gedge = dynamic_cast<vtkDiscreteModelEdge*>(thisModel->BuildModelEdge(v0, v1));
      // Don't show edges and vertex by default
      gedge->SetVisibility(bShow);
      v0->SetVisibility(bShow);
      v1->SetVisibility(bShow);

      // The map of <NewEdgeId, VertId1, VertId2>
      vtkIdType newEdgeVV[3] = { gedge->GetUniquePersistentId(), v0->GetUniquePersistentId(),
        v1->GetUniquePersistentId() };
      splitInfo.CreatedModelEdgeVertIDs->InsertNextTypedTuple(newEdgeVV);

      gedge->AddCellsToGeometry(edgeCells.GetPointer());
      newEdges[currentModelEdgeId] = gedge;
      // std::cout << "Create Model Edge: " << gedge->GetUniquePersistentId()
      //           << " (ID = " << currentModelEdgeId
      //           << ") From Vertices ("
      //           << v0->GetUniquePersistentId() << ","
      //           << v1->GetUniquePersistentId() << ")\n";
    }
    else
    {
      std::cout << "Could not build model edge!!\n";
    }
  }
}

void vtkDiscreteModelFace::SplitEdges(
  vtkDiscreteModelFace* newModelFace, FaceEdgeSplitInfo& splitInfo)
{
  vtkDiscreteModel* thisModel = vtkDiscreteModel::SafeDownCast(this->GetModel());

  // If we have a new face breaking existing edges,
  // that edge needs to be splitted
  const DiscreteMesh& mesh = thisModel->GetMesh();
  vtkDiscreteModel::ClassificationType& classificationInfo = thisModel->GetMeshClassification();
  vtkDiscreteModelFace* dface;

  std::vector<vtkModelEdge*> adjEdges;
  this->GetModelEdges(adjEdges);
  for (std::vector<vtkModelEdge*>::iterator eit = adjEdges.begin(); eit != adjEdges.end(); ++eit)
  {
    vtkDiscreteModelEdge* splitEdge = vtkDiscreteModelEdge::SafeDownCast(*eit);
    vtkPolyData* grid = vtkPolyData::SafeDownCast(splitEdge->GetGeometry());
    vtkCellArray* lines = grid->GetLines();
    vtkIdType npts, *pts;
    int nids;
    vtkNew<vtkIdList> cellIds;
    bool newFace1, newFace2;
    std::vector<vtkIdType> NewVertexList;
    for (lines->InitTraversal(); lines->GetNextCell(npts, pts);)
    {
      // process first point
      mesh.GetCellsUsingPoint(pts[0], cellIds.GetPointer());
      // Convert the facet Ids to model Ids
      nids = cellIds->GetNumberOfIds();
      newFace1 = false;
      for (int j = 0; j < nids; j++)
      {
        dface =
          dynamic_cast<vtkDiscreteModelFace*>(classificationInfo.GetEntity(cellIds->GetId(j)));
        newFace1 = (dface == newModelFace);
        if (newFace1)
        {
          break;
        }
      }
      // process second point
      mesh.GetCellsUsingPoint(pts[1], cellIds.GetPointer());
      // Convert the facet Ids to model Ids
      nids = cellIds->GetNumberOfIds();
      newFace2 = false;
      for (int j = 0; j < nids; j++)
      {
        dface =
          dynamic_cast<vtkDiscreteModelFace*>(classificationInfo.GetEntity(cellIds->GetId(j)));
        newFace2 = (dface == newModelFace);
        if (newFace2)
        {
          break;
        }
      }

      // if only one point in the cell (xor) is part of the new face,
      // we should split at the point that's part of the new face.
      if (newFace1 ^ newFace2)
      {
        NewVertexList.push_back(newFace1 ? pts[0] : pts[1]);
      }
    }
    for (std::vector<vtkIdType>::iterator pidIt = NewVertexList.begin();
         splitEdge && pidIt != NewVertexList.end(); ++pidIt)
    {
      vtkIdType newEdgeId = -1, newVertexId = -1;
      if (splitEdge->Split(*pidIt, newVertexId, newEdgeId))
      {
        vtkIdType splitEdgeVEdge[3] = { splitEdge->GetUniquePersistentId(), newVertexId,
          newEdgeId };
        splitInfo.SplitEdgeVertIds->InsertNextTypedTuple(splitEdgeVEdge);
        splitEdge = vtkDiscreteModelEdge::SafeDownCast(
          thisModel->GetModelEntity(vtkModelEdgeType, newEdgeId));
        // The map of <OldEdgeId, NewVertId, NewEdgId>
      }
      else
      {
        splitEdge = NULL;
      }
    }
  }
}

void vtkDiscreteModelFace::Serialize(vtkSerializer* ser)
{
  this->Superclass::Serialize(ser);
}

bool vtkDiscreteModelFace::Destroy()
{
  if (this->Superclass::Destroy())
  {
    this->RemoveAllAssociations(vtkDiscreteModelEntityGroupType);
    this->Modified();
    return true;
  }
  return false;
}

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