From b51661eb5c469b235da03b2dc3cc9c277f518cd5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jakub=20Klinkovsk=C3=BD?= <klinkovsky@mmg.fjfi.cvut.cz>
Date: Thu, 30 Dec 2021 12:34:32 +0100
Subject: [PATCH] Added function to deduplicate points in MeshBuilder
---
src/TNL/Meshes/MeshBuilder.h | 131 +++++++++++++++++++++++++++++++++--
1 file changed, 126 insertions(+), 5 deletions(-)
diff --git a/src/TNL/Meshes/MeshBuilder.h b/src/TNL/Meshes/MeshBuilder.h
index aa10a8a245..e3676484b4 100644
--- a/src/TNL/Meshes/MeshBuilder.h
+++ b/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;
--
GitLab