Refactoring mesh readers: created base class MeshReader (597022eb) · Commits · TNL / tnl-dev

src/Python/pytnl/tnl/Mesh.cpp

+9 −9

Original line number	Diff line number	Diff line
		@@ -14,17 +14,17 @@ void export_Meshes( py::module & m )

		py::class_< Reader >( m, "VTKReader" )
		.def(py::init<std::string>())
		.def("readMesh", &Reader::template readMesh< MeshOfEdges >)
		.def("readMesh", &Reader::template readMesh< MeshOfTriangles >)
		.def("readMesh", &Reader::template readMesh< MeshOfTetrahedrons >)
		// .def("readMesh", []( Reader& reader, const std::string& name, MeshOfEdges & mesh ) {
		// return reader.readMesh( name.c_str(), mesh );
		.def("loadMesh", &Reader::template loadMesh< MeshOfEdges >)
		.def("loadMesh", &Reader::template loadMesh< MeshOfTriangles >)
		.def("loadMesh", &Reader::template loadMesh< MeshOfTetrahedrons >)
		// .def("loadMesh", []( Reader& reader, const std::string& name, MeshOfEdges & mesh ) {
		// return reader.loadMesh( name.c_str(), mesh );
		// } )
		// .def("readMesh", []( Reader& reader, const std::string& name, MeshOfTriangles & mesh ) {
		// return reader.readMesh( name.c_str(), mesh );
		// .def("loadMesh", []( Reader& reader, const std::string& name, MeshOfTriangles & mesh ) {
		// return reader.loadMesh( name.c_str(), mesh );
		// } )
		// .def("readMesh", []( Reader& reader, const std::string& name, MeshOfTetrahedrons & mesh ) {
		// return reader.readMesh( name.c_str(), mesh );
		// .def("loadMesh", []( Reader& reader, const std::string& name, MeshOfTetrahedrons & mesh ) {
		// return reader.loadMesh( name.c_str(), mesh );
		// } )
		;
		}

src/TNL/Meshes/Readers/MeshReader.h

0 → 100644

+220 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		MeshReader.h - description
		-------------------
		begin : Apr 10, 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 <string>
		#include <vector>
		#include <TNL/variant.hpp>

		#include <TNL/Meshes/MeshBuilder.h>
		#include <TNL/Meshes/VTKTraits.h>

		namespace TNL {
		namespace Meshes {
		namespace Readers {

		struct MeshReaderError
		: public std::runtime_error
		{
		MeshReaderError( std::string readerName, std::string msg )
		: std::runtime_error( readerName + " error: " + msg )
		{}
		};

		class MeshReader
		{
		public:
		using VariantVector = mpark::variant< std::vector< std::int8_t >,
		std::vector< std::uint8_t >,
		std::vector< std::int16_t >,
		std::vector< std::uint16_t >,
		std::vector< std::int32_t >,
		std::vector< std::uint32_t >,
		std::vector< std::int64_t >,
		std::vector< std::uint64_t >,
		std::vector< float >,
		std::vector< double > >;

		MeshReader() = delete;

		MeshReader( const std::string& fileName )
		: fileName( fileName )
		{}

		virtual ~MeshReader() {}

		/**
		* \brief This method resets the reader to an empty state.
		*
		* In particular, implementations should call the \ref resetBase method to
		* reset the arrays holding the intermediate mesh representation.
		*/
		virtual void reset()
		{
		resetBase();
		}

		/**
		* \brief Main method responsible for reading the mesh file.
		*
		* The implementation has to set all protected attributes of this class such
		* that the mesh representation can be loaded into the mesh object by the
		* \ref loadMesh method.
		*/
		virtual void detectMesh() = 0;

		/**
		* \brief Method which loads the intermediate mesh representation into a
		* mesh object.
		*
		* When the method exits, the intermediate mesh representation is destroyed
		* to save memory. However, depending on the specific file format, the mesh
		* file may remain open so that the user can load additional data.
		*/
		template< typename MeshType >
		void loadMesh( MeshType& mesh )
		{
		// check that detectMesh has been called
		if( ! meshDetected )
		detectMesh();

		// check that the cell shape mathes
		const VTK::EntityShape meshCellShape = VTK::TopologyToEntityShape< typename MeshType::template EntityTraits< MeshType::getMeshDimension() >::EntityTopology >::shape;
		if( meshCellShape != cellShape )
		throw MeshReaderError( "MeshReader", "the mesh cell shape " + VTK::getShapeName(meshCellShape) + " does not match the shape "
		+ "of cells used in the file (" + VTK::getShapeName(cellShape) + ")" );

		using MeshBuilder = MeshBuilder< MeshType >;
		using IndexType = typename MeshType::GlobalIndexType;
		using PointType = typename MeshType::PointType;
		using CellSeedType = typename MeshBuilder::CellSeedType;

		MeshBuilder meshBuilder;
		meshBuilder.setPointsCount( NumberOfPoints );
		meshBuilder.setCellsCount( NumberOfCells );

		// assign points
		visit( [&meshBuilder](auto&& array) {
		PointType p;
		std::size_t i = 0;
		for( auto c : array ) {
		int dim = i++ % 3;
		if( dim >= PointType::getSize() )
		continue;
		p[dim] = c;
		if( dim == PointType::getSize() - 1 )
		meshBuilder.setPoint( (i - 1) / 3, p );
		}
		},
		pointsArray
		);

		// assign cells
		visit( [this, &meshBuilder](auto&& connectivity) {
		// let's just assume that the connectivity and offsets arrays have the same type...
		using mpark::get;
		const auto& offsets = get< std::decay_t<decltype(connectivity)> >( offsetsArray );
		std::size_t offsetStart = 0;
		for( std::size_t i = 0; i < NumberOfCells; i++ ) {
		CellSeedType& seed = meshBuilder.getCellSeed( i );
		const std::size_t offsetEnd = offsets[ i ];
		for( std::size_t o = offsetStart; o < offsetEnd; o++ )
		seed.setCornerId( o - offsetStart, connectivity[ o ] );
		offsetStart = offsetEnd;
		}
		},
		connectivityArray
		);

		// reset arrays since they are not needed anymore
		pointsArray = connectivityArray = offsetsArray = typesArray = {};

		if( ! meshBuilder.build( mesh ) )
		throw MeshReaderError( "VTKReader", "MeshBuilder failed" );
		}

		std::string
		getMeshType() const
		{
		return "Meshes::Mesh";
		}

		int
		getMeshDimension() const
		{
		return meshDimension;
		}

		int
		getWorldDimension() const
		{
		return worldDimension;
		}

		VTK::EntityShape
		getCellShape() const
		{
		return cellShape;
		}

		std::string
		getRealType() const
		{
		return pointsType;
		}

		std::string
		getGlobalIndexType() const
		{
		return connectivityType;
		}

		std::string
		getLocalIndexType() const
		{
		// not stored in any file format
		return "short int";
		}

		protected:
		// input file name
		std::string fileName;

		// indicator that detectMesh has been successfully called
		bool meshDetected = false;

		// attributes of the mesh
		std::size_t NumberOfPoints, NumberOfCells;
		int meshDimension, worldDimension;
		VTK::EntityShape cellShape = VTK::EntityShape::Vertex;

		// intermediate representation of the unstructured mesh (matches the VTU
		// file format, other formats have to be converted)
		VariantVector pointsArray, connectivityArray, offsetsArray, typesArray;
		// string representation of each array's value type
		std::string pointsType, connectivityType, offsetsType, typesType;

		void resetBase()
		{
		meshDetected = false;
		NumberOfPoints = NumberOfCells = 0;
		meshDimension = worldDimension = 0;
		cellShape = VTK::EntityShape::Vertex;
		pointsArray = connectivityArray = offsetsArray = typesArray = {};
		pointsType = connectivityType = offsetsType = typesType = "";
		}
		};

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

src/TNL/Meshes/Readers/NetgenReader.h

+94 −240

Original line number	Diff line number	Diff line
		@@ -10,296 +10,150 @@

		/***
		* Authors:
		* Oberhuber Tomas, tomas.oberhuber@fjfi.cvut.cz
		* Zabka Vitezslav, zabkav@gmail.com
		* Tomas Oberhuber
		* Vitezslav Zabka
		* Jakub Klinkovsky
		*/

		#pragma once

		#include <fstream>
		#include <istream>
		#include <sstream>

		#include <TNL/Meshes/MeshBuilder.h>
		#include <TNL/Meshes/VTKTraits.h>
		#include <TNL/Meshes/Readers/MeshReader.h>

		namespace TNL {
		namespace Meshes {
		namespace Readers {

		class NetgenReader
		: public MeshReader
		{
		public:
		NetgenReader() = delete;

		NetgenReader( const String& fileName )
		: fileName( fileName )
		NetgenReader( const std::string& fileName )
		: MeshReader( fileName )
		{}

		bool detectMesh()
		virtual void detectMesh() override
		{
		this->reset();
		reset();

		std::ifstream inputFile( fileName.getString() );
		std::ifstream inputFile( fileName );
		if( ! inputFile )
		{
		std::cerr << "I am not able to open the file " << fileName << "." << std::endl;
		return false;
		}
		throw MeshReaderError( "NetgenReader", "failed to open the file '" + fileName + "'." );

		std::string line;
		std::istringstream iss;

		/****
		* Skip whitespaces
		*/
		// skip whitespace
		inputFile >> std::ws;

		/****
		* Skip number of vertices
		*/
		if( ! inputFile )
		return false;
		getline( inputFile, line );
		iss.str( line );
		long int numberOfVertices;
		iss >> numberOfVertices;

		//cout << "There are " << numberOfVertices << " vertices." << std::endl;
		throw MeshReaderError( "NetgenReader", "unexpected error when reading the file '" + fileName + "'." );

		/****
		* Read the first vertex and compute number of components
		*/
		if( ! inputFile )
		return false;
		// read number of points
		getline( inputFile, line );
		iss.clear();
		iss.str( line );
		meshDimension = worldDimension = -1;
		while( iss )
		{
		double aux;
		iss >> aux;
		meshDimension = ++worldDimension;
		iss >> NumberOfPoints;

		// real type is not stored in Netgen files
		pointsType = "double";
		// global index type is not stored in Netgen files
		connectivityType = offsetsType = "std::int32_t";
		// only std::uint8_t makes sense for entity types
		typesType = "std::uint8_t";

		// arrays holding the data from the file
		std::vector< double > pointsArray;
		std::vector< std::int32_t > connectivityArray, offsetsArray;
		std::vector< std::uint8_t > typesArray;

		// read points
		worldDimension = 0;
		for( std::size_t pointIndex = 0; pointIndex < NumberOfPoints; pointIndex++ ) {
		if( ! inputFile ) {
		reset();
		throw MeshReaderError( "NetgenReader", "unable to read enough vertices, the file may be invalid or corrupted." );
		}

		/****
		* Skip vertices
		*/
		long int verticesRead( 1 );
		while( verticesRead < numberOfVertices )
		{
		getline( inputFile, line );
		if( ! inputFile )
		{
		std::cerr << "The mesh file " << fileName << " is probably corrupted, some vertices are missing." << std::endl;
		return false;
		}
		verticesRead++;
		}

		/****
		* Skip whitespaces
		*/
		inputFile >> std::ws;

		/****
		* Get number of cells
		*/
		long int numberOfCells;
		getline( inputFile, line );
		iss.clear();
		iss.str( line );
		iss >> numberOfCells;
		//cout << "There are " << numberOfCells << " cells." << std::endl;

		/****
		* Get number of vertices in a cell
		*/
		getline( inputFile, line );
		// read the coordinates and compute the world dimension
		iss.clear();
		iss.str( line );
		int verticesInCell = -2;
		while( iss )
		{
		int aux;
		for( int i = 0; i < 3; i++ ) {
		double aux;
		iss >> aux;
		verticesInCell++;
		if( ! iss ) {
		// the intermediate mesh representation uses the VTK convention - all points must have 3 coordinates
		aux = 0;
		}
		if( aux != 0.0 )
		worldDimension = std::max( worldDimension, i + 1 );
		pointsArray.push_back( aux );
		}
		}
		//cout << "There are " << verticesInCell << " vertices in cell ..." << std::endl;

		if( meshDimension == 1 && verticesInCell == 2 )
		// netgen supports only triangular and tetrahedral meshes
		meshDimension = worldDimension;
		if( meshDimension == 1 )
		cellShape = VTK::EntityShape::Line;
		else if( meshDimension == 2 ) {
		if( verticesInCell == 3 )
		else if( meshDimension == 2 )
		cellShape = VTK::EntityShape::Triangle;
		else if( verticesInCell == 4 )
		cellShape = VTK::EntityShape::Quad;
		}
		else if( meshDimension == 3 ) {
		if( verticesInCell == 4 )
		else if( meshDimension == 3 )
		cellShape = VTK::EntityShape::Tetra;
		else if( verticesInCell == 8 )
		cellShape = VTK::EntityShape::Hexahedron;
		}
		if( cellShape == VTK::EntityShape::Vertex ) {
		std::cerr << "Unknown cell topology: mesh dimension is " << meshDimension << ", number of vertices in cells is " << verticesInCell << "." << std::endl;
		return false;
		}

		meshDetected = true;
		return true;
		}
		else
		throw MeshReaderError( "NetgenReader", "unsupported mesh dimension: " + std::to_string(meshDimension) );

		template< typename MeshType >
		bool readMesh( MeshType& mesh )
		{
		// check that detectMesh has been called
		if( ! meshDetected )
		detectMesh();

		typedef typename MeshType::PointType PointType;
		typedef MeshBuilder< MeshType > MeshBuilder;

		const int dimension = PointType::getSize();

		// TODO: reuse inputFile from the detectMesh method
		std::ifstream inputFile( fileName.getString() );
		if( ! inputFile )
		{
		std::cerr << "I am not able to open the file " << fileName << "." << std::endl;
		return false;
		}

		MeshBuilder meshBuilder;
		std::string line;
		std::istringstream iss;

		/****
		* Skip white spaces
		*/
		// skip whitespace
		inputFile >> std::ws;

		/****
		* Read the number of vertices
		*/
		if( ! inputFile )
		return false;
		getline( inputFile, line );
		iss.str( line );
		typedef typename MeshType::Config::GlobalIndexType VertexIndexType;
		VertexIndexType pointsCount;
		iss >> pointsCount;
		meshBuilder.setPointsCount( pointsCount );
		throw MeshReaderError( "NetgenReader", "unexpected error when reading the file '" + fileName + "'." );

		for( VertexIndexType i = 0; i < pointsCount; i++ )
		{
		// read number of cells
		getline( inputFile, line );
		iss.clear();
		iss.str( line );
		PointType p;
		for( int d = 0; d < dimension; d++ )
		iss >> p[ d ];
		//cout << "Setting point number " << i << " of " << pointsCount << std::endl;
		meshBuilder.setPoint( i, p );
		//const PointType& point = mesh.getVertex( i ).getPoint();
		}
		iss >> NumberOfCells;

		/****
		* Skip white spaces
		*/
		inputFile >> std::ws;

		/****
		* Read number of cells
		*/
		typedef typename MeshType::Config::GlobalIndexType CellIndexType;
		if( ! inputFile )
		{
		std::cerr << "I cannot read the mesh cells." << std::endl;
		return false;
		// read cells
		for( std::size_t cellIndex = 0; cellIndex < NumberOfCells; cellIndex++ ) {
		if( ! inputFile ) {
		reset();
		throw MeshReaderError( "NetgenReader", "unable to read enough cells, the file may be invalid or corrupted."
		" (cellIndex = " + std::to_string(cellIndex) + ")" );
		}
		getline( inputFile, line );

		iss.clear();
		iss.str( line );
		CellIndexType numberOfCells = atoi( line.data() );
		//iss >> numberOfCells; // TODO: I do not know why this does not work
		meshBuilder.setCellsCount( numberOfCells );
		for( CellIndexType i = 0; i < numberOfCells; i++ )
		{
		getline( inputFile, line );
		iss.clear();
		iss.str( line );
		int subdomainIndex;
		iss >> subdomainIndex;
		//cout << "Setting cell number " << i << " of " << numberOfCells << std::endl;
		typedef typename MeshBuilder::CellSeedType CellSeedType;
		for( int cellVertex = 0; cellVertex < CellSeedType::getCornersCount(); cellVertex++ )
		{
		VertexIndexType vertexIdx;
		iss >> vertexIdx;
		meshBuilder.getCellSeed( i ).setCornerId( cellVertex, vertexIdx - 1 );
		}
		}
		meshBuilder.build( mesh );
		return true;
		// skip subdomain number
		int subdomain;
		iss >> subdomain;
		for( int v = 0; v <= meshDimension; v++ ) {
		std::size_t vid;
		iss >> vid;
		if( ! iss ) {
		reset();
		throw MeshReaderError( "NetgenReader", "unable to read enough cells, the file may be invalid or corrupted."
		" (cellIndex = " + std::to_string(cellIndex) + ", subvertex = " + std::to_string(v) + ")" );
		}

		String
		getMeshType() const
		{
		return "Meshes::Mesh";
		// convert point index from 1-based to 0-based
		connectivityArray.push_back( vid - 1 );
		}

		int getMeshDimension() const
		{
		return this->meshDimension;
		offsetsArray.push_back( connectivityArray.size() );
		}

		int
		getWorldDimension() const
		{
		return worldDimension;
		}
		// set cell types
		typesArray.resize( NumberOfCells, (std::uint8_t) cellShape );

		VTK::EntityShape
		getCellShape() const
		{
		return cellShape;
		}
		// set the arrays to the base class
		this->pointsArray = std::move(pointsArray);
		this->connectivityArray = std::move(connectivityArray);
		this->offsetsArray = std::move(offsetsArray);
		this->typesArray = std::move(typesArray);

		String
		getRealType() const
		{
		// not stored in the Netgen file
		return "float";
		}

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

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

		protected:
		String fileName;
		bool meshDetected = false;
		int meshDimension, worldDimension;
		VTK::EntityShape cellShape = VTK::EntityShape::Vertex;

		void reset()
		{
		meshDetected = false;
		meshDimension = worldDimension = 0;
		cellShape = VTK::EntityShape::Vertex;
		meshDetected = true;
		}
		};

src/TNL/Meshes/Readers/TNLReader.h

+9 −10

Original line number	Diff line number	Diff line
		@@ -59,8 +59,8 @@ public:
		std::cerr << "Unable to parse the mesh config type " << parsedMeshType[ 1 ] << "." << std::endl;
		return false;
		}
		if( parsedMeshConfig.size() != 7 ) {
		std::cerr << "The parsed mesh config type has wrong size (expected 7 elements):" << std::endl;
		if( parsedMeshConfig.size() != 6 ) {
		std::cerr << "The parsed mesh config type has wrong size (expected 6 elements):" << std::endl;
		std::cerr << "[ ";
		for( std::size_t i = 0; i < parsedMeshConfig.size() - 1; i++ )
		std::cerr << parsedMeshConfig[ i ] << ", ";
		@@ -75,15 +75,15 @@ public:
		globalIndexType = parsedMeshConfig[ 4 ];
		localIndexType = parsedMeshConfig[ 5 ];

		if( topology == "MeshEdgeTopology" )
		if( topology == "TNL::Meshes::Topologies::Edge" )
		cellShape = VTK::EntityShape::Line;
		else if( topology == "MeshTriangleTopology" )
		else if( topology == "TNL::Meshes::Topologies::Triangle" )
		cellShape = VTK::EntityShape::Triangle;
		else if( topology == "MeshQuadrilateralTopology" )
		else if( topology == "TNL::Meshes::Topologies::Quadrilateral" )
		cellShape = VTK::EntityShape::Quad;
		else if( topology == "MeshTetrahedronTopology" )
		else if( topology == "TNL::Meshes::Topologies::Tetrahedron" )
		cellShape = VTK::EntityShape::Tetra;
		else if( topology == "MeshHexahedronTopology" )
		else if( topology == "TNL::Meshes::Topologies::Hexahedron" )
		cellShape = VTK::EntityShape::Hexahedron;
		else {
		std::cerr << "Detected topology '" << topology << "' is not supported." << std::endl;
		@@ -99,11 +99,10 @@ public:
		}

		template< typename MeshType >
		bool
		readMesh( MeshType& mesh )
		void
		loadMesh( MeshType& mesh )
		{
		mesh.load( fileName );
		return true;
		}

		String

src/TNL/Meshes/Readers/VTKReader.h

+65 −167

File changed.

Preview size limit exceeded, changes collapsed.