qgstopologicalmesh.cpp

Go to the documentation of this file.

/***************************************************************************
  qgstopologicalmesh.cpp - QgsTopologicalMesh
 
 ---------------------
 begin                : 18.6.2021
 copyright            : (C) 2021 by Vincent Cloarec
 email                : vcloarec at gmail dot com
 ***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgis.h"
#include "qgstopologicalmesh.h"
#include "qgsmesheditor.h"
#include "qgsmessagelog.h"
#include "qgsgeometryutils.h"
#include "qgscircle.h"
 
#include <poly2tri.h>
#include <QSet>
#include <QQueue>
 
/*static*/ int QgsTopologicalMesh::vertexPositionInFace( const QgsMesh &mesh, int vertexIndex, int faceIndex )
{
  if ( faceIndex < 0 || faceIndex >= mesh.faceCount() )
    return -1;
 
  return QgsTopologicalMesh::vertexPositionInFace( vertexIndex, mesh.face( faceIndex ) );
}
 
 
static double crossProduct( const QgsMeshVertex &vc, const QgsMeshVertex &v1, const QgsMeshVertex &v2 )
{
  double ux1 = v1.x() - vc.x();
  double uy1 = v1.y() - vc.y();
  double vx1 = v2.x() - vc.x();
  double vy1 = v2.y() - vc.y();
 
  return ux1 * vy1 - uy1 * vx1;
}
 
static double crossProduct( int centralVertex, int vertex1, int vertex2, const QgsMesh &mesh )
{
  QgsMeshVertex vc = mesh.vertices.at( centralVertex );
  QgsMeshVertex v1 = mesh.vertices.at( vertex1 );
  QgsMeshVertex v2 = mesh.vertices.at( vertex2 );
 
  return crossProduct( vc, v1, v2 );
}
 
 
QgsMeshVertexCirculator::QgsMeshVertexCirculator( const QgsTopologicalMesh &topologicalMesh, int vertexIndex )
  : mFaces( topologicalMesh.mMesh->faces )
  , mFacesNeighborhood( topologicalMesh.mFacesNeighborhood )
  ,  mVertexIndex( vertexIndex )
{
  if ( vertexIndex >= 0 && vertexIndex < topologicalMesh.mMesh->vertexCount() )
  {
    mCurrentFace = topologicalMesh.mVertexToFace[vertexIndex];
    mIsValid = QgsTopologicalMesh::vertexPositionInFace( *topologicalMesh.mesh(), vertexIndex, mCurrentFace ) != -1;
  }
  else
  {
    mIsValid = false;
  }
 
  if ( mIsValid )
    mLastValidFace = mCurrentFace;
}
 
QgsMeshVertexCirculator::QgsMeshVertexCirculator( const QgsTopologicalMesh::TopologicalFaces &topologicalFaces, int faceIndex, int vertexIndex )
  : mFaces( topologicalFaces.mFaces )
  , mFacesNeighborhood( topologicalFaces.mFacesNeighborhood )
  , mVertexIndex( vertexIndex )
{
  const QgsMeshFace &face = topologicalFaces.mFaces.at( faceIndex );
  mIsValid = QgsTopologicalMesh::vertexPositionInFace( vertexIndex, face ) != -1;
 
  mCurrentFace = faceIndex;
  mLastValidFace = mCurrentFace;
}
 
QgsMeshVertexCirculator::QgsMeshVertexCirculator( const QgsTopologicalMesh::TopologicalFaces &topologicalFaces, int vertexIndex )
  : mFaces( topologicalFaces.mFaces )
  , mFacesNeighborhood( topologicalFaces.mFacesNeighborhood )
  , mVertexIndex( vertexIndex )
{
  if ( topologicalFaces.mVerticesToFace.contains( vertexIndex ) )
    mCurrentFace = topologicalFaces.mVerticesToFace.values( vertexIndex ).first();
  mLastValidFace = mCurrentFace;
  mIsValid = mCurrentFace != -1;
}
 
int QgsMeshVertexCirculator::turnCounterClockwise() const
{
  if ( mCurrentFace == -1 )
    mCurrentFace = mLastValidFace;
  else
  {
    int currentPos = positionInCurrentFace();
    Q_ASSERT( currentPos != -1 );
 
    const QgsMeshFace &currentFace = mFaces.at( mCurrentFace );
    int faceSize = currentFace.size();
    mLastValidFace = mCurrentFace;
    mCurrentFace = mFacesNeighborhood.at( mCurrentFace ).at( ( currentPos + faceSize - 1 ) % currentFace.count() );
  }
 
  return mCurrentFace;
}
 
int QgsMeshVertexCirculator::turnClockwise() const
{
  if ( mCurrentFace == -1 )
    mCurrentFace = mLastValidFace;
  else
  {
    int currentPos = positionInCurrentFace();
    Q_ASSERT( currentPos != -1 );
 
    const QgsMeshFace &currentFace = mFaces.at( mCurrentFace );
    int faceSize = currentFace.size();
    mLastValidFace = mCurrentFace;
    mCurrentFace = mFacesNeighborhood.at( mCurrentFace ).at( ( currentPos ) % faceSize );
  }
 
  return mCurrentFace;
}
 
int QgsMeshVertexCirculator::currentFaceIndex() const
{
  return mCurrentFace;
}
 
QgsMeshFace QgsMeshVertexCirculator::currentFace() const
{
  if ( mCurrentFace != -1 )
    return mFaces.at( mCurrentFace );
  else
    return QgsMeshFace();
}
 
bool QgsMeshVertexCirculator::goBoundaryClockwise() const
{
  if ( !isValid() )
    return false;
 
  if ( mCurrentFace == -1 )
    mCurrentFace = mLastValidFace;
 
  int firstFace = mCurrentFace;
 
  while ( turnClockwise() != -1 && currentFaceIndex() != firstFace ) {}
  if ( mCurrentFace == firstFace )
    return false; // a complete turn, so the vertex is not a boundary vertex, something went wrong
 
  return ( turnCounterClockwise() != -1 );
}
 
bool QgsMeshVertexCirculator::goBoundaryCounterClockwise() const
{
  if ( !isValid() )
    return false;
 
  if ( mCurrentFace == -1 )
    mCurrentFace = mLastValidFace;
 
  int firstFace = mCurrentFace;
 
  while ( turnCounterClockwise() != -1 && currentFaceIndex() != firstFace ) {}
  if ( mCurrentFace == firstFace )
    return false; // a complete turn, so the vertex is not a boundary vertex, something went wrong
 
  return ( turnClockwise() != -1 );
}
 
int QgsMeshVertexCirculator::oppositeVertexClockwise() const
{
  if ( mCurrentFace == -1 )
    return -1;
 
  const QgsMeshFace &face = currentFace();
 
  if ( face.isEmpty() )
    return -1;
 
  int vertexPosition = QgsTopologicalMesh::vertexPositionInFace( mVertexIndex, currentFace() );
 
  if ( vertexPosition == -1 )
    return -1;
 
  return face.at( ( vertexPosition + 1 ) % face.count() );
}
 
int QgsMeshVertexCirculator::oppositeVertexCounterClockwise() const
{
  if ( mCurrentFace == -1 )
    return -1;
 
  const QgsMeshFace &face = currentFace();
 
  if ( face.isEmpty() )
    return -1;
 
  int vertexPosition = QgsTopologicalMesh::vertexPositionInFace( mVertexIndex, currentFace() );
 
  if ( vertexPosition == -1 )
    return -1;
 
  return face.at( ( vertexPosition - 1 + face.count() ) % face.count() );
}
 
bool QgsMeshVertexCirculator::isValid() const
{
  return mIsValid;
}
 
QList<int> QgsMeshVertexCirculator::facesAround() const
{
  QList<int> ret;
  if ( !isValid() )
    return ret;
 
  if ( mCurrentFace != -1 )
    ret.append( mCurrentFace );
  while ( turnCounterClockwise() != ret.first() && currentFaceIndex() != -1 )
    ret.append( currentFaceIndex() );
 
 
  if ( currentFaceIndex() == -1 )  //we encounter a boundary, restart with other direction
  {
    ret.clear();
    if ( turnClockwise() == -1 )
      return ret;
    ret.append( currentFaceIndex() );
    while ( turnClockwise() != -1 )
    {
      ret.append( currentFaceIndex() );
    }
  }
 
  return ret;
}
 
int QgsMeshVertexCirculator::degree() const
{
  if ( mDegree != -1 )
    return mDegree;
 
  mDegree = 0;
  if ( !mIsValid )
    return mDegree;
 
  // if we are on the boundary, we count one more to take account of the circulator will
  // not cover the last vertex (the other vertex on boundary)
  if ( goBoundaryCounterClockwise() )
    mDegree = 2;
  else
    mDegree = 1;
 
  int firstFace = currentFaceIndex();
 
  while ( turnClockwise() != firstFace && currentFaceIndex() != -1 )
    ++mDegree;
 
  return mDegree;
}
 
 
int QgsMeshVertexCirculator::positionInCurrentFace() const
{
  if ( mCurrentFace < 0 || mCurrentFace >= mFaces.count() )
    return -1;
 
  return QgsTopologicalMesh::vertexPositionInFace( mVertexIndex, mFaces.at( mCurrentFace ) );
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::addFaces( const QgsTopologicalMesh::TopologicalFaces &topologicalFaces )
{
  Changes changes;
  changes.mFacesToAdd = topologicalFaces.mFaces;
  changes.mAddedFacesFirstIndex = mMesh->faceCount();
 
  changes.mFacesNeighborhoodToAdd.resize( changes.mFacesToAdd.count() );
  for ( int i = 0; i < changes.mFacesToAdd.count(); ++i )
    changes.mFacesNeighborhoodToAdd[i] = QVector<int>( changes.mFacesToAdd.at( i ).count(), -1 );
 
  for ( int boundary : topologicalFaces.mBoundaries )
  {
    //if the boundary is a free vertex in the destination mesh, no need to check
    if ( mVertexToFace.at( boundary ) == -1 )
      continue;
 
    const QList<int> &linkedFaces = topologicalFaces.mVerticesToFace.values( boundary );
    for ( int linkedFace : linkedFaces )
    {
      QgsMeshVertexCirculator newFacesCirculator( topologicalFaces, linkedFace, boundary );
      //search for face boundary on clockwise side of new faces
      newFacesCirculator.goBoundaryClockwise();
      int oppositeVertexForNewFace = newFacesCirculator.oppositeVertexClockwise();
      if ( mVertexToFace.at( oppositeVertexForNewFace ) == -1 )
        continue;
 
      QgsMeshVertexCirculator meshCirculator = vertexCirculator( boundary );
      meshCirculator.goBoundaryCounterClockwise();
      int oppositeVertexForMeshFace = meshCirculator.oppositeVertexCounterClockwise();
 
      const QgsMeshFace &newFaceBoundary = newFacesCirculator.currentFace();
      int boundaryPositionInNewFace = vertexPositionInFace( boundary, newFaceBoundary );
 
      if ( oppositeVertexForMeshFace != oppositeVertexForNewFace )
      {
        changes.mFacesNeighborhoodToAdd[newFacesCirculator.currentFaceIndex()][boundaryPositionInNewFace] = -1 ;
      }
      else
      {
        const QgsMeshFace &meshFaceBoundary = meshCirculator.currentFace();
        int boundaryPositionInMeshFace = vertexPositionInFace( meshCirculator.oppositeVertexCounterClockwise(), meshFaceBoundary );
 
        changes.mNeighborhoodChanges.append( std::array<int, 4>(
        {
          meshCirculator.currentFaceIndex(),
          boundaryPositionInMeshFace,
          -1,
          changes.addedFaceIndexInMesh( newFacesCirculator.currentFaceIndex() )
        } ) );
 
        changes.mFacesNeighborhoodToAdd[newFacesCirculator.currentFaceIndex()][boundaryPositionInNewFace] = meshCirculator.currentFaceIndex();
      }
    }
  }
 
  for ( int f = 0; f < changes.mFacesToAdd.count(); ++f )
    for ( int n = 0; n < changes.mFacesToAdd.at( f ).count(); ++n )
      if ( changes.mFacesNeighborhoodToAdd.at( f ).at( n ) == -1 )
        changes.mFacesNeighborhoodToAdd[f][n] = changes.addedFaceIndexInMesh( topologicalFaces.mFacesNeighborhood.at( f ).at( n ) );
 
  const QList<int> &verticesToFaceToChange = topologicalFaces.mVerticesToFace.keys();
  for ( const int vtc : verticesToFaceToChange )
    if ( mVertexToFace.at( vtc ) == -1 )
      changes.mVerticesToFaceChanges.append( {vtc,
                                              mVertexToFace.at( vtc ),
                                              changes.addedFaceIndexInMesh( topologicalFaces.mVerticesToFace.values( vtc ).first() )
                                             } );
 
  applyChanges( changes );
 
  return changes;
}
 
void QgsTopologicalMesh::applyChanges( const QgsTopologicalMesh::Changes &changes )
{
  int initialVerticesCount = mMesh->vertices.count();
  if ( !changes.mVerticesToAdd.empty() )
  {
    int newSize = mMesh->vertices.count() + changes.mVerticesToAdd.count();
    mMesh->vertices.resize( newSize );
    mVertexToFace.resize( newSize );
  }
 
  if ( !changes.mFacesToAdd.empty() )
  {
    int newSize = mMesh->faceCount() + changes.mFacesToAdd.count();
    mMesh->faces.resize( newSize );
    mFacesNeighborhood.resize( newSize );
  }
 
  for ( int i = 0; i < changes.mFacesToRemove.count(); ++i )
  {
    mMesh->faces[changes.removedFaceIndexInMesh( i )] = QgsMeshFace();
    mFacesNeighborhood[changes.removedFaceIndexInMesh( i )] = FaceNeighbors();//changes.facesNeighborhoodToRemove[i];
  }
 
  for ( int i = 0; i < changes.mVerticesToRemoveIndexes.count(); ++i )
  {
    int vertexIndex = changes.mVerticesToRemoveIndexes.at( i );
    if ( mVertexToFace.at( vertexIndex ) == -1 )
      dereferenceAsFreeVertex( vertexIndex );
    mMesh->vertices[vertexIndex] = QgsMeshVertex();
    mVertexToFace[vertexIndex] = -1;
  }
 
  for ( int i = 0; i < changes.mVerticesToAdd.count(); ++i )
  {
    mMesh->vertices[initialVerticesCount + i] = changes.mVerticesToAdd.at( i );
    mVertexToFace[initialVerticesCount + i] = changes.mVertexToFaceToAdd.at( i );
    if ( changes.mVertexToFaceToAdd.at( i ) == -1 )
      referenceAsFreeVertex( initialVerticesCount + i );
  }
 
  for ( int i = 0; i < changes.mFacesToAdd.count(); ++i )
  {
    mMesh->faces[changes.addedFaceIndexInMesh( i )] = changes.mFacesToAdd.at( i );
    mFacesNeighborhood[changes.addedFaceIndexInMesh( i )] = changes.mFacesNeighborhoodToAdd.at( i );
  }
 
  for ( const std::array<int, 4> neighborChange : std::as_const( changes.mNeighborhoodChanges ) )
  {
    const int faceIndex = neighborChange.at( 0 );
    const int positionInFace = neighborChange.at( 1 );
    const int valueToApply = neighborChange.at( 3 );
    mFacesNeighborhood[faceIndex][positionInFace] = valueToApply;
  }
 
  for ( const std::array<int, 3> vertexToFaceChange : std::as_const( changes.mVerticesToFaceChanges ) )
  {
    int vertexIndex = vertexToFaceChange.at( 0 );
    mVertexToFace[vertexToFaceChange.at( 0 )] = vertexToFaceChange.at( 2 );
 
    if ( vertexToFaceChange.at( 2 ) == -1 &&
         vertexToFaceChange.at( 1 ) != -1 &&
         !mMesh->vertices.at( vertexIndex ).isEmpty() )
      referenceAsFreeVertex( vertexIndex );
 
    if ( vertexToFaceChange.at( 1 ) == -1 && vertexToFaceChange.at( 2 ) != -1 )
      dereferenceAsFreeVertex( vertexIndex );
  }
 
  for ( int i = 0; i < changes.mChangeCoordinateVerticesIndexes.count(); ++i )
  {
    int vertexIndex = changes.mChangeCoordinateVerticesIndexes.at( i );
    if ( !changes.mNewZValues.isEmpty() )
      mMesh->vertices[vertexIndex].setZ( changes.mNewZValues.at( i ) );
    if ( !changes.mNewXYValues.isEmpty() )
    {
      const QgsPointXY &pt = changes.mNewXYValues.at( i );
      mMesh->vertices[vertexIndex].setX( pt.x() );
      mMesh->vertices[vertexIndex].setY( pt.y() );
    }
  }
}
 
void QgsTopologicalMesh::reverseChanges( const QgsTopologicalMesh::Changes &changes )
{
  for ( const std::array<int, 4> neighborChange : std::as_const( changes.mNeighborhoodChanges ) )
  {
    const int faceIndex = neighborChange.at( 0 );
    const int positionInFace = neighborChange.at( 1 );
    const int valueToApply = neighborChange.at( 2 );
    mFacesNeighborhood[faceIndex][positionInFace] = valueToApply;
  }
 
  for ( int i = 0; i < changes.mFacesToRemove.count(); ++i )
  {
    mMesh->faces[changes.removedFaceIndexInMesh( i )] = changes.mFacesToRemove.at( i );
    mFacesNeighborhood[changes.removedFaceIndexInMesh( i )] = changes.mFacesNeighborhoodToRemove.at( i );
  }
 
  for ( int i = 0; i < changes.mVerticesToRemoveIndexes.count(); ++i )
  {
    int vertexIndex = changes.mVerticesToRemoveIndexes.at( i );
    mMesh->vertices[vertexIndex] = changes.mRemovedVertices.at( i );
    mVertexToFace[vertexIndex] = changes.mVerticesToFaceRemoved.at( i );
    if ( mVertexToFace.at( vertexIndex ) == -1 )
      referenceAsFreeVertex( vertexIndex );
  }
 
  int verticesToFaceChangesCount = changes.mVerticesToFaceChanges.count();
  for ( int i = 0; i < verticesToFaceChangesCount; ++i )
  {
    const std::array<int, 3> vertexToFaceChange = changes.mVerticesToFaceChanges.at( verticesToFaceChangesCount - i - 1 );
    int vertexIndex = vertexToFaceChange.at( 0 );
    mVertexToFace[vertexIndex] = vertexToFaceChange.at( 1 );
 
    if ( vertexToFaceChange.at( 2 ) == -1 && vertexToFaceChange.at( 1 ) != -1 )
      dereferenceAsFreeVertex( vertexIndex );
 
    if ( vertexToFaceChange.at( 1 ) == -1 &&
         vertexToFaceChange.at( 2 ) != -1 &&
         !mMesh->vertex( vertexIndex ).isEmpty() )
      referenceAsFreeVertex( vertexIndex );
  }
 
  if ( !changes.mFacesToAdd.empty() )
  {
    int newSize = mMesh->faceCount() - changes.mFacesToAdd.count();
    mMesh->faces.resize( newSize );
    mFacesNeighborhood.resize( newSize );
  }
 
  if ( !changes.mVerticesToAdd.isEmpty() )
  {
    int newSize = mMesh->vertexCount() - changes.mVerticesToAdd.count();
 
    for ( int i = newSize; i < mMesh->vertexCount(); ++i )
      if ( mVertexToFace.at( i ) == -1 )
        dereferenceAsFreeVertex( i );
 
    mMesh->vertices.resize( newSize );
    mVertexToFace.resize( newSize );
  }
 
  for ( int i = 0; i < changes.mChangeCoordinateVerticesIndexes.count(); ++i )
  {
    int vertexIndex = changes.mChangeCoordinateVerticesIndexes.at( i );
    if ( !changes.mOldZValues.isEmpty() )
      mMesh->vertices[vertexIndex].setZ( changes.mOldZValues.at( i ) );
    if ( !changes.mOldXYValues.isEmpty() )
    {
      const QgsPointXY &pt = changes.mOldXYValues.at( i );
      mMesh->vertices[vertexIndex].setX( pt.x() );
      mMesh->vertices[vertexIndex].setY( pt.y() );
    }
  }
}
 
QgsMeshVertexCirculator QgsTopologicalMesh::vertexCirculator( int vertexIndex ) const
{
  return QgsMeshVertexCirculator( *this, vertexIndex );
}
 
QSet<int> QgsTopologicalMesh::concernedFacesBy( const QList<int> &faceIndexes ) const
{
  QSet<int> faces;
  for ( const int faceIndex : faceIndexes )
  {
    const QgsMeshFace &face = mMesh->face( faceIndex );
    for ( int i = 0; i < face.count(); ++i )
    {
      const QList<int> around =  facesAroundVertex( face.at( i ) );
      faces.unite( QSet< int >( around.begin(), around.end() ) );
    }
  }
  return faces;
}
 
void QgsTopologicalMesh::dereferenceAsFreeVertex( int vertexIndex )
{
  mFreeVertices.remove( vertexIndex );
}
 
void QgsTopologicalMesh::referenceAsFreeVertex( int vertexIndex )
{
  // QSet used to retrieve free vertex without going through all the vertices container.
  // But maybe later we could use more sophisticated reference (spatial index?), to retrieve free vertex in an extent
  mFreeVertices.insert( vertexIndex );
}
 
QgsMeshEditingError QgsTopologicalMesh::checkConsistency() const
{
  for ( int faceIndex = 0 ; faceIndex < mMesh->faces.count( ); ++faceIndex )
  {
    const QgsMeshFace &face = mMesh->faces.at( faceIndex );
    const FaceNeighbors &neighborhood = mFacesNeighborhood.at( faceIndex );
    if ( face.count() != neighborhood.count() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidFace, faceIndex );
    for ( int i = 0; i < face.count(); ++i )
    {
      int vertexIndex = face.at( i );
      // check if each vertices is linked to a face (not free vertex)
      if ( mVertexToFace.at( vertexIndex ) == -1 )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, vertexIndex );
 
      int neighborIndex = neighborhood.at( i );
      if ( neighborIndex != -1 )
      {
        const QgsMeshFace &neighborFace = mMesh->faces.at( neighborIndex );
        if ( neighborFace.isEmpty() )
          return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidFace, faceIndex );
        int neighborSize = neighborFace.size();
        const FaceNeighbors &neighborhoodOfNeighbor = mFacesNeighborhood.at( neighborIndex );
        int posInNeighbor = QgsTopologicalMesh::vertexPositionInFace( *mMesh, vertexIndex, neighborIndex );
        if ( neighborhoodOfNeighbor.isEmpty() || neighborhoodOfNeighbor.at( ( posInNeighbor + neighborSize - 1 ) % neighborSize ) != faceIndex )
          return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidFace, faceIndex );
      }
    }
  }
 
  for ( int vertexIndex = 0; vertexIndex < mMesh->vertexCount(); ++vertexIndex )
  {
    if ( mVertexToFace.at( vertexIndex ) != -1 )
    {
      if ( !mMesh->face( mVertexToFace.at( vertexIndex ) ).contains( vertexIndex ) )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, vertexIndex );
 
      if ( facesAroundVertex( vertexIndex ).count() == 0 )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, vertexIndex );
    }
  }
 
  return QgsMeshEditingError();
}
 
QgsMeshEditingError QgsTopologicalMesh::checkTopology( const QgsMesh &mesh, int maxVerticesPerFace )
{
  QgsMesh temp = mesh;
  QgsMeshEditingError error;
  createTopologicalMesh( &temp, maxVerticesPerFace, error );
  return error;
}
 
QgsMesh *QgsTopologicalMesh::mesh() const
{
  return mMesh;
}
 
int QgsTopologicalMesh::firstFaceLinked( int vertexIndex ) const
{
  if ( vertexIndex < 0 || vertexIndex >= mMesh->vertexCount() )
    return -1;
  return mVertexToFace.at( vertexIndex );
}
 
bool QgsTopologicalMesh::isVertexOnBoundary( int vertexIndex ) const
{
  QgsMeshVertexCirculator circulator = vertexCirculator( vertexIndex );
 
  if ( circulator.isValid() )
    return circulator.goBoundaryClockwise();
 
  return false;
}
 
bool QgsTopologicalMesh::isVertexFree( int vertexIndex ) const
{
  if ( vertexIndex < 0 || vertexIndex >= mMesh->vertexCount() )
    return false;
 
  if ( mMesh->vertices.at( vertexIndex ).isEmpty() )
    return false;
 
  return mVertexToFace.at( vertexIndex ) == -1;
}
 
QList<int> QgsTopologicalMesh::freeVerticesIndexes() const
{
  return QList<int>( mFreeVertices.begin(), mFreeVertices.end() );
}
 
QgsMeshEditingError QgsTopologicalMesh::checkTopologyOfVerticesAsFace( const QVector<QgsMeshVertex> &vertices, bool &clockwise )
{
  int size = vertices.size();
  int direction = 0;
  for ( int i = 0; i < size; ++i )
  {
    int iv0 =  i;
    int iv1 = ( i + 1 ) % size;
    int iv2 = ( i + 2 ) % size;
 
    const QgsMeshVertex &v0 = vertices.at( iv0 ) ;
    const QgsMeshVertex &v1 = vertices.at( iv1 ) ;
    const QgsMeshVertex &v2 = vertices.at( iv2 ) ;
 
    if ( v0.isEmpty() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, iv0 );
 
    if ( v1.isEmpty() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, iv1 );
 
    if ( v2.isEmpty() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, iv2 );
 
    double crossProd = crossProduct( v1, v0, v2 ); //if cross product>0, we have two edges clockwise
    if ( direction != 0 && crossProd * direction < 0 )   // We have a convex face or a (partially) flat face
    {
      clockwise = direction > 0;
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidFace, -1 );
    }
    else if ( crossProd == 0 )
    {
      clockwise = direction > 0;
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::FlatFace, -1 );
    }
    else if ( direction == 0 )
      direction = crossProd / std::fabs( crossProd );
  }
 
  clockwise = direction > 0;
 
  return QgsMeshEditingError( Qgis::MeshEditingErrorType::NoError, -1 );
}
 
QgsMeshEditingError QgsTopologicalMesh::counterClockwiseFaces( QgsMeshFace &face, QgsMesh *mesh )
{
  // First check the topology of the face, then put it counter clockwise if needed
  // If the indexes are not well ordered (edges intersect), invalid face
  int faceSize = face.count();
  if ( faceSize < 3 )
    return QgsMeshEditingError( Qgis::MeshEditingErrorType::FlatFace, -1 );
 
  QVector<QgsMeshVertex> vertices( face.size() );
 
  for ( int i = 0; i < faceSize; ++i )
  {
    int iv =  face[i];
    if ( iv < 0 || iv >= mesh->vertexCount() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, iv );
 
    vertices[i] = mesh->vertices.at( face[i] );
  }
 
  bool clockwise = false;
  QgsMeshEditingError error = QgsTopologicalMesh::checkTopologyOfVerticesAsFace( vertices, clockwise );
  if ( error != QgsMeshEditingError() )
    return error;
 
  if ( clockwise )// clockwise --> reverse the order of the index;
  {
    for ( int i = 0; i < faceSize / 2; ++i )
    {
      int temp = face[i];
      face[i] = face.at( faceSize - i - 1 );
      face[faceSize - i - 1] = temp;
    }
  }
 
  return QgsMeshEditingError( Qgis::MeshEditingErrorType::NoError, -1 );
}
 
void QgsTopologicalMesh::reindex()
{
  QVector<int> oldToNewIndex( mMesh->vertices.count(), -1 );
  int verticesTotalCount = mMesh->vertexCount();
  int oldIndex = 0;
  int newIndex = 0;
  while ( oldIndex < verticesTotalCount )
  {
    if ( mMesh->vertex( oldIndex ).isEmpty() )
    {
      oldIndex++;
    }
    else
    {
      oldToNewIndex[oldIndex] = newIndex;
      if ( oldIndex != newIndex )
        mMesh->vertices[newIndex] = mMesh->vertices[oldIndex];
      oldToNewIndex[oldIndex] = newIndex;
      oldIndex++;
      newIndex++;
    }
  }
 
  mMesh->vertices.resize( newIndex );
 
  oldIndex = 0;
  newIndex = 0;
  int facesTotalCount = mMesh->faceCount();
  while ( oldIndex < facesTotalCount )
  {
    if ( mMesh->face( oldIndex ).isEmpty() )
      oldIndex++;
    else
    {
      if ( oldIndex != newIndex )
        mMesh->faces[newIndex] = mMesh->faces[oldIndex];
      QgsMeshFace &face = mMesh->faces[newIndex];
      for ( int i = 0; i < face.count(); ++i )
        face[i] = oldToNewIndex[face.at( i )];
      newIndex++;
      oldIndex++;
    }
  }
 
  mMesh->faces.resize( newIndex );
 
  mVertexToFace.clear();
  mFacesNeighborhood.clear();
}
 
bool QgsTopologicalMesh::renumber()
{
  QVector<int> oldToNewVerticesIndexes;
  if ( !renumberVertices( oldToNewVerticesIndexes ) )
    return false;
 
 
  QVector<int> oldToNewFacesIndexes;
 
  if ( !renumberFaces( oldToNewFacesIndexes ) )
    return false;
 
  // If we are here, we can apply the renumbering
 
  QVector<QgsMeshVertex> tempVertices( mMesh->vertices.count() );
  for ( int i = 0; i < oldToNewVerticesIndexes.count(); ++i )
  {
    tempVertices[oldToNewVerticesIndexes.at( i )] = mMesh->vertex( i );
  }
  mMesh->vertices = tempVertices;
 
  QVector<QgsMeshFace> tempFaces( mMesh->faces.count() );
  for ( int i = 0; i < oldToNewFacesIndexes.count(); ++i )
  {
    tempFaces[oldToNewFacesIndexes.at( i )] = mMesh->face( i );
 
    QgsMeshFace &face = tempFaces[oldToNewFacesIndexes.at( i )];
 
    for ( int fi = 0; fi < face.count(); ++fi )
    {
      face[fi] = oldToNewVerticesIndexes.at( face.at( fi ) );
    }
  }
 
  mMesh->faces = tempFaces;
 
  return true;
 
}
 
bool QgsTopologicalMesh::renumberVertices( QVector<int> &oldToNewIndex ) const
{
  std::vector<QgsMeshVertexCirculator> circulators;
  circulators.reserve( mMesh->vertices.count() );
  int minDegree = std::numeric_limits<int>::max();
  int minDegreeVertex = -1;
 
  QQueue<int> queue;
  QSet<int> nonThreadedVertex;
  oldToNewIndex = QVector<int> ( mMesh->vertexCount(), -1 );
  for ( int i = 0; i < mMesh->vertexCount(); ++i )
  {
    circulators.emplace_back( *this, i );
    const QgsMeshVertexCirculator &circulator = circulators.back();
    if ( circulators.back().degree() < minDegree )
    {
      minDegreeVertex = circulators.size() - 1;
      minDegree = circulator.degree();
    }
    nonThreadedVertex.insert( i );
  }
 
  auto sortedNeighbor = [ = ]( QList<int> &neighbors, int index )
  {
    const QgsMeshVertexCirculator &circ = circulators.at( index );
 
    if ( !circ.isValid() )
      return;
 
    circ.goBoundaryCounterClockwise();
    neighbors.append( circ.oppositeVertexCounterClockwise() );
 
    int firstFace = circ.currentFaceIndex();
    do
    {
      int neighborIndex = circ.oppositeVertexClockwise();
      int degree = circulators.at( neighborIndex ).degree();
      QList<int>::Iterator it = neighbors.begin();
      while ( it != neighbors.end() )
      {
        if ( degree <= circulators.at( *it ).degree() )
        {
          neighbors.insert( it, neighborIndex );
          break;
        }
        it++;
      }
      if ( it == neighbors.end() )
        neighbors.append( neighborIndex );
    }
    while ( circ.turnClockwise() != firstFace && circ.currentFaceIndex() != -1 );
  };
 
  int newIndex = 0;
  int currentVertex = minDegreeVertex;
  nonThreadedVertex.remove( minDegreeVertex );
 
  while ( newIndex < mMesh->vertexCount() )
  {
    if ( oldToNewIndex[currentVertex] == -1 )
      oldToNewIndex[currentVertex] = newIndex++;
 
    if ( circulators.at( currentVertex ).isValid() )
    {
      QList<int> neighbors;
      sortedNeighbor( neighbors, currentVertex );
 
      for ( const int i : std::as_const( neighbors ) )
        if ( oldToNewIndex.at( i ) == -1 && nonThreadedVertex.contains( i ) )
        {
          queue.enqueue( i );
          nonThreadedVertex.remove( i );
        }
    }
 
    if ( queue.isEmpty() )
    {
      if ( nonThreadedVertex.isEmpty() && newIndex < mMesh->vertexCount() )
        return false;
 
      const QList<int> remainingVertex( nonThreadedVertex.constBegin(), nonThreadedVertex.constEnd() );
      int minRemainingDegree = std::numeric_limits<int>::max();
      int minRemainingVertex = -1;
      for ( const int i : remainingVertex )
      {
        int degree = circulators.at( i ).degree();
        if ( degree < minRemainingDegree )
        {
          minRemainingDegree = degree;
          minRemainingVertex = i;
        }
      }
      currentVertex = minRemainingVertex;
      nonThreadedVertex.remove( currentVertex );
    }
    else
    {
      currentVertex = queue.dequeue();
    }
  }
 
  return true;
}
 
bool QgsTopologicalMesh::renumberFaces( QVector<int> &oldToNewIndex ) const
{
  QQueue<int> queue;
  QSet<int> nonThreadedFaces;
 
  oldToNewIndex = QVector<int>( mMesh->faceCount(), -1 );
 
  QVector<int> faceDegrees( mMesh->faceCount(), 0 );
 
  int minDegree = std::numeric_limits<int>::max();
  int minDegreeFace = -1;
  for ( int faceIndex = 0; faceIndex < mMesh->faceCount(); ++faceIndex )
  {
    const FaceNeighbors &neighbors = mFacesNeighborhood.at( faceIndex );
 
    int degree = 0;
    for ( int n = 0; n < neighbors.size(); ++n )
    {
      if ( neighbors.at( n ) != -1 )
        degree++;
    }
 
    if ( degree < minDegree )
    {
      minDegree = degree;
      minDegreeFace = faceIndex;
    }
 
    faceDegrees[faceIndex] = degree;
    nonThreadedFaces.insert( faceIndex );
  }
 
  int newIndex = 0;
  int currentFace = minDegreeFace;
  nonThreadedFaces.remove( minDegreeFace );
 
  auto sortedNeighbor = [this, faceDegrees]( QList<int> &neighbors, int index )
  {
    const FaceNeighbors &neighborhood = mFacesNeighborhood.at( index );
 
    for ( int i = 0; i < neighborhood.count(); ++i )
    {
      int neighborIndex = neighborhood.at( i );
      if ( neighborIndex == -1 )
        continue;
 
      int degree = faceDegrees.at( neighborIndex );
      if ( neighbors.isEmpty() )
        neighbors.append( neighborIndex );
      else
      {
        QList<int>::Iterator it = neighbors.begin();
        while ( it != neighbors.end() )
        {
          if ( degree <= faceDegrees.at( *it ) )
          {
            neighbors.insert( it, neighborIndex );
            break;
          }
          it++;
        }
        if ( it == neighbors.end() )
          neighbors.append( neighborIndex );
      }
    }
  };
 
  while ( newIndex < mMesh->faceCount() )
  {
    if ( oldToNewIndex[currentFace] == -1 )
      oldToNewIndex[currentFace] = newIndex++;
 
    QList<int> neighbors;
    sortedNeighbor( neighbors, currentFace );
 
    for ( const int i : std::as_const( neighbors ) )
      if ( oldToNewIndex.at( i ) == -1 && nonThreadedFaces.contains( i ) )
      {
        queue.enqueue( i );
        nonThreadedFaces.remove( i );
      }
 
    if ( queue.isEmpty() )
    {
      if ( nonThreadedFaces.isEmpty() && newIndex < mMesh->faceCount() )
        return false;
 
      const QList<int> remainingFace( nonThreadedFaces.constBegin(), nonThreadedFaces.constEnd() );
      int minRemainingDegree = std::numeric_limits<int>::max();
      int minRemainingFace = -1;
      for ( const int i : remainingFace )
      {
        int degree = faceDegrees.at( i );
        if ( degree < minRemainingDegree )
        {
          minRemainingDegree = degree;
          minRemainingFace = i;
        }
      }
      currentFace = minRemainingFace;
      nonThreadedFaces.remove( currentFace );
    }
    else
    {
      currentFace = queue.dequeue();
    }
 
  }
 
  return true;
}
 
QVector<QgsMeshFace> QgsTopologicalMesh::Changes::addedFaces() const
{
  return mFacesToAdd;
}
 
QVector<QgsMeshFace> QgsTopologicalMesh::Changes::removedFaces() const
{
  return mFacesToRemove;
}
 
QList<int> QgsTopologicalMesh::Changes::removedFaceIndexes() const
{
  return mFaceIndexesToRemove;
}
 
QVector<QgsMeshVertex> QgsTopologicalMesh::Changes::addedVertices() const
{
  return mVerticesToAdd;
}
 
QList<int> QgsTopologicalMesh::Changes::changedCoordinatesVerticesIndexes() const
{
  return mChangeCoordinateVerticesIndexes;
}
 
QList<double> QgsTopologicalMesh::Changes::newVerticesZValues() const
{
  return mNewZValues;
}
 
QList<QgsPointXY> QgsTopologicalMesh::Changes::newVerticesXYValues() const
{
  return mNewXYValues;
}
 
QList<QgsPointXY> QgsTopologicalMesh::Changes::oldVerticesXYValues() const
{
  return mOldXYValues;
}
 
QList<int> QgsTopologicalMesh::Changes::nativeFacesIndexesGeometryChanged() const
{
  return mNativeFacesIndexesGeometryChanged;
}
 
bool QgsTopologicalMesh::Changes::isEmpty() const
{
  return ( mFaceIndexesToRemove.isEmpty() &&
           mFacesToAdd.isEmpty() &&
           mFacesNeighborhoodToAdd.isEmpty() &&
           mFacesToRemove.isEmpty() &&
           mFacesNeighborhoodToRemove.isEmpty() &&
           mNeighborhoodChanges.isEmpty() &&
           mVerticesToAdd.isEmpty() &&
           mVertexToFaceToAdd.isEmpty() &&
           mVerticesToRemoveIndexes.isEmpty() &&
           mRemovedVertices.isEmpty() &&
           mVerticesToFaceRemoved.isEmpty() &&
           mVerticesToFaceChanges.isEmpty() &&
           mChangeCoordinateVerticesIndexes.isEmpty() &&
           mNewZValues.isEmpty() &&
           mOldZValues.isEmpty() &&
           mNewXYValues.isEmpty() &&
           mOldXYValues.isEmpty() &&
           mNativeFacesIndexesGeometryChanged.isEmpty() );
}
 
QList<int> QgsTopologicalMesh::Changes::verticesToRemoveIndexes() const
{
  return mVerticesToRemoveIndexes;
}
 
int QgsTopologicalMesh::Changes::addedFaceIndexInMesh( int internalIndex ) const
{
  if ( internalIndex == -1 )
    return -1;
 
  return internalIndex + mAddedFacesFirstIndex;
}
 
int QgsTopologicalMesh::Changes::removedFaceIndexInMesh( int internalIndex ) const
{
  if ( internalIndex == -1 )
    return -1;
 
  return mFaceIndexesToRemove.at( internalIndex );
}
 
void QgsTopologicalMesh::Changes::clearChanges()
{
  mAddedFacesFirstIndex = 0;
  mFaceIndexesToRemove.clear();
  mFacesToAdd.clear();
  mFacesNeighborhoodToAdd.clear();
  mFacesToRemove.clear();
  mFacesNeighborhoodToRemove.clear();
  mNeighborhoodChanges.clear();
 
  mVerticesToAdd.clear();
  mVertexToFaceToAdd.clear();
  mVerticesToRemoveIndexes.clear();
  mRemovedVertices.clear();
  mVerticesToFaceRemoved.clear();
  mVerticesToFaceChanges.clear();
 
  mChangeCoordinateVerticesIndexes.clear();
  mNewZValues.clear();
  mOldZValues.clear();
  mNewXYValues.clear();
  mOldXYValues.clear();
  mNativeFacesIndexesGeometryChanged.clear();
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::addFreeVertex( const QgsMeshVertex &vertex )
{
  Changes changes;
  changes.mVerticesToAdd.append( vertex );
  changes.mVertexToFaceToAdd.append( -1 );
 
  mMesh->vertices.append( vertex );
  mVertexToFace.append( -1 );
  referenceAsFreeVertex( mMesh->vertices.count() - 1 );
 
  return changes;
}
 
// Returns the orientation of the polygon formed by mesh vertices, <0 counter clockwise; >0 clockwise
static double vertexPolygonOrientation( const QgsMesh &mesh, const QList<int> &vertexIndexes )
{
  if ( vertexIndexes.count() < 3 )
    return 0.0;
 
  int hullDomainVertexPos = -1;
  double xMin = std::numeric_limits<double>::max();
  double yMin = std::numeric_limits<double>::max();
  for ( int i = 0; i < vertexIndexes.count(); ++i )
  {
    const QgsMeshVertex &vertex = mesh.vertices.at( vertexIndexes.at( i ) );
    if ( xMin >= vertex.x() && yMin > vertex.y() )
    {
      hullDomainVertexPos = i;
      xMin = vertex.x();
      yMin = vertex.y();
    }
  }
 
  if ( hullDomainVertexPos >= 0 )
  {
    int iv1 = vertexIndexes.at( ( hullDomainVertexPos - 1 + vertexIndexes.count() ) % vertexIndexes.count() );
    int iv2 = vertexIndexes.at( ( hullDomainVertexPos + 1 ) % vertexIndexes.count() );
    int ivc = vertexIndexes.at( ( hullDomainVertexPos ) );
    double cp = crossProduct( ivc, iv1, iv2, mesh );
    return cp;
  }
 
  return 0.0;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::removeVertexFillHole( int vertexIndex )
{
  if ( vertexIndex >= mVertexToFace.count() )
    return Changes();
 
  if ( mVertexToFace.at( vertexIndex ) == -1 ) //it is a free vertex
  {
    Changes changes;
    changes.mRemovedVertices.append( mMesh->vertices.at( vertexIndex ) );
    changes.mVerticesToRemoveIndexes.append( vertexIndex );
    changes.mVerticesToFaceRemoved.append( -1 );
    dereferenceAsFreeVertex( vertexIndex );
    mMesh->vertices[vertexIndex] = QgsMeshVertex();
    return changes;
  }
 
  //search concerned faces
  QgsMeshVertexCirculator circulator = vertexCirculator( vertexIndex );
  circulator.goBoundaryClockwise();
  QList<int> boundariesVertexIndex;
  QList<int> associateFaceToBoundaries;
  QList<int> removedFacesIndexes;
  QSet<int> boundaryInGlobalMesh;
 
  do
  {
    removedFacesIndexes.append( circulator.currentFaceIndex() );
    boundariesVertexIndex.append( circulator.oppositeVertexClockwise() );
    Q_ASSERT( !mMesh->vertices.at( boundariesVertexIndex.last() ).isEmpty() );
    const QgsMeshFace &currentFace = circulator.currentFace();
    associateFaceToBoundaries.append( mFacesNeighborhood.at( circulator.currentFaceIndex() ).at(
                                        vertexPositionInFace( boundariesVertexIndex.last(), currentFace ) ) );
 
    if ( currentFace.count() > 3 ) // quad or more, need other vertices
    {
      int posInface = vertexPositionInFace( vertexIndex, currentFace );
      for ( int i = 2; i < currentFace.count() - 1; ++i )
      {
        boundariesVertexIndex.append( currentFace.at( ( posInface + i ) % currentFace.count() ) );
        Q_ASSERT( !mMesh->vertices.at( boundariesVertexIndex.last() ).isEmpty() );
        associateFaceToBoundaries.append( mFacesNeighborhood.at( circulator.currentFaceIndex() ).at(
                                            vertexPositionInFace( boundariesVertexIndex.last(), currentFace ) ) );
      }
    }
  }
  while ( circulator.turnCounterClockwise() != -1 && circulator.currentFaceIndex() != removedFacesIndexes.first() );
 
  bool boundaryFill = false;
  if ( circulator.currentFaceIndex() == -1 ) //we are on boundary of the mesh
  {
    boundaryFill = true;
    //we need to add last vertex/boundary faces that was not added because we are on mesh boundary
    circulator.goBoundaryCounterClockwise();
    int lastVertexIndex = circulator.oppositeVertexCounterClockwise();
    boundariesVertexIndex.append( lastVertexIndex );
 
    // but we can be on the case where the last vertex share an edge with the first point,
    // in the case, the associate face on boundarie is not -1, but the face on the other side of the edge
    QgsMeshVertexCirculator boundaryCirculator = vertexCirculator( lastVertexIndex );
    boundaryCirculator.goBoundaryCounterClockwise();
    if ( boundaryCirculator.oppositeVertexCounterClockwise() == boundariesVertexIndex.first() )
    {
      associateFaceToBoundaries.append( boundaryCirculator.currentFaceIndex() );
      boundaryFill = false; //we are not a boundary fill anymore
    }
    else
      associateFaceToBoundaries.append( -1 );
 
    for ( const int index : std::as_const( boundariesVertexIndex ) )
    {
      if ( isVertexOnBoundary( index ) )
        boundaryInGlobalMesh.insert( index );
    }
  }
 
  int currentVertexToFace = mVertexToFace.at( vertexIndex );
  // here, we use the method removeFaces that effectivly removes and then constructs changes
  Changes changes = removeFaces( removedFacesIndexes );
 
  QList<QList<int>> holes;
  QList<QList<int>> associateMeshFacesToHoles;
 
  bool cancelOperation = false;
 
  if ( boundaryFill )
  {
    // The hole is not a closed polygon, we need to close it, but the closing segment can intersect another segments/vertices.
    // In this case we consider as many polygons as necessary.
 
    int startPos = 0;
    int finalPos = boundariesVertexIndex.count() - 1;
    QList<int> uncoveredVertex;
 
    QList<int> partToCheck = boundariesVertexIndex.mid( startPos, finalPos - startPos + 1 );
    QList<int> associateFacePart = associateFaceToBoundaries.mid( startPos, finalPos - startPos + 1 );
    while ( startPos < finalPos && !partToCheck.isEmpty() )
    {
      // check if we intersect an edge between first and second
      int secondPos = partToCheck.count() - 1;
      const QgsPoint &closingSegmentExtremety1 = mMesh->vertex( partToCheck.at( 0 ) );
      const QgsPoint &closingSegmentExtremety2 = mMesh->vertex( partToCheck.last() );
      bool isEdgeIntersect = false;
      for ( int i = 1; i < secondPos - 1; ++i )
      {
        const QgsPoint &p1 = mMesh->vertex( partToCheck.at( i ) );
        const QgsPoint &p2 = mMesh->vertex( partToCheck.at( i + 1 ) );
        bool isLineIntersection;
        QgsPoint intersectPoint;
        isEdgeIntersect = QgsGeometryUtils::segmentIntersection( closingSegmentExtremety1, closingSegmentExtremety2, p1, p2, intersectPoint, isLineIntersection, 1e-8, true );
        if ( isEdgeIntersect )
          break;
      }
 
      int index = partToCheck.at( 0 );
      if ( boundaryInGlobalMesh.contains( index ) && index != boundariesVertexIndex.at( 0 ) )
      {
        cancelOperation = true;
        break;
      }
 
      // if uncovered vertex is a boundary vertex in the global mesh (except first that is always a boundary in the global mesh)
      // This operation will leads to a unique shared vertex that is not allowed, you have to cancel the operation
      if ( isEdgeIntersect || vertexPolygonOrientation( *mMesh, partToCheck ) >= 0 )
      {
        partToCheck.removeLast();
        associateFacePart.removeAt( associateFacePart.count() - 2 );
        if ( partToCheck.count() == 1 )
        {
          uncoveredVertex.append( index );
          startPos = startPos + 1;
          partToCheck = boundariesVertexIndex.mid( startPos, finalPos - startPos + 1 );
          associateFacePart = associateFaceToBoundaries.mid( startPos, finalPos - startPos + 1 );
        }
      }
      else
      {
        // store the well defined hole
        holes.append( partToCheck );
        associateMeshFacesToHoles.append( associateFacePart );
 
        startPos = startPos + partToCheck.count() - 1;
        uncoveredVertex.append( partToCheck.at( 0 ) );
        partToCheck = boundariesVertexIndex.mid( startPos, finalPos - startPos + 1 );
        associateFacePart = associateFaceToBoundaries.mid( startPos, finalPos - startPos + 1 );
      }
    }
  }
  else
  {
    holes.append( boundariesVertexIndex );
    associateMeshFacesToHoles.append( associateFaceToBoundaries );
  }
 
  if ( cancelOperation )
  {
    reverseChanges( changes );
    return Changes();
  }
 
  Q_ASSERT( holes.count() == associateMeshFacesToHoles.count() );
 
  changes.mRemovedVertices.append( mMesh->vertices.at( vertexIndex ) );
  changes.mVerticesToRemoveIndexes.append( vertexIndex );
  changes.mVerticesToFaceRemoved.append( currentVertexToFace );
  // these changes contain information that will lead to reference the removed vertex as free vertex when reverse/reapply
  dereferenceAsFreeVertex( vertexIndex );
  mMesh->vertices[vertexIndex] = QgsMeshVertex();
  mVertexToFace[vertexIndex] = -1;
 
  int oldFacesCount = mMesh->faceCount();
  for ( int h = 0; h < holes.count(); ++h )
  {
    const QList<int> &holeVertices = holes.at( h );
    const QList<int> &associateMeshFacesToHole = associateMeshFacesToHoles.at( h );
    QHash<p2t::Point *, int> mapPoly2TriPointToVertex;
    std::vector<p2t::Point *> holeToFill( holeVertices.count() );
    try
    {
      for ( int i = 0; i < holeVertices.count(); ++i )
      {
        const QgsMeshVertex &vertex = mMesh->vertex( holeVertices.at( i ) );
        holeToFill[i] = new p2t::Point( vertex.x(), vertex.y() );
        mapPoly2TriPointToVertex.insert( holeToFill[i], holeVertices.at( i ) );
      }
 
      std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( holeToFill ) );
 
      cdt->Triangulate();
      std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
      QVector<QgsMeshFace> newFaces( triangles.size() );
      for ( size_t i = 0; i < triangles.size(); ++i )
      {
        QgsMeshFace &face = newFaces[i];
        face.resize( 3 );
        for ( int j = 0; j < 3; j++ )
        {
          int vertInd = mapPoly2TriPointToVertex.value( triangles.at( i )->GetPoint( j ), -1 );
          if ( vertInd == -1 )
            throw std::exception();
          Q_ASSERT( !mMesh->vertices.at( vertInd ).isEmpty() );
          face[j] = vertInd;
        }
      }
 
      QgsMeshEditingError error;
      QgsTopologicalMesh::TopologicalFaces topologicalFaces = createNewTopologicalFaces( newFaces, false, error );
      if ( error.errorType != Qgis::MeshEditingErrorType::NoError )
        throw std::exception();
      int newFaceIndexStartIndex = mMesh->faceCount();
      QgsTopologicalMesh::Changes addChanges;
      addChanges.mAddedFacesFirstIndex = newFaceIndexStartIndex;
      addChanges.mFacesToAdd = topologicalFaces.meshFaces();
      addChanges.mFacesNeighborhoodToAdd = topologicalFaces.mFacesNeighborhood;
 
      // vertices to face changes
      const QList<int> &verticesToFaceToChange = topologicalFaces.mVerticesToFace.keys();
      for ( const int vtc : verticesToFaceToChange )
        if ( mVertexToFace.at( vtc ) == -1 )
          addChanges.mVerticesToFaceChanges.append( {vtc,
              mVertexToFace.at( vtc ),
              addChanges.addedFaceIndexInMesh( topologicalFaces.mVerticesToFace.values( vtc ).first() )
                                                    } );
 
 
      // reindex neighborhood for new faces
      for ( int i = 0; i < topologicalFaces.mFaces.count(); ++i )
      {
        FaceNeighbors &faceNeighbors = addChanges.mFacesNeighborhoodToAdd[i];
        faceNeighbors = topologicalFaces.mFacesNeighborhood.at( i );
        for ( int n = 0; n < faceNeighbors.count(); ++n )
        {
          if ( faceNeighbors.at( n ) != -1 )
            faceNeighbors[n] += newFaceIndexStartIndex; //reindex internal neighborhood
        }
      }
 
      // link neighborhood for boundaries of each side
      for ( int i = 0 ; i < holeVertices.count(); ++i )
      {
        int vertexHoleIndex = holeVertices.at( i );
        int meshFaceBoundaryIndex = associateMeshFacesToHole.at( i );
 
        const QgsMeshVertexCirculator circulator = QgsMeshVertexCirculator( topologicalFaces, vertexHoleIndex );
        circulator.goBoundaryClockwise();
        int newFaceBoundaryLocalIndex = circulator.currentFaceIndex();
        int newFaceBoundaryIndexInMesh = circulator.currentFaceIndex() + newFaceIndexStartIndex;
        const QgsMeshFace &newFace = circulator.currentFace();
        int positionInNewFaces = vertexPositionInFace( vertexHoleIndex, newFace );
 
        if ( meshFaceBoundaryIndex != -1 )
        {
          const QgsMeshFace meshFace = mMesh->face( meshFaceBoundaryIndex );
          int positionInMeshFaceBoundary = vertexPositionInFace( *mMesh, vertexHoleIndex, meshFaceBoundaryIndex );
          positionInMeshFaceBoundary = ( positionInMeshFaceBoundary - 1 + meshFace.count() ) % meshFace.count(); //take the position just before
 
          addChanges.mNeighborhoodChanges.append( {meshFaceBoundaryIndex, positionInMeshFaceBoundary, -1, newFaceBoundaryIndexInMesh} );
        }
 
        addChanges.mFacesNeighborhoodToAdd[newFaceBoundaryLocalIndex][positionInNewFaces] = meshFaceBoundaryIndex;
      }
 
      applyChanges( addChanges );
 
      changes.mFacesToAdd.append( addChanges.mFacesToAdd );
      changes.mFacesNeighborhoodToAdd.append( addChanges.mFacesNeighborhoodToAdd );
      //for each neighborhood change, check if a corresponding change already exist and merge them, if not just append
      for ( const std::array<int, 4> &neighborChangeToAdd : std::as_const( addChanges.mNeighborhoodChanges ) )
      {
        bool merged = false;
        for ( std::array<int, 4> &existingNeighborChange : changes.mNeighborhoodChanges )
        {
          if ( existingNeighborChange.at( 0 ) == neighborChangeToAdd.at( 0 ) &&
               existingNeighborChange.at( 1 ) == neighborChangeToAdd.at( 1 ) )
          {
            merged = true;
            Q_ASSERT( existingNeighborChange.at( 3 ) == neighborChangeToAdd.at( 2 ) );
            existingNeighborChange[3] = neighborChangeToAdd.at( 3 );
          }
        }
        if ( !merged )
          changes.mNeighborhoodChanges.append( neighborChangeToAdd );
      }
      //for each vertex to face change, check if a corresponding change already exist and merge them, if not just append
      for ( const std::array<int, 3> &verticesToFaceToAdd : std::as_const( addChanges.mVerticesToFaceChanges ) )
      {
        bool merged = false;
        for ( std::array<int, 3> &existingVerticesToFace : changes.mVerticesToFaceChanges )
        {
          if ( existingVerticesToFace.at( 0 ) == verticesToFaceToAdd.at( 0 ) )
          {
            merged = true;
            Q_ASSERT( existingVerticesToFace.at( 2 ) == verticesToFaceToAdd.at( 1 ) );
            existingVerticesToFace[2] = verticesToFaceToAdd.at( 2 );
          }
        }
        if ( !merged )
          changes.mVerticesToFaceChanges.append( verticesToFaceToAdd );
      }
 
      qDeleteAll( holeToFill );
    }
    catch ( ... )
    {
      qDeleteAll( holeToFill );
      QgsMessageLog::logMessage( QObject::tr( "Triangulation failed. Skipping hole…" ), QObject::tr( "Mesh Editing" ) );
    }
  }
  changes.mAddedFacesFirstIndex = oldFacesCount;
 
 
 
  changes.mAddedFacesFirstIndex = oldFacesCount;
 
  return changes;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::removeVertices( const QList<int> &vertices )
{
  QSet<int> facesIndex;
  //Search for associated faces
  for ( int vertexIndex : vertices )
  {
    const QList< int > faces = facesAroundVertex( vertexIndex );
    facesIndex.unite( QSet< int >( faces.begin(), faces.end() ) );
  }
 
  // remove the faces
  Changes changes = removeFaces( facesIndex.values() );
 
  // removes the vertices
  for ( int vertexIndex : vertices )
  {
    int currentVertexToFace = mVertexToFace.at( vertexIndex );
    // here, we use the method removeFaces that effectivly removes and then constructs changes
    changes.mRemovedVertices.append( mMesh->vertices.at( vertexIndex ) );
    changes.mVerticesToRemoveIndexes.append( vertexIndex );
    changes.mVerticesToFaceRemoved.append( currentVertexToFace );
    // these changes contain information that will lead to reference the removed vertex as free vertex when reverse/reapply
    dereferenceAsFreeVertex( vertexIndex );
    mMesh->vertices[vertexIndex] = QgsMeshVertex();
    mVertexToFace[vertexIndex] = -1;
  }
 
  return changes;
}
 
QgsMeshEditingError QgsTopologicalMesh::facesCanBeAdded( const QgsTopologicalMesh::TopologicalFaces &topologicFaces ) const
{
  QList<int> boundariesToCheckClockwiseInNewFaces = topologicFaces.mBoundaries;
  QList<std::array<int, 2>> boundariesToCheckCounterClockwiseInNewFaces; //couple boundary / associate face in topologicFaces.mVerticesToFace
  QList<int> uniqueSharedVertexBoundary;
 
 
  // Go through the boundary and search for opposite boundary vertex clockwise in new faces and compare
  // with boundary opposite vertices on the mesh
  // If, in the mesh, the opposite vertex counter clockwise is not the same , another check will be done
  // later with counter clockwise in new faces
  // If, in the mesh, the opposite vertex clockwise is the same, this is an error
  while ( !boundariesToCheckClockwiseInNewFaces.isEmpty() )
  {
    int boundary = boundariesToCheckClockwiseInNewFaces.takeLast();
 
    const QList<int> &linkedFaces = topologicFaces.mVerticesToFace.values( boundary );
 
    for ( int const linkedFace : linkedFaces )
    {
 
      //if the boundary is a free vertex in the destination mesh, no need to check
      if ( mVertexToFace.at( boundary ) == -1 )
        continue;
 
      QgsMeshVertexCirculator newFacescirculator( topologicFaces, linkedFace, boundary );
      QgsMeshVertexCirculator meshCirculator = vertexCirculator( boundary );
 
      if ( !newFacescirculator.isValid() )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, boundary );
 
      //Get the opposite vertex on the clockwise side with new faces block
      if ( !newFacescirculator.goBoundaryClockwise() )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, boundary );
 
      int oppositeVertexInNewFaces = newFacescirculator.oppositeVertexClockwise();
 
      if ( !meshCirculator.goBoundaryCounterClockwise() )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, boundary );
 
      int oppositeVertexCCWInMesh = meshCirculator.oppositeVertexCounterClockwise();
 
      if ( oppositeVertexCCWInMesh == oppositeVertexInNewFaces ) //this boundary is OK, continue wit next one
        continue;
      else
      {
        //to avoid manifold face that could pass through the check, compare not only the boundary edges but also with the opposite internal edge in new face
        const QgsMeshFace &newFaceOnBoundary = newFacescirculator.currentFace();
        int faceSize = newFaceOnBoundary.size();
        int posInNewFace = vertexPositionInFace( boundary, newFaceOnBoundary );
        int previousVertexIndex = ( posInNewFace + faceSize - 1 ) % faceSize;
        if ( newFaceOnBoundary.at( previousVertexIndex ) == oppositeVertexCCWInMesh )
          return QgsMeshEditingError( Qgis::MeshEditingErrorType::ManifoldFace, newFacescirculator.currentFaceIndex() );
      }
 
      meshCirculator.goBoundaryClockwise();
 
      int oppositeVertexCWInMesh = meshCirculator.oppositeVertexClockwise();
 
      if ( oppositeVertexCWInMesh == oppositeVertexInNewFaces )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::ManifoldFace, newFacescirculator.currentFaceIndex() );
 
      //if we are here we need more checks
      boundariesToCheckCounterClockwiseInNewFaces.append( {boundary, linkedFace} );
    }
  }
 
  // Check now with opposite boundary vertex counter clockwise in new faces
  while ( !boundariesToCheckCounterClockwiseInNewFaces.isEmpty() )
  {
    std::array<int, 2> boundaryLinkedface = boundariesToCheckCounterClockwiseInNewFaces.takeLast();
    int boundary = boundaryLinkedface.at( 0 );
    int linkedFace = boundaryLinkedface.at( 1 );
 
    QgsMeshVertexCirculator newFacescirculator( topologicFaces, linkedFace, boundary );
    QgsMeshVertexCirculator meshCirculator = vertexCirculator( boundary );
 
    if ( !newFacescirculator.goBoundaryCounterClockwise() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, boundary );
 
    int oppositeVertexInNewFaces = newFacescirculator.oppositeVertexCounterClockwise();
 
    if ( !meshCirculator.goBoundaryClockwise() )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, boundary );
 
    int oppositeVertexCWInMesh = meshCirculator.oppositeVertexClockwise();
 
    if ( oppositeVertexCWInMesh == oppositeVertexInNewFaces ) //this boundary is OK, continue with next one
      continue;
 
    meshCirculator.goBoundaryCounterClockwise();
 
    int oppositeVertexCCWInMesh = meshCirculator.oppositeVertexCounterClockwise();
 
    if ( oppositeVertexCCWInMesh == oppositeVertexInNewFaces )
      return QgsMeshEditingError( Qgis::MeshEditingErrorType::ManifoldFace, newFacescirculator.currentFaceIndex() );
 
    uniqueSharedVertexBoundary.append( boundary );
  }
 
  if ( !uniqueSharedVertexBoundary.isEmpty() )
    return QgsMeshEditingError( Qgis::MeshEditingErrorType::UniqueSharedVertex, uniqueSharedVertexBoundary.first() );
 
  // Check if internal vertices of new faces block are free in the mesh
  QSet<int> boundaryVertices( topologicFaces.mBoundaries.constBegin(), topologicFaces.mBoundaries.constEnd() );
  for ( const QgsMeshFace &newFace : std::as_const( topologicFaces.mFaces ) )
  {
    for ( const int vertexIndex : newFace )
    {
      if ( boundaryVertices.contains( vertexIndex ) )
        continue;
      if ( mVertexToFace.at( vertexIndex ) != -1 )
        return QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidVertex, vertexIndex );
    }
  }
 
  return QgsMeshEditingError();
}
 
void QgsTopologicalMesh::TopologicalFaces::clear()
{
  mFaces.clear();
  mFacesNeighborhood.clear();
  mVerticesToFace.clear();
  mBoundaries.clear();
}
 
QVector<QgsTopologicalMesh::FaceNeighbors> QgsTopologicalMesh::TopologicalFaces::facesNeighborhood() const
{
  return mFacesNeighborhood;
}
 
int QgsTopologicalMesh::TopologicalFaces::vertexToFace( int vertexIndex ) const
{
  if ( mVerticesToFace.contains( vertexIndex ) )
    return mVerticesToFace.values( vertexIndex ).at( 0 );
 
  return -1;
}
 
QgsTopologicalMesh QgsTopologicalMesh::createTopologicalMesh( QgsMesh *mesh, int maxVerticesPerFace, QgsMeshEditingError &error )
{
  QgsTopologicalMesh topologicMesh;
  topologicMesh.mMesh = mesh;
  topologicMesh.mVertexToFace = QVector<int>( mesh->vertexCount(), -1 );
  topologicMesh.mMaximumVerticesPerFace = maxVerticesPerFace;
  error.errorType = Qgis::MeshEditingErrorType::NoError;
 
  for ( int i = 0; i < mesh->faceCount(); ++i )
  {
    if ( mesh->face( i ).isEmpty() )
      continue;
    if ( maxVerticesPerFace != 0 && mesh->face( i ).count() > maxVerticesPerFace )
    {
      error = QgsMeshEditingError( Qgis::MeshEditingErrorType::InvalidFace, i );
      break;
    }
 
    error = counterClockwiseFaces( mesh->faces[i], mesh );
    if ( error.errorType != Qgis::MeshEditingErrorType::NoError )
    {
      if ( error.errorType == Qgis::MeshEditingErrorType::InvalidFace || error.errorType == Qgis::MeshEditingErrorType::FlatFace )
        error.elementIndex = i;
      break;
    }
  }
 
  if ( error.errorType == Qgis::MeshEditingErrorType::NoError )
  {
    TopologicalFaces subMesh = createTopologicalFaces( mesh->faces, &topologicMesh.mVertexToFace, error, false );
    topologicMesh.mFacesNeighborhood = subMesh.mFacesNeighborhood;
 
    for ( int i = 0; i < topologicMesh.mMesh->vertexCount(); ++i )
    {
      if ( topologicMesh.mVertexToFace.at( i ) == -1 )
        topologicMesh.mFreeVertices.insert( i );
    }
  }
 
  return topologicMesh;
}
 
QgsTopologicalMesh::TopologicalFaces QgsTopologicalMesh::createNewTopologicalFaces( const QVector<QgsMeshFace> &faces,  bool uniqueSharedVertexAllowed, QgsMeshEditingError &error )
{
  return createTopologicalFaces( faces, nullptr, error, uniqueSharedVertexAllowed );
}
 
 
QgsTopologicalMesh::TopologicalFaces QgsTopologicalMesh::createTopologicalFaces(
  const QVector<QgsMeshFace> &faces,
  QVector<int> *globalVertexToFace,
  QgsMeshEditingError &error,
  bool allowUniqueSharedVertex )
{
  int facesCount = faces.count();
  QVector<FaceNeighbors> faceTopologies;
  QMultiHash<int, int> verticesToFace;
 
  error = QgsMeshEditingError();
  TopologicalFaces ret;
 
  // Contains for each vertex a map (opposite vertex # edge) --> face index
  // when turning counter clockwise, if v1 first vertex and v2 second one, [v1][v2]--> neighbor face
  QMap<int, QMap<int, int>> verticesToNeighbor;
 
  for ( int faceIndex = 0; faceIndex < facesCount; ++faceIndex )
  {
    const QgsMeshFace &face = faces.at( faceIndex );
    int faceSize = face.count();
    //Fill vertices to neighbor faces
    for ( int i = 0; i < faceSize; ++i )
    {
      int v1 = face[i % faceSize];
      int v2 = face[( i + 1 ) % faceSize];
      if ( verticesToNeighbor[v2].contains( v1 ) )
      {
        error = QgsMeshEditingError( Qgis::MeshEditingErrorType::ManifoldFace, faceIndex );
        return ret;
      }
      else
        verticesToNeighbor[v2].insert( v1, faceIndex );
    }
  }
 
  faceTopologies = QVector<FaceNeighbors>( faces.count() );
 
  QSet<int> boundaryVertices;
 
  for ( int faceIndex = 0; faceIndex < facesCount; ++faceIndex )
  {
    const QgsMeshFace &face = faces.at( faceIndex );
    int faceSize = face.size();
    FaceNeighbors &faceTopology = faceTopologies[faceIndex];
    faceTopology.resize( faceSize );
 
    for ( int i = 0; i < faceSize; ++i )
    {
      int v1 = face.at( i );
      int v2 = face.at( ( i + 1 ) % faceSize );
 
      if ( globalVertexToFace )
      {
        if ( ( *globalVertexToFace )[v1] == -1 )
          ( *globalVertexToFace )[v1] = faceIndex ;
      }
      else
      {
        if ( allowUniqueSharedVertex || !verticesToFace.contains( v1 ) )
          verticesToFace.insert( v1, faceIndex ) ;
      }
 
      QMap<int, int> &edges = verticesToNeighbor[v1];
      if ( edges.contains( v2 ) )
        faceTopology[i] = edges.value( v2 );
      else
      {
        faceTopology[i] = -1;
 
        if ( !allowUniqueSharedVertex )
        {
          if ( boundaryVertices.contains( v1 ) )
          {
            error = QgsMeshEditingError( Qgis::MeshEditingErrorType::UniqueSharedVertex, v1 ); // if a vertices is more than one time in the boundary, that means faces share only one vertices
            return ret;
          }
        }
        boundaryVertices.insert( v1 );
      }
    }
  }
 
  ret.mFaces = faces;
  ret.mFacesNeighborhood = faceTopologies;
  ret.mBoundaries = boundaryVertices.values();
  ret.mVerticesToFace = verticesToFace;
  return ret;
}
 
QVector<int> QgsTopologicalMesh::neighborsOfFace( int faceIndex ) const
{
  return mFacesNeighborhood.at( faceIndex );
}
 
QList<int> QgsTopologicalMesh::facesAroundVertex( int vertexIndex ) const
{
  QgsMeshVertexCirculator circ = vertexCirculator( vertexIndex );
 
  return circ.facesAround();
}
 
QgsMeshEditingError QgsTopologicalMesh::facesCanBeRemoved( const QList<int> &facesIndexes )
{
  QSet<int> removedFaces( facesIndexes.begin(), facesIndexes.end() );
  QSet<int> concernedFaces = concernedFacesBy( facesIndexes );
 
  for ( const int f : std::as_const( removedFaces ) )
    concernedFaces.remove( f );
 
  QVector<QgsMeshFace> remainingFaces;
  remainingFaces.reserve( concernedFaces.count() );
  for ( const int f : std::as_const( concernedFaces ) )
    remainingFaces.append( mMesh->face( f ) );
 
  QgsMeshEditingError error;
  createTopologicalFaces( remainingFaces, nullptr, error, false );
 
  return error;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::removeFaces( const QList<int> &facesIndexesToRemove )
{
  Changes changes;
  changes.mFaceIndexesToRemove = facesIndexesToRemove;
  changes.mFacesToRemove.resize( facesIndexesToRemove.count() );
  changes.mFacesNeighborhoodToRemove.resize( facesIndexesToRemove.count() );
 
  QSet<int> indexSet( facesIndexesToRemove.begin(), facesIndexesToRemove.end() );
  QSet<int> threatedVertex;
 
  for ( int i = 0; i < facesIndexesToRemove.count(); ++i )
  {
    const int faceIndex = facesIndexesToRemove.at( i );
    const QgsMeshFace &face = mMesh->face( faceIndex );
    changes.mFacesToRemove[i] = face;
    const FaceNeighbors &neighborhood = mFacesNeighborhood.at( faceIndex );
    changes.mFacesNeighborhoodToRemove[i] = neighborhood;
    for ( int j = 0; j < face.count(); ++j )
    {
      //change the neighborhood for each neighbor face
      int neighborIndex = neighborhood.at( j );
      if ( neighborIndex != -1 && !indexSet.contains( neighborIndex ) )
      {
        int positionInNeighbor = mFacesNeighborhood.at( neighborIndex ).indexOf( faceIndex );
        changes.mNeighborhoodChanges.append( {neighborIndex, positionInNeighbor, faceIndex, -1} );
      }
 
      //change vertexToFace
      int vertexIndex = face.at( j );
      if ( !threatedVertex.contains( vertexIndex ) && indexSet.contains( mVertexToFace.at( vertexIndex ) ) )
      {
        int oldValue = mVertexToFace.at( vertexIndex );
        //look for another face linked to this vertex
        int refValue = -1;
        if ( neighborIndex != -1 && !indexSet.contains( neighborIndex ) ) //if exist, simpler to take it
          refValue = neighborIndex;
        else
        {
          QList<int> aroundFaces = facesAroundVertex( vertexIndex );
          aroundFaces.removeOne( faceIndex );
          if ( !aroundFaces.isEmpty() )
          {
            while ( !aroundFaces.isEmpty() && refValue == -1 )
            {
              if ( !indexSet.contains( aroundFaces.first() ) )
                refValue = aroundFaces.first();
              else
                aroundFaces.removeFirst();
            }
          }
        }
        changes.mVerticesToFaceChanges.append( {vertexIndex, oldValue, refValue} );
        threatedVertex.insert( vertexIndex );
      }
    }
  }
 
  applyChanges( changes );
 
  return changes;
}
 
bool QgsTopologicalMesh::eitherSideFacesAndVertices( int vertexIndex1,
    int vertexIndex2,
    int &face1,
    int &face2,
    int &neighborVertex1InFace1,
    int &neighborVertex1InFace2,
    int &neighborVertex2inFace1,
    int &neighborVertex2inFace2 ) const
{
  QgsMeshVertexCirculator circulator1 = vertexCirculator( vertexIndex1 );
  QgsMeshVertexCirculator circulator2 = vertexCirculator( vertexIndex2 );
 
  circulator1.goBoundaryClockwise();
  int firstFace1 = circulator1.currentFaceIndex();
  circulator2.goBoundaryClockwise();
  int firstFace2 = circulator2.currentFaceIndex();
 
  if ( circulator1.oppositeVertexCounterClockwise() != vertexIndex2 )
    while ( circulator1.turnCounterClockwise() != -1 &&
            circulator1.currentFaceIndex() != firstFace1 &&
            circulator1.oppositeVertexCounterClockwise() != vertexIndex2 ) {}
 
  if ( circulator2.oppositeVertexCounterClockwise() != vertexIndex1 )
    while ( circulator2.turnCounterClockwise() != -1 &&
            circulator2.currentFaceIndex() != firstFace2 &&
            circulator2.oppositeVertexCounterClockwise() != vertexIndex1 ) {}
 
  if ( circulator1.oppositeVertexCounterClockwise() != vertexIndex2
       || circulator2.oppositeVertexCounterClockwise() != vertexIndex1 )
    return false;
 
  face1 = circulator1.currentFaceIndex();
  face2 = circulator2.currentFaceIndex();
 
  neighborVertex1InFace1 = circulator1.oppositeVertexClockwise();
  circulator1.turnCounterClockwise();
  neighborVertex1InFace2 = circulator1.oppositeVertexCounterClockwise();
 
  neighborVertex2inFace2 = circulator2.oppositeVertexClockwise();
  circulator2.turnCounterClockwise();
  neighborVertex2inFace1 = circulator2.oppositeVertexCounterClockwise();
 
  return true;
}
 
bool QgsTopologicalMesh::edgeCanBeFlipped( int vertexIndex1, int vertexIndex2 ) const
{
  int faceIndex1;
  int faceIndex2;
  int oppositeVertexFace1;
  int oppositeVertexFace2;
  int supposedOppositeVertexFace1;
  int supposedoppositeVertexFace2;
 
  bool result = eitherSideFacesAndVertices(
                  vertexIndex1,
                  vertexIndex2,
                  faceIndex1,
                  faceIndex2,
                  oppositeVertexFace1,
                  supposedoppositeVertexFace2,
                  supposedOppositeVertexFace1,
                  oppositeVertexFace2 );
 
  if ( !result ||
       faceIndex1 < 0 ||
       faceIndex2 < 0 ||
       oppositeVertexFace1 < 0 ||
       oppositeVertexFace2 < 0 ||
       supposedOppositeVertexFace1 != oppositeVertexFace1 ||
       supposedoppositeVertexFace2 != oppositeVertexFace2 )
    return false;
 
  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
 
 
  if ( face1.count() != 3 || face2.count() != 3 )
    return false;
 
  double crossProduct1 = crossProduct( vertexIndex1, oppositeVertexFace1, oppositeVertexFace2, *mMesh );
  double crossProduct2 = crossProduct( vertexIndex2, oppositeVertexFace1, oppositeVertexFace2, *mMesh );
 
  return crossProduct1 * crossProduct2 < 0;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::flipEdge( int vertexIndex1, int vertexIndex2 )
{
  int faceIndex1;
  int faceIndex2;
  int oppositeVertexFace1;
  int oppositeVertexFace2;
  int supposedOppositeVertexFace1;
  int supposedoppositeVertexFace2;
 
  bool result = eitherSideFacesAndVertices(
                  vertexIndex1,
                  vertexIndex2,
                  faceIndex1,
                  faceIndex2,
                  oppositeVertexFace1,
                  supposedoppositeVertexFace2,
                  supposedOppositeVertexFace1,
                  oppositeVertexFace2 );
 
  if ( !result ||
       faceIndex1 < 0 ||
       faceIndex2 < 0 ||
       oppositeVertexFace1 < 0 ||
       oppositeVertexFace2 < 0 ||
       supposedOppositeVertexFace1 != oppositeVertexFace1 ||
       supposedoppositeVertexFace2 != oppositeVertexFace2 )
    return Changes();
 
 
  Changes changes;
 
  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
 
  Q_ASSERT( face1.count() == 3 );
  Q_ASSERT( face2.count() == 3 );
 
  int pos1 = vertexPositionInFace( vertexIndex1, face1 );
  int pos2 = vertexPositionInFace( vertexIndex2, face2 );
 
  int neighborFace1 = mFacesNeighborhood.at( faceIndex1 ).at( pos1 );
  int posInNeighbor1 = QgsTopologicalMesh::vertexPositionInFace( *mMesh, oppositeVertexFace1, neighborFace1 );
  int neighborFace2 = mFacesNeighborhood.at( faceIndex1 ).at( ( pos1 + 1 ) % 3 );
  int posInNeighbor2 = QgsTopologicalMesh::vertexPositionInFace( *mMesh, vertexIndex2, neighborFace2 );
  int neighborFace3 = mFacesNeighborhood.at( faceIndex2 ).at( pos2 );
  int posInNeighbor3 = QgsTopologicalMesh::vertexPositionInFace( *mMesh, oppositeVertexFace2, neighborFace3 );
  int neighborFace4 = mFacesNeighborhood.at( faceIndex2 ).at( ( pos2 + 1 ) % 3 );
  int posInNeighbor4 = QgsTopologicalMesh::vertexPositionInFace( *mMesh, vertexIndex1, neighborFace4 );
 
  changes.mFaceIndexesToRemove.append( faceIndex1 );
  changes.mFaceIndexesToRemove.append( faceIndex2 );
  changes.mFacesToRemove.append( face1 );
  changes.mFacesToRemove.append( face2 );
  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex1 ) );
  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex2 ) );
  int startIndex = mMesh->faceCount();
  changes.mAddedFacesFirstIndex = startIndex;
  changes.mFacesToAdd.append( {oppositeVertexFace1, oppositeVertexFace2, vertexIndex1} );
  changes.mFacesToAdd.append( {oppositeVertexFace2, oppositeVertexFace1, vertexIndex2} );
  changes.mFacesNeighborhoodToAdd.append( {startIndex + 1,
                                          mFacesNeighborhood.at( faceIndex2 ).at( ( pos2 + 1 ) % 3 ),
                                          mFacesNeighborhood.at( faceIndex1 ).at( pos1 )
                                          } );
  changes.mFacesNeighborhoodToAdd.append( {startIndex,
                                          mFacesNeighborhood.at( faceIndex1 ).at( ( pos1 + 1 ) % 3 ),
                                          mFacesNeighborhood.at( faceIndex2 ).at( pos2 )
                                          } );
 
  if ( neighborFace1 >= 0 )
    changes.mNeighborhoodChanges.append( {neighborFace1, posInNeighbor1, faceIndex1, startIndex} );
  if ( neighborFace2 >= 0 )
    changes.mNeighborhoodChanges.append( {neighborFace2, posInNeighbor2, faceIndex1, startIndex + 1} );
  if ( neighborFace3 >= 0 )
    changes.mNeighborhoodChanges.append( {neighborFace3, posInNeighbor3, faceIndex2, startIndex + 1} );
  if ( neighborFace4 >= 0 )
    changes.mNeighborhoodChanges.append( {neighborFace4, posInNeighbor4, faceIndex2, startIndex} );
 
 
  if ( mVertexToFace.at( vertexIndex1 ) == faceIndex1 || mVertexToFace.at( vertexIndex1 ) == faceIndex2 )
    changes.mVerticesToFaceChanges.append( {vertexIndex1,  mVertexToFace.at( vertexIndex1 ), startIndex} );
  if ( mVertexToFace.at( vertexIndex2 ) == faceIndex1 || mVertexToFace.at( vertexIndex2 ) == faceIndex2 )
    changes.mVerticesToFaceChanges.append( {vertexIndex2,  mVertexToFace.at( vertexIndex2 ), startIndex + 1} );
 
  if ( mVertexToFace.at( oppositeVertexFace1 ) == faceIndex1 )
    changes.mVerticesToFaceChanges.append( {oppositeVertexFace1,  faceIndex1, startIndex} );
 
  if ( mVertexToFace.at( oppositeVertexFace2 ) == faceIndex2 )
    changes.mVerticesToFaceChanges.append( {oppositeVertexFace2,  faceIndex2, startIndex + 1} );
 
  applyChanges( changes );
 
  return changes;
}
 
bool QgsTopologicalMesh::canBeMerged( int vertexIndex1, int vertexIndex2 ) const
{
  int faceIndex1;
  int faceIndex2;
  int neighborVertex1InFace1;
  int neighborVertex1InFace2;
  int neighborVertex2inFace1;
  int neighborVertex2inFace2;
 
  bool result = eitherSideFacesAndVertices(
                  vertexIndex1,
                  vertexIndex2,
                  faceIndex1,
                  faceIndex2,
                  neighborVertex1InFace1,
                  neighborVertex1InFace2,
                  neighborVertex2inFace1,
                  neighborVertex2inFace2 );
 
  if ( !result ||
       faceIndex1 < 0 ||
       faceIndex2 < 0 )
    return false;
 
  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
 
  if ( face1.count() + face2.count() - 2 > mMaximumVerticesPerFace )
    return false;
 
  QgsMeshVertex v1 = mMesh->vertices.at( vertexIndex1 );
  QgsMeshVertex v2 = mMesh->vertices.at( vertexIndex2 );
  QgsMeshVertex nv11 = mMesh->vertices.at( neighborVertex1InFace1 );
  QgsMeshVertex nv12 = mMesh->vertices.at( neighborVertex1InFace2 );
  QgsMeshVertex nv21 = mMesh->vertices.at( neighborVertex2inFace1 );
  QgsMeshVertex nv22 = mMesh->vertices.at( neighborVertex2inFace2 );
 
  double crossProduct1 = crossProduct( vertexIndex1, neighborVertex1InFace1, neighborVertex1InFace2, *mMesh );
  double crossProduct2 = crossProduct( vertexIndex2, neighborVertex2inFace1, neighborVertex2inFace2, *mMesh );
 
  return crossProduct1 * crossProduct2 < 0;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::merge( int vertexIndex1, int vertexIndex2 )
{
  int faceIndex1;
  int faceIndex2;
  int neighborVertex1InFace1;
  int neighborVertex1InFace2;
  int neighborVertex2inFace1;
  int neighborVertex2inFace2;
 
  bool result = eitherSideFacesAndVertices(
                  vertexIndex1,
                  vertexIndex2,
                  faceIndex1,
                  faceIndex2,
                  neighborVertex1InFace1,
                  neighborVertex1InFace2,
                  neighborVertex2inFace1,
                  neighborVertex2inFace2 );
 
  if ( !result ||
       faceIndex1 < 0 ||
       faceIndex2 < 0 )
    return Changes();
 
  Changes changes;
 
  const QgsMeshFace &face1 = mMesh->face( faceIndex1 );
  const QgsMeshFace &face2 = mMesh->face( faceIndex2 );
  int faceSize1 = face1.count();
  int faceSize2 = face2.count();
 
  int pos1 = QgsTopologicalMesh::vertexPositionInFace( vertexIndex1, face1 );
  int pos2 = QgsTopologicalMesh::vertexPositionInFace( vertexIndex2, face2 );
 
  changes.mFaceIndexesToRemove.append( faceIndex1 );
  changes.mFaceIndexesToRemove.append( faceIndex2 );
  changes.mFacesToRemove.append( face1 );
  changes.mFacesToRemove.append( face2 );
  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex1 ) );
  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex2 ) );
  int startIndex = mMesh->faceCount();
  changes.mAddedFacesFirstIndex = startIndex;
 
  QgsMeshFace newface;
  FaceNeighbors newNeighborhood;
  for ( int i = 0; i < faceSize1 - 1; ++i )
  {
    int currentPos = ( pos1 + i ) % faceSize1;
    newface.append( face1.at( currentPos ) ); //add vertex of face1
 
    int currentNeighbor = mFacesNeighborhood.at( faceIndex1 ).at( currentPos );
    newNeighborhood.append( currentNeighbor );
 
    if ( currentNeighbor != -1 )
    {
      int currentPosInNeighbor = QgsTopologicalMesh::vertexPositionInFace( *mMesh, face1.at( ( currentPos + 1 ) % faceSize1 ), currentNeighbor );
      changes.mNeighborhoodChanges.append( {currentNeighbor, currentPosInNeighbor, faceIndex1, startIndex} );
    }
  }
  for ( int i = 0; i < faceSize2 - 1; ++i )
  {
    int currentPos = ( pos2 + i ) % faceSize2;
    newface.append( face2.at( currentPos ) ); //add vertex of face2
 
    int currentNeighbor = mFacesNeighborhood.at( faceIndex2 ).at( currentPos );
    newNeighborhood.append( currentNeighbor );
 
    if ( currentNeighbor != -1 )
    {
      int currentPosInNeighbor = QgsTopologicalMesh::vertexPositionInFace( *mMesh, face2.at( ( currentPos + 1 ) % faceSize2 ), currentNeighbor );
      changes.mNeighborhoodChanges.append( {currentNeighbor, currentPosInNeighbor, faceIndex2, startIndex} );
    }
  }
 
  for ( int i = 0; i < faceSize1; ++i )
    if ( mVertexToFace.at( face1.at( i ) ) == faceIndex1 )
      changes.mVerticesToFaceChanges.append( {face1.at( i ), faceIndex1, startIndex} );
 
  for ( int i = 0; i < faceSize2; ++i )
    if ( mVertexToFace.at( face2.at( i ) ) == faceIndex2 )
      changes.mVerticesToFaceChanges.append( {face2.at( i ), faceIndex2, startIndex} );
 
  changes.mFacesToAdd.append( newface );
  changes.mFacesNeighborhoodToAdd.append( newNeighborhood );
 
  applyChanges( changes );
 
  return changes;
}
 
bool QgsTopologicalMesh::canBeSplit( int faceIndex ) const
{
  const QgsMeshFace face = mMesh->face( faceIndex );
 
  return face.count() == 4;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::splitFace( int faceIndex )
{
  //search for the spliited angle (greater angle)
  const QgsMeshFace &face = mMesh->face( faceIndex );
  int faceSize = face.count();
 
  Q_ASSERT( faceSize == 4 );
 
  double maxAngle = 0;
  int splitVertexPos = -1;
  for ( int i = 0; i < faceSize; ++i )
  {
    QgsVector vect1( mMesh->vertex( face.at( i ) ) - mMesh->vertex( face.at( ( i + 1 ) % faceSize ) ) );
    QgsVector vect2( mMesh->vertex( face.at( ( i + 2 ) % faceSize ) ) - mMesh->vertex( face.at( ( i + 1 ) % faceSize ) ) );
 
    double angle = std::abs( vect1.angle( vect2 ) );
    angle = std::min( angle, 2.0 * M_PI - angle );
    if ( angle > maxAngle )
    {
      maxAngle = angle;
      splitVertexPos = ( i + 1 ) % faceSize;
    }
  }
 
  Changes changes;
  if ( splitVertexPos == -1 )
    return changes;
 
  const QgsMeshFace newFace1 = {face.at( splitVertexPos ),
                                face.at( ( splitVertexPos + 1 ) % faceSize ),
                                face.at( ( splitVertexPos + 2 ) % faceSize )
                               };
 
  const QgsMeshFace newFace2 = {face.at( splitVertexPos ),
                                face.at( ( splitVertexPos + 2 ) % faceSize ),
                                face.at( ( splitVertexPos + 3 ) % faceSize )
                               };
 
  QVector<int> neighborIndex( faceSize );
  QVector<int> posInNeighbor( faceSize );
 
  for ( int i = 0; i < faceSize; ++i )
  {
    neighborIndex[i] = mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + i ) % faceSize );
    posInNeighbor[i] = QgsTopologicalMesh::vertexPositionInFace( *mMesh,  face.at( ( splitVertexPos + i + 1 ) % faceSize ), neighborIndex[i] );
  }
 
  changes.mFaceIndexesToRemove.append( faceIndex );
  changes.mFacesToRemove.append( face );
  changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( faceIndex ) );
  int startIndex = mMesh->faceCount();
  changes.mAddedFacesFirstIndex = startIndex;
  changes.mFacesToAdd.append( newFace1 );
  changes.mFacesToAdd.append( newFace2 );
 
  changes.mFacesNeighborhoodToAdd.append( {mFacesNeighborhood.at( faceIndex ).at( splitVertexPos ),
                                          mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + 1 ) % faceSize ),
                                          startIndex + 1
                                          } );
  changes.mFacesNeighborhoodToAdd.append( {startIndex,
                                          mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + 2 ) % faceSize ),
                                          mFacesNeighborhood.at( faceIndex ).at( ( splitVertexPos + 3 ) % faceSize )
                                          } );
 
  for ( int i = 0; i < faceSize; ++i )
  {
    if ( neighborIndex[i] >= 0 )
      changes.mNeighborhoodChanges.append( {neighborIndex[i], posInNeighbor[i], faceIndex, startIndex + int( i / 2 )} );
 
    int vertexIndex = face.at( ( splitVertexPos + i ) % faceSize );
    if ( mVertexToFace.at( vertexIndex ) == faceIndex )
      changes.mVerticesToFaceChanges.append( {vertexIndex, faceIndex, startIndex + int( i / 2 )} );
  }
 
  applyChanges( changes );
 
  return changes;
}
 
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::addVertexInFace( int includingFaceIndex, const QgsMeshVertex &vertex )
{
  Changes changes;
  changes.mVerticesToAdd.append( vertex );
  changes.mVertexToFaceToAdd.append( -1 );
 
  mMesh->vertices.append( vertex );
  mVertexToFace.append( -1 );
  changes.mAddedFacesFirstIndex = mMesh->faceCount();
 
  const QgsMeshFace includingFace = mMesh->face( includingFaceIndex );
  const FaceNeighbors includingFaceNeighborhood = mFacesNeighborhood.at( includingFaceIndex );
  int includingFaceSize = includingFace.count();
 
  for ( int i = 0; i < includingFaceSize; ++i )
  {
    // add a new face
    QgsMeshFace face( 3 );
    face[0] = mMesh->vertexCount() - 1;
    face[1] = includingFace.at( i );
    face[2] = includingFace.at( ( i + 1 ) % includingFaceSize );
    mMesh->faces.append( face );
    changes.mFacesToAdd.append( face );
 
    int currentVertexIndex = includingFace.at( i );
    if ( mVertexToFace.at( currentVertexIndex ) == includingFaceIndex )
    {
      int newFaceIndex = mMesh->faceCount() - 1;
      mVertexToFace[currentVertexIndex] = newFaceIndex;
      changes.mVerticesToFaceChanges.append( {currentVertexIndex, includingFaceIndex, newFaceIndex} );
    }
 
    int includingFaceNeighbor = includingFaceNeighborhood.at( i );
    FaceNeighbors neighbors(
    {
      changes.mAddedFacesFirstIndex + ( i + includingFaceSize - 1 ) % includingFaceSize,
      includingFaceNeighbor,
      changes.mAddedFacesFirstIndex + ( i + includingFaceSize + 1 ) % includingFaceSize
    } );
    mFacesNeighborhood.append( neighbors );
    changes.mFacesNeighborhoodToAdd.append( neighbors );
 
    if ( includingFaceNeighbor != -1 )
    {
      int indexInNeighbor = QgsTopologicalMesh::vertexPositionInFace( *mMesh, includingFace.at( ( i + 1 ) % includingFaceSize ), includingFaceNeighbor );
      int oldValue = mFacesNeighborhood[includingFaceNeighbor][indexInNeighbor];
      mFacesNeighborhood[includingFaceNeighbor][indexInNeighbor] = changes.mAddedFacesFirstIndex + i;
      changes.mNeighborhoodChanges.append( {includingFaceNeighbor, indexInNeighbor, oldValue, changes.mAddedFacesFirstIndex + i} );
    }
  }
 
  changes.mFacesToRemove.append( includingFace );
  changes.mFaceIndexesToRemove.append( includingFaceIndex );
  changes.mFacesNeighborhoodToRemove.append( includingFaceNeighborhood );
 
  mFacesNeighborhood[includingFaceIndex] = FaceNeighbors();
  mMesh->faces[includingFaceIndex] = QgsMeshFace();
  mVertexToFace[mVertexToFace.count() - 1] = mMesh->faceCount() - 1;
  changes.mVertexToFaceToAdd[changes.mVertexToFaceToAdd.count() - 1] = mMesh->faceCount() - 1 ;
 
  return changes;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::insertVertexInFacesEdge( int faceIndex, int position, const QgsMeshVertex &vertexToInsert )
{
  const QgsMeshFace face1 = mMesh->face( faceIndex );
 
  Changes changes;
  changes.mVerticesToAdd.append( vertexToInsert );
  changes.mAddedFacesFirstIndex = mMesh->faceCount();
 
  // triangulate the first face
  int newVertexPositionInFace1 = position + 1;
 
  auto triangulate = [this, &changes]( int removedFaceIndex, const QgsMeshVertex & newVertex, int newVertexPosition, QVector<int> &edgeFacesIndexes )->bool
  {
    const QgsMeshFace &initialFace = mMesh->face( removedFaceIndex );
    changes.mFacesToRemove.append( initialFace );
    changes.mFaceIndexesToRemove.append( removedFaceIndex );
    changes.mFacesNeighborhoodToRemove.append( mFacesNeighborhood.at( removedFaceIndex ) );
    const int addedVertexIndex = mMesh->vertexCount();
 
    int faceStartGlobalIndex = mMesh->faceCount() + changes.mFacesToAdd.count();
    int localStartIndex = changes.mFacesToAdd.count();
 
    QVector<int> newBoundary = initialFace;
    newBoundary.insert( newVertexPosition, addedVertexIndex );
 
    try
    {
      QHash<p2t::Point *, int> mapPoly2TriPointToVertex;
      std::vector<p2t::Point *> faceToFill( newBoundary.count() );
      for ( int i = 0; i < newBoundary.count(); ++i )
      {
        QgsMeshVertex vert;
 
        if ( newBoundary.at( i ) == addedVertexIndex )
          vert = newVertex;
        else
          vert = mMesh->vertex( newBoundary.at( i ) );
 
        faceToFill[i] = new p2t::Point( vert.x(), vert.y() );
        mapPoly2TriPointToVertex.insert( faceToFill[i], newBoundary.at( i ) );
      }
 
      std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( faceToFill ) );
      cdt->Triangulate();
      std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
      QVector<QgsMeshFace> newFaces( triangles.size() );
      for ( size_t i = 0; i < triangles.size(); ++i )
      {
        QgsMeshFace &face = newFaces[i];
        face.resize( 3 );
        for ( int j = 0; j < 3; j++ )
        {
          int vertInd = mapPoly2TriPointToVertex.value( triangles.at( i )->GetPoint( j ), -1 );
          if ( vertInd == -1 )
            throw std::exception();
          face[j] = vertInd;
        }
      }
 
      QgsMeshEditingError error;
      QgsTopologicalMesh::TopologicalFaces topologicalFaces = createNewTopologicalFaces( newFaces, false, error );
      if ( error.errorType != Qgis::MeshEditingErrorType::NoError )
        throw std::exception();
 
      changes.mFacesToAdd.append( topologicalFaces.meshFaces() );
      changes.mFacesNeighborhoodToAdd.append( topologicalFaces.mFacesNeighborhood );
 
      // vertices to face changes
      const QList<int> &verticesToFaceToChange = topologicalFaces.mVerticesToFace.keys();
      for ( const int vtc : verticesToFaceToChange )
        if ( vtc != addedVertexIndex && mVertexToFace.at( vtc ) == removedFaceIndex )
          changes.mVerticesToFaceChanges.append(
        {
          vtc,
          removedFaceIndex,
          topologicalFaces.vertexToFace( vtc ) + faceStartGlobalIndex
        } );
 
      // reindex neighborhood for new faces
      for ( int i = 0; i < topologicalFaces.mFaces.count(); ++i )
      {
        FaceNeighbors &faceNeighbors = changes.mFacesNeighborhoodToAdd[i + localStartIndex];
        for ( int n = 0; n < faceNeighbors.count(); ++n )
        {
          if ( faceNeighbors.at( n ) != -1 )
            faceNeighbors[n] += faceStartGlobalIndex; //reindex internal neighborhood
        }
      }
 
      edgeFacesIndexes.resize( 2 );
      // link neighborhood for boundaries of each side
      for ( int i = 0 ; i < newBoundary.count(); ++i )
      {
        int vertexIndex = newBoundary.at( i );
        QgsMeshVertexCirculator circulator = QgsMeshVertexCirculator( topologicalFaces, vertexIndex );
        circulator.goBoundaryClockwise();
        int newFaceBoundaryLocalIndex = localStartIndex + circulator.currentFaceIndex();
 
        int newFaceBoundaryIndexInMesh = faceStartGlobalIndex;
 
        int meshFaceBoundaryIndex;
        if ( i == newVertexPosition )
        {
          meshFaceBoundaryIndex = -1; //face that are on the opposite side of the edge, filled later
          edgeFacesIndexes[0] =  newFaceBoundaryLocalIndex;
        }
        else if ( i == ( newVertexPosition + newBoundary.count() - 1 ) % newBoundary.count() )
        {
          meshFaceBoundaryIndex = -1; //face that are on the opposite side of the edge, filled later
          edgeFacesIndexes[1] =  newFaceBoundaryLocalIndex;
        }
        else
          meshFaceBoundaryIndex = mFacesNeighborhood.at( removedFaceIndex ).at( QgsTopologicalMesh::vertexPositionInFace( vertexIndex, initialFace ) );
 
        const QgsMeshFace &newFace = circulator.currentFace();
        int positionInNewFaces = QgsTopologicalMesh::vertexPositionInFace( vertexIndex, newFace );
 
        if ( meshFaceBoundaryIndex != -1 )
        {
          const QgsMeshFace meshFace = mMesh->face( meshFaceBoundaryIndex );
          int positionInMeshFaceBoundary = QgsTopologicalMesh::vertexPositionInFace( *mMesh, vertexIndex, meshFaceBoundaryIndex );
          positionInMeshFaceBoundary = ( positionInMeshFaceBoundary - 1 + meshFace.count() ) % meshFace.count(); //take the position just before
 
          changes.mNeighborhoodChanges.append( {meshFaceBoundaryIndex,
                                                positionInMeshFaceBoundary,
                                                removedFaceIndex,
                                                newFaceBoundaryIndexInMesh +
                                                circulator.currentFaceIndex()
                                               } );
        }
 
        changes.mFacesNeighborhoodToAdd[newFaceBoundaryLocalIndex][positionInNewFaces] = meshFaceBoundaryIndex;
      }
 
      qDeleteAll( faceToFill );
    }
    catch ( ... )
    {
      return false;
    }
 
    return true;
  };
 
  QVector<int> edgeFacesIndexes;
  if ( !triangulate( faceIndex, vertexToInsert, newVertexPositionInFace1, edgeFacesIndexes ) )
    return Changes();
 
  changes.mVertexToFaceToAdd.append( edgeFacesIndexes.at( 0 ) + changes.mAddedFacesFirstIndex );
 
  int addedVertexIndex = mMesh->vertexCount();
 
  //if exist, triangulate the second face if exists
  int face2Index = mFacesNeighborhood.at( faceIndex ).at( position );
  if ( face2Index != -1 )
  {
    const QgsMeshFace &face2 = mMesh->face( face2Index );
    int vertexPositionInFace2 = vertexPositionInFace( face1.at( position ), face2 );
    QVector<int> edgeFacesIndexesFace2;
    if ( !triangulate( face2Index, vertexToInsert, vertexPositionInFace2, edgeFacesIndexesFace2 ) )
      return Changes();
 
    //link neighborhood with other side
    const QgsMeshFace &firstFaceSide1 = changes.mFacesToAdd.at( edgeFacesIndexes.at( 0 ) );
    int pos1InFaceSide1 = vertexPositionInFace( addedVertexIndex, firstFaceSide1 );
 
    const QgsMeshFace &secondFaceSide1 = changes.mFacesToAdd.at( edgeFacesIndexes.at( 1 ) );
    int pos2InFaceSide1 = vertexPositionInFace( addedVertexIndex, secondFaceSide1 );
    pos2InFaceSide1 = ( pos2InFaceSide1 + secondFaceSide1.size() - 1 ) % secondFaceSide1.size();
 
    const QgsMeshFace &firstFaceSide2 = changes.mFacesToAdd.at( edgeFacesIndexesFace2.at( 0 ) );
    int pos1InFaceSide2 = vertexPositionInFace( addedVertexIndex, firstFaceSide2 );
 
    const QgsMeshFace &secondFaceSide2 = changes.mFacesToAdd.at( edgeFacesIndexesFace2.at( 1 ) );
    int pos2InFaceSide2 = vertexPositionInFace( addedVertexIndex, secondFaceSide2 );
    pos2InFaceSide2 = ( pos2InFaceSide2 + secondFaceSide1.size() - 1 ) % secondFaceSide1.size();
 
    changes.mFacesNeighborhoodToAdd[edgeFacesIndexes.at( 0 )][pos1InFaceSide1] = edgeFacesIndexesFace2.at( 1 ) + changes.mAddedFacesFirstIndex;
    changes.mFacesNeighborhoodToAdd[edgeFacesIndexes.at( 1 )][pos2InFaceSide1] = edgeFacesIndexesFace2.at( 0 ) + changes.mAddedFacesFirstIndex;
    changes.mFacesNeighborhoodToAdd[edgeFacesIndexesFace2.at( 0 )][pos1InFaceSide2] = edgeFacesIndexes.at( 1 ) + changes.mAddedFacesFirstIndex;
    changes.mFacesNeighborhoodToAdd[edgeFacesIndexesFace2.at( 1 )][pos2InFaceSide2] = edgeFacesIndexes.at( 0 ) + changes.mAddedFacesFirstIndex;
  }
 
  applyChanges( changes );
  return changes;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::changeZValue( const QList<int> &verticesIndexes, const QList<double> &newValues )
{
  Q_ASSERT( verticesIndexes.count() == newValues.count() );
  Changes changes;
  changes.mChangeCoordinateVerticesIndexes.reserve( verticesIndexes.count() );
  changes.mNewZValues.reserve( verticesIndexes.count() );
  changes.mOldZValues.reserve( verticesIndexes.count() );
  for ( int i = 0; i < verticesIndexes.count(); ++i )
  {
    changes.mChangeCoordinateVerticesIndexes.append( verticesIndexes.at( i ) );
    changes.mOldZValues.append( mMesh->vertices.at( verticesIndexes.at( i ) ).z() );
    changes.mNewZValues.append( newValues.at( i ) );
  }
 
  applyChanges( changes );
 
  return changes;
}
 
QgsTopologicalMesh::Changes QgsTopologicalMesh::changeXYValue( const QList<int> &verticesIndexes, const QList<QgsPointXY> &newValues )
{
  Q_ASSERT( verticesIndexes.count() == newValues.count() );
  Changes changes;
  changes.mChangeCoordinateVerticesIndexes.reserve( verticesIndexes.count() );
  changes.mNewXYValues.reserve( verticesIndexes.count() );
  changes.mOldXYValues.reserve( verticesIndexes.count() );
  QSet<int> concernedFace;
  for ( int i = 0; i < verticesIndexes.count(); ++i )
  {
    changes.mChangeCoordinateVerticesIndexes.append( verticesIndexes.at( i ) );
    changes.mOldXYValues.append( mMesh->vertices.at( verticesIndexes.at( i ) ) );
    changes.mNewXYValues.append( newValues.at( i ) );
    const QList< int > faces = facesAroundVertex( verticesIndexes.at( i ) );
    concernedFace.unite( QSet< int>( faces.begin(), faces.end() ) );
  }
 
  changes.mNativeFacesIndexesGeometryChanged = concernedFace.values();
 
  applyChanges( changes );
 
  return changes;
}
 
bool QgsTopologicalMesh::delaunayConditionForEdge( int vertexIndex1, int vertexIndex2 )
{
  int faceIndex1;
  int faceIndex2;
  int oppositeVertexFace1;
  int oppositeVertexFace2;
  int supposedOppositeVertexFace1;
  int supposedoppositeVertexFace2;
 
  bool result =  eitherSideFacesAndVertices(
                   vertexIndex1,
                   vertexIndex2,
                   faceIndex1,
                   faceIndex2,
                   oppositeVertexFace1,
                   supposedoppositeVertexFace2,
                   supposedOppositeVertexFace1,
                   oppositeVertexFace2 );
 
  if ( ! result )
    return false;
 
  const QgsMeshFace face1 = mMesh->face( faceIndex1 );
  const QgsMeshFace face2 = mMesh->face( faceIndex2 );
 
  QgsCircle circle = QgsCircle::from3Points( mMesh->vertex( face1.at( 0 ) ),
                     mMesh->vertex( face1.at( 1 ) ),
                     mMesh->vertex( face1.at( 2 ) ) );
  bool circle1ContainsPoint = circle.contains( mMesh->vertex( supposedoppositeVertexFace2 ) );
 
  circle = QgsCircle::from3Points( mMesh->vertex( face2.at( 0 ) ),
                                   mMesh->vertex( face2.at( 1 ) ),
                                   mMesh->vertex( face2.at( 2 ) ) );
  bool circle2ContainsPoint = circle.contains( mMesh->vertex( supposedOppositeVertexFace1 ) );
 
  return !( circle1ContainsPoint || circle2ContainsPoint );
}