Commit 9076cf0f authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

fixed conversion of non-regular c3t3 to TNL mesh

parent 20f9200d
Loading
Loading
Loading
Loading
Loading
+28 −12
Original line number Diff line number Diff line
@@ -49,24 +49,40 @@ cgal_mesh_to_tnl( const C3T3& c3t3 )
	using MeshType = TetrahedralMesh< typename C3T3::Triangulation::Geom_traits::FT >;
	using PointType = typename MeshType::PointType;

	// NOTE: this function cannot work with the underlying triangulation alone:
	// - c3t3's triangulation contains even points outside the domain (on the
	//   bounding ball).
	// - Moreover, the triangulation may be non-regular after the remeshing
	//   algorithm, i.e. the empty-sphere property may be violated (there may
	//   additional points in the circumsphere of a cell). Especially on highly
	//   non-convex domains, the point may be actually inside the domain, so it
	//   is not counted in c3t3.number_of_far_points() and its in_dimension()
	//   may be > -1.
	// - 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<typename C3T3::Triangulation::Vertex_handle, std::size_t> V;
	std::size_t idx = 0;
	for( auto cit = c3t3.cells_in_complex_begin();
		cit != c3t3.cells_in_complex_end();
		cit++ )
	{
		for (int i = 0; i < 4; i++)
			if (V.count(cit->vertex(i)) == 0)
				V[cit->vertex(i)] = idx++;
	}

	// initialize TNL mesh builder
	TNL::Meshes::MeshBuilder< MeshType > builder;
	// NOTE: c3t3's triangulation contains even points outside the domain (on
	// the bounding ball), which is why this function must work with c3t3 even
	// though working with just the triangulation would be easier
	builder.setPointsCount(c3t3.triangulation().number_of_vertices() - c3t3.number_of_far_points());
	builder.setPointsCount(idx);
	builder.setCellsCount(c3t3.number_of_cells());

	// index mapping for vertices in the complex
	std::unordered_map<typename C3T3::Triangulation::Vertex_handle, std::size_t> V;

	// set points
	std::size_t idx = 0;
	for( auto vit : c3t3.triangulation().finite_vertex_handles() )
	{
		if(vit->in_dimension() <= -1) continue;
		V[vit] = idx;
		builder.setPoint(idx++, {vit->point()[0], vit->point()[1], vit->point()[2]});
		if (V.count(vit) == 0)
			continue;
		builder.setPoint(V[vit], {vit->point()[0], vit->point()[1], vit->point()[2]});
	}

	// set cell seeds