Merge remote-tracking branch 'origin/v0.8.0' into ifcviewer-wgpu

src/ifcgeom/kernels/cgal/CgalConversionResult.cpp

@@ -11,11 +11,10 @@
 
 using IfcGeom::OpaqueNumber;
 using IfcGeom::OpaqueCoordinate;
-using IfcGeom::NumberNativeDouble;
 using IfcGeom::ConversionResultShape;
 
 #ifdef IFOPSH_SIMPLE_KERNEL
-#define NumberType NumberNativeDouble
+#define NumberType OpaqueNumber
 #else
 using ifcopenshell::geometry::NumberEpeck;
 #define NumberType NumberEpeck
@@ -30,6 +29,174 @@ typedef Polyhedron::Facet_const_handle Facet_const_handle;
 typedef Polyhedron::Halfedge_around_facet_const_circulator Halfedge_around_facet_circulator;
 
 namespace {
+	cgal_placement_t make_transform(const ifcopenshell::geometry::taxonomy::matrix4& place) {
+		const auto& m = place.ccomponents();
+		return cgal_placement_t(
+			m(0, 0), m(0, 1), m(0, 2), m(0, 3),
+			m(1, 0), m(1, 1), m(1, 2), m(1, 3),
+			m(2, 0), m(2, 1), m(2, 2), m(2, 3));
+	}
+
+	OpaqueCoordinate<3> opaque_point(const cgal_point_t& p) {
+		return OpaqueCoordinate<3>(
+			NumberType(p.cartesian(0)),
+			NumberType(p.cartesian(1)),
+			NumberType(p.cartesian(2))
+		);
+	}
+
+	typename Kernel_::FT max_abs3(const typename Kernel_::FT& a, const typename Kernel_::FT& b, const typename Kernel_::FT& c) {
+		std::array<typename Kernel_::FT, 3> abc{ a, b, c };
+		auto minel = std::min_element(abc.begin(), abc.end());
+		auto maxel = std::max_element(abc.begin(), abc.end());
+		return ((-*minel) > *maxel) ? (-*minel) : *maxel;
+	}
+
+	OpaqueCoordinate<3> opaque_axis(const cgal_vector_t& v) {
+		auto maxval = max_abs3(v.x(), v.y(), v.z());
+		if (maxval == 0) {
+			throw std::runtime_error("Invalid shape type");
+		}
+		return OpaqueCoordinate<3>(
+			NumberType(v.x() / maxval),
+			NumberType(v.y() / maxval),
+			NumberType(v.z() / maxval)
+		);
+	}
+
+	OpaqueCoordinate<4> opaque_plane(const cgal_plane_t& p) {
+		auto maxval = max_abs3(p.a(), p.b(), p.c());
+		if (maxval == 0) {
+			throw std::runtime_error("Invalid shape type");
+		}
+		return OpaqueCoordinate<4>(
+			NumberType(p.a() / maxval),
+			NumberType(p.b() / maxval),
+			NumberType(p.c() / maxval),
+			NumberType(p.d() / maxval)
+		);
+	}
+
+	cgal_plane_t plane_from_opaque(const OpaqueCoordinate<4>& p) {
+#ifdef IFOPSH_SIMPLE_KERNEL
+		return cgal_plane_t(
+			p.get(0).to_double(),
+			p.get(1).to_double(),
+			p.get(2).to_double(),
+			p.get(3).to_double()
+		);
+#else
+		return cgal_plane_t(
+			p.get(0).value_as<CGAL::Epeck::FT>(),
+			p.get(1).value_as<CGAL::Epeck::FT>(),
+			p.get(2).value_as<CGAL::Epeck::FT>(),
+			p.get(3).value_as<CGAL::Epeck::FT>()
+		);
+#endif
+	}
+
+	void insert_normalized_plane_map(plane_map<Kernel_>& mp, const OpaqueCoordinate<4>& from, const OpaqueCoordinate<4>& to) {
+		mp.insert({
+			normalized_plane_for_map<Kernel_>(plane_from_opaque(from)),
+			normalized_plane_for_map<Kernel_>(plane_from_opaque(to))
+		});
+	}
+
+	void apply_normalized_plane_map(const plane_map<Kernel_>& mp, std::list<cgal_plane_t>& planes) {
+		for (auto& plane : planes) {
+			auto it = mp.find(normalized_plane_for_map<Kernel_>(plane));
+			if (it != mp.end()) {
+				plane = it->second;
+			}
+		}
+	}
+
+	cgal_vector_t wire_normal(const cgal_wire_t& wire) {
+		typename Kernel_::FT a(0), b(0), c(0);
+		if (wire.size() < 3) {
+			return cgal_vector_t(a, b, c);
+		}
+		for (std::size_t i = 0; i < wire.size(); ++i) {
+			const auto& curr = wire[i];
+			const auto& next = wire[(i + 1) % wire.size()];
+			a += (curr.y() - next.y()) * (curr.z() + next.z());
+			b += (curr.z() - next.z()) * (curr.x() + next.x());
+			c += (curr.x() - next.x()) * (curr.y() + next.y());
+		}
+		return cgal_vector_t(a, b, c);
+	}
+
+	cgal_point_t wire_centroid(const cgal_wire_t& wire) {
+		if (wire.empty()) {
+			throw std::runtime_error("Invalid shape type");
+		}
+		std::array<Kernel_::FT, 3> p{ Kernel_::FT(0), Kernel_::FT(0), Kernel_::FT(0) };
+		for (const auto& point : wire) {
+			for (int i = 0; i < 3; ++i) {
+				p[i] += point.cartesian(i);
+			}
+		}
+		Kernel_::FT n(wire.size());
+		return cgal_point_t(p[0] / n, p[1] / n, p[2] / n);
+	}
+
+	Kernel_::FT wire_length(const cgal_wire_t& wire) {
+		Kernel_::FT len(0);
+		if (wire.size() < 2) {
+			return len;
+		}
+		for (std::size_t i = 1; i < wire.size(); ++i) {
+			len += CGAL::approximate_sqrt(CGAL::Segment_3<Kernel_>(wire[i - 1], wire[i]).squared_length());
+		}
+		if (wire.size() > 2) {
+			len += CGAL::approximate_sqrt(CGAL::Segment_3<Kernel_>(wire.back(), wire.front()).squared_length());
+		}
+		return len;
+	}
+
+	Kernel_::FT wire_area(const cgal_wire_t& wire) {
+		Kernel_::FT area(0);
+		if (wire.size() < 3) {
+			return area;
+		}
+		const auto& origin = wire.front();
+		for (std::size_t i = 1; i + 1 < wire.size(); ++i) {
+			auto v1 = wire[i] - origin;
+			auto v2 = wire[i + 1] - origin;
+			area += CGAL::approximate_sqrt(CGAL::cross_product(v1, v2).squared_length()) / Kernel_::FT(2);
+		}
+		return area;
+	}
+
+	cgal_wire_t moved_wire(const cgal_wire_t& wire, const cgal_placement_t& trsf) {
+		cgal_wire_t result;
+		result.reserve(wire.size());
+		for (const auto& point : wire) {
+			result.push_back(point.transform(trsf));
+		}
+		return result;
+	}
+
+	void write_off_point(std::stringstream& sstream, const cgal_point_t& point) {
+		sstream << "OFF\n1 0 0\n";
+		sstream << point.x() << " " << point.y() << " " << point.z() << "\n";
+	}
+
+	void write_off_wire(std::stringstream& sstream, const cgal_wire_t& wire) {
+		const bool face = wire.size() >= 3;
+		sstream << "OFF\n" << wire.size() << " " << (face ? 1 : 0) << " 0\n";
+		for (const auto& point : wire) {
+			sstream << point.x() << " " << point.y() << " " << point.z() << "\n";
+		}
+		if (face) {
+			sstream << wire.size();
+			for (std::size_t i = 0; i < wire.size(); ++i) {
+				sstream << " " << i;
+			}
+			sstream << "\n";
+		}
+	}
+
 	template <typename Facet>
 	CGAL::Direction_3<Kernel_> newell(Facet& face) {
 		typename Kernel_::FT a(0), b(0), c(0);
@@ -99,20 +266,21 @@ namespace {
 	}
 }
 
-ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t& shape, bool convex) {
+ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t& shape, bool convex, Logger& logger) {
 	shape_ = shape;
 	convex_tag_ = convex;
+	auto& poly = std::get<cgal_shape_t>(*shape_);
 
 	std::set<cgal_shape_t::Facet_handle> faces_to_remove;
 
-	for (const auto& face : CGAL::faces(*shape_)) {
+	for (const auto& face : CGAL::faces(poly)) {
 		auto V = newell(*face).to_vector();
 		CGAL::Plane_3<Kernel_> plane(CGAL::Point_3<Kernel_>(), V);
 		auto b1 = plane.base1();
 		auto b2 = plane.base2();
 
 		if (V.squared_length() == 0) {
-			logger::warning("Removed face due to self-intersections");
+			logger.Warning("GEO", 62, "Removed face due to self-intersections");
 			faces_to_remove.insert(face);
 			continue;
 		}
@@ -127,46 +295,54 @@ ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t& shape, bool con
 
 		std::vector<CGAL::Point_2<Kernel_>> ps;
 			
-		for (auto& he1 : CGAL::halfedges_around_face(face->halfedge(), *shape_)) {
+		for (auto& he1 : CGAL::halfedges_around_face(face->halfedge(), poly)) {
 			const auto& source = he1->vertex()->point();
 			ps.push_back(transform_point(source));
 		}
 
 		if (!CGAL::Polygon_2<Kernel_>(ps.begin(), ps.end()).is_simple()) {
-			logger::warning("Removed face due to self-intersections");
+			logger.Warning("GEO", 63, "Removed face due to self-intersections");
 			faces_to_remove.insert(face);
 		}
 	}
 
 	{
 		for (auto& face : faces_to_remove) {
-			CGAL::Euler::remove_face(face->halfedge(), *shape_);
+			CGAL::Euler::remove_face(face->halfedge(), poly);
 		}
 	}
+}
 
-	if (shape.size_of_facets() != 1) {
-		// the size_of_facets() == 1 check is for handling the specical case of
-		// storing a single point in a polyhedron as a degenerate triangle
-		// 
-		// @todo come up with a proper variant for storing lower dimensional entities
-		
-		// @todo we don't have access to settings here so we don't know whether we should triangulate
-		// remove_degenerate_faces() is also called in the triangulate() call below though...
-		// CGAL::Polygon_mesh_processing::triangulate_faces(*shape_);
-		// CGAL::Polygon_mesh_processing::remove_degenerate_faces(*shape_);
+ifcopenshell::geometry::CgalShape::CgalShape(const cgal_point_t& point, bool convex) {
+	shape_ = point;
+	convex_tag_ = convex;
+}
+
+ifcopenshell::geometry::CgalShape::CgalShape(const cgal_wire_t& wire, bool convex) {
+	shape_ = wire;
+	convex_tag_ = convex;
+}
+
+const cgal_shape_t& ifcopenshell::geometry::CgalShape::poly() const {
+#ifndef IFOPSH_SIMPLE_KERNEL
+	to_poly();
+#endif
+	if (!shape_ || !std::holds_alternative<cgal_shape_t>(*shape_)) {
+		throw std::runtime_error("Invalid shape type");
 	}
+	return std::get<cgal_shape_t>(*shape_);
 }
 
 #ifndef IFOPSH_SIMPLE_KERNEL
 void ifcopenshell::geometry::CgalShape::to_poly() const {
 	if (!shape_) {
-		shape_.emplace();
-
-		convert_to_polyhedron(*nef_, *shape_);
-		if (shape_->size_of_vertices() > 0) {
+		cgal_shape_t poly;
+		convert_to_polyhedron(*nef_, poly, std::numeric_limits<std::size_t>::max());
+		if (poly.size_of_vertices() > 0) {
 			// @todo why is this necessary? we have the mark of the volumes?
-			CGAL::Polygon_mesh_processing::orient_to_bound_a_volume(*shape_);
+			CGAL::Polygon_mesh_processing::orient_to_bound_a_volume(poly);
 		}
+		shape_ = poly;
 		
 		// nef_->convert_to_polyhedron(*shape_);
 	}
@@ -174,18 +350,24 @@ void ifcopenshell::geometry::CgalShape::to_poly() const {
 
 void ifcopenshell::geometry::CgalShape::to_nef() const {
 	if (!nef_) {
+		auto shp = poly();
 		if (!convex_tag_) {
-			if (CGAL::Polygon_mesh_processing::does_self_intersect(*shape_)) {
+			CGAL::Polygon_mesh_processing::triangulate_faces(shp);
+			if (CGAL::Polygon_mesh_processing::does_self_intersect(shp)) {
 				throw std::runtime_error("Self-intersections detected, unable to proceed");
 			}
 		}
-		nef_ = utils::create_nef_polyhedron(*shape_);
+		nef_ = utils::create_nef_polyhedron(shp);
 	}
 }
 #endif
 
-void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const {
-	const bool all_triangles = std::all_of(shape_->facets_begin(), shape_->facets_end(), [](auto f) { return f.is_triangle(); });
+void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id, Logger& logger) const {
+	if (is_point() || is_wire()) {
+		return;
+	}
+	const auto& base_shape = poly();
+	const bool all_triangles = std::all_of(base_shape.facets_begin(), base_shape.facets_end(), [](auto f) { return f.is_triangle(); });
 	const bool has_iden_transform = place.is_identity();
 
 	std::unique_ptr<cgal_shape_t> shape_copy_holder;
@@ -193,10 +375,10 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett
 
 	if (!all_triangles || !has_iden_transform) {
 		// A copy is made when triangulate_faces() is required or when vertex positions need be transformed
-		shape_copy_holder.reset(new cgal_shape_t(*this));
+		shape_copy_holder.reset(new cgal_shape_t(base_shape));
 		shape_to_use = shape_copy_holder.get();
 	} else {
-		shape_to_use = &*shape_;
+		shape_to_use = const_cast<cgal_shape_t*>(&base_shape);
 	}
 
 	const bool setting_use_original_edges = settings.get<ifcopenshell::geometry::settings::CgalEmitOriginalEdges>().get();
@@ -233,7 +415,7 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett
 
 	if (!all_triangles) {
 		if (!shape_to_use->is_valid()) {
-			logger::message(logger::LOG_ERROR, "Invalid Polyhedron_3 in object (before triangulation)");
+			logger.Message(Logger::LOG_ERROR, "GEO", 64, "Invalid Polyhedron_3 in object (before triangulation)");
 			return;
 		}
 
@@ -241,19 +423,19 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett
 		try {
 			success = CGAL::Polygon_mesh_processing::triangulate_faces(*shape_to_use);
 		} catch (...) {
-			logger::message(logger::LOG_ERROR, "Triangulation crashed");
+			logger.Message(Logger::LOG_ERROR, "GEO", 65, "Triangulation crashed");
 			return;
 		}
 
 		CGAL::Polygon_mesh_processing::remove_degenerate_faces(*shape_to_use);
 
 		if (!success) {
-			logger::message(logger::LOG_ERROR, "Triangulation failed");
+			logger.Message(Logger::LOG_ERROR, "GEO", 66, "Triangulation failed");
 			return;
 		}
 
 		if (!shape_to_use->is_valid()) {
-			logger::message(logger::LOG_ERROR, "Invalid Polyhedron_3 in object (after triangulation)");
+			logger.Message(Logger::LOG_ERROR, "GEO", 67, "Invalid Polyhedron_3 in object (after triangulation)");
 			return;
 		}
 	}
@@ -282,7 +464,7 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett
 	try {
 		CGAL::Polygon_mesh_processing::compute_face_normals(*shape_to_use, face_normals_map);
 	} catch (...) {
-		logger::message(logger::LOG_ERROR, "Face normal calculation failed");
+		logger.Message(Logger::LOG_ERROR, "GEO", 68, "Face normal calculation failed");
 		return;
 	}
 
@@ -403,25 +585,33 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett
 }
 
 void ifcopenshell::geometry::CgalShape::Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string& r) const {
-	cgal_shape_t s = *this;
-
-	if (!place.is_identity()) {
-		const auto& m = place.ccomponents();
-
-		// @todo check
-		const cgal_placement_t trsf(
-			m(0, 0), m(0, 1), m(0, 2), m(0, 3),
-			m(1, 0), m(1, 1), m(1, 2), m(1, 3),
-			m(2, 0), m(2, 1), m(2, 2), m(2, 3));
-
-		// Apply transformation
-		for (auto &vertex : s.vertex_handles()) {
-			vertex->point() = vertex->point().transform(trsf);
-		}
-	}
-
 	std::stringstream sstream;
-	sstream << s;
+	if (is_point()) {
+		auto p = point();
+		if (!place.is_identity()) {
+			p = p.transform(make_transform(place));
+		}
+		write_off_point(sstream, p);
+	} else if (is_wire()) {
+		auto w = wire();
+		if (!place.is_identity()) {
+			w = moved_wire(w, make_transform(place));
+		}
+		write_off_wire(sstream, w);
+	} else {
+		cgal_shape_t s = poly();
+
+		if (!place.is_identity()) {
+			const auto trsf = make_transform(place);
+
+			// Apply transformation
+			for (auto &vertex : s.vertex_handles()) {
+				vertex->point() = vertex->point().transform(trsf);
+			}
+		}
+
+		sstream << s;
+	}
 	r = sstream.str();
 }
 
@@ -432,12 +622,26 @@ double ifcopenshell::geometry::CgalShape::bounding_box(void *& b) const {
 		b = new CGAL::Bbox_3;
 	}
 	auto& bb = (*((CGAL::Bbox_3*)b));
-	bb += CGAL::Polygon_mesh_processing::bbox(static_cast<cgal_shape_t>(*this));
+	if (is_point()) {
+		bb += point().bbox();
+	} else if (is_wire()) {
+		for (const auto& point : wire()) {
+			bb += point.bbox();
+		}
+	} else {
+		bb += CGAL::Polygon_mesh_processing::bbox(poly());
+	}
 	return (bb.xmax() - bb.xmin()) * (bb.ymax() - bb.ymin()) * (bb.zmax() - bb.zmin());
 }
 
 int ifcopenshell::geometry::CgalShape::num_vertices() const {
-	return (int) static_cast<cgal_shape_t>(*this).size_of_vertices();
+	if (is_point()) {
+		return 1;
+	}
+	if (is_wire()) {
+		return (int) wire().size();
+	}
+	return (int) poly().size_of_vertices();
 }
 
 void ifcopenshell::geometry::CgalShape::set_box(void * b) {
@@ -448,10 +652,13 @@ void ifcopenshell::geometry::CgalShape::set_box(void * b) {
 }
 
 int ifcopenshell::geometry::CgalShape::surface_genus() const {
-	to_poly();
-	auto nv = shape_->size_of_vertices();
-	auto ne = shape_->size_of_halfedges() / 2;
-	auto nf = shape_->size_of_facets();
+	if (is_point() || is_wire()) {
+		return 0;
+	}
+	const auto& shp = poly();
+	auto nv = shp.size_of_vertices();
+	auto ne = shp.size_of_halfedges() / 2;
+	auto nf = shp.size_of_facets();
 
 	auto euler = nv - ne + nf;
 	auto genus = (2 - euler) / 2;
@@ -461,14 +668,22 @@ int ifcopenshell::geometry::CgalShape::surface_genus() const {
 
 bool ifcopenshell::geometry::CgalShape::is_manifold() const {
 	// @todo ?
-	to_poly();
-	return shape_->is_valid();
+	return (is_point() || is_wire()) ? true : poly().is_valid();
 }
 
 int ifcopenshell::geometry::CgalShape::num_edges() const
 {
-	to_poly();
-	return (int) shape_->size_of_halfedges() / 2;
+	if (is_point()) {
+		return 0;
+	}
+	if (is_wire()) {
+		const auto n = wire().size();
+		if (n < 2) {
+			return 0;
+		}
+		return (int)(n == 2 ? 1 : n);
+	}
+	return (int) poly().size_of_halfedges() / 2;
 }
 
 int ifcopenshell::geometry::CgalShape::num_faces() const
@@ -479,61 +694,84 @@ int ifcopenshell::geometry::CgalShape::num_faces() const
 	} else
 #endif
 	if (shape_) {
-		return (int) shape_->size_of_facets();
+		if (is_poly()) {
+			return (int) poly().size_of_facets();
+		}
+		if (is_wire() && wire().size() >= 3) {
+			return 1;
+		}
+		return 0;
 	} else {
 		return 0;
 	}
 }
 
-OpaqueNumber* ifcopenshell::geometry::CgalShape::CgalShape::length()
+OpaqueNumber ifcopenshell::geometry::CgalShape::CgalShape::length()
 {
-	to_poly();
 	Kernel_::FT len = 0;
-	for (auto it = shape_->edges_begin(); it != shape_->edges_end(); ++it) {
-		len += CGAL::approximate_sqrt(CGAL::Segment_3<Kernel_>(
-			it->vertex()->point(),
-			it->next()->vertex()->point()
-		).squared_length());
+	if (is_wire()) {
+		len = wire_length(wire());
+	} else if (!is_point()) {
+		const auto& shp = poly();
+		for (auto it = shp.edges_begin(); it != shp.edges_end(); ++it) {
+			len += CGAL::approximate_sqrt(CGAL::Segment_3<Kernel_>(
+				it->vertex()->point(),
+				it->opposite()->vertex()->point()
+			).squared_length());
+		}
 	}
-	return new NumberType(len);
+	return NumberType(len);
 }
 
-OpaqueNumber* ifcopenshell::geometry::CgalShape::area()
+OpaqueNumber ifcopenshell::geometry::CgalShape::area()
 {
-	to_poly();
-	auto s = *shape_;
+	if (is_wire()) {
+		return NumberType(wire_area(wire()));
+	}
+	if (is_point()) {
+		return NumberType(Kernel_::FT(0));
+	}
+	auto s = poly();
 	CGAL::Polygon_mesh_processing::triangulate_faces(s);
-	return new NumberType(CGAL::Polygon_mesh_processing::area(s));
+	return NumberType(CGAL::Polygon_mesh_processing::area(s));
 }
 
-OpaqueNumber* ifcopenshell::geometry::CgalShape::volume()
+OpaqueNumber ifcopenshell::geometry::CgalShape::volume()
 {
-	to_poly();
-	auto s = *shape_;
+	if (is_point() || is_wire()) {
+		return NumberType(Kernel_::FT(0));
+	}
+	auto s = poly();
 	CGAL::Polygon_mesh_processing::triangulate_faces(s);
-	return new NumberType(CGAL::Polygon_mesh_processing::volume(s));
+	return NumberType(CGAL::Polygon_mesh_processing::volume(s));
 }
 
 OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::position()
 {
-	to_poly();
-	if (shape_->size_of_facets() == 1) {
+	if (is_point()) {
+		return opaque_point(point());
+	}
+	if (is_wire()) {
+		return opaque_point(wire_centroid(wire()));
+	}
+	const auto& shp = poly();
+	if (shp.size_of_facets() == 1) {
 		// return centroid;
 		// CGAL::Vector_3<Kernel_> p;
-		std::array<Kernel_::FT, 3> p;
-		for (auto it = shape_->points_begin(); it != shape_->points_end(); ++it) {
+		std::array<Kernel_::FT, 3> p{ Kernel_::FT(0), Kernel_::FT(0), Kernel_::FT(0) };
+		for (auto it = shp.points_begin(); it != shp.points_end(); ++it) {
 			for (int i = 0; i < 3; ++i) {
 				p[i] += it->cartesian(i);
 			}			
 		}
-		Kernel_::FT N(std::distance(shape_->points_begin(), shape_->points_end()));
+		Kernel_::FT N(std::distance(shp.points_begin(), shp.points_end()));
 		for (int i = 0; i < 3; ++i) {
 			p[i] /= N;
 		}
 		return OpaqueCoordinate<3>(
-			new NumberType(p[0]),
-			new NumberType(p[1]),
-			new NumberType(p[2])
+			NumberType(p[0]),
+			NumberType(p[1]),
+			NumberType(p[2])
 		);
 	} else {
 		throw std::runtime_error("Invalid shape type");
@@ -542,19 +780,19 @@ OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::position()
 
 OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::axis()
 {
-	to_poly();
-	if (shape_->size_of_facets() == 1) {
-		auto pl = Plane_equation()(*shape_->facets_begin());
-		std::array<typename Kernel_::FT, 3> abc{ pl.a(), pl.b(), pl.c() };
-		auto minel = std::min_element(abc.begin(), abc.end());
-		auto maxel = std::max_element(abc.begin(), abc.end());
-		auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
-
-		return OpaqueCoordinate<3>(
-			new NumberType(pl.a() / maxval),
-			new NumberType(pl.b() / maxval),
-			new NumberType(pl.c() / maxval)
-		);
+	if (is_wire()) {
+		if (wire().size() == 2) {
+			return opaque_axis(wire()[1] - wire()[0]);
+		}
+		if (wire().size() >= 3) {
+			return opaque_axis(wire_normal(wire()));
+		}
+		throw std::runtime_error("Invalid shape type");
+	}
+	auto shp = poly();
+	if (shp.size_of_facets() == 1) {
+		auto pl = Plane_equation()(*shp.facets_begin());
+		return opaque_axis(cgal_vector_t(pl.a(), pl.b(), pl.c()));
 	} else {
 		throw std::runtime_error("Invalid shape type");
 	}
@@ -562,6 +800,14 @@ OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::axis()
 
 OpaqueCoordinate<4> ifcopenshell::geometry::CgalShape::plane_equation()
 {
+	if (is_wire() && wire().size() >= 3) {
+		auto normal = wire_normal(wire());
+		return opaque_plane(cgal_plane_t(wire().front(), CGAL::Direction_3<Kernel_>(normal)));
+	}
+	auto shp = poly();
+	if (shp.size_of_facets() == 1) {
+		return opaque_plane(Plane_equation()(*shp.facets_begin()));
+	}
 	throw std::runtime_error("Invalid shape type");
 }
 
@@ -612,75 +858,71 @@ ConversionResultShape * ifcopenshell::geometry::CgalShape::box()
 
 ConversionResultShape* ifcopenshell::geometry::CgalShape::wrap_in_compound()
 {
-	return new CgalShape(poly(), convex_tag_);
+	return clone();
 }
 
 std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShape::vertices()
 {
-	// @todo this is ridiculous
-	to_poly();
 	std::vector<ConversionResultShape*> result;
-	for (auto& p : shape_->points()) {
-		std::vector<cgal_point_t> ps = {
-			p, p, p
-		};
-
-		std::vector<std::vector<size_t>> ids(1);
-		ids.front().push_back(0);
-		ids.front().push_back(1);
-		ids.front().push_back(2);
-
-		cgal_shape_t poly;
-		CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(ps, ids, poly);
-
-		result.push_back(new CgalShape(poly));
+	if (is_point()) {
+		result.push_back(new CgalShape(point()));
+		return result;
+	}
+	if (is_wire()) {
+		for (const auto& p : wire()) {
+			result.push_back(new CgalShape(p));
+		}
+		return result;
+	}
+	for (const auto& p : poly().points()) {
+		result.push_back(new CgalShape(p));
 	}
 	return result;
 }
 
 std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShape::edges()
 {
-	// @todo this is ridiculous
-	to_poly();
 	std::vector<ConversionResultShape*> result;
-	for (auto& ed : shape_->edges()) {
-		std::vector<cgal_point_t> ps = {
-			ed.vertex()->point(),
-			ed.vertex()->point(),
-			ed.next()->vertex()->point()
-		};
-
-		std::vector<std::vector<size_t>> ids(1);
-		ids.front().push_back(0);
-		ids.front().push_back(1);
-		ids.front().push_back(2);
-
-		cgal_shape_t poly;
-		CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(ps, ids, poly);
-
-		result.push_back(new CgalShape(poly));
+	if (is_point()) {
+		return result;
+	}
+	if (is_wire()) {
+		const auto& w = wire();
+		for (std::size_t i = 1; i < w.size(); ++i) {
+			result.push_back(new CgalShape(cgal_wire_t{ w[i - 1], w[i] }));
+		}
+		if (w.size() > 2) {
+			result.push_back(new CgalShape(cgal_wire_t{ w.back(), w.front() }));
+		}
+		return result;
+	}
+	for (auto ed : poly().edges()) {
+		result.push_back(new CgalShape(cgal_wire_t{ ed.vertex()->point(), ed.opposite()->vertex()->point() }));
 	}
 	return result;
 }
 
 std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShape::facets()
 {
-	to_poly();
 	std::vector<ConversionResultShape*> result;
-	for (auto &face : faces(*shape_)) {
+	if (is_point()) {
+		return result;
+	}
+	if (is_wire()) {
+		if (wire().size() >= 3) {
+			result.push_back(new CgalShape(wire()));
+		}
+		return result;
+	}
+	for (auto face : faces(poly())) {
 		std::vector<cgal_point_t> ps;
-		std::vector<std::vector<size_t>> ids(1);
 
 		auto it = face->facet_begin();
 		do {
 			ps.push_back(it->vertex()->point());
-			ids.front().push_back(ids.front().size());
 		} while (++it != face->facet_begin());
 
-		cgal_shape_t poly;
-		CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(ps, ids, poly);
-
-		result.push_back(new CgalShape(poly));
+		result.push_back(new CgalShape(ps));
 	}
 	return result;
 }
@@ -756,31 +998,31 @@ std::pair<OpaqueCoordinate<3>, OpaqueCoordinate<3>> ifcopenshell::geometry::Cgal
 
 ConversionResultShape* ifcopenshell::geometry::CgalShape::moved(ifcopenshell::geometry::taxonomy::matrix4::ptr place) const
 {
-	cgal_shape_t s = *this;
+	if (place->is_identity()) {
+		return clone();
+	}
 
-	if (!place->is_identity()) {
-		const auto& m = place->ccomponents();
+	const auto trsf = make_transform(*place);
+	if (is_point()) {
+		return new CgalShape(point().transform(trsf), convex_tag_);
+	}
+	if (is_wire()) {
+		return new CgalShape(moved_wire(wire(), trsf), convex_tag_);
+	}
 
-		// @todo check
-		const cgal_placement_t trsf(
-			m(0, 0), m(0, 1), m(0, 2), m(0, 3),
-			m(1, 0), m(1, 1), m(1, 2), m(1, 3),
-			m(2, 0), m(2, 1), m(2, 2), m(2, 3));
-
-		// Apply transformation
-		for (auto &vertex : s.vertex_handles()) {
-			vertex->point() = vertex->point().transform(trsf);
-		}
+	cgal_shape_t s = poly();
+	for (auto &vertex : s.vertex_handles()) {
+		vertex->point() = vertex->point().transform(trsf);
 	}
 
 	return new CgalShape(s, convex_tag_);
 }
 
-void ifcopenshell::geometry::CgalShape::map(OpaqueCoordinate<4>&, OpaqueCoordinate<4>&) {
+std::size_t ifcopenshell::geometry::CgalShape::map(OpaqueCoordinate<4>&, OpaqueCoordinate<4>&) {
 	throw std::runtime_error("Not implemented");
 }
 
-void ifcopenshell::geometry::CgalShape::map(const std::vector<OpaqueCoordinate<4>>&, const std::vector<OpaqueCoordinate<4>>&) {
+std::size_t ifcopenshell::geometry::CgalShape::map(const std::vector<OpaqueCoordinate<4>>&, const std::vector<OpaqueCoordinate<4>>&) {
 	throw std::runtime_error("Not implemented");
 }
 
@@ -791,7 +1033,7 @@ bool ifcopenshell::geometry::CgalShape::surface_area_along_direction(double tol,
 
 #ifndef IFOPSH_SIMPLE_KERNEL
 
-void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const {
+void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id, Logger& logger) const {
 	throw std::runtime_error("Not implemented");
 }
 
@@ -825,17 +1067,17 @@ int ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::num_faces() const
 	throw std::runtime_error("Not implemented");
 }
 
-OpaqueNumber* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::CgalShapeHalfSpaceDecomposition::length()
+OpaqueNumber ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::CgalShapeHalfSpaceDecomposition::length()
 {
 	throw std::runtime_error("Not implemented");
 }
 
-OpaqueNumber* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::area()
+OpaqueNumber ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::area()
 {
 	throw std::runtime_error("Not implemented");
 }
 
-OpaqueNumber* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::volume()
+OpaqueNumber ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::volume()
 {
 	throw std::runtime_error("Not implemented");
 }
@@ -845,9 +1087,9 @@ OpaqueCoordinate<3> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::pos
 	if (planes_.size() == 1) {
 		auto xyz = CGAL::ORIGIN + planes_.front().d() * CGAL::Vector_3<Kernel_>(planes_.front().a(), planes_.front().b(), planes_.front().c());
 		return OpaqueCoordinate<3>(
-			new NumberType(xyz.cartesian(0)),
-			new NumberType(xyz.cartesian(1)),
-			new NumberType(xyz.cartesian(2))
+			NumberType(xyz.cartesian(0)),
+			NumberType(xyz.cartesian(1)),
+			NumberType(xyz.cartesian(2))
 		);
 	} else {
 		throw std::runtime_error("Invalid shape type");
@@ -862,9 +1104,9 @@ OpaqueCoordinate<3> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::axi
 		auto maxel = std::max_element(abc.begin(), abc.end());
 		auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
 		return OpaqueCoordinate<3>(
-			new NumberType(planes_.front().a() / maxval),
-			new NumberType(planes_.front().b() / maxval),
-			new NumberType(planes_.front().c() / maxval)
+			NumberType(planes_.front().a() / maxval),
+			NumberType(planes_.front().b() / maxval),
+			NumberType(planes_.front().c() / maxval)
 		);
 	} else {
 		throw std::runtime_error("Invalid shape type");
@@ -879,10 +1121,10 @@ OpaqueCoordinate<4> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::pla
 		auto maxel = std::max_element(abc.begin(), abc.end());
 		auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
 		return OpaqueCoordinate<4>(
-			new NumberType(planes_.front().a() / maxval),
-			new NumberType(planes_.front().b() / maxval),
-			new NumberType(planes_.front().c() / maxval),
-			new NumberType(planes_.front().d() / maxval)
+			NumberType(planes_.front().a() / maxval),
+			NumberType(planes_.front().b() / maxval),
+			NumberType(planes_.front().c() / maxval),
+			NumberType(planes_.front().d() / maxval)
 		);
 	} else {
 		throw std::runtime_error("Invalid shape type");
@@ -957,27 +1199,17 @@ ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::
 	throw std::runtime_error("Not implemented");
 }
 
-void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to) {
+std::size_t ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to) {
 	plane_map<Kernel_> mp;
-	mp.insert({
-		CGAL::Plane_3<Kernel_>(
-			static_cast<NumberEpeck*>(from.get(0))->value(),
-			static_cast<NumberEpeck*>(from.get(1))->value(),
-			static_cast<NumberEpeck*>(from.get(2))->value(),
-			static_cast<NumberEpeck*>(from.get(3))->value()
-		),
-		CGAL::Plane_3<Kernel_>(
-			static_cast<NumberEpeck*>(to.get(0))->value(),
-			static_cast<NumberEpeck*>(to.get(1))->value(),
-			static_cast<NumberEpeck*>(to.get(2))->value(),
-			static_cast<NumberEpeck*>(to.get(3))->value()
-		)
-	});
-	auto nw = shape_->map(mp);
+	insert_normalized_plane_map(mp, from, to);
+	std::size_t mutated = 0;
+	auto nw = shape_->map(mp, mutated);
 	shape_ = std::move(nw);
+	apply_normalized_plane_map(mp, planes_);
+	return mutated;
 }
 
-void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::map(const std::vector<OpaqueCoordinate<4>>& froms, const std::vector<OpaqueCoordinate<4>>& tos) {
+std::size_t ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::map(const std::vector<OpaqueCoordinate<4>>& froms, const std::vector<OpaqueCoordinate<4>>& tos) {
 	plane_map<Kernel_> mp;
 	if (froms.size() != tos.size()) {
 		throw std::runtime_error("Expected equal size");
@@ -987,23 +1219,13 @@ void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::map(const std::vec
 	for (; it < froms.end(); ++it, ++jt) {
 		auto& from = *it;
 		auto& to = *jt;
-		mp.insert({
-			CGAL::Plane_3<Kernel_>(
-				static_cast<NumberEpeck*>(from.get(0))->value(),
-				static_cast<NumberEpeck*>(from.get(1))->value(),
-				static_cast<NumberEpeck*>(from.get(2))->value(),
-				static_cast<NumberEpeck*>(from.get(3))->value()
-			),
-			CGAL::Plane_3<Kernel_>(
-				static_cast<NumberEpeck*>(to.get(0))->value(),
-				static_cast<NumberEpeck*>(to.get(1))->value(),
-				static_cast<NumberEpeck*>(to.get(2))->value(),
-				static_cast<NumberEpeck*>(to.get(3))->value()
-			)
-			});
+		insert_normalized_plane_map(mp, from, to);
 	}
-	auto nw = shape_->map(mp);
+	std::size_t mutated = 0;
+	auto nw = shape_->map(mp, mutated);
 	shape_ = std::move(nw);
+	apply_normalized_plane_map(mp, planes_);
+	return mutated;
 }
 
 
@@ -1012,4 +1234,4 @@ ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::
 	throw std::runtime_error("Not implemented");
 }
 
-#endif
+#endif

src/ifcgeom/kernels/cgal/CgalConversionResult.h

@@ -39,6 +39,8 @@
 #include <CGAL/Polygon_mesh_processing/compute_normal.h>
 #include <CGAL/Polygon_mesh_processing/self_intersections.h>
 
+#include <variant>
+
 #ifdef IFOPSH_SIMPLE_KERNEL
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
@@ -93,99 +95,107 @@ namespace ifcopenshell { namespace geometry {
 	using IfcGeom::OpaqueCoordinate;
 	using IfcGeom::OpaqueNumber;
 
-	using IfcGeom::add_;
-	using IfcGeom::subtract_;
-	using IfcGeom::multiply_;
-	using IfcGeom::divide_;
-	using IfcGeom::equals_;
-	using IfcGeom::less_than_;
-	using IfcGeom::negate_;
-
 #ifndef IFOPSH_SIMPLE_KERNEL
 	class IFC_GEOMLIBRARY_API NumberEpeck : public OpaqueNumber {
 	private:
-		CGAL::Epeck::FT value_;
+		struct Model : OpaqueNumber::NumberConcept {
+			CGAL::Epeck::FT value;
 
-		template <CGAL::Epeck::FT(*Fn)(CGAL::Epeck::FT, CGAL::Epeck::FT)>
-		OpaqueNumber* binary_op(OpaqueNumber* other) const {
-			auto nnd = dynamic_cast<NumberEpeck*>(other);
-			if (nnd) {
-				return new NumberEpeck(Fn(value_, nnd->value_));
-			} else {
-				return nullptr;
+			Model(const CGAL::Epeck::FT& v)
+				: value(v) {}
+
+			static const Model& as_same(const NumberConcept& other) {
+				auto same = dynamic_cast<const Model*>(&other);
+				if (same == nullptr) {
+					throw std::runtime_error("Incompatible opaque number types");
+				}
+				return *same;
 			}
-		}
 
-		template <bool(*Fn)(CGAL::Epeck::FT, CGAL::Epeck::FT)>
-		bool binary_op_bool(OpaqueNumber* other) const {
-			auto nnd = dynamic_cast<NumberEpeck*>(other);
-			if (nnd) {
-				return Fn(value_, nnd->value_);
-			} else {
-				return false;
+			virtual double to_double() const {
+				return CGAL::to_double(value);
 			}
-		}
 
-		template <CGAL::Epeck::FT(*Fn)(CGAL::Epeck::FT)>
-		OpaqueNumber* unary_op() const {
-			return new NumberEpeck(Fn(value_));
-		}
+			virtual std::string to_string() const {
+				std::stringstream ss;
+				ss << value.exact();
+				return ss.str();
+			}
+
+			virtual std::shared_ptr<const NumberConcept> add(const NumberConcept& other) const {
+				return std::make_shared<Model>(value + as_same(other).value);
+			}
+
+			virtual std::shared_ptr<const NumberConcept> subtract(const NumberConcept& other) const {
+				return std::make_shared<Model>(value - as_same(other).value);
+			}
+
+			virtual std::shared_ptr<const NumberConcept> multiply(const NumberConcept& other) const {
+				return std::make_shared<Model>(value * as_same(other).value);
+			}
+
+			virtual std::shared_ptr<const NumberConcept> divide(const NumberConcept& other) const {
+				return std::make_shared<Model>(value / as_same(other).value);
+			}
+
+			virtual std::shared_ptr<const NumberConcept> negate() const {
+				return std::make_shared<Model>(-value);
+			}
+
+			virtual std::shared_ptr<const NumberConcept> from_double(double v) const {
+				return std::make_shared<Model>(CGAL::Epeck::FT(v));
+			}
+
+			virtual std::shared_ptr<const NumberConcept> from_int(int v) const {
+                return std::make_shared<Model>(CGAL::Epeck::FT(v));
+            }
+
+			virtual bool equals(const NumberConcept& other) const {
+				return value == as_same(other).value;
+			}
+
+			virtual bool less_than(const NumberConcept& other) const {
+				return value < as_same(other).value;
+			}
+
+			virtual const std::type_info& type() const {
+				return typeid(CGAL::Epeck::FT);
+			}
+
+			virtual const void* value_ptr() const {
+				return &value;
+			}
+		};
+
 	public:
 		NumberEpeck(const CGAL::Epeck::FT& v)
-			: value_(v) {}
-
-		virtual ~NumberEpeck() { }
-
-		virtual double to_double() const {
-			return CGAL::to_double(value_);
-		}
-
-		virtual std::string to_string() const {
-			std::stringstream ss;
-			ss << value_.exact();
-			return ss.str();
-		}
+			: OpaqueNumber(std::make_shared<Model>(v)) {}
 
 		const CGAL::Epeck::FT& value() const {
-			return value_;
-		}
-
-		virtual OpaqueNumber* operator+(OpaqueNumber* other) const {
-			return binary_op<add_<CGAL::Epeck::FT>>(other);
-		}
-		virtual OpaqueNumber* operator-(OpaqueNumber* other) const {
-			return binary_op<subtract_<CGAL::Epeck::FT>>(other);
-		}
-		virtual OpaqueNumber* operator*(OpaqueNumber* other) const {
-			return binary_op<multiply_<CGAL::Epeck::FT>>(other);
-		}
-		virtual OpaqueNumber* operator/(OpaqueNumber* other) const {
-			return binary_op<divide_<CGAL::Epeck::FT>>(other);
-		}
-		virtual bool operator==(OpaqueNumber* other) const {
-			return binary_op_bool<equals_<CGAL::Epeck::FT>>(other);
-		}
-		virtual bool operator<(OpaqueNumber* other) const {
-			return binary_op_bool<less_than_<CGAL::Epeck::FT>>(other);
-		}
-		virtual OpaqueNumber* operator-() const {
-			return unary_op<negate_<CGAL::Epeck::FT>>();
-		}
-		virtual OpaqueNumber* clone() const {
-			return new NumberEpeck(value_);
+			return value_as<CGAL::Epeck::FT>();
 		}
 	};
 #endif
 
 	class IFC_GEOMLIBRARY_API CgalShape : public IfcGeom::ConversionResultShape {
 	private:
+		typedef std::variant<cgal_shape_t, cgal_point_t, cgal_wire_t> cgal_shape_storage_t;
+
 		bool convex_tag_ = false;
-		mutable std::optional<cgal_shape_t> shape_;
+		mutable std::optional<cgal_shape_storage_t> shape_;
 #ifndef IFOPSH_SIMPLE_KERNEL
 		mutable std::optional<CGAL::Nef_polyhedron_3<Kernel_>> nef_;
 #endif
-    public:
-		CgalShape(const cgal_shape_t& shape, bool convex = false);
+      public:
+#ifdef IFOPSH_SIMPLE_KERNEL
+        std::string type() const override { return "CgalSimpleShape"; }
+#else
+        std::string type() const override { return "CgalShape"; }
+#endif
+
+		CgalShape(const cgal_shape_t& shape, bool convex = false, Logger& logger = Logger::Root());
+		CgalShape(const cgal_point_t& point, bool convex = false);
+		CgalShape(const cgal_wire_t& wire, bool convex = false);
 
 #ifndef IFOPSH_SIMPLE_KERNEL
 		CgalShape(const CGAL::Nef_polyhedron_3<Kernel_>& shape, bool convex = false) {
@@ -206,8 +216,6 @@ namespace ifcopenshell { namespace geometry {
 		void to_poly() const {}
 #endif
 
-		operator const cgal_shape_t& () const { to_poly();  return *shape_; }
-		const cgal_shape_t& poly() const { to_poly(); return *shape_; }
 		virtual std::string_view backend_id() const {
 #ifdef IFOPSH_SIMPLE_KERNEL
 			return "cgal-simple";
@@ -215,12 +223,29 @@ namespace ifcopenshell { namespace geometry {
 			return "cgal";
 #endif
 		}
+		operator const cgal_shape_t& () const { return poly(); }
+		const cgal_shape_t& poly() const;
+		bool is_poly() const { return shape_ && std::holds_alternative<cgal_shape_t>(*shape_); }
+		bool is_point() const { return shape_ && std::holds_alternative<cgal_point_t>(*shape_); }
+		bool is_wire() const { return shape_ && std::holds_alternative<cgal_wire_t>(*shape_); }
+		const cgal_point_t& point() const { return std::get<cgal_point_t>(*shape_); }
+		const cgal_wire_t& wire() const { return std::get<cgal_wire_t>(*shape_); }
 
-		virtual void Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const;
+		virtual void Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id, Logger& logger = Logger::Root()) const;
 		virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const;
 
 		virtual IfcGeom::ConversionResultShape* clone() const {
-			return new CgalShape(*shape_);
+			if (shape_) {
+				return std::visit([this](const auto& value) -> IfcGeom::ConversionResultShape* {
+					return new CgalShape(value, convex_tag_);
+				}, *shape_);
+			}
+#ifndef IFOPSH_SIMPLE_KERNEL
+			if (nef_) {
+				return new CgalShape(*nef_, convex_tag_);
+			}
+#endif
+			return nullptr;
 		}
 
 		virtual bool is_manifold() const;
@@ -239,9 +264,9 @@ namespace ifcopenshell { namespace geometry {
 		// @todo this must be something with a virtual dtor so that we can delete it.
 		virtual std::pair<OpaqueCoordinate<3>, OpaqueCoordinate<3>> bounding_box() const;
 
-		virtual OpaqueNumber* length();
-		virtual OpaqueNumber* area();
-		virtual OpaqueNumber* volume();
+		virtual OpaqueNumber length();
+		virtual OpaqueNumber area();
+		virtual OpaqueNumber volume();
 
 		virtual OpaqueCoordinate<3> position();
 		virtual OpaqueCoordinate<3> axis();
@@ -262,8 +287,8 @@ namespace ifcopenshell { namespace geometry {
 		virtual ConversionResultShape* intersect(ConversionResultShape*);
 		virtual ConversionResultShape* concat(ConversionResultShape*);
 
-		virtual void map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
-		virtual void map(const std::vector<OpaqueCoordinate<4>>& from, const std::vector<OpaqueCoordinate<4>>& to);
+		virtual std::size_t map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
+		virtual std::size_t map(const std::vector<OpaqueCoordinate<4>>& from, const std::vector<OpaqueCoordinate<4>>& to);
 		virtual ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const;
 
 		virtual bool surface_area_along_direction(double tol, const ifcopenshell::geometry::taxonomy::matrix4::ptr&, double& along_x, double& along_y, double& along_z) const;
@@ -279,6 +304,8 @@ namespace ifcopenshell { namespace geometry {
 		std::list<CGAL::Plane_3<Kernel_>> planes_;
 
 	public:
+        std::string type() const override { return "CgalShapeHalfSpaceDecomposition"; }
+
 		CgalShapeHalfSpaceDecomposition(const CGAL::Nef_polyhedron_3<Kernel_>& shape, bool is_convex) {
 			if (is_convex) {
 				shape_ = std::move(build_halfspace_tree_is_decomposed(shape, planes_));
@@ -299,7 +326,7 @@ namespace ifcopenshell { namespace geometry {
 #endif
 		}
 
-		virtual void Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const;
+		virtual void Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id, Logger& logger = Logger::Root()) const;
 		virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const;
 
 		virtual int surface_genus() const;
@@ -315,9 +342,9 @@ namespace ifcopenshell { namespace geometry {
 		virtual std::pair<OpaqueCoordinate<3>, OpaqueCoordinate<3>> bounding_box() const;
 		virtual void set_box(void* b);
 
-		virtual OpaqueNumber* length();
-		virtual OpaqueNumber* area();
-		virtual OpaqueNumber* volume();
+		virtual OpaqueNumber length();
+		virtual OpaqueNumber area();
+		virtual OpaqueNumber volume();
 
 		virtual OpaqueCoordinate<3> position();
 		virtual OpaqueCoordinate<3> axis();
@@ -340,8 +367,8 @@ namespace ifcopenshell { namespace geometry {
 			return nullptr;
 		}
 
-		virtual void map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
-		virtual void map(const std::vector<OpaqueCoordinate<4>>& from, const std::vector<OpaqueCoordinate<4>>& to);
+		virtual std::size_t map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
+		virtual std::size_t map(const std::vector<OpaqueCoordinate<4>>& from, const std::vector<OpaqueCoordinate<4>>& to);
 		virtual ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const;
 
 		virtual bool surface_area_along_direction(double tol, const ifcopenshell::geometry::taxonomy::matrix4::ptr&, double& along_x, double& along_y, double& along_z) const {

src/ifcgeom/kernels/cgal/CgalKernel.cpp

@@ -46,18 +46,19 @@ namespace {
 	struct PolyhedronBuilder : public CGAL::Modifier_base<CGAL::Polyhedron_3<Kernel_>::HalfedgeDS> {
 	private:
 		std::list<cgal_face_t> *face_list;
+      Logger& logger_;
 	public:
 		std::optional<cgal_shape_t> from_soup;
-		PolyhedronBuilder(std::list<cgal_face_t> *face_list);
+		PolyhedronBuilder(std::list<cgal_face_t> *face_list, Logger& logger = Logger::Root());
 		void operator()(CGAL::Polyhedron_3<Kernel_>::HalfedgeDS &hds);
 	};
 }
 
-CGAL::Polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_polyhedron(std::list<cgal_face_t> &face_list, bool stitch_borders) {
+CGAL::Polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_polyhedron(std::list<cgal_face_t> &face_list, bool stitch_borders, Logger& logger) {
 
 	// Naive creation
 	CGAL::Polyhedron_3<Kernel_> polyhedron;
-	PolyhedronBuilder builder(&face_list);
+	PolyhedronBuilder builder(&face_list, logger);
 	polyhedron.delegate(builder);
 	if (builder.from_soup) {
 		polyhedron = *builder.from_soup;
@@ -76,7 +77,7 @@ CGAL::Polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_polyhedron(std
 
 	polyhedron.normalize_border();
 	if (!polyhedron.is_valid(false, 1)) {
-		logger::message(logger::LOG_ERROR, "create_polyhedron: Polyhedron not valid!");
+		logger.Message(Logger::LOG_ERROR, "GEO", 69, "create_polyhedron: Polyhedron not valid!");
 		//    std::ofstream fresult;
 		//    fresult.open("/Users/ken/Desktop/invalid.off");
 		//    fresult << polyhedron << std::endl;
@@ -90,31 +91,31 @@ CGAL::Polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_polyhedron(std
 }
 
 #ifndef IFOPSH_SIMPLE_KERNEL
-CGAL::Polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_polyhedron(const CGAL::Nef_polyhedron_3<Kernel_>& nef_polyhedron) {
+CGAL::Polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_polyhedron(const CGAL::Nef_polyhedron_3<Kernel_>& nef_polyhedron, Logger& logger) {
 	if (nef_polyhedron.is_simple()) {
 		try {
 			CGAL::Polyhedron_3<Kernel_> polyhedron;
 			nef_polyhedron.convert_to_polyhedron(polyhedron);
 			return polyhedron;
 		} catch (...) {
-			logger::message(logger::LOG_ERROR, "Conversion from Nef to polyhedron failed!");
+			logger.Message(Logger::LOG_ERROR, "GEO", 70, "Conversion from Nef to polyhedron failed!");
 			return CGAL::Polyhedron_3<Kernel_>();
 		}
 	} else {
-		logger::message(logger::LOG_ERROR, "Nef polyhedron not simple: cannot create polyhedron!");
+        logger.Message(Logger::LOG_ERROR, "GEO", 71, "Nef polyhedron not simple: cannot create polyhedron!");
 		return CGAL::Polyhedron_3<Kernel_>();
 	}
 }
 
-CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhedron(std::list<cgal_face_t> &face_list) {
-	CGAL::Polyhedron_3<Kernel_> polyhedron = create_polyhedron(face_list);
+CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhedron(std::list<cgal_face_t> &face_list, Logger& logger) {
+	CGAL::Polyhedron_3<Kernel_> polyhedron = create_polyhedron(face_list, true, logger);
 	if (polyhedron.is_closed()) {
 		try {
 			if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) {
 				CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
 			}
 		} catch (CGAL::Failure_exception& e) {
-			logger::message(logger::LOG_ERROR, e);
+            logger.Message(Logger::LOG_ERROR, "GEO", 72, e);
 		}
 	}
 	CGAL::Polygon_mesh_processing::triangulate_faces(polyhedron);
@@ -122,12 +123,12 @@ CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhe
 	try {
 		nef_polyhedron = CGAL::Nef_polyhedron_3<Kernel_>(polyhedron);
 	} catch (...) {
-		logger::message(logger::LOG_ERROR, "Conversion to Nef polyhedron failed!");
+        logger.Message(Logger::LOG_ERROR, "GEO", 73, "Conversion to Nef polyhedron failed!");
 	}
 	return nef_polyhedron;
 }
 
-CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhedron(CGAL::Polyhedron_3<Kernel_> &polyhedron) {
+CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhedron(CGAL::Polyhedron_3<Kernel_> &polyhedron, Logger& logger) {
 	// @todo needed?
 	polyhedron.normalize_border();
 
@@ -137,7 +138,7 @@ CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhe
 				CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron);
 			}
 		} catch (CGAL::Failure_exception& e) {
-			logger::message(logger::LOG_ERROR, e);
+			logger.Message(Logger::LOG_ERROR, "GEO", 74, e);
 		}
 	}
 
@@ -148,11 +149,11 @@ CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhe
 		try {
 			nef_polyhedron = CGAL::Nef_polyhedron_3<Kernel_>(polyhedron);
 		} catch (...) {
-			logger::message(logger::LOG_ERROR, "Conversion to Nef polyhedron failed!");
+            logger.Message(Logger::LOG_ERROR, "GEO", 75, "Conversion to Nef polyhedron failed!");
 		}
 		return nef_polyhedron;
 	} else {
-		logger::message(logger::LOG_ERROR, "Polyhedron not valid: cannot create Nef polyhedron!");
+		logger.Message(Logger::LOG_ERROR, "GEO", 76, "Polyhedron not valid: cannot create Nef polyhedron!");
 		return CGAL::Nef_polyhedron_3<Kernel_>();
 	}
 }
@@ -161,13 +162,13 @@ CGAL::Nef_polyhedron_3<Kernel_> ifcopenshell::geometry::utils::create_nef_polyhe
 bool CgalKernel::convert(const taxonomy::shell::ptr l, cgal_shape_t& shape) {
 	for (auto& f : l->children) {
 		if (f->basis && f->basis->kind() != taxonomy::PLANE) {
-			logger::error("CGAL Kernel: Non-planar faces not supported at the moment");
+			logger().Error("UNS", 3, "CGAL Kernel: Non-planar faces not supported at the moment");
 			throw not_supported_error();
 		}
 		for (auto& w : f->children) {
 			for (auto& e : w->children) {
 				if (e->basis && e->basis->kind() == taxonomy::BSPLINE_CURVE) {
-					logger::error("CGAL Kernel: B-spline edge curves not supported at the moment");
+					logger().Error("UNS", 4, "CGAL Kernel: B-spline edge curves not supported at the moment");
 					throw not_supported_error();
 				}
 			}
@@ -196,9 +197,9 @@ bool CgalKernel::convert(const taxonomy::shell::ptr l, cgal_shape_t& shape) {
 		double volume = diag(0) * diag(1) * diag(2);
 		// @todo volume van be zero also..
 		double density = num_points / volume;
-		logger::notice("Density " + boost::lexical_cast<std::string>(density), l->instance);
+		logger().Notice("GEO", 77, "Density " + boost::lexical_cast<std::string>(density), l->instance);
 		if (density > 5000) {
-			logger::notice("Substituted element with " + boost::lexical_cast<std::string>(density) + " vertices / m3 with a bounding box");
+			logger().Notice("GEO", 78, "Substituted element with " + boost::lexical_cast<std::string>(density) + " vertices / m3 with a bounding box");
 			CGAL::Point_3<Kernel_> lower(minmax.first(0), minmax.first(1), minmax.first(2));
 			CGAL::Point_3<Kernel_> upper(minmax.second(0), minmax.second(1), minmax.second(2));
 			shape = utils::create_cube(lower, upper);
@@ -215,10 +216,10 @@ bool CgalKernel::convert(const taxonomy::shell::ptr l, cgal_shape_t& shape) {
 
 		if (!success) {
 			if (this->partial_success_is_success) {
-				logger::message(logger::LOG_WARNING, "Failed to convert face, skipping:", f->instance);
+				logger().message(logger::LOG_WARNING, "Failed to convert face, skipping:", f->instance);
                 continue;
 			} else {
-				logger::message(logger::LOG_ERROR, "Failed to convert face:", f->instance);
+				logger().message(logger::LOG_ERROR, "Failed to convert face:", f->instance);
                 return false;
             }			
 		}
@@ -241,7 +242,7 @@ bool CgalKernel::convert(const taxonomy::face::ptr face, std::list<cgal_face_t>&
 	}
 
 	if (face->children.size() > 1 && num_outer_bounds > 1 && face->children.size() != num_outer_bounds) {
-		logger::message(logger::LOG_ERROR, "Invalid configuration of boundaries for:", face->instance);
+		logger().Message(Logger::LOG_ERROR, "GEO", 80, "Invalid configuration of boundaries for:", face->instance);
 		return false;
 	}
 
@@ -254,7 +255,7 @@ bool CgalKernel::convert(const taxonomy::face::ptr face, std::list<cgal_face_t>&
 
 		cgal_wire_t wire;
 		if (!convert(bound, wire)) {
-			logger::message(logger::LOG_ERROR, "Failed to process face boundary loop", bound->instance);
+            logger().Message(Logger::LOG_ERROR, "GEO", 81, "Failed to process face boundary loop", bound->instance);
 			return false;
 		}
 
@@ -708,7 +709,7 @@ bool CgalKernel::convert(const taxonomy::loop::ptr loop, cgal_wire_t& result) {
 		if (d < 1.e-5) {
 			points.erase(points.end() - 1);
 		} else {
-			logger::warning("Loop not closed", loop->instance);
+            logger().Warning("GEO", 82, "Loop not closed", loop->instance);
 		}
 	}
 
@@ -722,7 +723,7 @@ bool CgalKernel::convert(const taxonomy::loop::ptr loop, cgal_wire_t& result) {
 	// A loop should consist of at least three vertices
 	std::size_t original_count = polygon.size();
 	if (original_count < 3) {
-		logger::warning("Not enough edges for:", loop->instance);
+        logger().Warning("GEO", 83, "Not enough edges for:", loop->instance);
 		return false;
 	}
 
@@ -733,14 +734,14 @@ bool CgalKernel::convert(const taxonomy::loop::ptr loop, cgal_wire_t& result) {
 	std::size_t count = polygon.size();
 	if (original_count - count != 0) {
 		std::stringstream ss; ss << (original_count - count) << " edges removed for:";
-		logger::warning(ss.str(), loop->instance);
+        logger().Warning("GEO", 84, ss.str(), loop->instance);
 	}
 
 	{
 		std::set<cgal_point_t> visited_points;
 		for (auto& p : polygon) {
 			if (visited_points.find(p) != visited_points.end()) {
-				logger::error("Skipping self-intersecting loop", loop->instance);
+                logger().Error("GEO", 85, "Skipping self-intersecting loop", loop->instance);
 				// @todo signal somehow that occt kernel might be able to solve this
 				// @todo implement cycle detection using Arrangement_2, but that only works in exact kernel
 				return false;
@@ -762,7 +763,7 @@ bool CgalKernel::convert(const taxonomy::loop::ptr loop, cgal_wire_t& result) {
 	}
 
 	if (do_segments_intersect(segments)) {
-		logger::message(logger::LOG_WARNING, "Skipping self-intersecting loop", loop->instance);
+        logger().Message(Logger::LOG_WARNING, "GEO", 86, "Skipping self-intersecting loop", loop->instance);
 		return false;
 	}
 
@@ -790,7 +791,7 @@ bool CgalKernel::convert(const taxonomy::loop::ptr loop, cgal_wire_t& result) {
 	*/
 
 	if (count < 3) {
-		logger::message(logger::LOG_ERROR, "Not enough edges for:", loop->instance);
+        logger().Message(Logger::LOG_ERROR, "GEO", 87, "Not enough edges for:", loop->instance);
 		return false;
 	}
 
@@ -824,7 +825,7 @@ bool CgalKernel::convert_impl(const taxonomy::shell::ptr shell, ConversionResult
 
 bool CgalKernel::convert_impl(const taxonomy::solid::ptr solid, ConversionResults& results) {
 	if (solid->children.size() > 1) {
-		logger::error("Multiple shells in solid not supported at the moment");
+        logger().Error("UNS", 5, "Multiple shells in solid not supported at the moment");
 		return false;
 	}
 	cgal_shape_t shape;
@@ -969,7 +970,7 @@ bool ifcopenshell::geometry::kernels::CgalKernel::convert_openings(const express
 		try {
 			a.convert_to_polyhedron(a_poly);
 		} catch (...) {
-			logger::message(logger::LOG_ERROR, "Could not convert from Nef:", entity);
+            logger().Message(Logger::LOG_ERROR, "GEO", 88, "Could not convert from Nef:", entity);
 			return false;
 		}
 
@@ -1194,7 +1195,7 @@ bool CgalKernel::process_extrusion(const cgal_face_t& bottom_face, taxonomy::dir
 bool CgalKernel::convert(const taxonomy::extrusion::ptr extrusion, cgal_shape_t &shape) {
 	const double& height = extrusion->depth;
 	if (height < settings_.get<settings::Precision>().get()) {
-		logger::message(logger::LOG_ERROR, "Non-positive extrusion height encountered for:", extrusion->instance);
+        logger().Message(Logger::LOG_ERROR, "GEO", 89, "Non-positive extrusion height encountered for:", extrusion->instance);
 		return false;
 	}
 
@@ -1330,13 +1331,13 @@ bool CgalKernel::preprocess_boolean_operand(const express::Base& log_reference,
 	cgal_shape_t shape = shape_const;
 
 	if (!shape.is_valid()) {
-		logger::message(logger::LOG_ERROR, "Conversion to Nef will fail. Invalid geometry:", log_reference);
+        logger().Message(Logger::LOG_ERROR, "GEO", 90, "Conversion to Nef will fail. Invalid geometry:", log_reference);
 		return false;
 	}
 
 	if (!shape.is_closed()) {
 		// TODO: There can be substractions to remove parts of non-volumetric objects. Maybe iterate over all faces of an entity and put them in a Nef_polyhedron_3 through Boolean union? Highly inefficient but maybe desirable...
-		logger::message(logger::LOG_ERROR, "Subtraction of openings not supported for non-closed geometry:", log_reference);
+        logger().Message(Logger::LOG_ERROR, "UNS", 6, "Subtraction of openings not supported for non-closed geometry:", log_reference);
 		return false;
 	}
 
@@ -1345,18 +1346,18 @@ bool CgalKernel::preprocess_boolean_operand(const express::Base& log_reference,
 	try {
 		success = CGAL::Polygon_mesh_processing::triangulate_faces(shape);
 	} catch (CGAL::Failure_exception& e) {
-		logger::notice(e);
-		logger::message(logger::LOG_ERROR, "Triangulation of geometry crashed:", log_reference);
+        logger().Notice("GEO", 91, e);
+        logger().Message(Logger::LOG_ERROR, "GEO", 92, "Triangulation of geometry crashed:", log_reference);
 		return false;
 	}
 
 	if (!success) {
-		logger::message(logger::LOG_ERROR, "Triangulation of geometry failed:", log_reference);
+        logger().Message(Logger::LOG_ERROR, "GEO", 93, "Triangulation of geometry failed:", log_reference);
 		return false;
 	}
 
 	if (CGAL::Polygon_mesh_processing::does_self_intersect(shape)) {
-		logger::message(logger::LOG_ERROR, "Conversion to Nef will fail. Self-intersecting geometry:", log_reference);
+		logger().Message(Logger::LOG_ERROR, "GEO", 94, "Conversion to Nef will fail. Self-intersecting geometry:", log_reference);
 		return false;
 	}
 
@@ -1428,8 +1429,8 @@ bool CgalKernel::preprocess_boolean_operand(const express::Base& log_reference,
 	try {
 		result = CGAL::Nef_polyhedron_3<Kernel_>(shape);
 	} catch (CGAL::Failure_exception& e) {
-		logger::notice(e);
-		logger::message(logger::LOG_ERROR, "Could not convert geometry to Nef:", log_reference);
+        logger().Notice("GEO", 95, e);
+        logger().Message(Logger::LOG_ERROR, "GEO", 96, "Could not convert geometry to Nef:", log_reference);
 		return false;
 	}
 
@@ -1501,8 +1502,8 @@ bool CgalKernel::preprocess_boolean_operand(const express::Base& log_reference,
 			// @todo don't dilate in 3 dimensions but only in the XY plane, orthogonal to wall axis.
 			result = CGAL::minkowski_sum_3(result, precision_cube_);
 		} catch (CGAL::Failure_exception& e) {
-			logger::notice(e);
-			logger::message(logger::LOG_ERROR, "Could not dilate boolean operand", log_reference);
+            logger().Notice("GEO", 97, e);
+            logger().Message(Logger::LOG_ERROR, "GEO", 98, "Could not dilate boolean operand", log_reference);
 			return false;
 		}
 	}
@@ -1527,8 +1528,8 @@ bool CgalKernel::preprocess_boolean_operand(const express::Base& log_reference,
 		cgal_shape_t convert_back;
 		result.convert_to_polyhedron(convert_back);
 	} catch (CGAL::Failure_exception& e) {
-		logger::notice(e);
-		logger::message(logger::LOG_WARNING, "Final conversion will likely fail. Could not convert geometry from Nef:", log_reference);
+        logger().Notice("GEO", 99, e);
+        logger().Message(Logger::LOG_WARNING, "GEO", 100, "Final conversion will likely fail. Could not convert geometry from Nef:", log_reference);
 	}
 
 	return true;
@@ -1850,7 +1851,7 @@ bool CgalKernel::convert_impl(const taxonomy::boolean_result::ptr br, Conversion
 			// even-odd fill rule will result in incorrect results.
 			// See for example the Duplex model roof.
 
-			logger::notice("Holes are not disjoint");
+			logger().Notice("GEO", 101, "Holes are not disjoint");
 
 			CGAL::Polygon_set_2<Kernel_> result;
 			auto it = loops.begin();
@@ -1903,7 +1904,7 @@ bool CgalKernel::convert_impl(const taxonomy::boolean_result::ptr br, Conversion
 			);
 		});
 
-		logger::notice("Processed boolean operation as 2d arrangement");
+		logger().Notice("GEO", 102, "Processed boolean operation as 2d arrangement");
 
 		return true;
 
@@ -1989,7 +1990,7 @@ bool CgalKernel::convert_impl(const taxonomy::boolean_result::ptr br, Conversion
 					ps.push_back({ p.x(), p.y() });
 				}
 				if (!ps.is_simple()) {
-					logger::warning("Polygonal boundary not simple", face->children[0]->instance);
+                    logger().Warning("GEO", 103, "Polygonal boundary not simple", face->children[0]->instance);
 					continue;
 				}
 
@@ -2135,7 +2136,7 @@ bool CgalKernel::convert_impl(const taxonomy::boolean_result::ptr br, Conversion
 	try {
 		a.convert_to_polyhedron(a_poly);
 	} catch (...) {
-		logger::message(logger::LOG_ERROR, "Could not convert geometry with openings from Nef:", br->instance);
+        logger().Message(Logger::LOG_ERROR, "GEO", 104, "Could not convert geometry with openings from Nef:", br->instance);
 		return false;
 	}
 
@@ -2150,8 +2151,7 @@ bool CgalKernel::convert_impl(const taxonomy::boolean_result::ptr br, Conversion
 #endif
 }
 
-PolyhedronBuilder::PolyhedronBuilder(std::list<cgal_face_t>* face_list) {
-	this->face_list = face_list;
+PolyhedronBuilder::PolyhedronBuilder(std::list<cgal_face_t>* face_list, Logger& logger) : face_list(face_list), logger_(logger) {
 }
 
 #include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
@@ -2228,6 +2228,7 @@ void PolyhedronBuilder::operator()(CGAL::Polyhedron_3<Kernel_>::HalfedgeDS &hds)
 						// the Aff_transformation_3 stored in place to convert the 2d
 						// coords back to 3d.
 						logger::warning("Ignoring triangulated facet with novel point likely due to self-intersections");
+                        logger_.Warning("GEO", 105, "Ignoring triangulated facet with novel point likely due to self-intersections");
 						facet_vertices.erase(facet_vertices.end() - 1);
 						break;
 					}

src/ifcgeom/kernels/cgal/CgalKernel.h

@@ -20,6 +20,8 @@
 #ifndef CGAL_KERNEL_H
 #define CGAL_KERNEL_H
 
+#include "../../../ifcparse/IfcLogger.h"
+
 /*
 #ifdef NO_CACHE
 
@@ -58,12 +60,12 @@ namespace ifcopenshell {
 		namespace utils {
 			IFC_GEOMLIBRARY_API CGAL::Polyhedron_3<Kernel_> create_cube(double d);
 			IFC_GEOMLIBRARY_API CGAL::Polyhedron_3<Kernel_> create_cube(const Kernel_::Point_3& lower, const Kernel_::Point_3& upper);
-			IFC_GEOMLIBRARY_API CGAL::Polyhedron_3<Kernel_> create_polyhedron(std::list<cgal_face_t> &face_list, bool stitch_borders = false);
+			IFC_GEOMLIBRARY_API CGAL::Polyhedron_3<Kernel_> create_polyhedron(std::list<cgal_face_t> &face_list, bool stitch_borders = false, Logger& logger = Logger::Root());
 
 #ifndef IFOPSH_SIMPLE_KERNEL
-			IFC_GEOMLIBRARY_API CGAL::Polyhedron_3<Kernel_> create_polyhedron(const CGAL::Nef_polyhedron_3<Kernel_> &nef_polyhedron);
-			IFC_GEOMLIBRARY_API CGAL::Nef_polyhedron_3<Kernel_> create_nef_polyhedron(std::list<cgal_face_t> &face_list);
-			IFC_GEOMLIBRARY_API CGAL::Nef_polyhedron_3<Kernel_> create_nef_polyhedron(CGAL::Polyhedron_3<Kernel_> &polyhedron);
+			IFC_GEOMLIBRARY_API CGAL::Polyhedron_3<Kernel_> create_polyhedron(const CGAL::Nef_polyhedron_3<Kernel_> &nef_polyhedron, Logger& logger = Logger::Root());
+			IFC_GEOMLIBRARY_API CGAL::Nef_polyhedron_3<Kernel_> create_nef_polyhedron(std::list<cgal_face_t> &face_list, Logger& logger = Logger::Root());
+			IFC_GEOMLIBRARY_API CGAL::Nef_polyhedron_3<Kernel_> create_nef_polyhedron(CGAL::Polyhedron_3<Kernel_> &polyhedron, Logger& logger = Logger::Root());
 #endif
 		}
 
@@ -91,12 +93,12 @@ namespace ifcopenshell {
 #endif
 			public:
 
-				CgalKernel(const Settings& settings)
-					: AbstractKernel("cgal", settings)
+				CgalKernel(const Settings& settings, Logger& logger = Logger::Root())
+					: AbstractKernel("cgal", settings, logger)
 				{}
 
-				virtual AbstractKernel* clone() const {
-					return new CgalKernel(settings());
+				virtual AbstractKernel* clone(Logger& logger) const {
+					return new CgalKernel(settings(), logger);
 				}
 
 				virtual bool supports_boolean_operations() const {
@@ -133,4 +135,4 @@ namespace ifcopenshell {
 		}
 	}
 }
-#endif
+#endif

src/ifcgeom/kernels/cgal/nef_to_halfspace_tree.h

@@ -46,6 +46,7 @@
 #include <boost/iterator/transform_iterator.hpp>
 #include <boost/graph/copy.hpp>
 
+#include <cstddef>
 #include <list>
 #include <queue>
 #include <memory>
@@ -116,6 +117,23 @@ template <typename Kernel>
 using plane_map = std::map<typename Kernel::Plane_3, typename Kernel::Plane_3, PlaneLess<Kernel>>;
 // using plane_map = std::unordered_map<typename Kernel::Plane_3, typename Kernel::Plane_3, PlaneHash<Kernel>>;
 
+template <typename Kernel>
+typename Kernel::Plane_3 normalized_plane_for_map(const typename Kernel::Plane_3& plane) {
+	std::array<typename Kernel::FT, 3> abc{ plane.a(), plane.b(), plane.c() };
+	auto minel = std::min_element(abc.begin(), abc.end());
+	auto maxel = std::max_element(abc.begin(), abc.end());
+	auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
+	if (maxval == 0) {
+		return plane;
+	}
+	return typename Kernel::Plane_3(
+		plane.a() / maxval,
+		plane.b() / maxval,
+		plane.c() / maxval,
+		plane.d() / maxval
+	);
+}
+
 // Lexicographic comparator for CGAL Point_d (operator< is deleted in CGAL 6.x)
 struct Point_d_4d_Less {
 	using Point_d = CGAL::Epick_d<CGAL::Dimension_tag<4>>::Point_d;
@@ -264,7 +282,11 @@ class halfspace_tree {
 public:
 	virtual CGAL::Nef_polyhedron_3<Kernel> evaluate() const = 0;
 	virtual void accumulate(std::list<typename Kernel::Plane_3>&) const = 0;
-	virtual std::unique_ptr<halfspace_tree> map(const plane_map<Kernel>&) const = 0;
+	std::unique_ptr<halfspace_tree> map(const plane_map<Kernel>& m) const {
+		std::size_t ignored = 0;
+		return map(m, ignored);
+	}
+	virtual std::unique_ptr<halfspace_tree> map(const plane_map<Kernel>&, std::size_t& mutated) const = 0;
 	virtual std::string dump(int level = 0) const = 0;
 	virtual tree_type kind() const = 0;
 	virtual void merge(CGAL::Nef_polyhedron_3<Kernel>&) const = 0;
@@ -366,10 +388,10 @@ public:
 			op->accumulate(points);
 		}
 	}
-	virtual std::unique_ptr<halfspace_tree<Kernel>> map(const plane_map<Kernel>& m) const {
+	virtual std::unique_ptr<halfspace_tree<Kernel>> map(const plane_map<Kernel>& m, std::size_t& mutated) const {
 		decltype(operands_) mapped;
 		for (auto& op : operands_) {
-			mapped.emplace_back(op->map(m));
+			mapped.emplace_back(op->map(m, mutated));
 		}
 		return std::unique_ptr<halfspace_tree<Kernel>>(new halfspace_tree_nary_branch(operation_, std::move(mapped)));
 	}
@@ -485,21 +507,12 @@ public:
 	virtual void accumulate(std::list<typename Kernel::Plane_3>& points) const {
 		points.push_back(plane_);
 	}
-	virtual std::unique_ptr<halfspace_tree<Kernel>> map(const plane_map<Kernel>& m) const {
-
-		std::array<typename Kernel::FT, 3> abc{ plane_.a(), plane_.b(), plane_.c() };
-		auto minel = std::min_element(abc.begin(), abc.end());
-		auto maxel = std::max_element(abc.begin(), abc.end());
-		auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
-		CGAL::Plane_3<Kernel> pp(
-			plane_.a() / maxval,
-			plane_.b() / maxval,
-			plane_.c() / maxval,
-			plane_.d() / maxval
-		);
+	virtual std::unique_ptr<halfspace_tree<Kernel>> map(const plane_map<Kernel>& m, std::size_t& mutated) const {
+		CGAL::Plane_3<Kernel> pp = normalized_plane_for_map<Kernel>(plane_);
 
 		auto it = m.find(pp);
 		if (it != m.end()) {
+			++mutated;
 			return std::unique_ptr<halfspace_tree<Kernel>>(new halfspace_tree_plane(it->second));
 		} else {
 			return std::unique_ptr<halfspace_tree<Kernel>>(new halfspace_tree_plane(plane_));
@@ -1319,20 +1332,27 @@ size_t edge_contract(Graph<Kernel>& G) {
 // For some reason gives better results then Nef_polyhedron_3.convert_to_polyhedron() in some cases
 template <typename Kernel>
 bool convert_to_polyhedron(const CGAL::Nef_polyhedron_3<Kernel>& a, CGAL::Polyhedron_3<Kernel>& b, size_t volume_index=0) {
+	const bool all_volumes = volume_index == std::numeric_limits<size_t>::max();
 	size_t v = 0;
+	Polysoup_builder<Kernel> vis;
 	for (auto it = a.volumes_begin(); it != a.volumes_end(); ++it) {
 		if (!it->mark()) {
 			continue;
 		}
 		for (auto jt = it->shells_begin(); jt != it->shells_end(); ++jt) {
-			if (v++ == volume_index) {
-				Polysoup_builder<Kernel> vis;
+			if (v++ == volume_index || all_volumes) {
 				a.visit_shell_objects(typename CGAL::Nef_polyhedron_3<Kernel>::SFace_const_handle(jt), vis);
-				vis.build(b);
-				return true;
+				if (!all_volumes) {
+					vis.build(b);
+					return true;
+				}
 			}
 		}
 	}
+	if (all_volumes && v > 0) {
+		vis.build(b);
+		return true;
+	}
 	return false;
 }
 
@@ -1358,4 +1378,4 @@ bool write_to_obj(const CGAL::Nef_polyhedron_3<Kernel>& a, std::ostream& ofs, si
 	return volume_index == std::numeric_limits<size_t>::max();
 }
 
-#endif
+#endif