updated polyhedral generator to use the VTUWriter from TNL (a7998566) · Commits · Jakub Klinkovský / meshing

cinolib-mesh-to-tnl.h

0 → 100644

+134 −0

Original line number	Diff line number	Diff line
		#pragma once

		#include <TNL/Meshes/MeshBuilder.h>
		#include <TNL/Meshes/Mesh.h>
		#include <TNL/Meshes/DefaultConfig.h>
		#include <TNL/Meshes/Topologies/Polyhedron.h>

		#include <cinolib/meshes/polyhedralmesh.h>

		template< typename Cell,
		int SpaceDimension = Cell::dimension,
		typename Real = double,
		typename GlobalIndex = int,
		typename LocalIndex = GlobalIndex >
		struct MinimalConfig
		: public TNL::Meshes::DefaultConfig< Cell, SpaceDimension, Real, GlobalIndex, LocalIndex >
		{
		static constexpr bool subentityStorage( int entityDimension, int subentityDimension )
		{
		if( subentityDimension == 0 && entityDimension >= Cell::dimension - 1 )
		return true;
		if( subentityDimension == Cell::dimension - 1 && entityDimension == Cell::dimension )
		return true;
		return false;
		}

		static constexpr bool superentityStorage( int entityDimension, int superentityDimension )
		{
		return false;
		}

		static constexpr bool entityTagsStorage( int entityDimension )
		{
		return false;
		}

		static constexpr bool dualGraphStorage()
		{
		return false;
		}
		};

		template< typename Real >
		using PolyhedralMesh = TNL::Meshes::Mesh< MinimalConfig<
		TNL::Meshes::Topologies::Polyhedron,
		3,
		Real,
		std::size_t,
		short int > >;

		template< typename CinoPolyMesh >
		PolyhedralMesh< double >
		cinolib_mesh_to_tnl( const CinoPolyMesh& m_poly )
		{
		using MeshType = PolyhedralMesh< double >;
		using PointType = typename MeshType::PointType;
		using NeighborCountsArray = typename TNL::Meshes::MeshBuilder< MeshType >::NeighborCountsArray;

		// NOTE:
		// - Cinolib's polyhedral 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.my_get_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();

		// initialize TNL mesh builder
		TNL::Meshes::MeshBuilder< MeshType > builder;
		builder.setEntitiesCount(nv, nc, nf);

		// set points
		idx = 0;
		for( const cinolib::vec3d & v : m_poly.my_get_verts() )
		{
		if (V.count(idx) == 0)
		{
		idx++;
		continue;
		}
		builder.setPoint(V[idx], {v.x(), v.y(), v.z()});
		idx++;
		}

		// set face corners counts
		NeighborCountsArray faceCorners( nf );
		idx = 0;
		for( const std::vector<uint> & f : m_poly.my_get_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() )
		{
		auto seed = builder.getFaceSeed(idx++);
		for( std::size_t i = 0; i < f.size(); i++ )
		seed.setCornerId(i, V[f[i]]);
		}

		// set cell corners counts
		NeighborCountsArray cellCorners( nc );
		idx = 0;
		for( const std::vector<uint> & c : m_poly.my_get_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() )
		{
		auto seed = builder.getCellSeed(idx++);
		for( std::size_t i = 0; i < c.size(); i++ )
		seed.setCornerId(i, c[i]);
		}

		// build the mesh
		MeshType mesh;
		if (!builder.build(mesh))
		throw std::runtime_error("TNL's MeshBuilder failed");
		return mesh;
		}

polyhedral-generator.cpp

+21 −69

Original line number	Diff line number	Diff line
		@@ -3,8 +3,12 @@
		#include <sstream>

		#include <cinolib/meshes/tetmesh.h>
		#include <cinolib/meshes/polyhedralmesh.h>
		#include <cinolib/dual_mesh.h>

		#include <TNL/Meshes/Writers/VTUWriter.h>
		#include "cinolib-mesh-to-tnl.h"

		#include "lyra/lyra.hpp"

		template <typename T>
		@@ -48,79 +52,23 @@ struct Parameters
		}
		};

		void write_FPMA(const char * filename,
		const std::vector<cinolib::vec3d> & verts,
		const std::vector<std::vector<uint>> & faces,
		const std::vector<std::vector<uint>> & polys,
		const std::vector<std::vector<bool>> & polys_winding)
		{
		setlocale(LC_NUMERIC, "en_US.UTF-8"); // makes sure "." is the decimal separator

		FILE *fp = fopen(filename, "w");

		if(!fp)
		{
		std::cerr << "ERROR : " << __FILE__ << ", line " << __LINE__ << " : save_HEDRA() : couldn't write output file " << filename << std::endl;
		exit(-1);
		}

		uint nv = verts.size();
		uint nf = faces.size();
		uint np = polys.size();

		// write vertices
		fprintf(fp, "%d\n", nv);

		// http://stackoverflow.com/questions/16839658/printf-width-specifier-to-maintain-precision-of-floating-point-value
		for(const cinolib::vec3d & v : verts) fprintf(fp, "%.17g %.17g %.17g\n", v.x(), v.y(), v.z());

		fprintf(fp, "\n");

		// write faces
		fprintf(fp, "%d\n", nf);

		for(const std::vector<uint> & f : faces)
		{
		fprintf(fp, "%d ", static_cast<int>(f.size()));
		for(uint vid : f) fprintf(fp, "%d ", vid);
		fprintf(fp, "\n");
		}

		fprintf(fp, "\n");

		// write cells
		fprintf(fp, "%d\n", np);

		for(uint pid=0; pid<np; ++pid)
		{
		fprintf(fp, "%d ", static_cast<int>(polys.at(pid).size()));
		for(uint off=0; off<polys.at(pid).size(); ++off)
		fprintf(fp, "%d ", polys.at(pid).at(off));
		fprintf(fp, "\n");
		}

		fprintf(fp, "\n");

		fclose(fp);
		}

		class MyPolyhedralmesh : public cinolib::Polyhedralmesh<>
		{
		public:
		// hack to expose attributes that we need for FPMA export
		auto my_get_verts()
		// hack to expose attributes that we need for custom output
		auto my_get_verts() const
		{
		return this->verts;
		}
		auto my_get_faces()
		auto my_get_faces() const
		{
		return this->faces;
		}
		auto my_get_polys()
		auto my_get_polys() const
		{
		return this->polys;
		}
		auto my_get_polys_face_winding()
		auto my_get_polys_face_winding() const
		{
		return this->polys_face_winding;
		}
		@@ -147,9 +95,13 @@ void generateDual(
		MyPolyhedralmesh m_poly;
		cinolib::dual_mesh(m_tet, m_poly, args.detect_features);

		// write the polyhedral mesh
		//m_poly.save("output.hedra"); // only .hedra is supported for output
		write_FPMA(output_file.string().c_str(), m_poly.my_get_verts(), m_poly.my_get_faces(), m_poly.my_get_polys(), m_poly.my_get_polys_face_winding());
		// Output to VTU using the writer from TNL
		using MeshType = PolyhedralMesh<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<3>(tnlMesh);
		}

		int main(int argc, char * argv[])
		@@ -171,7 +123,7 @@ int main(int argc, char * argv[])
		("Output path. If it is an existing directory, the stem of "
		"input_file (i.e., base file name without suffix) is appended to "
		"output_path. Otherwise, output_path denotes the output file and "
		"it must end with the \".fpma\" suffix." + fmt_default(output_path));
		"it must end with the \".vtu\" suffix." + fmt_default(output_path));

		auto result = cli.parse({ argc, argv });

		@@ -196,10 +148,10 @@ int main(int argc, char * argv[])

		// check if output_path is an existing directory
		if (std::filesystem::is_directory(output_path))
		output_path /= (input_file.stem().string() + ".fpma");
		// ensure that output_path ends with ".fpma"
		else if (output_path.extension() != ".fpma") {
		std::cerr << "Error: output_path \"" + output_path.string() + "\" does not end with \".fpma\"." << std::endl;
		output_path /= (input_file.stem().string() + ".vtu");
		// ensure that output_path ends with ".vtu"
		else if (output_path.extension() != ".vtu") {
		std::cerr << "Error: output_path \"" + output_path.string() + "\" does not end with \".vtu\"." << std::endl;
		return EXIT_FAILURE;
		}