Added function to deduplicate points in MeshBuilder

src/TNL/Meshes/MeshBuilder.h

@@ -16,6 +16,9 @@
 
 #pragma once
 
+#include <numeric>   // std::iota
+#include <vector>
+
 #include <TNL/Containers/Vector.h>
 #include <TNL/Meshes/Topologies/Polyhedron.h>
 
@@ -31,14 +34,16 @@ public:
    using GlobalIndexType    = typename MeshTraitsType::GlobalIndexType;
    using LocalIndexType     = typename MeshTraitsType::LocalIndexType;
    using PointType          = typename MeshTraitsType::PointType;
+   using PointArrayType     = typename MeshTraitsType::PointArrayType;
+   using BoolVector         = Containers::Vector< bool, Devices::Host, GlobalIndexType >;
    using CellTopology       = typename MeshTraitsType::CellTopology;
    using CellSeedMatrixType = typename MeshTraitsType::CellSeedMatrixType;
    using CellSeedType       = typename CellSeedMatrixType::EntitySeedMatrixSeed;
    using FaceSeedMatrixType = typename MeshTraitsType::FaceSeedMatrixType;
    using FaceSeedType       = typename FaceSeedMatrixType::EntitySeedMatrixSeed;
    using NeighborCountsArray = typename MeshTraitsType::NeighborCountsArray;
-   
-   void setEntitiesCount( const GlobalIndexType& points, 
+
+   void setEntitiesCount( const GlobalIndexType& points,
                           const GlobalIndexType& cells = 0,
                           const GlobalIndexType& faces = 0 )
    {
@@ -109,20 +114,136 @@ public:
       return this->cellSeeds.getSeed( index );
    }
 
+   void deduplicatePoints( const double numericalThreshold = 1e-9 )
+   {
+      // prepare vector with an identity permutationutation
+      std::vector< GlobalIndexType > permutation( points.getSize() );
+      std::iota( permutation.begin(), permutation.end(), (GlobalIndexType) 0 );
+
+      // workaround for lexicographical sorting
+      // FIXME: https://mmg-gitlab.fjfi.cvut.cz/gitlab/tnl/tnl-dev/-/issues/79
+      auto lexless = [numericalThreshold, this] (const GlobalIndexType& a, const GlobalIndexType& b) -> bool
+      {
+         const PointType& left = this->points[ a ];
+         const PointType& right = this->points[ b ];
+         for( LocalIndexType i = 0; i < PointType::getSize(); i++ )
+            if( TNL::abs( left[ i ] - right[ i ] ) > numericalThreshold )
+               return left[ i ] < right[ i ];
+         return false;
+      };
+
+      // sort points in lexicographical order
+      std::stable_sort( permutation.begin(), permutation.end(), lexless );
+
+      // old -> new index mapping for points
+      std::vector< GlobalIndexType > points_perm_to_new( points.getSize() );
+
+      // find duplicate points
+      GlobalIndexType uniquePointsCount = 0;
+      // first index is unique
+      points_perm_to_new[ permutation[ 0 ] ] = uniquePointsCount++;
+      for( GlobalIndexType i = 1; i < points.getSize(); i++ ) {
+         const PointType& curr = points[ permutation[ i ] ];
+         const PointType& prev = points[ permutation[ i - 1 ] ];
+         if( maxNorm( curr - prev ) > numericalThreshold )
+            // unique point
+            points_perm_to_new[ permutation[ i ] ] = uniquePointsCount++;
+         else
+            // duplicate point - use previous index
+            points_perm_to_new[ permutation[ i ] ] = uniquePointsCount - 1;
+      }
+
+      // if all points are unique, we are done
+      if( uniquePointsCount == points.getSize() )
+         return;
+
+      std::cout << "Found " << points.getSize() - uniquePointsCount << " duplicate points (total " << points.getSize() << ", unique " << uniquePointsCount << ")" << std::endl;
+
+      // copy this->points and this->pointsSet, drop duplicate points
+      // (trying to do this in-place is not worth it, since even Array::reallocate
+      // needs to allocate a temporary array and copy the elements)
+      PointArrayType newPoints( uniquePointsCount );
+      BoolVector newPointsSet( uniquePointsCount );
+      std::vector< GlobalIndexType > points_old_to_new( points.getSize() );
+      // TODO: this can almost be parallelized, except we have multiple writes for the duplicate points
+      for( std::size_t i = 0; i < points_perm_to_new.size(); i++ ) {
+         const GlobalIndexType oldIndex = permutation[ i ];
+         const GlobalIndexType newIndex = points_perm_to_new[ oldIndex ];
+         newPoints[ newIndex ] = points[ oldIndex ];
+         newPointsSet[ newIndex ] = pointsSet[ oldIndex ];
+         points_old_to_new[ oldIndex ] = newIndex;
+      }
+      points = std::move( newPoints );
+      pointsSet = std::move( newPointsSet );
+      // reset permutation and points_perm_to_new - we need just points_old_to_new further on
+      permutation.clear();
+      permutation.shrink_to_fit();
+      points_perm_to_new.clear();
+      points_perm_to_new.shrink_to_fit();
+
+      auto remap_matrix = [uniquePointsCount, &points_old_to_new] ( auto& seeds )
+      {
+         // TODO: optimize out useless copy and counting corner counts - it would be sufficient to just remap column indices and then shrink the columns count
+         std::remove_reference_t< decltype( seeds ) > newSeeds;
+         newSeeds.setDimensions( seeds.getMatrix().getRows(), uniquePointsCount );
+         newSeeds.setEntityCornersCounts( seeds.getEntityCornerCounts() );
+         // TODO: parallelize (we have the IndexPermutationApplier)
+         for( GlobalIndexType i = 0; i < newSeeds.getEntitiesCount(); i++ ) {
+            const auto oldSeed = seeds.getSeed( i );
+            auto seed = newSeeds.getSeed( i );
+            for( LocalIndexType j = 0; j < seed.getCornersCount(); j++ ) {
+               const GlobalIndexType newIndex = points_old_to_new[ oldSeed.getCornerId( j ) ];
+               seed.setCornerId( j, newIndex );
+            }
+         }
+         seeds = std::move( newSeeds );
+      };
+
+      // remap points in this->faceSeeds or this->cellSeeds
+      if( faceSeeds.empty() )
+         remap_matrix( cellSeeds );
+      else
+         remap_matrix( faceSeeds );
+   }
+
    bool build( MeshType& mesh )
    {
       if( ! this->validate() )
          return false;
+      pointsSet.reset();
       mesh.init( this->points, this->faceSeeds, this->cellSeeds );
       return true;
    }
 
 private:
-   using PointArrayType     = typename MeshTraitsType::PointArrayType;
-   using BoolVector         = Containers::Vector< bool, Devices::Host, GlobalIndexType >;
-
    bool validate() const
    {
+      // verify that matrix dimensions are consistent with points
+      if( faceSeeds.empty() ) {
+         // no face seeds - cell seeds refer to points
+         if( cellSeeds.getMatrix().getColumns() != points.getSize() ) {
+            std::cerr << "Mesh builder error: Inconsistent size of the cellSeeds matrix (it has "
+                      << cellSeeds.getMatrix().getColumns() << " columns, but there are " << points.getSize()
+                      << " points)." << std::endl;
+            return false;
+         }
+      }
+      else {
+         // cell seeds refer to faces and face seeds refer to points
+         if( cellSeeds.getMatrix().getColumns() != faceSeeds.getMatrix().getRows() ) {
+            std::cerr << "Mesh builder error: Inconsistent size of the cellSeeds matrix (it has "
+                      << cellSeeds.getMatrix().getColumns() << " columns, but there are " << faceSeeds.getMatrix().getRows()
+                      << " faces)." << std::endl;
+            return false;
+         }
+         if( faceSeeds.getMatrix().getColumns() != points.getSize() ) {
+            std::cerr << "Mesh builder error: Inconsistent size of the faceSeeds matrix (it has "
+                      << faceSeeds.getMatrix().getColumns() << " columns, but there are " << points.getSize()
+                      << " points)." << std::endl;
+            return false;
+         }
+      }
+
       if( min( pointsSet ) != true ) {
          std::cerr << "Mesh builder error: Not all points were set." << std::endl;
          return false;