Added VTKReader for testing (using the VTK library) (268acecd) · Commits · TNL / tnl-dev

src/TNL/Meshes/Readers/CMakeLists.txt

+1 −0

Original line number	Diff line number	Diff line
		SET( headers NetgenReader.h
		TNLReader.h
		VTKEntityType.h
		VTKReader.h
		)

		INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/TNL/Meshes/Readers )

src/TNL/Meshes/Readers/VTKReader.h

0 → 100644

+339 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		VTKReader.h - description
		-------------------
		begin : Nov 6, 2016
		copyright : (C) 2016 by Tomas Oberhuber et al.
		email : tomas.oberhuber@fjfi.cvut.cz
		***************************************************************************/

		/* See Copyright Notice in tnl/Copyright */

		#pragma once

		#include <fstream>
		#include <istream>
		#include <vector>
		#include <map>
		#include <unordered_map>
		#include <type_traits>

		#include <TNL/Containers/StaticVector.h>
		#include <TNL/Meshes/MeshBuilder.h>
		#include <TNL/Meshes/Topologies/MeshSubtopology.h>
		#include <TNL/Meshes/Readers/VTKEntityType.h>

		#ifdef HAVE_VTK
		#include <vtk/vtkSmartPointer.h>
		#include <vtk/vtkCell.h>
		#include <vtk/vtkUnstructuredGrid.h>
		#include <vtk/vtkUnstructuredGridReader.h>
		#endif

		namespace TNL {
		namespace Meshes {
		namespace Readers {

		// types used in VTK
		using VTKRealType = double;
		#ifdef HAVE_VTK
		using VTKIndexType = vtkIdType;
		#else
		using VTKIndexType = int;
		#endif
		// wrapper type for physical coordinates to preserve my sanity
		using VTKPointType = Containers::StaticVector< 3, VTKRealType >;

		template< typename Index = int >
		class VTKReader
		{
		public:
		using IndexType = Index;

		bool
		detectMesh( const String& fileName )
		{
		resetAll();
		if( ! loadVTKFile( fileName ) ) {
		resetAll();
		return false;
		}
		reset();
		return true;
		}

		template< typename MeshType >
		bool
		readMesh( const String& fileName, MeshType& mesh )
		{
		static_assert( std::is_same< IndexType, typename MeshType::GlobalIndexType >::value, "VTKReader::IndexType and MeshType::GlobalIndexType must be the same type." );

		resetAll();
		if( ! loadVTKFile( fileName ) ) {
		resetAll();
		return false;
		}

		TNL_ASSERT( this->worldDimension == MeshType::Config::worldDimension, );
		TNL_ASSERT( this->meshDimension == MeshType::Config::meshDimension, );
		const int subvertices = MeshSubtopology< typename MeshType::Config::CellTopology, 0 >::count;
		TNL_ASSERT( this->verticesInEntities.at( this->meshDimension ) == subvertices, );

		using MeshBuilder = MeshBuilder< MeshType >;
		using VertexIndexType = typename MeshType::MeshTraitsType::template EntityTraits< 0 >::GlobalIndexType;
		using CellIndexType = typename MeshType::MeshTraitsType::template EntityTraits< MeshType::Config::meshDimension >::GlobalIndexType;

		const VertexIndexType numberOfPoints = this->pointsData.size();
		const CellIndexType numberOfCells = this->entityIdMappings.at( this->meshDimension ).size();

		MeshBuilder meshBuilder;

		if( ! meshBuilder.setPointsCount( numberOfPoints ) ) {
		std::cerr << "Failed to allocate memory for " << numberOfPoints << " vertices." << std::endl;
		return false;
		}

		for( VertexIndexType i = 0; i < numberOfPoints; i++ ) {
		typename MeshType::PointType p;
		for( int j = 0; j < p.size; j++ )
		p[ j ] = this->pointsData.at( i )[ j ];
		meshBuilder.setPoint( i, p );
		}

		if( ! meshBuilder.setCellsCount( numberOfCells ) ) {
		std::cerr << "Failed to allocate memory for " << numberOfCells << " cells." << std::endl;
		return false;
		}

		const auto& cellIdMap = this->entityIdMappings.at( this->meshDimension );
		for( CellIndexType i = 0; i < numberOfCells; i++ ) {
		const VTKIndexType vtkCellIndex = cellIdMap.at( i );
		const auto& vtkCellSeeds = this->entitySeeds.at( vtkCellIndex );
		using CellSeedType = typename MeshBuilder::CellSeedType;
		TNL_ASSERT( CellSeedType::getCornersCount() == vtkCellSeeds.size(),
		std::cerr << "CellSeedType::getCornersCount() = " << CellSeedType::getCornersCount()
		<< ", vtkCellseeds.size() = " << vtkCellSeeds.size()
		<< std::endl; );
		CellSeedType& seed = meshBuilder.getCellSeed( i );
		for( int v = 0; v < CellSeedType::getCornersCount(); v++ ) {
		seed.setCornerId( v, vtkCellSeeds[ v ] );
		}
		}

		// drop all data from the VTK file since it's not needed anymore
		this->resetAll();

		return meshBuilder.build( mesh );
		}

		String
		getMeshType() const
		{
		// we can read only the UNSTRUCTURED_GRID dataset
		return "Meshes::Mesh";
		}

		int
		getWorldDimension() const
		{
		return this->worldDimension;
		}

		int
		getMeshDimension() const
		{
		return this->meshDimension;
		}

		VTKEntityType
		getCellVTKType() const
		{
		return this->entityTypes.at( this->meshDimension );
		}

		// int
		// getVerticesInCell() const
		// {
		// return this->verticesInEntities.at( this->getMeshDimension() );
		// }

		// VTKEntityType
		// getEntityType( int entityDimension ) const
		// {
		// return this->entityTypes.at( entityDimension );
		// }

		String
		getRealType() const
		{
		// TODO: how to extract it from the VTK object?
		// return "float";
		return "double";
		}

		String
		getGlobalIndexType() const
		{
		// not stored in the VTK file
		return "int";
		}

		String
		getLocalIndexType() const
		{
		// not stored in the VTK file
		return "short int";
		}

		String
		getIdType() const
		{
		// not stored in the VTK file
		return "int";
		}

		protected:
		int worldDimension = 0;
		int meshDimension = 0;

		// maps vertex indices to physical coordinates
		std::unordered_map< VTKIndexType, VTKPointType > pointsData;

		// maps entity dimension to the number of vertices in the entity
		std::map< int, int > verticesInEntities;

		// type for mapping TNL entity indices to VTK cell indices
		using TNL2VTKindexmap = std::unordered_map< IndexType, VTKIndexType >;
		// maps dimension to maps of TNL entity IDs to corresponding VTK cell IDs
		std::unordered_map< int, TNL2VTKindexmap > entityIdMappings;

		// maps VTK cell indices to entity seeds (set of indices of subvertices)
		std::unordered_map< VTKIndexType, std::vector< VTKIndexType > > entitySeeds;

		// maps dimension to VTK type of the entity with given dimension
		std::unordered_map< int, VTKEntityType > entityTypes;

		void reset()
		{
		pointsData.clear();
		verticesInEntities.clear();
		entityIdMappings.clear();
		entitySeeds.clear();
		}

		void resetAll()
		{
		reset();
		worldDimension = 0;
		meshDimension = 0;
		entityTypes.clear();
		}

		bool loadVTKFile( const String& fileName )
		{
		#ifdef HAVE_VTK
		vtkSmartPointer< vtkUnstructuredGridReader > vtkReader
		= vtkSmartPointer<vtkUnstructuredGridReader>::New();
		vtkReader->SetFileName( fileName.getString() );
		vtkReader->Update();

		if( ! vtkReader->IsFileUnstructuredGrid() ) {
		std::cerr << "The file '" << fileName << "' is not a VTK Legacy file with an UNSTRUCTURED_GRID dataset type." << std::endl;
		return false;
		}

		auto& vtkMesh = *vtkReader->GetOutput();

		/* To determine the mesh type, we need:
		* - world dimension (default 1D, if some 2nd coordinate is non-zero -> 2D, some 3rd coordinate non-zero -> 3D)
		* - mesh dimension (highest dimension of the entities present in mesh)
		* - check that it is homogeneous - all entities of the same dimension should have the same number of vertices and type
		* - types of entities for each dimension
		*
		* To initialize the mesh, we need:
		* - points data (vertex indices and world coordinates)
		* - cell seeds (indices of vertices composing each cell)
		*/
		const VTKIndexType numberOfPoints = vtkMesh.GetNumberOfPoints();
		const VTKIndexType numberOfCells = vtkMesh.GetNumberOfCells();

		if( numberOfPoints == 0 ) {
		std::cerr << "There are no points data in the file '" << fileName << "'." << std::endl;
		return false;
		}

		if( numberOfCells == 0 ) {
		std::cerr << "There are no cells data in the file '" << fileName << "'." << std::endl;
		return false;
		}

		for( VTKIndexType i = 0; i < numberOfPoints; i++ ) {
		VTKRealType* p = vtkMesh.GetPoint( i );

		// get world dimension
		for( int j = 0; j < 3; j++ )
		if( p[ j ] != 0.0 )
		this->worldDimension = std::max( this->worldDimension, j + 1 );

		// copy points data
		this->pointsData[ i ] = VTKPointType( p[ 0 ], p[ 1 ], p[ 2 ] );
		}

		for( VTKIndexType i = 0; i < numberOfCells; i++ ) {
		vtkCell* cell = vtkMesh.GetCell( i );
		const int dimension = cell->GetCellDimension();
		const int points = cell->GetNumberOfPoints();
		const VTKEntityType type = (VTKEntityType) cell->GetCellType();

		// number of vertices in entities
		if( this->verticesInEntities.find( dimension ) == this->verticesInEntities.cend() )
		this->verticesInEntities[ dimension ] = points;
		else if( this->verticesInEntities[ dimension ] != points ) {
		std::cerr << "Mixed unstructured meshes are not supported. There are elements of dimension " << dimension
		<< " with " << this->verticesInEntities[ dimension ] << " vertices and with " << points
		<< " vertices. The number of vertices per entity must be constant." << std::endl;
		this->reset();
		return false;
		}

		// entity types
		if( this->entityTypes.find( dimension ) == this->entityTypes.cend() )
		this->entityTypes.emplace( std::make_pair( dimension, type ) );
		else if( this->entityTypes[ dimension ] != type ) {
		std::cerr << "Mixed unstructured meshes are not supported. There are elements of dimension " << dimension
		<< " with type " << this->entityTypes[ dimension ] << " and " << type
		<< ". The type of all entities with the same dimension must be the same." << std::endl;
		this->reset();
		return false;
		}

		// mapping between TNL and VTK indices
		auto& map = this->entityIdMappings[ dimension ];
		map[ map.size() ] = i;

		// copy seed
		auto& seed = this->entitySeeds[ i ];
		for( int j = 0; j < points; j++ )
		seed.push_back( cell->GetPointId( j ) );
		}

		// std::map is sorted, so the last key is the maximum
		this->meshDimension = this->verticesInEntities.rbegin()->first;

		if( this->meshDimension > this->worldDimension ) {
		std::cerr << "Invalid mesh: world dimension is " << this->worldDimension
		<< ", but mesh dimension is " << this->meshDimension << "." << std::endl;
		this->reset();
		return false;
		}

		return true;
		#else
		std::cerr << "The VTKReader needs to be compiled with the VTK library." << std::endl;
		return false;
		#endif
		}
		};

		} // namespace Readers
		} // namespace Meshes
		} // namespace TNL

src/TNL/Meshes/TypeResolver/TypeResolver_impl.h

+18 −5

Original line number	Diff line number	Diff line
		@@ -16,6 +16,7 @@
		#include <TNL/Meshes/TypeResolver/TypeResolver.h>
		#include <TNL/Meshes/Readers/TNLReader.h>
		#include <TNL/Meshes/Readers/NetgenReader.h>
		#include <TNL/Meshes/Readers/VTKReader.h>
		#include <TNL/Meshes/TypeResolver/GridTypeResolver.h>
		#include <TNL/Meshes/TypeResolver/MeshTypeResolver.h>

		@@ -56,7 +57,7 @@ template< typename ConfigTag,
		bool resolveMeshType( const String& fileName_,
		ProblemSetterArgs&&... problemSetterArgs )
		{
		std::cout << "Loading mesh from file " << fileName_ << " ..." << std::endl;
		std::cout << "Detecting mesh from file " << fileName_ << " ..." << std::endl;
		std::string fileName( fileName_.getString() );
		if( ends_with( fileName, ".tnl" ) ) {
		Readers::TNLReader reader;
		@@ -88,11 +89,23 @@ bool resolveMeshType( const String& fileName_,
		return false;
		}
		}
		// TODO
		// else if( ends_with( fileName, ".vtk" ) ) {
		// }
		else if( ends_with( fileName, ".vtk" ) ) {
		// FIXME: The VTK files don't store the global index, local index and id types.
		// The reader has some defaults, but they might be disabled by the BuildConfigTags - in
		// this case we should use the first enabled type.
		Readers::VTKReader<> reader;
		if( ! reader.detectMesh( fileName_ ) )
		return false;
		if( reader.getMeshType() == "Meshes::Mesh" )
		return MeshTypeResolver< decltype(reader), ConfigTag, Device, ProblemSetter, ProblemSetterArgs... >::
		run( reader, std::forward<ProblemSetterArgs>(problemSetterArgs)... );
		else {
		std::cerr << "The mesh type " << reader.getMeshType() << " is not supported in the VTK reader." << std::endl;
		return false;
		}
		}
		else {
		std::cerr << "File '" << fileName << "' has unknown extension. Supported extensions are '.tnl' and '.ng'." << std::endl;
		std::cerr << "File '" << fileName << "' has unknown extension. Supported extensions are '.tnl', '.vtk' and '.ng'." << std::endl;
		return false;
		}
		}