Commit 656e8e1f authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

generalized polyhedral-generator to handle also 2D polygonal meshes

parent 1a19c42c

cinolib-mesh-to-tnl.h

+97 −19
Original line number Diff line number Diff line
@@ -3,9 +3,11 @@ 
#include <TNL/Meshes/MeshBuilder.h>

#include <TNL/Meshes/Mesh.h>

#include <TNL/Meshes/DefaultConfig.h>

#include <TNL/Meshes/Topologies/Polygon.h>

#include <TNL/Meshes/Topologies/Polyhedron.h>



#include <cinolib/meshes/polyhedralmesh.h>

#include <cinolib/meshes/abstract_polyhedralmesh.h>

#include <cinolib/meshes/abstract_polygonmesh.h>



template< typename Cell,

          int SpaceDimension = Cell::dimension,
@@ -17,9 +19,10 @@ struct MinimalConfig 
{

	static constexpr bool subentityStorage( int entityDimension, int subentityDimension )

	{

		if( subentityDimension == 0 && entityDimension >= Cell::dimension - 1 )

		constexpr int D = Cell::dimension;

		if( subentityDimension == 0 && entityDimension >= D - 1 )

			return true;

		if( subentityDimension == Cell::dimension - 1 && entityDimension == Cell::dimension )

		if( D == 3 && subentityDimension == D - 1 && entityDimension == D )

			return true;

		return false;

	}
@@ -40,6 +43,14 @@ struct MinimalConfig 
	}

};



template< typename Real >

using PolygonalMesh = TNL::Meshes::Mesh< MinimalConfig<

						TNL::Meshes::Topologies::Polygon,

						2,

						Real,

						std::size_t,

						short int > >;



template< typename Real >

using PolyhedralMesh = TNL::Meshes::Mesh< MinimalConfig<

						TNL::Meshes::Topologies::Polyhedron,
@@ -48,9 +59,76 @@ using PolyhedralMesh = TNL::Meshes::Mesh< MinimalConfig< 
						std::size_t,

						short int > >;



template< typename CinoPolyMesh >

template< typename _M, typename _V, typename _E, typename _P >

PolygonalMesh< double >

cinolib_mesh_to_tnl( const cinolib::AbstractPolygonMesh<_M, _V, _E, _P>& m_poly )

{

	using MeshType = PolygonalMesh< double >;

	using PointType = typename MeshType::PointType;

	using NeighborCountsArray = typename TNL::Meshes::MeshBuilder< MeshType >::NeighborCountsArray;



	// NOTE:

	// - Cinolib's polygonal mesh may contain some points which are not used

	//   for cells (e.g. when the domain features are not preserved).

	// - Consequently, the mesh builder may fail with the error: Some points

	//   were not used for cells.

	// Hence, we need to first find the points that are actually used for cells:

	std::unordered_map<std::size_t, std::size_t> V;

	std::size_t idx = 0;

	for( const std::vector<uint> & f : m_poly.vector_polys() )

	{

		for( std::size_t i = 0; i < f.size(); i++)

			if (V.count(f[i]) == 0)

				V[f[i]] = idx++;

	}



	const std::size_t nv = idx;  // m_poly.vector_verts().size();

	const std::size_t nc = m_poly.vector_polys().size();



	// initialize TNL mesh builder

	TNL::Meshes::MeshBuilder< MeshType > builder;

	builder.setEntitiesCount(nv, nc);



	// set points

	idx = 0;

	for( const cinolib::vec3d & v : m_poly.vector_verts() )

	{

		if (V.count(idx) == 0)

		{

			idx++;

			continue;

		}

		builder.setPoint(V[idx], {v.x(), v.y(), v.z()});

		idx++;

	}



	// set cell corners counts

	NeighborCountsArray cellCorners( nc );

	idx = 0;

	for( const std::vector<uint> & c : m_poly.vector_polys() )

		cellCorners[idx++] = c.size();

	builder.setCellCornersCounts(cellCorners);

	cellCorners.reset();



	// set cell seeds

	idx = 0;

	for( const std::vector<uint> & c : m_poly.vector_polys() )

	{

		auto seed = builder.getCellSeed(idx++);

		for( std::size_t i = 0; i < c.size(); i++ )

			seed.setCornerId(i, V[c[i]]);

	}



	// build the mesh

	MeshType mesh;

	if (!builder.build(mesh))

		throw std::runtime_error("TNL's MeshBuilder failed");

	return mesh;

}



template< typename _M, typename _V, typename _E, typename _F, typename _P >

PolyhedralMesh< double >

cinolib_mesh_to_tnl( const CinoPolyMesh& m_poly )

cinolib_mesh_to_tnl( const cinolib::AbstractPolyhedralMesh<_M, _V, _E, _F, _P>& m_poly )

{

	using MeshType = PolyhedralMesh< double >;

	using PointType = typename MeshType::PointType;
@@ -64,16 +142,16 @@ cinolib_mesh_to_tnl( const CinoPolyMesh& m_poly ) 
	// Hence, we need to first find the points that are actually used for cells:

	std::unordered_map<std::size_t, std::size_t> V;

	std::size_t idx = 0;

	for( const std::vector<uint> & f : m_poly.my_get_faces() )

	for( const std::vector<uint> & f : m_poly.vector_faces() )

	{

		for( std::size_t i = 0; i < f.size(); i++)

			if (V.count(f[i]) == 0)

				V[f[i]] = idx++;

	}



	const std::size_t nv = idx;  // m_poly.my_get_verts().size();

	const std::size_t nf = m_poly.my_get_faces().size();

	const std::size_t nc = m_poly.my_get_polys().size();

	const std::size_t nv = idx;  // m_poly.vector_verts().size();

	const std::size_t nf = m_poly.vector_faces().size();

	const std::size_t nc = m_poly.vector_polys().size();



	// initialize TNL mesh builder

	TNL::Meshes::MeshBuilder< MeshType > builder;
@@ -81,7 +159,7 @@ cinolib_mesh_to_tnl( const CinoPolyMesh& m_poly ) 


	// set points

	idx = 0;

	for( const cinolib::vec3d & v : m_poly.my_get_verts() )

	for( const cinolib::vec3d & v : m_poly.vector_verts() )

	{

		if (V.count(idx) == 0)

		{
@@ -95,14 +173,14 @@ cinolib_mesh_to_tnl( const CinoPolyMesh& m_poly ) 
	// set face corners counts

	NeighborCountsArray faceCorners( nf );

	idx = 0;

	for( const std::vector<uint> & f : m_poly.my_get_faces() )

	for( const std::vector<uint> & f : m_poly.vector_faces() )

		faceCorners[idx++] = f.size();

	builder.setFaceCornersCounts(faceCorners);

	faceCorners.reset();



	// set face seeds

	idx = 0;

	for( const std::vector<uint> & f : m_poly.my_get_faces() )

	for( const std::vector<uint> & f : m_poly.vector_faces() )

	{

		auto seed = builder.getFaceSeed(idx++);

		for( std::size_t i = 0; i < f.size(); i++ )
@@ -112,14 +190,14 @@ cinolib_mesh_to_tnl( const CinoPolyMesh& m_poly ) 
	// set cell corners counts

	NeighborCountsArray cellCorners( nc );

	idx = 0;

	for( const std::vector<uint> & c : m_poly.my_get_polys() )

	for( const std::vector<uint> & c : m_poly.vector_polys() )

		cellCorners[idx++] = c.size();

	builder.setCellCornersCounts(cellCorners);

	cellCorners.reset();



	// set cell seeds

	idx = 0;

	for( const std::vector<uint> & c : m_poly.my_get_polys() )

	for( const std::vector<uint> & c : m_poly.vector_polys() )

	{

		auto seed = builder.getCellSeed(idx++);

		for( std::size_t i = 0; i < c.size(); i++ )

polyhedral-generator.cpp

+52 −43
Original line number Diff line number Diff line
@@ -3,6 +3,7 @@ 
#include <sstream>



#include <cinolib/meshes/tetmesh.h>

#include <cinolib/meshes/polygonmesh.h>

#include <cinolib/meshes/polyhedralmesh.h>

#include <cinolib/dual_mesh.h>
@@ -52,37 +53,17 @@ struct Parameters 
	}

};



class MyPolyhedralmesh : public cinolib::Polyhedralmesh<>

{

public:

	// hack to expose attributes that we need for custom output

	auto my_get_verts() const

	{

		return this->verts;

	}

	auto my_get_faces() const

	{

		return this->faces;

	}

	auto my_get_polys() const

	{

		return this->polys;

	}

	auto my_get_polys_face_winding() const

	{

		return this->polys_face_winding;

	}

};



void generateDual(

		const std::filesystem::path & input_file,

		const std::filesystem::path & output_file,

		const Parameters & args)

{

	// load the input mesh

	// try to load the input as tetrahedral mesh

	cinolib::Tetmesh<> m_tet;

	m_tet.load(input_file.string().c_str());



	if( m_tet.num_polys() > 0 )

	{

		// detect feature edges

		if (args.detect_features) {

			// convert degrees to radians
@@ -91,8 +72,7 @@ void generateDual( 
		}



		// make polyhedral (voronoi) mesh

	//cinolib::Polyhedralmesh<> m_poly;

	MyPolyhedralmesh m_poly;

		cinolib::Polyhedralmesh<> m_poly;

		cinolib::dual_mesh(m_tet, m_poly, args.detect_features);



		// Output to VTU using the writer from TNL
@@ -103,6 +83,33 @@ void generateDual( 
		MeshWriter writer(out);

		writer.template writeEntities<3>(tnlMesh);

	}

	else

	{

		// try to load the input as polygonal mesh

		std::cout << "Input mesh does not contain any 3D cell, loading as 2D polygonal mesh." << std::endl;

		cinolib::Polygonmesh<> m_2d;

		m_2d.load(input_file.string().c_str());



		// detect feature edges

		if (args.detect_features) {

			// convert degrees to radians

			const double angle_threshold = args.features_angle_bound * 3.1415926535897932384626433 / 180.;

			m_tet.edge_mark_sharp_creases(angle_threshold);

		}



		// make polygonal (voronoi) mesh

		cinolib::Polygonmesh<> m_poly;

		cinolib::dual_mesh(m_2d, m_poly, args.detect_features);



		// Output to VTU using the writer from TNL

		using MeshType = PolygonalMesh<double>;

		const MeshType tnlMesh = cinolib_mesh_to_tnl(m_poly);

		std::ofstream out(output_file);

		using MeshWriter = TNL::Meshes::Writers::VTUWriter<MeshType>;

		MeshWriter writer(out);

		writer.template writeEntities<2>(tnlMesh);

	}

}



int main(int argc, char * argv[])

{
@@ -112,13 +119,15 @@ int main(int argc, char * argv[]) 
	bool show_help = false;



	auto cli = lyra::cli();

	cli |= lyra::help(show_help).description("polyhedral mesh generator based on the cinolib library");

	cli |= lyra::help(show_help).description("polygonal and polyhedral mesh generator based on the cinolib library");



	args.set_cli_options(cli);



	// positional arguments (should be generally defined after optional arguments)

	cli |= lyra::arg(input_file, "input_file")

			("Input file of a tetrahedral mesh (supported formats: VTK, VTU, TET, MESH).").required();

			("Input file of 2D triangular (or general polygonal) mesh, or a "

			 "3D tetrahedral mesh. Supported formats: for 2D: OFF, OBJ, STL; "

			 "for 3D: VTK, VTU, TET, MESH.").required();

	cli |= lyra::arg(output_path, "output_path")

			("Output path. If it is an existing directory, the stem of "

			 "input_file (i.e., base file name without suffix) is appended to "