Added mesh layer for the dual graph storage (56652119) · Commits · TNL / tnl-dev

src/TNL/Meshes/DefaultConfig.h

+16 −0

Original line number	Diff line number	Diff line
		@@ -87,6 +87,22 @@ struct DefaultConfig
		( EntityTopology::dimension >= meshDimension - 1 \|\| subentityStorage( FaceTopology(), EntityTopology::dimension ) );
		//return false;
		}

		/****
		* Storage of the dual graph.
		*
		* If enabled, links from vertices to cells must be stored.
		*/
		static constexpr bool dualGraphStorage()
		{
		return true;
		}

		/****
		* Cells must have at least this number of common vertices to be considered
		* as neighbors in the dual graph.
		*/
		static constexpr int dualGraphMinCommonVertices = meshDimension;
		};

		} // namespace Meshes

src/TNL/Meshes/Mesh.h

+9 −0

Original line number	Diff line number	Diff line
		@@ -165,6 +165,15 @@ class Mesh
		__cuda_callable__
		GlobalIndexType getSuperentityIndex( const GlobalIndexType entityIndex, const LocalIndexType superentityIndex ) const;

		/**
		* Cell neighbors - access the dual graph
		*/
		__cuda_callable__
		LocalIndexType getCellNeighborsCount( const GlobalIndexType cellIndex ) const;

		__cuda_callable__
		GlobalIndexType getCellNeighborIndex( const GlobalIndexType cellIndex, const LocalIndexType neighborIndex ) const;


		/*
		* The permutations follow the definition used in the Metis library: Let M

src/TNL/Meshes/Mesh.hpp

+41 −5

Original line number	Diff line number	Diff line
		@@ -29,10 +29,13 @@ MeshInitializableBase< MeshConfig, Device, MeshType >::
		init( typename MeshTraitsType::PointArrayType& points,
		typename MeshTraitsType::CellSeedArrayType& cellSeeds )
		{
		MeshType* mesh = static_cast< MeshType* >( this );
		Initializer< typename MeshType::Config > initializer;
		initializer.createMesh( points, cellSeeds, static_cast<MeshType>(this) );
		initializer.createMesh( points, cellSeeds, *mesh );
		// init boundary tags
		static_cast< BoundaryTags::LayerFamily< MeshConfig, Device, MeshType >* >( static_cast< MeshType* >( this ) )->initLayer();
		static_cast< BoundaryTags::LayerFamily< MeshConfig, Device, MeshType >* >( mesh )->initLayer();
		// init dual graph
		mesh->initializeDualGraph( *mesh );
		}


		@@ -205,6 +208,30 @@ getSuperentityIndex( const GlobalIndexType entityIndex, const LocalIndexType sup
		return this->template getSuperentityStorageNetwork< EntityDimension, SuperentityDimension >().getValue( entityIndex, superentityIndex );
		}

		template< typename MeshConfig, typename Device >
		__cuda_callable__
		typename Mesh< MeshConfig, Device >::LocalIndexType
		Mesh< MeshConfig, Device >::
		getCellNeighborsCount( const GlobalIndexType cellIndex ) const
		{
		static_assert( MeshConfig::dualGraphStorage(),
		"You try to access the dual graph which is disabled in the mesh configuration." );
		return this->getNeighborCounts()[ cellIndex ];
		}

		template< typename MeshConfig, typename Device >
		__cuda_callable__
		typename Mesh< MeshConfig, Device >::GlobalIndexType
		Mesh< MeshConfig, Device >::
		getCellNeighborIndex( const GlobalIndexType cellIndex, const LocalIndexType neighborIndex ) const
		{
		static_assert( MeshConfig::dualGraphStorage(),
		"You try to access the dual graph which is disabled in the mesh configuration." );
		TNL_ASSERT_GE( neighborIndex, 0, "Invalid cell neighbor index." );
		TNL_ASSERT_LT( neighborIndex, getCellNeighborsCount( cellIndex ), "Invalid cell neighbor index." );
		return this->getDualGraph().getValues( cellIndex )[ neighborIndex ];
		}


		template< typename MeshConfig, typename Device >
		template< int Dimension >
		@@ -253,9 +280,18 @@ void
		Mesh< MeshConfig, Device >::
		load( File& file )
		{
		// loading via host is necessary for the initialization of the dual graph
		if( std::is_same< Device, Devices::Cuda >::value ) {
		Mesh< MeshConfig, Devices::Host > hostMesh;
		hostMesh.load( file );
		*this = hostMesh;
		}
		else {
		Object::load( file );
		StorageBaseType::load( file );
		BoundaryTagsLayerFamily::load( file );
		this->initializeDualGraph( *this );
		}
		}

		template< typename MeshConfig, typename Device >

src/TNL/Meshes/MeshDetails/IndexPermutationApplier.h

+35 −24

Original line number	Diff line number	Diff line
		@@ -120,49 +120,55 @@ private:
		template< int Superdimension >
		using SuperentitiesWorker = IndexPermutationApplierSuperentitiesWorker< Superdimension >;

		public:
		static void permutePoints( Mesh& mesh,
		const GlobalIndexVector& perm,
		const GlobalIndexVector& iperm )
		template< typename Mesh_, std::enable_if_t< Mesh_::Config::dualGraphStorage(), bool > = true >
		static void permuteDualGraph( Mesh_& mesh, const GlobalIndexVector& perm, const GlobalIndexVector& iperm )
		{
		using IndexType = typename Mesh::GlobalIndexType;
		using DeviceType = typename Mesh::DeviceType;
		using PointArrayType = typename Mesh::MeshTraitsType::PointArrayType;
		permuteArray( mesh.getNeighborCounts(), perm );
		auto& graph = mesh.getDualGraph();
		Containers::Multimaps::permuteMultimapKeys( graph, perm );
		Containers::Multimaps::permuteMultimapValues( graph, iperm );
		}

		const IndexType pointsCount = mesh.template getEntitiesCount< 0 >();
		template< typename Mesh_, std::enable_if_t< ! Mesh_::Config::dualGraphStorage(), bool > = true >
		static void permuteDualGraph( Mesh_& mesh, const GlobalIndexVector& perm, const GlobalIndexVector& iperm ) {}

		PointArrayType points;
		points.setSize( pointsCount );
		public:
		template< typename Array >
		static void permuteArray( Array& array, const GlobalIndexVector& perm )
		{
		using IndexType = typename Array::IndexType;
		using DeviceType = typename Array::DeviceType;

		Array buffer( array.getSize() );

		// kernel to copy entities to new array, applying the permutation
		auto kernel1 = [] __cuda_callable__
		( IndexType i,
		const Mesh* mesh,
		typename PointArrayType::ValueType* pointsArray,
		const typename Array::ValueType* array,
		typename Array::ValueType* buffer,
		const IndexType* perm )
		{
		pointsArray[ i ] = mesh->getPoint( perm[ i ] );
		buffer[ i ] = array[ perm[ i ] ];
		};

		// kernel to copy permuted entities back to the mesh
		auto kernel2 = [] __cuda_callable__
		( IndexType i,
		Mesh* mesh,
		const typename PointArrayType::ValueType* pointsArray )
		typename Array::ValueType* array,
		const typename Array::ValueType* buffer )
		{
		mesh->getPoint( i ) = pointsArray[ i ];
		array[ i ] = buffer[ i ];
		};

		Pointers::DevicePointer< Mesh > meshPointer( mesh );
		Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, pointsCount,
		Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, array.getSize(),
		kernel1,
		&meshPointer.template getData< DeviceType >(),
		points.getData(),
		array.getData(),
		buffer.getData(),
		perm.getData() );
		Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, pointsCount,
		Algorithms::ParallelFor< DeviceType >::exec( (IndexType) 0, array.getSize(),
		kernel2,
		&meshPointer.template modifyData< DeviceType >(),
		points.getData() );
		array.getData(),
		buffer.getData() );
		}

		static void exec( Mesh& mesh,
		@@ -173,7 +179,7 @@ public:
		using DeviceType = typename Mesh::DeviceType;

		if( Dimension == 0 )
		permutePoints( mesh, perm, iperm );
		permuteArray( mesh.getPoints(), perm );

		// permute superentities storage
		Algorithms::TemplateStaticFor< int, 0, Dimension, SubentitiesStorageWorker >::execHost( mesh, perm );
		@@ -186,6 +192,11 @@ public:

		// update subentity indices from the superentities
		Algorithms::TemplateStaticFor< int, Dimension + 1, Mesh::getMeshDimension() + 1, SuperentitiesWorker >::execHost( mesh, iperm );

		if( Dimension == Mesh::getMeshDimension() ) {
		// permute dual graph
		permuteDualGraph( mesh, perm, iperm );
		}
		}
		};

src/TNL/Meshes/MeshDetails/layers/DualGraphLayer.h

0 → 100644

+178 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		DualGraphLayer.h - description
		-------------------
		begin : Apr 4, 2020
		copyright : (C) 2020 by Tomas Oberhuber et al.
		email : tomas.oberhuber@fjfi.cvut.cz
		***************************************************************************/

		/* See Copyright Notice in tnl/Copyright */

		// Implemented by: Jakub Klinkovský

		#pragma once

		#include <TNL/Meshes/MeshDetails/traits/MeshTraits.h>

		namespace TNL {
		namespace Meshes {

		template< typename MeshConfig,
		typename Device,
		bool enabled = MeshConfig::dualGraphStorage() >
		class DualGraphLayer
		{
		public:
		using MeshTraitsType = MeshTraits< MeshConfig, Device >;
		using GlobalIndexType = typename MeshTraitsType::GlobalIndexType;
		using LocalIndexType = typename MeshTraitsType::LocalIndexType;
		using DualGraph = typename MeshTraitsType::DualGraph;
		using NeighborCountsArray = typename DualGraph::ValuesAllocationVectorType;

		DualGraphLayer() = default;

		explicit DualGraphLayer( const DualGraphLayer& ) = default;

		template< typename Device_ >
		DualGraphLayer( const DualGraphLayer< MeshConfig, Device_ >& other )
		{
		operator=( other );
		}

		DualGraphLayer& operator=( const DualGraphLayer& ) = default;

		template< typename Device_ >
		DualGraphLayer& operator=( const DualGraphLayer< MeshConfig, Device_ >& other )
		{
		neighborCounts = other.getNeighborCounts();
		graph.setLike( other.getDualGraph() );
		graph = other.getDualGraph();
		return *this;
		}

		bool operator==( const DualGraphLayer& other ) const
		{
		return neighborCounts == other.getNeighborCounts() &&
		graph == other.getDualGraph();
		}

		__cuda_callable__
		const NeighborCountsArray& getNeighborCounts() const
		{
		return neighborCounts;
		}

		__cuda_callable__
		NeighborCountsArray& getNeighborCounts()
		{
		return neighborCounts;
		}

		__cuda_callable__
		const DualGraph& getDualGraph() const
		{
		return graph;
		}

		__cuda_callable__
		DualGraph& getDualGraph()
		{
		return graph;
		}

		// algorithm inspired by the CreateGraphDual function in METIS
		template< typename Mesh >
		void initializeDualGraph( const Mesh& mesh,
		// when this parameter is <= 0, it will be replaced with MeshConfig::dualGraphMinCommonVertices
		LocalIndexType minCommon = 0 )
		{
		static_assert( std::is_same< MeshConfig, typename Mesh::Config >::value,
		"mismatched MeshConfig type" );
		static_assert( MeshConfig::superentityStorage( typename Mesh::Vertex::EntityTopology{}, Mesh::getMeshDimension() ),
		"The dual graph cannot be initialized when links from vertices to cells are not stored in the mesh." );
		static_assert( MeshConfig::dualGraphMinCommonVertices >= 1,
		"MeshConfig error: dualGraphMinCommonVertices must be at least 1." );
		if( minCommon <= 0 )
		minCommon = MeshConfig::dualGraphMinCommonVertices;

		const GlobalIndexType cellsCount = mesh.template getEntitiesCount< Mesh::getMeshDimension() >();

		// allocate row lengths vector
		neighborCounts.setSize( cellsCount );

		// allocate working arrays
		using GlobalIndexArray = Containers::Array< GlobalIndexType, Devices::Sequential, GlobalIndexType >;
		using LocalIndexArray = Containers::Array< LocalIndexType, Devices::Sequential, GlobalIndexType >;
		GlobalIndexArray neighbors( cellsCount );
		LocalIndexArray marker( cellsCount );
		marker.setValue( 0 );

		auto findNeighbors = [&] ( const GlobalIndexType k )
		{
		const LocalIndexType subvertices = mesh.template getSubentitiesCount< Mesh::getMeshDimension(), 0 >( k );

		// find all elements that share at least one vertex with k
		LocalIndexType counter = 0;
		for( LocalIndexType v = 0; v < subvertices; v++ ) {
		const GlobalIndexType gv = mesh.template getSubentityIndex< Mesh::getMeshDimension(), 0 >( k, v );
		const LocalIndexType supercells = mesh.template getSuperentitiesCount< 0, Mesh::getMeshDimension() >( gv );
		for( LocalIndexType sc = 0; sc < supercells; sc++ ) {
		const GlobalIndexType nk = mesh.template getSuperentityIndex< 0, Mesh::getMeshDimension() >( gv, sc );
		if( marker[ nk ] == 0 )
		neighbors[ counter++ ] = nk;
		marker[ nk ]++;
		}
		}

		// mark k to be removed from the neighbor list in the next step
		marker[ k ] = 0;

		// compact the list to contain only those with at least minCommon vertices
		LocalIndexType compacted = 0;
		for( LocalIndexType i = 0; i < counter; i++ ) {
		const GlobalIndexType nk = neighbors[ i ];
		const LocalIndexType neighborSubvertices = mesh.template getSubentitiesCount< Mesh::getMeshDimension(), 0 >( nk );
		const LocalIndexType overlap = marker[ nk ];
		if( overlap >= minCommon \|\| overlap >= subvertices - 1 \|\| overlap >= neighborSubvertices - 1 )
		neighbors[ compacted++ ] = nk;
		marker[ nk ] = 0;
		}

		return compacted;
		};

		// count neighbors of each cell
		for( GlobalIndexType k = 0; k < cellsCount; k++ )
		neighborCounts[ k ] = findNeighbors( k );

		// allocate adjacency matrix
		graph.setKeysRange( cellsCount );
		graph.allocate( neighborCounts );

		// fill in neighbor indices
		for( GlobalIndexType k = 0; k < cellsCount; k++) {
		auto adjncy = graph.getValues( k );
		const LocalIndexType nnbrs = findNeighbors( k );
		for( LocalIndexType j = 0; j < nnbrs; j++)
		adjncy[ j ] = neighbors[ j ];
		}
		}

		protected:
		NeighborCountsArray neighborCounts;
		DualGraph graph;
		};

		template< typename MeshConfig,
		typename Device >
		class DualGraphLayer< MeshConfig, Device, false >
		{
		public:
		template< typename Mesh >
		void initializeDualGraph( const Mesh& mesh,
		int minCommon = 0 )
		{}
		};

		} // namespace Meshes
		} // namespace TNL