Added mesh ordering algorithms from tnl-mhfem (8e340584) · Commits · TNL / tnl-dev

src/TNL/Meshes/MeshOrdering.h

0 → 100644

+50 −0

Original line number	Diff line number	Diff line
		// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
		//
		// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
		//
		// SPDX-License-Identifier: MIT

		// Implemented by: Jakub Klinkovský

		#pragma once

		#include <TNL/Meshes/Mesh.h>

		namespace TNL {
		namespace Meshes {

		// general implementation covering grids
		template< typename Mesh, typename Ordering >
		struct MeshOrdering
		{
		void
		reorder( Mesh& mesh )
		{}
		};

		// reordering makes sense only for unstructured meshes
		template< typename MeshConfig, typename Device, typename Ordering >
		struct MeshOrdering< TNL::Meshes::Mesh< MeshConfig, Device >, Ordering >
		{
		using Mesh = TNL::Meshes::Mesh< MeshConfig, Device >;

		void
		reorder( Mesh& mesh )
		{
		Algorithms::staticFor< int, 0, Mesh::getMeshDimension() >(
		[ & ]( auto i )
		{
		// make sure to reorder cells first
		constexpr int dim = Mesh::getMeshDimension() - i;
		using PermutationArray = typename Mesh::GlobalIndexArray;
		PermutationArray perm;
		PermutationArray iperm;
		using EntityType = typename Mesh::template EntityType< dim >;
		Ordering::template getPermutations< EntityType >( mesh, perm, iperm );
		mesh.template reorderEntities< dim >( perm, iperm );
		} );
		}
		};

		} // namespace Meshes
		} // namespace TNL

src/TNL/Meshes/MeshOrdering/CuthillMcKeeOrdering.h

0 → 100644

+163 −0

Original line number	Diff line number	Diff line
		// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
		//
		// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
		//
		// SPDX-License-Identifier: MIT

		// Implemented by: Jakub Klinkovský

		#pragma once

		#include <type_traits>

		#include <TNL/Containers/Array.h>

		namespace TNL {
		namespace Meshes {

		template< bool reverse = true >
		struct CuthillMcKeeOrdering
		{
		template< typename MeshEntity, typename Mesh, typename PermutationArray >
		static void
		getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		static_assert( std::is_same< typename Mesh::DeviceType, TNL::Devices::Host >::value, "" );
		static_assert( std::is_same< typename PermutationArray::DeviceType, TNL::Devices::Host >::value, "" );

		// The reverse Cuthill-McKee ordering is implemented only for cells,
		// other entities are ordered from the current order of cells exactly
		// as if the mesh was re-initialized from the cell seeds.
		Wrapper< MeshEntity, Mesh >::getPermutations( mesh, perm, iperm );
		}

		private:
		template< typename MeshEntity, typename Mesh, bool is_cell = MeshEntity::getEntityDimension() == Mesh::getMeshDimension() >
		struct Wrapper
		{
		template< typename PermutationArray >
		static void
		getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		reorderCells( mesh, perm, iperm );
		}
		};

		template< typename MeshEntity, typename Mesh >
		struct Wrapper< MeshEntity, Mesh, false >
		{
		template< typename PermutationArray >
		static void
		getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		reorderEntities< MeshEntity >( mesh, perm, iperm );
		}
		};

		template< typename Mesh, typename PermutationArray >
		static void
		reorderCells( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		using IndexType = typename Mesh::GlobalIndexType;
		const IndexType numberOfCells = mesh.template getEntitiesCount< typename Mesh::Cell >();

		// allocate permutation vectors
		perm.setSize( numberOfCells );
		iperm.setSize( numberOfCells );

		// vector view for the neighbor counts
		const auto neighborCounts = mesh.getNeighborCounts().getConstView();

		// worker array - marker for inserted elements
		TNL::Containers::Array< bool, TNL::Devices::Host, IndexType > marker( numberOfCells );
		marker.setValue( false );
		// worker vector for collecting neighbors
		std::vector< IndexType > neighbors;

		// comparator functor
		auto comparator = [ & ]( IndexType a, IndexType b )
		{
		return neighborCounts[ a ] < neighborCounts[ b ];
		};

		// counter for assigning indices
		IndexType permIndex = 0;

		// modifier for the reversed variant
		auto mod = [ numberOfCells ]( IndexType i )
		{
		if( reverse )
		return numberOfCells - 1 - i;
		else
		return i;
		};

		// start with a peripheral node
		const IndexType peripheral = argMin( neighborCounts ).first;
		perm[ mod( permIndex ) ] = peripheral;
		iperm[ peripheral ] = mod( permIndex );
		permIndex++;
		marker[ peripheral ] = true;

		// Cuthill--McKee
		IndexType i = 0;
		while( permIndex < numberOfCells ) {
		const IndexType k = perm[ mod( i ) ];
		// collect all neighbors which were not marked yet
		const IndexType count = neighborCounts[ k ];
		for( IndexType n = 0; n < count; n++ ) {
		const IndexType nk = mesh.getCellNeighborIndex( k, n );
		if( marker[ nk ] == false ) {
		neighbors.push_back( nk );
		marker[ nk ] = true;
		}
		}
		// sort collected neighbors with ascending neighbors count
		std::sort( neighbors.begin(), neighbors.end(), comparator );
		// assign an index to the neighbors in this order
		for( auto nk : neighbors ) {
		perm[ mod( permIndex ) ] = nk;
		iperm[ nk ] = mod( permIndex );
		permIndex++;
		}
		// next iteration
		i++;
		neighbors.clear();
		}
		}

		template< typename MeshEntity, typename Mesh, typename PermutationArray >
		static void
		reorderEntities( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		using IndexType = typename Mesh::GlobalIndexType;
		const IndexType numberOfEntities = mesh.template getEntitiesCount< MeshEntity >();
		const IndexType numberOfCells = mesh.template getEntitiesCount< typename Mesh::Cell >();

		// allocate permutation vectors
		perm.setSize( numberOfEntities );
		iperm.setSize( numberOfEntities );

		// worker array - marker for numbered entities
		TNL::Containers::Array< bool, TNL::Devices::Host, IndexType > marker( numberOfEntities );
		marker.setValue( false );

		IndexType permIndex = 0;
		for( IndexType K = 0; K < numberOfCells; K++ ) {
		const auto& cell = mesh.template getEntity< Mesh::getMeshDimension() >( K );
		for( typename Mesh::LocalIndexType e = 0; e < cell.template getSubentitiesCount< MeshEntity::getEntityDimension() >();
		e++ ) {
		const auto E = cell.template getSubentityIndex< MeshEntity::getEntityDimension() >( e );
		if( marker[ E ] == false ) {
		marker[ E ] = true;
		perm[ permIndex ] = E;
		iperm[ E ] = permIndex;
		permIndex++;
		}
		}
		}
		}
		};

		} // namespace Meshes
		} // namespace TNL

src/TNL/Meshes/MeshOrdering/HilbertOrdering.h

0 → 100644

+108 −0

Original line number	Diff line number	Diff line
		// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
		//
		// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
		//
		// SPDX-License-Identifier: MIT

		// Implemented by: Jakub Klinkovský

		#pragma once

		#ifdef HAVE_CGAL

		#include <TNL/Meshes/Geometry/getEntityCenter.h>

		// Documentation:
		// - https://doc.cgal.org/latest/Spatial_sorting/index.html
		#include <CGAL/hilbert_sort.h>

		namespace TNL {
		namespace Meshes {

		struct HilbertOrdering
		{
		template< typename MeshEntity, typename Mesh, typename PermutationArray >
		static void
		getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		static_assert( std::is_same< typename Mesh::DeviceType, TNL::Devices::Host >::value, "" );
		static_assert( std::is_same< typename PermutationArray::DeviceType, TNL::Devices::Host >::value, "" );
		using GlobalIndexType = typename Mesh::GlobalIndexType;
		using PointType = typename Mesh::PointType;
		using RealType = typename Mesh::RealType;

		// wrappers for CGAL
		using pair = std::pair< PointType, GlobalIndexType >;
		struct Compute_d
		{
		RealType
		operator()( const pair& p, int d ) const
		{
		return p.first[ d ];
		}
		};
		struct Less_d
		{
		bool
		operator()( const pair& p, const pair& q, int d ) const
		{
		return p.first[ d ] < q.first[ d ];
		}
		};
		struct Point_dim_d
		{
		int
		operator()( const pair& p ) const
		{
		return p.first.getSize();
		}
		};
		struct SortingTraits
		{
		using Point_d = pair;
		using Compute_coordinate_d = Compute_d;
		using Less_coordinate_d = Less_d;
		using Point_dimension_d = Point_dim_d;
		Compute_coordinate_d
		compute_coordinate_d_object() const
		{
		return {};
		}
		Less_coordinate_d
		less_coordinate_d_object() const
		{
		return {};
		}
		Point_dimension_d
		point_dimension_d_object() const
		{
		return {};
		}
		};

		// create a vector with entity centers and initial indices
		std::vector< pair > points;
		const GlobalIndexType numberOfEntities = mesh.template getEntitiesCount< MeshEntity >();
		for( GlobalIndexType i = 0; i < numberOfEntities; i++ ) {
		const auto& entity = mesh.template getEntity< MeshEntity >( i );
		const auto center = getEntityCenter( mesh, entity );
		points.push_back( std::make_pair( center, i ) );
		}

		// sort the points
		CGAL::hilbert_sort( points.begin(), points.end(), SortingTraits(), CGAL::Hilbert_sort_middle_policy() );

		// build the permutations
		perm.setSize( numberOfEntities );
		iperm.setSize( numberOfEntities );
		for( GlobalIndexType i = 0; i < numberOfEntities; i++ ) {
		perm[ i ] = points[ i ].second;
		iperm[ points[ i ].second ] = i;
		}
		}
		};

		} // namespace Meshes
		} // namespace TNL

		#endif // HAVE_CGAL

src/TNL/Meshes/MeshOrdering/KdTreeOrdering.h

0 → 100644

+95 −0

Original line number	Diff line number	Diff line
		// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
		//
		// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
		//
		// SPDX-License-Identifier: MIT

		// Implemented by: Jakub Klinkovský

		#pragma once

		#ifdef HAVE_CGAL

		#include <TNL/Meshes/Geometry/getEntityCenter.h>

		// Documentation:
		// - https://doc.cgal.org/latest/Spatial_searching/index.html
		// - https://doc.cgal.org/latest/Spatial_searching/group__PkgSpatialSearchingDRef.html
		// - https://doc.cgal.org/latest/Spatial_searching/classCGAL_1_1Kd__tree.html
		#include <CGAL/Search_traits.h>
		#include <CGAL/Search_traits_adapter.h>
		#include <CGAL/Kd_tree.h>

		namespace TNL {
		namespace Meshes {

		struct KdTreeOrdering
		{
		template< typename MeshEntity, typename Mesh, typename PermutationArray >
		static void
		getPermutations( const Mesh& mesh, PermutationArray& perm, PermutationArray& iperm )
		{
		static_assert( std::is_same< typename Mesh::DeviceType, TNL::Devices::Host >::value, "" );
		static_assert( std::is_same< typename PermutationArray::DeviceType, TNL::Devices::Host >::value, "" );
		using GlobalIndexType = typename Mesh::GlobalIndexType;
		using PointType = typename Mesh::PointType;
		using RealType = typename Mesh::RealType;

		const GlobalIndexType numberOfEntities = mesh.template getEntitiesCount< MeshEntity >();

		// allocate permutation vectors
		perm.setSize( numberOfEntities );
		iperm.setSize( numberOfEntities );

		// wrapper for CGAL
		struct Construct_coord_iterator
		{
		using result_type = const RealType*;
		const RealType*
		operator()( const PointType& p ) const
		{
		return static_cast< const RealType* >( &p[ 0 ] );
		}
		const RealType*
		operator()( const PointType& p, int ) const
		{
		return static_cast< const RealType* >( &p[ 0 ] + p.getSize() );
		}
		};

		// CGAL types
		using DimTag = CGAL::Dimension_tag< PointType::getSize() >;
		using Point_and_index = std::tuple< PointType, GlobalIndexType >;
		using Traits_base = CGAL::Search_traits< RealType, PointType, const RealType*, Construct_coord_iterator, DimTag >;
		using TreeTraits =
		CGAL::Search_traits_adapter< Point_and_index, CGAL::Nth_of_tuple_property_map< 0, Point_and_index >, Traits_base >;
		// Note: splitters affect the quality of the ordering
		// available splitters are: https://doc.cgal.org/latest/Spatial_searching/classSplitter.html
		using Splitter = CGAL::Sliding_fair< TreeTraits >;
		using Tree = CGAL::Kd_tree< TreeTraits, Splitter >;

		// build a k-d tree
		Tree tree;
		for( GlobalIndexType i = 0; i < numberOfEntities; i++ ) {
		const auto& entity = mesh.template getEntity< MeshEntity >( i );
		const auto center = getEntityCenter( mesh, entity );
		tree.insert( std::make_tuple( center, i ) );
		}
		tree.build();

		GlobalIndexType permIndex = 0;

		// in-order traversal of the k-d tree
		for( auto iter : tree ) {
		const GlobalIndexType i = std::get< 1 >( iter );
		perm[ permIndex ] = i;
		iperm[ i ] = permIndex;
		permIndex++;
		}
		}
		};

		} // namespace Meshes
		} // namespace TNL

		#endif // HAVE_CGAL

src/Tools/CMakeLists.txt

+22 −8

Original line number	Diff line number	Diff line
		@@ -10,6 +10,8 @@ set( MESH_TOOLS
		tnl-planar-correct-mesh
		tnl-refine-mesh
		tnl-game-of-life
		tnl-reorder-mesh
		tnl-plot-mesh-ordering
		)
		set( CUDA_MESH_TOOLS )
		if( BUILD_CUDA )
		@@ -57,6 +59,18 @@ if( METIS_FOUND )
		install( TARGETS tnl-decompose-mesh DESTINATION bin )
		endif()

		set( CGAL_DATA_DIR "/" ) # to disable stupid warning: https://github.com/CGAL/cgal/pull/6649
		find_package( CGAL )
		if( CGAL_FOUND )
		option( CGAL_DO_NOT_WARN_ABOUT_CMAKE_BUILD_TYPE "disable warning" ON )

		target_compile_definitions( tnl-reorder-mesh PUBLIC "-DHAVE_CGAL" )
		target_link_libraries( tnl-reorder-mesh CGAL::CGAL )

		target_compile_definitions( tnl-plot-mesh-ordering PUBLIC "-DHAVE_CGAL" )
		target_link_libraries( tnl-plot-mesh-ordering CGAL::CGAL )
		endif()


		add_executable( tnl-grid-setup tnl-grid-setup.cpp )