diff --git a/src/ifcconvert/validation_utils.h b/src/ifcconvert/validation_utils.h
index ecd119df76..3c0eb8b222 100644
--- a/src/ifcconvert/validation_utils.h
+++ b/src/ifcconvert/validation_utils.h
@@ -476,7 +476,7 @@ struct intersection_validator {
 			std::wcout << sss.c_str() << std::endl;
 
 			for (auto& g : geom_object->geometry()) {
-				auto s = ((ifcopenshell::geometry::CgalShape*) g.Shape())->shape();
+				cgal_shape_t s = *std::static_pointer_cast<ifcopenshell::geometry::CgalShape>(g.Shape());
 				const auto& m = g.Placement()->ccomponents();
 				const auto& n = geom_object->transformation().data()->ccomponents();
 
diff --git a/src/ifcgeom/ConversionResult.cpp b/src/ifcgeom/ConversionResult.cpp
new file mode 100644
index 0000000000..b8f657e22a
--- /dev/null
+++ b/src/ifcgeom/ConversionResult.cpp
@@ -0,0 +1,10 @@
+#include "ConversionResult.h"
+#include "IfcGeomRepresentation.h"
+
+IfcGeom::Representation::Triangulation * IfcGeom::ConversionResultShape::Triangulate(const IfcGeom::IteratorSettings & settings) const
+{
+	auto t = IfcGeom::Representation::Triangulation::empty(settings);
+	static ifcopenshell::geometry::taxonomy::matrix4 iden;
+	Triangulate(settings, iden, t, -1);
+	return t;
+}
diff --git a/src/ifcgeom/ConversionResult.h b/src/ifcgeom/ConversionResult.h
index a21f7e14a4..411f9e0845 100644
--- a/src/ifcgeom/ConversionResult.h
+++ b/src/ifcgeom/ConversionResult.h
@@ -24,6 +24,7 @@
 #include "../ifcgeom/IteratorSettings.h"
 #include "../ifcgeom/taxonomy.h"
 
+#include <memory>
 #include <vector>
 
 namespace IfcGeom {	
@@ -32,18 +33,195 @@ namespace IfcGeom {
 		class IFC_GEOM_API Triangulation;
 	}
 
+	template <typename T>
+	constexpr T add_(T a, T b) {
+		return a + b;
+	}
+
+	template <typename T>
+	constexpr T subtract_(T a, T b) {
+		return a - b;
+	}
+
+	template <typename T>
+	constexpr T multiply_(T a, T b) {
+		return a * b;
+	}
+
+	template <typename T>
+	constexpr T divide_(T a, T b) {
+		return a / b;
+	}
+
+	template <typename T>
+	constexpr bool equals_(T a, T b) {
+		return a == b;
+	}
+
+	template <typename T>
+	constexpr bool less_than_(T a, T b) {
+		return a < b;
+	}
+
+	template <typename T>
+	constexpr T negate_(T a) {
+		return -a;
+	}
+
+	class IFC_GEOM_API OpaqueNumber {
+	public:
+		virtual double to_double() const = 0;
+		virtual ~OpaqueNumber() {}
+
+		virtual OpaqueNumber* operator+(OpaqueNumber* other) const = 0;
+		virtual OpaqueNumber* operator-(OpaqueNumber* other) const = 0;
+		virtual OpaqueNumber* operator*(OpaqueNumber* other) const = 0;
+		virtual OpaqueNumber* operator/(OpaqueNumber* other) const = 0;
+		virtual bool operator==(OpaqueNumber* other) const = 0;
+		virtual bool operator<(OpaqueNumber* other) const = 0;
+		virtual OpaqueNumber* operator-() const = 0;
+	};
+
+	// @todo this can simply be a template class, to remove the need for the NumberEpeck in CGAL kernel.
+	class IFC_GEOM_API NumberNativeDouble : public OpaqueNumber {
+	private:
+		double value_;
+
+		template <double (*Fn)(double, double)>
+		OpaqueNumber* binary_op(OpaqueNumber* other) const {
+			auto nnd = dynamic_cast<NumberNativeDouble*>(other);
+			if (nnd) {
+				return new NumberNativeDouble(Fn(value_, nnd->value_));
+			} else {
+				return nullptr;
+			}
+		}
+
+		template <bool(*Fn)(double, double)>
+		bool binary_op_bool(OpaqueNumber* other) const {
+			auto nnd = dynamic_cast<NumberNativeDouble*>(other);
+			if (nnd) {
+				return new NumberNativeDouble(Fn(value_, nnd->value_));
+			} else {
+				return nullptr;
+			}
+		}
+
+		template <double(*Fn)(double)>
+		OpaqueNumber* unary_op() const {
+			return new NumberNativeDouble(Fn(value_));
+		}
+	public:
+		NumberNativeDouble(double v)
+			: value_(v) {}
+
+		virtual double to_double() const {
+			return value_;
+		}
+
+		virtual OpaqueNumber* operator+(OpaqueNumber* other) const {
+			return binary_op<add_<double>>(other);
+		}
+
+		virtual OpaqueNumber* operator-(OpaqueNumber* other) const {
+			return binary_op<subtract_<double>>(other);
+		}
+		virtual OpaqueNumber* operator*(OpaqueNumber* other) const {
+			return binary_op<multiply_<double>>(other);
+		}
+		virtual OpaqueNumber* operator/(OpaqueNumber* other) const {
+			return binary_op<divide_<double>>(other);
+		}
+		virtual bool operator==(OpaqueNumber* other) const {
+			return binary_op_bool<equals_<double>>(other);
+		}
+		virtual bool operator<(OpaqueNumber* other) const {
+			return binary_op_bool<less_than_<double>>(other);
+		}
+		virtual OpaqueNumber* operator-() const {
+			return unary_op<negate_<double>>();
+		}
+	};
+
+	template <size_t N>
+	struct IFC_GEOM_API OpaqueCoordinate {
+		std::array<std::shared_ptr<OpaqueNumber>, N> values;
+
+		template <typename... Args>
+		OpaqueCoordinate(Args... args) {
+			static_assert(sizeof...(args) == N, "Incorrect number of arguments provided");
+			init_<0>(args...);
+		}
+
+		OpaqueCoordinate() {
+			for (auto it = values.begin(); it != values.end(); ++it) {
+				*it = nullptr;
+			}
+		}
+
+		std::shared_ptr<OpaqueNumber> get(size_t i) {
+			if (i >= N) {
+				return nullptr;
+			}
+			return values[i];
+		}
+
+		void set(size_t i, std::shared_ptr<OpaqueNumber> n) {
+			if (i < N) {
+				values[i] = n;
+			}
+		}
+	private:
+		template <size_t Index, typename... Args>
+		void init_(std::shared_ptr<OpaqueNumber> value, Args... args) {
+			values[Index] = value;
+			if constexpr (Index + 1 < N) {
+				init_<Index + 1>(args...);
+			}
+		}
+	};
+
 	class IFC_GEOM_API ConversionResultShape {
 	public:
-		virtual void Triangulate(const IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, Representation::Triangulation* t, int surface_style_id) const = 0;
-
+		virtual void Triangulate(const IfcGeom::IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, Representation::Triangulation* t, int surface_style_id) const = 0;
+		IfcGeom::Representation::Triangulation* Triangulate(const IfcGeom::IteratorSettings& settings) const;
 		virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const = 0;
-		virtual ConversionResultShape* clone() const = 0;
+				
 		virtual int surface_genus() const = 0;
 		virtual bool is_manifold() const = 0;
-		// @todo this must be something with a virtual dtor so that we can delete it.
-		virtual double bounding_box(void*& b) const = 0;
+		
 		virtual int num_vertices() const = 0;
+		virtual int num_edges() const = 0;
+		virtual int num_faces() const = 0;
+			
+		// @todo choose one prototype
+		virtual double bounding_box(void*&) const = 0;
+		// @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 = 0;
 		virtual void set_box(void* b) = 0;
+		
+		virtual std::shared_ptr<OpaqueNumber> length() = 0;
+		virtual std::shared_ptr<OpaqueNumber> area() = 0;
+		virtual std::shared_ptr<OpaqueNumber> volume() = 0;
+
+		virtual OpaqueCoordinate<3> position() = 0;
+		virtual OpaqueCoordinate<3> axis() = 0;
+		virtual OpaqueCoordinate<4> plane_equation() = 0;
+
+		virtual std::vector<ConversionResultShape*> convex_decomposition() = 0;
+		virtual ConversionResultShape* halfspaces() = 0;
+		virtual ConversionResultShape* box() = 0;
+		virtual ConversionResultShape* solid() = 0;
+		virtual std::vector<ConversionResultShape*> edges() = 0;
+		virtual std::vector<ConversionResultShape*> facets() = 0;
+
+		virtual ConversionResultShape* add(ConversionResultShape*) = 0;
+		virtual ConversionResultShape* subtract(ConversionResultShape*) = 0;
+		virtual ConversionResultShape* intersect(ConversionResultShape*) = 0;
+
+		virtual void map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to) = 0;
+		virtual ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const = 0;
+
 		virtual ~ConversionResultShape() {}
 	};
 
@@ -51,7 +229,7 @@ namespace IfcGeom {
 	private:
 		int id;
 		ifcopenshell::geometry::taxonomy::matrix4::ptr placement_;
-		ConversionResultShape* shape_;
+		std::shared_ptr<ConversionResultShape> shape_;
 		ifcopenshell::geometry::taxonomy::style::ptr style_;
 	public:
 		ConversionResult(int id, ifcopenshell::geometry::taxonomy::matrix4::ptr placement, ConversionResultShape* shape, ifcopenshell::geometry::taxonomy::style::ptr style)
@@ -74,7 +252,7 @@ namespace IfcGeom {
 			// @todo verify order
 			placement_->components() = trsf->ccomponents() * placement_->ccomponents();
 		}
-		ConversionResultShape* Shape() const { return shape_; }
+		std::shared_ptr<ConversionResultShape> Shape() const { return shape_; }
 		ifcopenshell::geometry::taxonomy::matrix4::ptr Placement() const { return placement_; }
 		bool hasStyle() const { return !!style_; }
 		const ifcopenshell::geometry::taxonomy::style& Style() const { return *style_; }
diff --git a/src/ifcgeom/IfcGeomRepresentation.cpp b/src/ifcgeom/IfcGeomRepresentation.cpp
index 235d7f11f8..0df4f6ea6b 100644
--- a/src/ifcgeom/IfcGeomRepresentation.cpp
+++ b/src/ifcgeom/IfcGeomRepresentation.cpp
@@ -171,7 +171,7 @@ IfcGeom::Representation::Serialization::Serialization(const BRep& brep)
 	}
 
 	if (brep.begin() != brep.end()) {
-		if (dynamic_cast<ifcopenshell::geometry::OpenCascadeShape*>(brep.begin()->Shape())) {
+		if (std::dynamic_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(brep.begin()->Shape())) {
 			ConversionResultShape* shape = brep.as_compound();
 			ifcopenshell::geometry::taxonomy::matrix4 identity;
 			shape->Serialize(identity, brep_data_);
@@ -198,7 +198,7 @@ IfcGeom::ConversionResultShape* IfcGeom::Representation::BRep::as_compound(bool
 	builder.MakeCompound(compound);
 
 	for (auto it = begin(); it != end(); ++it) {
-		const TopoDS_Shape& s = *(ifcopenshell::geometry::OpenCascadeShape*)it->Shape();
+		const TopoDS_Shape& s = *std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape());
 
 		// @todo, check
 		gp_GTrsf trsf;
@@ -237,7 +237,7 @@ bool IfcGeom::Representation::BRep::calculate_surface_area(double& area) const {
 
 		for (IfcGeom::ConversionResults::const_iterator it = begin(); it != end(); ++it) {
 			GProp_GProps prop;
-			BRepGProp::SurfaceProperties(*(ifcopenshell::geometry::OpenCascadeShape*)it->Shape(), prop);
+			BRepGProp::SurfaceProperties(*std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape()), prop);
 			area += prop.Mass();
 		}
 
@@ -257,9 +257,9 @@ bool IfcGeom::Representation::BRep::calculate_volume(double& volume) const {
 		volume = 0.;
 
 		for (IfcGeom::ConversionResults::const_iterator it = begin(); it != end(); ++it) {
-			if (util::is_manifold(*(ifcopenshell::geometry::OpenCascadeShape*)it->Shape())) {
+			if (util::is_manifold(*std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape()))) {
 				GProp_GProps prop;
-				BRepGProp::VolumeProperties(*(ifcopenshell::geometry::OpenCascadeShape*)it->Shape(), prop);
+				BRepGProp::VolumeProperties(*std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape()), prop);
 				volume += prop.Mass();
 			} else {
 				return false;
@@ -299,9 +299,9 @@ bool IfcGeom::Representation::BRep::calculate_projected_surface_area(const ifcop
 
 		for (IfcGeom::ConversionResults::const_iterator it = begin(); it != end(); ++it) {
 			double x, y, z;
-			surface_area_along_direction(settings().deflection_tolerance(), *(ifcopenshell::geometry::OpenCascadeShape*)it->Shape(), ax, x, y, z);
+			surface_area_along_direction(settings().deflection_tolerance(), *std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape()), ax, x, y, z);
 
-			if (util::is_manifold(*(ifcopenshell::geometry::OpenCascadeShape*)it->Shape())) {
+			if (util::is_manifold(*std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape()))) {
 				x /= 2.;
 				y /= 2.;
 				z /= 2.;
@@ -414,3 +414,11 @@ void IfcGeom::Representation::Triangulation::addEdge(int n1, int n2, std::map<st
 	else edgecount[e] ++;
 	edges_temp.push_back(e);
 }
+
+const IfcGeom::ConversionResultShape* IfcGeom::Representation::BRep::item(int i) const {
+	if (i >= 0 && i < shapes_.size()) {
+		return shapes_[i].Shape()->moved(shapes_[i].Placement());
+	} else {
+		return nullptr;
+	}
+}
diff --git a/src/ifcgeom/IfcGeomRepresentation.h b/src/ifcgeom/IfcGeomRepresentation.h
index f53551e4b4..0de11dda83 100644
--- a/src/ifcgeom/IfcGeomRepresentation.h
+++ b/src/ifcgeom/IfcGeomRepresentation.h
@@ -64,6 +64,9 @@ namespace IfcGeom {
 			bool calculate_volume(double&) const;
 			bool calculate_surface_area(double&) const;
 			bool calculate_projected_surface_area(const ifcopenshell::geometry::taxonomy::matrix4& ax, double& along_x, double& along_y, double& along_z) const;
+
+			int size() const { return shapes_.size(); }
+			const IfcGeom::ConversionResultShape* item(int i) const;
 		};
 
 		class IFC_GEOM_API Serialization : public Representation  {
@@ -105,6 +108,11 @@ namespace IfcGeom {
 			size_t weld_offset_;
 			VertexKeyMap welds;
 
+			Triangulation(IfcGeom::IteratorSettings settings)
+				: Representation(IfcGeom::ElementSettings{ settings, 1., "" })
+				, weld_offset_(0)
+				{}
+
 		public:
 			const std::string& id() const { return id_; }
 			const std::vector<double>& verts() const { return _verts; }
@@ -146,6 +154,8 @@ namespace IfcGeom {
             /// @todo Very simple impl. Assumes that input vertices and normals match 1:1.
 			static std::vector<double> box_project_uvs(const std::vector<double> &vertices, const std::vector<double> &normals);
 
+			static Triangulation* empty(IfcGeom::IteratorSettings settings) { return new Triangulation(settings); }
+
 			/// Welds vertices that belong to different faces
 			int addVertex(int material_index, double X, double Y, double Z);
 
diff --git a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp
index e7bbfcf6b6..b62c0ff1dd 100644
--- a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp
+++ b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp
@@ -6,10 +6,22 @@
 #include "../../../ifcparse/IfcLogger.h"
 #include "../../../ifcgeom/IfcGeomRepresentation.h"
 
+using IfcGeom::OpaqueNumber;
+using IfcGeom::OpaqueCoordinate;
+using IfcGeom::NumberNativeDouble;
+using IfcGeom::ConversionResultShape;
+
+#ifdef IFOPSH_SIMPLE_KERNEL
+#define NumberType NumberNativeDouble
+#else
+using ifcopenshell::geometry::NumberEpeck;
+#define NumberType NumberEpeck
+#endif
+
 void ifcopenshell::geometry::CgalShape::Triangulate(const IfcGeom::IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int surface_style_id) const {
 	// Copy is made because triangulate_faces() obviously does not accept a const argument
 	// ... also becuase of transforming the vertex positions, right?
-	cgal_shape_t s = shape_;
+	cgal_shape_t s = *this;
 
 	if (!place.is_identity()) {
 		const auto& m = place.ccomponents();
@@ -71,6 +83,11 @@ void ifcopenshell::geometry::CgalShape::Triangulate(const IfcGeom::IteratorSetti
 		return;
 	}
 
+	// We do welding here in addition to in the triangulation item, because
+	// CGAL does not have a concept of vertices with identity like OCCT has.
+	typedef std::tuple<Kernel_::FT, Kernel_::FT, Kernel_::FT, Kernel_::FT, Kernel_::FT, Kernel_::FT> postion_normal;
+	std::map<postion_normal, size_t> welds;
+
 	int num_faces = 0, num_vertices = 0;
 	for (auto &face : faces(s)) {
 		if (!face->is_triangle()) {
@@ -81,20 +98,38 @@ void ifcopenshell::geometry::CgalShape::Triangulate(const IfcGeom::IteratorSetti
 		int vertexidx[3];
 		int i = 0;
 		do {
-			vertexidx[i++] = t->addVertex(surface_style_id,
-				CGAL::to_double(current_halfedge->vertex()->point().cartesian(0)),
-				CGAL::to_double(current_halfedge->vertex()->point().cartesian(1)),
-				CGAL::to_double(current_halfedge->vertex()->point().cartesian(2)));
+			postion_normal pn = {
+				current_halfedge->vertex()->point().cartesian(0),
+				current_halfedge->vertex()->point().cartesian(1),
+				current_halfedge->vertex()->point().cartesian(2),
+				face_normals_map[face].cartesian(0),
+				face_normals_map[face].cartesian(1),
+				face_normals_map[face].cartesian(2)
+			};
 
-			double nx = 0.;
-			double ny = 0.;
-			double nz = 1.;
+			// @todo normalzie based on largest component?
 
-			nx = CGAL::to_double(face_normals_map[face].cartesian(0));
-			ny = CGAL::to_double(face_normals_map[face].cartesian(1));
-			nz = CGAL::to_double(face_normals_map[face].cartesian(2));
-		
-			t->addNormal(nx, ny, nz);
+			size_t vidx;
+			auto it = welds.find(pn);
+			if (it == welds.end()) {
+				vidx = t->addVertex(
+					surface_style_id,
+					CGAL::to_double(current_halfedge->vertex()->point().cartesian(0)),
+					CGAL::to_double(current_halfedge->vertex()->point().cartesian(1)),
+					CGAL::to_double(current_halfedge->vertex()->point().cartesian(2))
+				);
+				welds.insert({ pn, vidx });
+
+				auto nx = CGAL::to_double(face_normals_map[face].cartesian(0));
+				auto ny = CGAL::to_double(face_normals_map[face].cartesian(1));
+				auto nz = CGAL::to_double(face_normals_map[face].cartesian(2));
+
+				t->addNormal(nx, ny, nz);
+			} else {
+				vidx = it->second;
+			}
+
+			vertexidx[i++] = vidx;
 
 			++num_vertices;
 			++current_halfedge;
@@ -108,7 +143,7 @@ void ifcopenshell::geometry::CgalShape::Triangulate(const IfcGeom::IteratorSetti
 }
 
 void ifcopenshell::geometry::CgalShape::Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string& r) const {
-	cgal_shape_t s = shape_;
+	cgal_shape_t s = *this;
 
 	if (!place.is_identity()) {
 		const auto& m = place.ccomponents();
@@ -137,12 +172,12 @@ 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(shape_);
+	bb += CGAL::Polygon_mesh_processing::bbox(static_cast<cgal_shape_t>(*this));
 	return (bb.xmax() - bb.xmin()) * (bb.ymax() - bb.ymin()) * (bb.zmax() - bb.zmin());
 }
 
 int ifcopenshell::geometry::CgalShape::num_vertices() const {
-	return shape_.size_of_vertices();
+	return static_cast<cgal_shape_t>(*this).size_of_vertices();
 }
 
 void ifcopenshell::geometry::CgalShape::set_box(void * b) {
@@ -153,5 +188,361 @@ void ifcopenshell::geometry::CgalShape::set_box(void * b) {
 }
 
 int ifcopenshell::geometry::CgalShape::surface_genus() const {
+	to_poly();
+	int nv = shape_->size_of_vertices();
+	int ne = shape_->size_of_halfedges() / 2;
+	int nf = shape_->size_of_facets();
+
+	const int euler = nv - ne + nf;
+	const int genus = (2 - euler) / 2;
+
+	return genus;
+}
+
+bool ifcopenshell::geometry::CgalShape::is_manifold() const {
+	// @todo ?
+	to_poly();
+	return shape_->is_valid();
+}
+
+int ifcopenshell::geometry::CgalShape::num_edges() const
+{
+	to_poly();
+	return shape_->size_of_halfedges() / 2;
+}
+
+int ifcopenshell::geometry::CgalShape::num_faces() const
+{
+	to_poly();
+	return shape_->size_of_facets();
+}
+
+std::shared_ptr<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());
+	}
+	return std::make_shared<NumberType>(len);
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::CgalShape::area()
+{
+	to_poly();
+	return std::make_shared<NumberType>(CGAL::Polygon_mesh_processing::area(*shape_));
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::CgalShape::volume()
+{
+	to_poly();
+	return std::make_shared<NumberType>(CGAL::Polygon_mesh_processing::volume(*shape_));
+}
+
+OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::position()
+{
+	throw std::runtime_error("Invalid shape type");
+}
+
+OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::axis()
+{
+	throw std::runtime_error("Invalid shape type");
+}
+
+OpaqueCoordinate<4> ifcopenshell::geometry::CgalShape::plane_equation()
+{
+	throw std::runtime_error("Invalid shape type");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShape::convex_decomposition()
+{
+#ifdef IFOPSH_SIMPLE_KERNEL
+	throw std::runtime_error("Not implemented");
+#else
+	std::vector<ConversionResultShape*> result;
+	auto copy = nef();
+	CGAL::convex_decomposition_3(copy);
+	// the first volume is the outer volume, which is
+	// ignored in the decomposition
+	auto ci = ++copy.volumes_begin();
+	int NN = 0;
+
+	for (; ci != copy.volumes_end(); ++ci, ++NN) {
+		if (ci->mark()) {
+			// @todo couldn't get it to work with the multiple volumes of a complex decomposition
+			// directly, so for now we need to isolate the individual volumes.
+			CGAL::Polyhedron_3<Kernel_> P;
+			copy.convert_inner_shell_to_polyhedron(ci->shells_begin(), P);
+			result.push_back(new CgalShape(P));
+		}
+	}
+	return result;
+#endif
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShape::halfspaces()
+{
+#ifdef IFOPSH_SIMPLE_KERNEL
+	throw std::runtime_error("Not implemented");
+#else
+	return new CgalShapeHalfSpaceDecomposition(nef());
+#endif
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShape::solid()
+{
 	throw std::runtime_error("Not implemented");
 }
+
+ConversionResultShape * ifcopenshell::geometry::CgalShape::box()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShape::edges()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShape::facets()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShape::add(ConversionResultShape* other)
+{
+#ifdef IFOPSH_SIMPLE_KERNEL
+	throw std::runtime_error("Not implemented");
+#else
+	return new CgalShape(this->nef() + ((CgalShape*)other)->nef());
+#endif
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShape::subtract(ConversionResultShape* other)
+{
+#ifdef IFOPSH_SIMPLE_KERNEL
+	throw std::runtime_error("Not implemented");
+#else
+	return new CgalShape(this->nef() - ((CgalShape*)other)->nef());
+#endif
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShape::intersect(ConversionResultShape* other)
+{
+#ifdef IFOPSH_SIMPLE_KERNEL
+	throw std::runtime_error("Not implemented");
+#else
+	return new CgalShape(this->nef() * ((CgalShape*)other)->nef());
+#endif
+}
+
+std::pair<OpaqueCoordinate<3>, OpaqueCoordinate<3>> ifcopenshell::geometry::CgalShape::bounding_box() const
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShape::moved(ifcopenshell::geometry::taxonomy::matrix4::ptr place) 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 : vertices(s)) {
+			vertex->point() = vertex->point().transform(trsf);
+		}
+	}
+
+	return new CgalShape(s);
+}
+
+void ifcopenshell::geometry::CgalShape::map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to) {
+	throw std::runtime_error("Not implemented");
+}
+
+#ifndef IFOPSH_SIMPLE_KERNEL
+
+void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::Triangulate(const IfcGeom::IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int surface_style_id) const {
+	throw std::runtime_error("Not implemented");
+}
+
+void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string& r) const {
+	throw std::runtime_error("Not implemented");
+}
+
+int ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::num_vertices() const {
+	throw std::runtime_error("Not implemented");
+}
+
+void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::set_box(void * b) {
+	throw std::runtime_error("Not implemented");
+}
+
+int ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::surface_genus() const {
+	throw std::runtime_error("Not implemented");
+}
+
+bool ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::is_manifold() const {
+	throw std::runtime_error("Not implemented");
+}
+
+int ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::num_edges() const
+{
+	throw std::runtime_error("Not implemented");
+}
+
+int ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::num_faces() const
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::CgalShapeHalfSpaceDecomposition::length()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::area()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::volume()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+OpaqueCoordinate<3> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::position()
+{
+	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>(
+			std::make_shared<NumberType>(xyz.cartesian(0)),
+			std::make_shared<NumberType>(xyz.cartesian(1)),
+			std::make_shared<NumberType>(xyz.cartesian(2))
+		);
+	} else {
+		throw std::runtime_error("Invalid shape type");
+	}
+}
+
+OpaqueCoordinate<3> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::axis()
+{
+	if (planes_.size() == 1) {
+		return OpaqueCoordinate<3>(
+			std::make_shared<NumberType>(planes_.front().a()),
+			std::make_shared<NumberType>(planes_.front().b()),
+			std::make_shared<NumberType>(planes_.front().c())
+		);
+	} else {
+		throw std::runtime_error("Invalid shape type");
+	}
+}
+
+OpaqueCoordinate<4> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::plane_equation()
+{
+	if (planes_.size() == 1) {
+		return OpaqueCoordinate<4>(
+			std::make_shared<NumberType>(planes_.front().a()),
+			std::make_shared<NumberType>(planes_.front().b()),
+			std::make_shared<NumberType>(planes_.front().c()),
+			std::make_shared<NumberType>(planes_.front().d())
+		);
+	} else {
+		throw std::runtime_error("Invalid shape type");
+	}
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::convex_decomposition()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::halfspaces()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::solid()
+{
+	return new CgalShape(shape_->evaluate());
+}
+
+ConversionResultShape * ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::box()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::edges()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::facets()
+{
+	std::vector<ConversionResultShape*> res;
+	for (auto& p : planes_) {
+		res.push_back(new CgalShapeHalfSpaceDecomposition(p));
+	}
+	return res;
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::add(ConversionResultShape* other)
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::subtract(ConversionResultShape* other)
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::intersect(ConversionResultShape* other)
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::pair<OpaqueCoordinate<3>, OpaqueCoordinate<3>> ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::bounding_box() const
+{
+	throw std::runtime_error("Not implemented");
+}
+
+double ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::bounding_box(void *& b) const {
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const
+{
+	throw std::runtime_error("Not implemented");
+}
+
+void ifcopenshell::geometry::CgalShapeHalfSpaceDecomposition::map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to) {
+	plane_map<Kernel_> mp;
+	mp.insert({
+		CGAL::Plane_3<Kernel_>(
+			std::static_pointer_cast<NumberEpeck>(from.values[0])->value(),
+			std::static_pointer_cast<NumberEpeck>(from.values[1])->value(),
+			std::static_pointer_cast<NumberEpeck>(from.values[2])->value(),
+			std::static_pointer_cast<NumberEpeck>(from.values[3])->value()
+		),
+		CGAL::Plane_3<Kernel_>(
+			std::static_pointer_cast<NumberEpeck>(to.values[0])->value(),
+			std::static_pointer_cast<NumberEpeck>(to.values[1])->value(),
+			std::static_pointer_cast<NumberEpeck>(to.values[2])->value(),
+			std::static_pointer_cast<NumberEpeck>(to.values[3])->value()
+		)
+	});
+	auto nw = shape_->map(mp);
+	shape_ = std::move(nw);
+}
+
+#endif
\ No newline at end of file
diff --git a/src/ifcgeom/kernels/cgal/CgalConversionResult.h b/src/ifcgeom/kernels/cgal/CgalConversionResult.h
index b02d259452..bda2dd90af 100644
--- a/src/ifcgeom/kernels/cgal/CgalConversionResult.h
+++ b/src/ifcgeom/kernels/cgal/CgalConversionResult.h
@@ -24,6 +24,8 @@
 
 #undef Handle
 
+#include "../../../ifcgeom/kernels/cgal/nef_to_halfspace_tree.h"
+
 #define CGAL_NO_DEPRECATED_CODE
 
 #include <boost/property_map/property_map.hpp>
@@ -86,25 +88,131 @@ typedef boost::graph_traits<CGAL::Polyhedron_3<Kernel_>>::face_descriptor cgal_f
 
 namespace ifcopenshell { namespace geometry {
 
-	class CgalShape : public IfcGeom::ConversionResultShape {
-    public:
-		CgalShape(const cgal_shape_t& shape)
-			: shape_(shape) {}
+	using IfcGeom::OpaqueCoordinate;
+	using IfcGeom::OpaqueNumber;
 
-		const cgal_shape_t& shape() const { return shape_; }
-		operator const cgal_shape_t& () { return shape_; }
+	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_GEOM_API NumberEpeck : public OpaqueNumber {
+	private:
+		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;
+			}
+		}
+
+		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 new NumberEpeck(Fn(value_, nnd->value_));
+			} else {
+				return nullptr;
+			}
+		}
+
+		template <CGAL::Epeck::FT(*Fn)(CGAL::Epeck::FT)>
+		OpaqueNumber* unary_op() const {
+			return new NumberEpeck(Fn(value_));
+		}
+	public:
+		NumberEpeck(const CGAL::Epeck::FT& v)
+			: value_(v) {}
+
+		virtual double to_double() const {
+			return CGAL::to_double(value_);
+		}
+
+		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>>();
+		}
+	};
+#endif
+
+	class CgalShape : public IfcGeom::ConversionResultShape {
+	private:
+		mutable boost::optional<cgal_shape_t> shape_;
+#ifndef IFOPSH_SIMPLE_KERNEL
+		mutable boost::optional<CGAL::Nef_polyhedron_3<Kernel_>> nef_;
+#endif
+    public:
+		CgalShape(const cgal_shape_t& shape) {
+			shape_ = shape;
+		}
+
+#ifndef IFOPSH_SIMPLE_KERNEL
+		CgalShape(const CGAL::Nef_polyhedron_3<Kernel_>& shape) {
+			nef_ = shape;
+		}
+#endif
+
+#ifndef IFOPSH_SIMPLE_KERNEL
+		void to_poly() const {
+			if (!shape_) {
+				shape_.emplace();
+				nef_->convert_to_polyhedron(*shape_);
+			}
+		}
+
+		void to_nef() const {
+			if (!nef_) {
+				nef_.emplace(*shape_);
+			}
+		}
+
+		operator const CGAL::Nef_polyhedron_3<Kernel_>& () const { to_nef(); return *nef_; }
+		const CGAL::Nef_polyhedron_3<Kernel_>& nef() const { to_nef(); return *nef_; }
+#else
+		// noop on simple kernel
+		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 void Triangulate(const IfcGeom::IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int surface_style_id) const;
-		
 		virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const;
 
 		virtual IfcGeom::ConversionResultShape* clone() const {
-			return new CgalShape(shape_);
+			return new CgalShape(*shape_);
 		}
 
-		virtual bool is_manifold() const {
-			throw std::runtime_error("Not implemented");
-		}
+		virtual bool is_manifold() const;
 
 		virtual double bounding_box(void*&) const;
 
@@ -114,10 +222,91 @@ namespace ifcopenshell { namespace geometry {
 
 		virtual int surface_genus() const;
 
-    private:
-        cgal_shape_t shape_;
+		virtual int num_edges() const;
+		virtual int num_faces() const;
+
+		// @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 std::shared_ptr<OpaqueNumber> length();
+		virtual std::shared_ptr<OpaqueNumber> area();
+		virtual std::shared_ptr<OpaqueNumber> volume();
+
+		virtual OpaqueCoordinate<3> position();
+		virtual OpaqueCoordinate<3> axis();
+		virtual OpaqueCoordinate<4> plane_equation();
+
+		virtual std::vector<ConversionResultShape*> convex_decomposition();
+		virtual ConversionResultShape* halfspaces();
+		virtual ConversionResultShape* solid();
+		virtual ConversionResultShape* box();
+		virtual std::vector<ConversionResultShape*> edges();
+		virtual std::vector<ConversionResultShape*> facets();
+
+		virtual ConversionResultShape* add(ConversionResultShape*);
+		virtual ConversionResultShape* subtract(ConversionResultShape*);
+		virtual ConversionResultShape* intersect(ConversionResultShape*);
+
+		virtual void map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
+		virtual ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const;
 	};
-    
+
+#ifndef IFOPSH_SIMPLE_KERNEL
+	class CgalShapeHalfSpaceDecomposition : public IfcGeom::ConversionResultShape {
+	private:
+		std::unique_ptr<halfspace_tree<Kernel_>> shape_;
+		std::list<CGAL::Plane_3<Kernel_>> planes_;
+
+	public:
+		CgalShapeHalfSpaceDecomposition(const CGAL::Nef_polyhedron_3<Kernel_>& shape) {
+			auto shape_copy = shape;
+			shape_ = std::move(build_halfspace_tree_decomposed(shape_copy, planes_));
+		}
+
+		CgalShapeHalfSpaceDecomposition(const CGAL::Plane_3<Kernel_>& shape) {
+			shape_.reset(new halfspace_tree_plane<Kernel_>(shape));
+			planes_.push_back(shape);
+		}
+
+		virtual void Triangulate(const IfcGeom::IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int surface_style_id) const;
+		virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const;
+
+		virtual int surface_genus() const;
+		virtual bool is_manifold() const;
+
+		virtual int num_vertices() const;
+		virtual int num_edges() const;
+		virtual int num_faces() const;
+
+		virtual double bounding_box(void*&) const;
+
+		// @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 void set_box(void* b);
+
+		virtual std::shared_ptr<OpaqueNumber> length();
+		virtual std::shared_ptr<OpaqueNumber> area();
+		virtual std::shared_ptr<OpaqueNumber> volume();
+
+		virtual OpaqueCoordinate<3> position();
+		virtual OpaqueCoordinate<3> axis();
+		virtual OpaqueCoordinate<4> plane_equation();
+
+		virtual std::vector<ConversionResultShape*> convex_decomposition();
+		virtual ConversionResultShape* halfspaces();
+		virtual ConversionResultShape* solid();
+		virtual ConversionResultShape* box();
+		virtual std::vector<ConversionResultShape*> edges();
+		virtual std::vector<ConversionResultShape*> facets();
+
+		virtual ConversionResultShape* add(ConversionResultShape*);
+		virtual ConversionResultShape* subtract(ConversionResultShape*);
+		virtual ConversionResultShape* intersect(ConversionResultShape*);
+
+		virtual void map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
+		virtual ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const;
+	};
+#endif
 }}
 
 #endif
diff --git a/src/ifcgeom/kernels/cgal/CgalKernel.cpp b/src/ifcgeom/kernels/cgal/CgalKernel.cpp
index ea4c3fad04..39511bb08d 100644
--- a/src/ifcgeom/kernels/cgal/CgalKernel.cpp
+++ b/src/ifcgeom/kernels/cgal/CgalKernel.cpp
@@ -23,7 +23,6 @@
 
 #include "../../../ifcparse/IfcLogger.h"
 #include "../../../ifcgeom/kernels/cgal/CgalConversionResult.h"
-#include "../../../ifcgeom/kernels/cgal/nef_to_halfspace_tree.h"
 
 #include <CGAL/minkowski_sum_3.h>
 #include <CGAL/exceptions.h>
@@ -817,7 +816,7 @@ bool ifcopenshell::geometry::kernels::CgalKernel::convert_openings(const IfcUtil
 	std::list<CGAL::Nef_polyhedron_3<Kernel_>> first_operands_nef, second_operands_nef;
 
 	for (auto& shp : entity_shapes) {
-		auto entity_shape = ((CgalShape*)shp.Shape())->shape();
+		cgal_shape_t entity_shape = *std::static_pointer_cast<CgalShape>(shp.Shape());
 		const auto& m = shp.Placement()->ccomponents();
 		if (!m.isIdentity()) {
 			cgal_placement_t trsf;
@@ -848,7 +847,7 @@ bool ifcopenshell::geometry::kernels::CgalKernel::convert_openings(const IfcUtil
 		AbstractKernel::convert(op.first, opening_shapes);
 
 		for (unsigned int i = 0; i < opening_shapes.size(); ++i) {
-			auto entity_shape_unlocated = ((CgalShape*)opening_shapes[i].Shape())->shape();
+			cgal_shape_t entity_shape_unlocated = *std::static_pointer_cast<CgalShape>(opening_shapes[i].Shape());
 			cgal_shape_t entity_shape(entity_shape_unlocated);
 			auto gtrsf = opening_shapes[i].Placement();
 			// @todo check
@@ -865,9 +864,7 @@ bool ifcopenshell::geometry::kernels::CgalKernel::convert_openings(const IfcUtil
 				continue;
 			}
 
-			auto graph = build_facet_edge_graph(nef);
-			auto tree = build_halfspace_tree(graph, nef);
-			tree->accumulate(all_operand_planes);
+			auto tree = build_halfspace_tree_decomposed(nef, all_operand_planes);
 
 			second_operand_instances.push_back(op.first->instance->as<IfcUtil::IfcBaseClass>());
 			second_operands.push_back(entity_shape);
@@ -1350,39 +1347,64 @@ bool CgalKernel::preprocess_boolean_operand(const IfcUtil::IfcBaseClass* log_ref
 
 	if (proc == PP_SNAP_PLANES_TO_FIRST_OPERAND) {
 		std::list<Kernel_::Plane_3> planes_fixed;
+		std::list<Kernel_::Plane_3> temp;
 		for (auto& nef : first_operands_nef) {
 			// @todo eliminate this copy (= to remove const)
 			auto nef_copy = nef;
-			auto graph = build_facet_edge_graph(nef);
-			auto tree = build_halfspace_tree(graph, nef_copy);
-			tree->accumulate(planes_fixed);
+			auto tree = build_halfspace_tree_decomposed(nef_copy, planes_fixed);  
 		}
 		{
 			// @nb we snap internally as well...
 			// @todo we can probably eliminate an evaluate() here
 			{
-				auto graph = build_facet_edge_graph(result);
 				// @todo is it deterministic enough so that rebuilding the same tree is identical/compatible?
-				auto tree = build_halfspace_tree(graph, result);
+				auto tree = build_halfspace_tree_decomposed(result, temp);
 				auto pmap = snap_halfspaces(all_operand_planes, 1.e-5);
 				result = tree->map(pmap)->evaluate();
 			}
 			{
 				std::list<Kernel_::Plane_3> planes;
-				auto graph = build_facet_edge_graph(result);
-				auto tree = build_halfspace_tree(graph, result);
-				tree->accumulate(planes);
+				auto tree = build_halfspace_tree_decomposed(result, planes);
 				auto pmap = snap_halfspaces_2(planes_fixed, planes, 1.e-5);
 				result = tree->map(pmap)->evaluate();
 			}
 		}
 	} else if (proc == PP_UNIFY_PLANES_INTERNALLY) {
 		std::list<Kernel_::Plane_3> planes;
-		auto graph = build_facet_edge_graph(result);
-		auto tree = build_halfspace_tree(graph, result);
-		tree->accumulate(planes);
-		auto pmap = snap_halfspaces(planes, 1.e-4);
-		result = tree->map(pmap)->evaluate();
+		auto tree = build_halfspace_tree_decomposed(result, planes);
+		auto pmap = snap_halfspaces(planes, 1.e-6);
+		std::wcout << tree->dump().c_str() << std::endl;
+		auto mapped = tree->map(pmap);
+		std::wcout << mapped->dump().c_str() << std::endl;
+
+		{
+			static int i = 1;
+			auto x = (halfspace_tree_nary_branch<Kernel_>*)&*tree;
+			int j = 0;
+			for (auto& a : x->operands_) {
+				auto A = a->evaluate();
+				cgal_shape_t p;
+				A.convert_to_Polyhedron(p);
+				std::string fn = "debug-orig-" + std::to_string(i) + "-" + std::to_string(j++) + ".off";
+				std::ofstream(fn.c_str()) << p;
+			}
+			i++;
+		}
+		{
+			static int i = 1;
+			auto x = (halfspace_tree_nary_branch<Kernel_>*)&*mapped;
+			int j = 0;
+			for (auto& a : x->operands_) {
+				auto A = a->evaluate();
+				cgal_shape_t p;
+				A.convert_to_Polyhedron(p);
+				std::string fn = "debug-mapped-" + std::to_string(i) + "-" + std::to_string(j++) + ".off";
+				std::ofstream(fn.c_str()) << p;
+			}
+			i++;
+		}
+
+		result = mapped->evaluate();
 	}
 
 	if (proc == PP_MINKOWSKY_DILATE) {
@@ -1397,6 +1419,22 @@ bool CgalKernel::preprocess_boolean_operand(const IfcUtil::IfcBaseClass* log_ref
 		}
 	}
 
+	
+
+	/*
+	{
+		size_t vi = 0;
+		static int i = 1;
+		std::string fn = "debug-" + std::to_string(i) + "-" + std::to_string(vi) + ".off";
+		auto ofs = std::make_unique<std::ofstream>(fn.c_str());
+		while (write_to_obj(result, *ofs, vi++)) {
+			fn = "debug-" + std::to_string(i++) + "-" + std::to_string(vi) + ".off";
+			ofs = std::make_unique<std::ofstream>(fn.c_str());
+		}
+		i += 1;
+	}
+	*/
+
 	try {
 		cgal_shape_t convert_back;
 		result.convert_to_polyhedron(convert_back);
@@ -1891,7 +1929,7 @@ bool CgalKernel::convert_impl(const taxonomy::boolean_result::ptr br, Conversion
 		}
 
 		for (auto it = cr.begin(); it != cr.end(); ++it) {
-			const cgal_shape_t& entity_shape_unlocated(((CgalShape*)it->Shape())->shape());
+			cgal_shape_t entity_shape_unlocated = *std::static_pointer_cast<CgalShape>(it->Shape());
 			cgal_shape_t entity_shape(entity_shape_unlocated);
 			if (!it->Placement()->is_identity()) {
 				cgal_placement_t trsf;
diff --git a/src/ifcgeom/kernels/cgal/nef_to_halfspace_tree.h b/src/ifcgeom/kernels/cgal/nef_to_halfspace_tree.h
index a915638c12..6e1dd61c70 100644
--- a/src/ifcgeom/kernels/cgal/nef_to_halfspace_tree.h
+++ b/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 <list>
 #include <queue>
 #include <memory>
 #include <functional>
@@ -253,6 +254,7 @@ public:
 	virtual CGAL::Nef_polyhedron_3<Kernel> evaluate(int level = 0) 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;
+	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;
 };
@@ -260,7 +262,7 @@ public:
 // Halfspace tree component as n-ary operands
 template <typename Kernel>
 class halfspace_tree_nary_branch : public halfspace_tree<Kernel> {
-private:
+public:
 	halfspace_operation operation_;
 	std::list<std::unique_ptr<halfspace_tree<Kernel>>> operands_;
 
@@ -277,10 +279,18 @@ public:
 		: operation_(operation)
 		, operands_(std::move(operands))
 	{}
-	virtual CGAL::Nef_polyhedron_3<Kernel> evaluate(int level) const {
+	std::string dump(int level) const {
 		static const char* const ops[] = { "union", "subtraction", "intersection" };
-		// std::cout << std::string(level * 2, ' ') << ops[operation_] << " (" << std::endl;
 
+		std::ostringstream ss;
+		ss << std::string(level * 2, ' ') << ops[operation_] << " (" << std::endl;
+		for (auto& op : operands_) {
+			ss << op->dump(level + 1);
+		}
+		ss << std::string(level * 2, ' ') << ")" << std::endl;
+		return ss.str();
+	}
+	virtual CGAL::Nef_polyhedron_3<Kernel> evaluate(int level) const {
 		CGAL::Nef_polyhedron_3<Kernel> result;
 
 		if (operation_ == OP_SUBTRACTION) {
@@ -337,9 +347,6 @@ public:
 			return r;
 			*/
 		}
-
-		// std::cout << std::string(level * 2, ' ') << ")" << std::endl;
-
 		return result;
 	}
 	virtual void accumulate(std::list<typename Kernel::Plane_3>& points) const {
@@ -438,9 +445,12 @@ public:
 	halfspace_tree_plane(const typename Kernel::Plane_3& plane)
 		: plane_(plane)
 	{}
+	std::string dump(int level) const {
+		std::ostringstream ss;
+		ss << std::string(level * 2, ' ') << "p " << std::setprecision(15) << plane_ << std::endl;
+		return ss.str();
+	}
 	virtual CGAL::Nef_polyhedron_3<Kernel> evaluate(int level) const {
-		// std::cout << std::string(level * 2, ' ') << "p " << plane_ << std::endl;
-
 		if constexpr(CGAL::Is_extended_kernel<Kernel>::value_type::value) {
 			throw std::runtime_error("Not implemented yet");
 			// typename Kernel::Plane_3 plane(plane_.a().exact(), plane_.b().exact(), plane_.c().exact(), plane_.d().exact());
@@ -464,7 +474,19 @@ public:
 		points.push_back(plane_);
 	}
 	virtual std::unique_ptr<halfspace_tree<Kernel>> map(const std::map<typename Kernel::Plane_3, typename Kernel::Plane_3, PlaneLess<Kernel>>& m) const {
-		auto it = m.find(plane_);
+
+		std::array<Kernel::FT, 3> abcd = { {plane_.a(), plane_.b(), plane_.c()} };
+		auto minel = std::min_element(abcd.begin(), abcd.end());
+		auto maxel = std::max_element(abcd.begin(), abcd.end());
+		auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
+		CGAL::Plane_3<Kernel> pp(
+			plane_.a() / maxval,
+			plane_.b() / maxval,
+			plane_.c() / maxval,
+			plane_.d() / maxval
+		);
+
+		auto it = m.find(pp);
 		if (it != m.end()) {
 			return std::unique_ptr<halfspace_tree<Kernel>>(new halfspace_tree_plane(it->second));
 		} else {
@@ -960,41 +982,22 @@ std::unique_ptr<halfspace_tree<TreeKernel>> build_halfspace_tree(Graph<Kernel>&
 	return std::move(tree_0);
 }
 
+// Visitor for Nef_polyhedron_3 shells to convert facets to Polyhedron_3
+// using the Polygon_mesh_processing package and Polygon_triangulation_decomposition_2
+// in case of facets with inner bounds.
 template <typename Kernel>
-size_t edge_contract(Graph<Kernel>& G) {
-	size_t n = 0;
-	typename boost::graph_traits<Graph<Kernel>>::edge_iterator ei, ei_end;
-	bool has_contracted = true;
-	while (has_contracted) {
-		has_contracted = false;
-		for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) {
-			auto srcid = boost::source(*ei, G);
-			auto tgtid = boost::target(*ei, G);
-			auto a = G[srcid].facet->plane().orthogonal_vector();
-			auto b = G[tgtid].facet->plane().orthogonal_vector();
-			// std::cout << srcid << " -- " << tgtid << std::endl << G[srcid].facet->plane() << std::endl << G[tgtid].facet->plane() << std::endl << CGAL::approximate_angle(a, b) << std::endl;
-			if (CGAL::approximate_angle(a, b) < 0.1) {
-				for (auto oe : boost::make_iterator_range(boost::out_edges(tgtid, G))) {
-					auto tt = boost::target(oe, G);
-					if (srcid != tt) { // Avoid self-loop
-						bool exists = boost::edge(srcid, tt, G).second;
-						if (!exists) {
-							boost::add_edge(srcid, tt, G);
-						}						
-					}
-				}
-				++n;
-				has_contracted = true;
-
-				boost::clear_vertex(tgtid, G);
-				boost::remove_vertex(tgtid , G);
-
-				break;
-			}
-		}
+struct Halffacet_collector {
+	std::set<typename CGAL::Nef_polyhedron_3<Kernel>::Halffacet_const_handle> facets;
+	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::Vertex_const_handle) {}
+	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::Halfedge_const_handle) {}
+	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::Halffacet_const_handle h) {
+		facets.insert(h);
 	}
-	return n;
-}
+	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::SHalfedge_const_handle) {}
+	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::SHalfloop_const_handle) {}
+	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::SFace_const_handle) {}
+};
+
 
 
 // Visitor for Nef_polyhedron_3 shells to convert facets to Polyhedron_3
@@ -1009,6 +1012,10 @@ public:
 	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::Vertex_const_handle) {}
 	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::Halfedge_const_handle) {}
 	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::Halffacet_const_handle h) {
+
+		auto verts_backup = verts;
+		auto facets_backup = facets;
+
 		boost::optional<CGAL::Polygon_with_holes_2<Kernel>> pwh;
 		auto nf = std::distance(h->facet_cycles_begin(), h->facet_cycles_end());
 		for (auto fc = h->facet_cycles_begin(); fc != h->facet_cycles_end(); ++fc) {
@@ -1025,7 +1032,7 @@ public:
 			if (nf == 1) {
 				facets.emplace_back();
 			}
-			
+
 			CGAL_For_all(hc, hc_end) {
 				auto p = hc->source()->center_vertex()->point();
 				if (nf == 1) {
@@ -1063,6 +1070,46 @@ public:
 				}
 			}
 		}
+
+		// Eliminate small slivers that create topologic connections between interior and exterior.
+		// Use connected components on polyhedron to eliminate.
+
+		std::vector<const CGAL::Point_3<Kernel>*> verts_vector(verts.size());
+		for (auto& p : verts) {
+			verts_vector[p.second] = &p.first;
+		}
+
+		// We need to normalize because we want to compare to a real-world area value
+		auto b1 = h->plane().base1();
+		double l = std::sqrt(CGAL::to_double(b1.squared_length()));
+		b1 /= l;
+
+		auto b2 = h->plane().base2();
+		l = std::sqrt(CGAL::to_double(b2.squared_length()));
+		b2 /= l;
+
+		// check for size of emitted facet, rollback if too insignificant
+		typename Kernel::FT area = 0;
+		for (auto it = facets.begin() + facets_backup.size(); it != facets.end(); ++it) {
+			CGAL::Polygon_2<Kernel> loop;
+			for (auto& i : *it) {
+				const auto& p = *verts_vector[i];
+				auto v = p - h->plane().point();
+				// std::cout << "v " << v << std::endl;
+				CGAL::Point_2<Kernel> uv(v * b1, v * b2);
+				// std::cout << "uv " << uv << std::endl;
+				loop.push_back(uv);
+			}
+			area += loop.area();
+		}
+
+		// auto pl = h->plane();
+		// l = std::sqrt(CGAL::to_double(pl.orthogonal_vector().squared_length()));
+		// std::cout << pl.a() / l << " " << pl.b() / l << " " << pl.c() / l << " " << pl.d() / l << ": " << area << std::endl;
+		if (area < 1.e-5) {
+			verts = verts_backup;
+			facets = facets_backup;
+		}
 	}
 	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::SHalfedge_const_handle) {}
 	void visit(typename CGAL::Nef_polyhedron_3<Kernel>::SHalfloop_const_handle) {}
@@ -1074,8 +1121,141 @@ public:
 		}
 		CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(verts_vector, facets, P);
 	}
+	void write(std::ostream& ofs) {
+		std::vector<CGAL::Point_3<Kernel>> verts_vector(verts.size());
+		for (auto& p : verts) {
+			verts_vector[p.second] = p.first;
+		}
+		for (auto& v : verts_vector) {
+			ofs << "v " << v.cartesian(0) << " " << v.cartesian(1) << " " << v.cartesian(2) << "\n";
+		}
+		for (auto& idxs : facets) {
+			ofs << "i";
+			for (auto& i : idxs) {
+				ofs << " " << i;
+			}
+			ofs << "\n";
+		}
+		ofs << std::flush;
+	}
 };
 
+template <typename Kernel, typename TreeKernel = Kernel>
+std::unique_ptr<halfspace_tree<TreeKernel>> build_halfspace_tree_decomposed(CGAL::Nef_polyhedron_3<Kernel>& poly_, std::list<CGAL::Plane_3<Kernel>>& planes) {
+	std::unique_ptr<halfspace_tree<TreeKernel>> tree;
+	std::list<std::unique_ptr<halfspace_tree<TreeKernel>>> root_expression;
+
+	// Don't add intermediate facets created by decomposition
+	for (auto it = poly_.halffacets_begin(); it != poly_.halffacets_end(); ++it) {
+		// but below is converted to and from vanilla polyhedron, which recomputes planes (but is still exact).
+		// Therefore, we do need to have some uniformization step, which in this case is divide by larged a,b or c
+		// component.
+		if (it->incident_volume()->mark()) {
+			std::array<Kernel::FT, 3> abcd = { {it->plane().a(), it->plane().b(), it->plane().c()} };
+			auto minel = std::min_element(abcd.begin(), abcd.end());
+			auto maxel = std::max_element(abcd.begin(), abcd.end());
+			auto maxval = ((-*minel) > *maxel) ? (-*minel) : *maxel;
+			planes.push_back(Kernel::Plane_3(
+				it->plane().a() / maxval,
+				it->plane().b() / maxval,
+				it->plane().c() / maxval,
+				it->plane().d() / maxval
+			));
+		}
+	}
+
+	auto poly = poly_;
+	CGAL::convex_decomposition_3(poly);
+	// the first volume is the outer volume, which is
+	// ignored in the decomposition
+	auto ci = ++poly.volumes_begin();
+	int NN = 0;
+
+	for (; ci != poly.volumes_end(); ++ci, ++NN) {
+		std::list<std::unique_ptr<halfspace_tree<TreeKernel>>> sub_expression;
+
+		if (ci->mark()) {
+			Halffacet_collector<Kernel> vis;
+			poly.visit_shell_objects(CGAL::Nef_polyhedron_3<Kernel>::SFace_const_handle(ci->shells_begin()), vis);
+			for (auto& f : vis.facets) {
+				sub_expression.emplace_back(new halfspace_tree_plane<TreeKernel>(f->plane()));
+			}
+
+			/*
+			// @todo couldn't get it to work with the multiple volumes of a complex decomposition
+			// directly, so for now we need to isolate the individual volumes.
+			CGAL::Polyhedron_3<Kernel> P;
+			poly.convert_inner_shell_to_polyhedron(ci->shells_begin(), P);
+			CGAL::Nef_polyhedron_3<Kernel> Pnef(P);			
+
+			for (auto it = Pnef.halffacets_begin(); it != Pnef.halffacets_end(); ++it) {
+				if (it->incident_volume()->mark()) {
+					sub_expression.emplace_back(new halfspace_tree_plane<TreeKernel>(it->plane()));
+				}
+			}
+
+			if (sub_expression.size() != vis.facets.size()) {
+
+				Polysoup_builder<Kernel> vis2;
+				poly.visit_shell_objects(CGAL::Nef_polyhedron_3<Kernel>::SFace_const_handle(ci->shells_begin()), vis2);
+
+				{
+					std::ofstream("debug-nef-decomp.off") << P;
+				}
+				{
+					CGAL::Polyhedron_3<Kernel> temp;
+					vis2.build(temp);
+					std::ofstream("debug-converted-poly.off") << temp;
+				}
+
+				// throw std::runtime_error("Unexpected");
+			}
+			*/
+		}
+
+		root_expression.emplace_back(new halfspace_tree_nary_branch<TreeKernel>(OP_INTERSECTION, std::move(sub_expression)));
+	}
+
+	tree.reset(new halfspace_tree_nary_branch<TreeKernel>(OP_UNION, std::move(root_expression)));
+	return tree;
+}
+
+template <typename Kernel>
+size_t edge_contract(Graph<Kernel>& G) {
+	size_t n = 0;
+	typename boost::graph_traits<Graph<Kernel>>::edge_iterator ei, ei_end;
+	bool has_contracted = true;
+	while (has_contracted) {
+		has_contracted = false;
+		for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) {
+			auto srcid = boost::source(*ei, G);
+			auto tgtid = boost::target(*ei, G);
+			auto a = G[srcid].facet->plane().orthogonal_vector();
+			auto b = G[tgtid].facet->plane().orthogonal_vector();
+			// std::cout << srcid << " -- " << tgtid << std::endl << G[srcid].facet->plane() << std::endl << G[tgtid].facet->plane() << std::endl << CGAL::approximate_angle(a, b) << std::endl;
+			if (CGAL::approximate_angle(a, b) < 0.1) {
+				for (auto oe : boost::make_iterator_range(boost::out_edges(tgtid, G))) {
+					auto tt = boost::target(oe, G);
+					if (srcid != tt) { // Avoid self-loop
+						bool exists = boost::edge(srcid, tt, G).second;
+						if (!exists) {
+							boost::add_edge(srcid, tt, G);
+						}						
+					}
+				}
+				++n;
+				has_contracted = true;
+
+				boost::clear_vertex(tgtid, G);
+				boost::remove_vertex(tgtid , G);
+
+				break;
+			}
+		}
+	}
+	return n;
+}
+
 // 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) {
@@ -1096,4 +1276,23 @@ bool convert_to_polyhedron(const CGAL::Nef_polyhedron_3<Kernel>& a, CGAL::Polyhe
 	return false;
 }
 
+template <typename Kernel>
+bool write_to_obj(const CGAL::Nef_polyhedron_3<Kernel>& a, std::ofstream& ofs, size_t volume_index = 0) {
+	size_t v = 0;
+	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;
+				a.visit_shell_objects(CGAL::Nef_polyhedron_3<Kernel>::SFace_const_handle(jt), vis);
+				vis.write(ofs);
+				return true;
+			}
+		}
+	}
+	return false;
+}
+
 #endif
\ No newline at end of file
diff --git a/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp b/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp
index 7b3127cddf..3a9df92d0f 100644
--- a/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp
+++ b/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp
@@ -1,13 +1,22 @@
-#include "OpenCascadeConversionResult.h"
+#include <map>
+
+#include <TopoDS.hxx>
+#include <TopExp.hxx>
+#include <BRepGProp.hxx>
+#include <GProp_GProps.hxx>
+#include <Geom_SphericalSurface.hxx>
+
+#include "OpenCascadeConversionResult.h"
 
 #include "../../../ifcparse/IfcLogger.h"
 #include "../../../ifcgeom/IfcGeomRepresentation.h"
 #include "base_utils.h"
+#include "boolean_utils.h"
 
-#include <TopoDS.hxx>
-#include <Geom_SphericalSurface.hxx>
-
-#include <map>
+using IfcGeom::OpaqueNumber;
+using IfcGeom::OpaqueCoordinate;
+using IfcGeom::NumberNativeDouble;
+using IfcGeom::ConversionResultShape;
 
 namespace {
 	// We bypass the conversion to gp_GTrsf, because it does not work
@@ -226,9 +235,190 @@ void ifcopenshell::geometry::OpenCascadeShape::Serialize(const ifcopenshell::geo
 }
 
 int ifcopenshell::geometry::OpenCascadeShape::surface_genus() const {
-	throw std::runtime_error("Not implemented");
+	return IfcGeom::util::surface_genus(shape_);
 }
 
 bool ifcopenshell::geometry::OpenCascadeShape::is_manifold() const {
+	return IfcGeom::util::is_manifold(shape_);
+}
+
+int ifcopenshell::geometry::OpenCascadeShape::num_vertices() const
+{
+	return IfcGeom::util::count(shape_, TopAbs_VERTEX);
+}
+
+int ifcopenshell::geometry::OpenCascadeShape::num_edges() const
+{
+	return IfcGeom::util::count(shape_, TopAbs_EDGE);
+}
+
+int ifcopenshell::geometry::OpenCascadeShape::num_faces() const
+{
+	return IfcGeom::util::count(shape_, TopAbs_FACE);
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::OpenCascadeShape::OpenCascadeShape::length()
+{
+	GProp_GProps prop;
+	BRepGProp::LinearProperties(shape_, prop);
+	double l = prop.Mass();
+	return std::make_shared<NumberNativeDouble>(l);
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::OpenCascadeShape::area()
+{
+	GProp_GProps prop;
+	BRepGProp::SurfaceProperties(shape_, prop);
+	double l = prop.Mass();
+	return std::make_shared<NumberNativeDouble>(l);
+}
+
+std::shared_ptr<OpaqueNumber> ifcopenshell::geometry::OpenCascadeShape::volume()
+{
+	GProp_GProps prop;
+	BRepGProp::VolumeProperties(shape_, prop);
+	double l = prop.Mass();
+	return std::make_shared<NumberNativeDouble>(l);
+}
+
+#include <Geom_Plane.hxx>
+
+OpaqueCoordinate<3> ifcopenshell::geometry::OpenCascadeShape::position()
+{
+	if (shape_.ShapeType() == TopAbs_FACE) {
+		auto surf = BRep_Tool::Surface(TopoDS::Face(shape_));
+		auto plane = Handle(Geom_Plane)::DownCast(surf);
+		if (plane) {
+			auto loc = plane->Location();
+			return OpaqueCoordinate<3>(
+				std::make_shared<NumberNativeDouble>(loc.X()),
+				std::make_shared<NumberNativeDouble>(loc.Y()),
+				std::make_shared<NumberNativeDouble>(loc.Z())
+			);
+		}
+	}
+	throw std::runtime_error("Invalid shape type");
+}
+
+OpaqueCoordinate<3> ifcopenshell::geometry::OpenCascadeShape::axis()
+{
+	if (shape_.ShapeType() == TopAbs_FACE) {
+		auto surf = BRep_Tool::Surface(TopoDS::Face(shape_));
+		auto plane = Handle(Geom_Plane)::DownCast(surf);
+		if (plane) {
+			auto dir = plane->Axis().Direction();
+			return OpaqueCoordinate<3>(
+				std::make_shared<NumberNativeDouble>(dir.X()),
+				std::make_shared<NumberNativeDouble>(dir.Y()),
+				std::make_shared<NumberNativeDouble>(dir.Z())
+			);
+		}
+	}
+	throw std::runtime_error("Invalid shape type");
+}
+
+OpaqueCoordinate<4> ifcopenshell::geometry::OpenCascadeShape::plane_equation()
+{
+	if (shape_.ShapeType() == TopAbs_FACE) {
+		auto surf = BRep_Tool::Surface(TopoDS::Face(shape_));
+		auto plane = Handle(Geom_Plane)::DownCast(surf);
+		if (plane) {
+			double a, b, c, d;
+			plane->Pln().Coefficients(a, b, c, d);
+			return OpaqueCoordinate<4>(
+				std::make_shared<NumberNativeDouble>(a),
+				std::make_shared<NumberNativeDouble>(b),
+				std::make_shared<NumberNativeDouble>(c),
+				std::make_shared<NumberNativeDouble>(d)
+			);
+		}
+	}
+	throw std::runtime_error("Invalid shape type");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::OpenCascadeShape::convex_decomposition()
+{
 	throw std::runtime_error("Not implemented");
-}
\ No newline at end of file
+}
+
+ConversionResultShape * ifcopenshell::geometry::OpenCascadeShape::halfspaces()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::OpenCascadeShape::solid()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape * ifcopenshell::geometry::OpenCascadeShape::box()
+{
+	throw std::runtime_error("Not implemented");
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::OpenCascadeShape::edges()
+{
+	TopTools_IndexedMapOfShape map;
+	TopExp::MapShapes(shape_, TopAbs_EDGE, map);
+	std::vector<ConversionResultShape*> vec;
+	for (int i = 1; i <= map.Extent(); ++i) {
+		vec.push_back(new OpenCascadeShape(map.FindKey(i)));
+	}
+	return vec;
+}
+
+std::vector<ConversionResultShape*> ifcopenshell::geometry::OpenCascadeShape::facets()
+{
+	TopTools_IndexedMapOfShape map;
+	TopExp::MapShapes(shape_, TopAbs_FACE, map);
+	std::vector<ConversionResultShape*> vec;
+	for (int i = 1; i <= map.Extent(); ++i) {
+		vec.push_back(new OpenCascadeShape(map.FindKey(i)));
+	}
+	return vec;
+}
+
+namespace {
+	ConversionResultShape* boolean_op(BOPAlgo_Operation op, const TopoDS_Shape& shape_, const TopoDS_Shape& other_shape) {
+		IfcGeom::util::boolean_settings st;
+		st.attempt_2d = true;
+		st.debug = false;
+		st.precision = 1.e-5;
+
+		TopoDS_Shape result;
+		if (IfcGeom::util::boolean_operation(st, shape_, other_shape, op, result)) {
+			return new ifcopenshell::geometry::OpenCascadeShape(result);
+		} else {
+			throw std::runtime_error("Failed to process boolean operation");
+		}
+	}
+}
+
+ConversionResultShape* ifcopenshell::geometry::OpenCascadeShape::add(ConversionResultShape* other)
+{
+	return boolean_op(BOPAlgo_FUSE, shape_, ((ifcopenshell::geometry::OpenCascadeShape*)other)->shape_);
+}
+
+ConversionResultShape* ifcopenshell::geometry::OpenCascadeShape::subtract(ConversionResultShape* other)
+{
+	return boolean_op(BOPAlgo_CUT, shape_, ((ifcopenshell::geometry::OpenCascadeShape*)other)->shape_);
+}
+
+ConversionResultShape* ifcopenshell::geometry::OpenCascadeShape::intersect(ConversionResultShape* other)
+{
+	return boolean_op(BOPAlgo_COMMON, shape_, ((ifcopenshell::geometry::OpenCascadeShape*)other)->shape_);
+}
+
+std::pair<OpaqueCoordinate<3>, OpaqueCoordinate<3>> ifcopenshell::geometry::OpenCascadeShape::bounding_box() const
+{
+	throw std::runtime_error("Not implemented");
+}
+
+ConversionResultShape* ifcopenshell::geometry::OpenCascadeShape::moved(ifcopenshell::geometry::taxonomy::matrix4::ptr t) const
+{
+	return new OpenCascadeShape(IfcGeom::util::apply_transformation(shape_, *t));
+}
+
+void ifcopenshell::geometry::OpenCascadeShape::map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to) {
+	throw std::runtime_error("Not implemented");
+}
diff --git a/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.h b/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.h
index 198ae07393..4b9600f25b 100644
--- a/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.h
+++ b/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.h
@@ -39,6 +39,9 @@
 namespace ifcopenshell {
 	namespace geometry {
 
+		using IfcGeom::OpaqueCoordinate;
+		using IfcGeom::OpaqueNumber;
+
 		class OpenCascadeShape : public IfcGeom::ConversionResultShape {
 		public:
 			OpenCascadeShape(const TopoDS_Shape& shape)
@@ -48,28 +51,51 @@ namespace ifcopenshell {
 			operator const TopoDS_Shape& () { return shape_; }
 
 			virtual void Triangulate(const IfcGeom::IteratorSettings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int surface_style_id) const;
-
 			virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const;
 
 			virtual IfcGeom::ConversionResultShape* clone() const {
 				return new OpenCascadeShape(shape_);
 			}
 
-			virtual bool is_manifold() const;
-
 			virtual double bounding_box(void*&) const {
 				throw std::runtime_error("Not implemented");
 			}
 
-			virtual int num_vertices() const {
-				throw std::runtime_error("Not implemented");
-			}
-
 			virtual void set_box(void*) {
 				throw std::runtime_error("Not implemented");
 			}
 
 			virtual int surface_genus() const;
+			virtual bool is_manifold() const;
+
+			virtual int num_vertices() const;
+			virtual int num_edges() const;
+			virtual int num_faces() const;
+
+			// @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 std::shared_ptr<OpaqueNumber> length();
+			virtual std::shared_ptr<OpaqueNumber> area();
+			virtual std::shared_ptr<OpaqueNumber> volume();
+
+			virtual OpaqueCoordinate<3> position();
+			virtual OpaqueCoordinate<3> axis();
+			virtual OpaqueCoordinate<4> plane_equation();
+
+			virtual std::vector<ConversionResultShape*> convex_decomposition();
+			virtual ConversionResultShape* halfspaces();
+			virtual ConversionResultShape* solid();
+			virtual ConversionResultShape* box();
+			virtual std::vector<ConversionResultShape*> edges();
+			virtual std::vector<ConversionResultShape*> facets();
+
+			virtual ConversionResultShape* add(ConversionResultShape*);
+			virtual ConversionResultShape* subtract(ConversionResultShape*);
+			virtual ConversionResultShape* intersect(ConversionResultShape*);
+
+			virtual void map(OpaqueCoordinate<4>& from, OpaqueCoordinate<4>& to);
+			virtual ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const;
 		private:
 			TopoDS_Shape shape_;
 		};
diff --git a/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.cpp b/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.cpp
index b41c5c40d9..4b770a3dbd 100644
--- a/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.cpp
+++ b/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.cpp
@@ -247,7 +247,7 @@ bool IfcGeom::OpenCascadeKernel::convert_openings(const IfcUtil::IfcBaseEntity*
 
 		for (unsigned int i = 0; i < opening_shapes.size(); ++i) {
 			TopoDS_Shape opening_shape_solid;
-			auto opening_shape_i = ((OpenCascadeShape*)opening_shapes[i].Shape())->shape();
+			auto opening_shape_i = std::static_pointer_cast<OpenCascadeShape>(opening_shapes[i].Shape())->shape();
 			const TopoDS_Shape& opening_shape_unlocated = util::ensure_fit_for_subtraction(opening_shape_i, opening_shape_solid, conv_settings_.getValue(ConversionSettings::GV_PRECISION));
 
 			auto gtrsf = opening_shapes[i].Placement();
@@ -277,7 +277,7 @@ bool IfcGeom::OpenCascadeKernel::convert_openings(const IfcUtil::IfcBaseEntity*
 
 		std::list<TopoDS_Shape> parts;
 
-		auto it3_shape = ((OpenCascadeShape*)it3->Shape())->shape();
+		auto it3_shape = std::static_pointer_cast<OpenCascadeShape>(it3->Shape())->shape();
 
 		bool is_multiple = it3_shape.ShapeType() == TopAbs_COMPOUND && TopoDS_Iterator(it3_shape).More() && util::is_nested_compound_of_solid(it3_shape);
 
diff --git a/src/ifcgeom/kernels/opencascade/base_utils.cpp b/src/ifcgeom/kernels/opencascade/base_utils.cpp
index 896bee9c8c..ea0e19f862 100644
--- a/src/ifcgeom/kernels/opencascade/base_utils.cpp
+++ b/src/ifcgeom/kernels/opencascade/base_utils.cpp
@@ -769,7 +769,7 @@ bool IfcGeom::util::flatten_shape_list(const IfcGeom::ConversionResults& shapes,
 
 	for (IfcGeom::ConversionResults::const_iterator it = shapes.begin(); it != shapes.end(); ++it) {
 		TopoDS_Shape merged;
-		const TopoDS_Shape& s = ((ifcopenshell::geometry::OpenCascadeShape*)it->Shape())->shape();
+		const TopoDS_Shape& s = std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape())->shape();
 		if (fuse) {
 			util::ensure_fit_for_subtraction(s, merged, tol);
 		} else {
diff --git a/src/ifcgeom/kernels/opencascade/boolean_result.cpp b/src/ifcgeom/kernels/opencascade/boolean_result.cpp
index 365419b500..9764316ee1 100644
--- a/src/ifcgeom/kernels/opencascade/boolean_result.cpp
+++ b/src/ifcgeom/kernels/opencascade/boolean_result.cpp
@@ -82,7 +82,7 @@ bool OpenCascadeKernel::convert_impl(const taxonomy::boolean_result::ptr br, Con
 		} else {
 
 			for (auto& r : cr) {
-				auto S = ((OpenCascadeShape*)r.Shape())->shape();
+				auto S = std::static_pointer_cast<OpenCascadeShape>(r.Shape())->shape();
 				gp_GTrsf trsf;
 				convert(r.Placement(), trsf);
 				// @todo it really confuses me why I cannot use Moved() here instead
diff --git a/src/ifcgeom/kernels/opencascade/layerset.cpp b/src/ifcgeom/kernels/opencascade/layerset.cpp
index 9f53f2fe38..ceaaee4805 100644
--- a/src/ifcgeom/kernels/opencascade/layerset.cpp
+++ b/src/ifcgeom/kernels/opencascade/layerset.cpp
@@ -248,7 +248,7 @@ bool IfcGeom::util::apply_folded_layerset(const ConversionResults& items, const
 
 		for (ConversionResults::const_iterator it = items.begin(); it != items.end(); ++it) {
 			TopoDS_Shape a, b;
-			if (split_solid_by_shell(((OpenCascadeShape*)it->Shape())->shape(), shells.First(), a, b, tol)) {
+			if (split_solid_by_shell(std::static_pointer_cast<OpenCascadeShape>(it->Shape())->shape(), shells.First(), a, b, tol)) {
 				result.push_back(ConversionResult(it->ItemId(), it->Placement(), new OpenCascadeShape(b), (!!styles[0] ? styles[0] : it->StylePtr())));
 				result.push_back(ConversionResult(it->ItemId(), it->Placement(), new OpenCascadeShape(a), (!!styles[1] ? styles[1] : it->StylePtr())));
 			} else {
@@ -262,7 +262,7 @@ bool IfcGeom::util::apply_folded_layerset(const ConversionResults& items, const
 
 		for (ConversionResults::const_iterator it = items.begin(); it != items.end(); ++it) {
 
-			const TopoDS_Shape& s = ((OpenCascadeShape*)it->Shape())->shape();
+			const TopoDS_Shape& s = std::static_pointer_cast<OpenCascadeShape>(it->Shape())->shape();
 			TopoDS_Solid sld;
 			ensure_fit_for_subtraction(s, sld, tol);
 
@@ -291,7 +291,7 @@ bool IfcGeom::util::apply_layerset(const ConversionResults& items, const std::ve
 
 		for (ConversionResults::const_iterator it = items.begin(); it != items.end(); ++it) {
 			TopoDS_Shape a, b;
-			if (split_solid_by_surface(((OpenCascadeShape*)it->Shape())->shape(), surfaces[1], a, b, tol)) {
+			if (split_solid_by_surface(std::static_pointer_cast<OpenCascadeShape>(it->Shape())->shape(), surfaces[1], a, b, tol)) {
 				result.push_back(ConversionResult(it->ItemId(), it->Placement(),new OpenCascadeShape(b), (!!styles[0] ? styles[0] : it->StylePtr())));
 				result.push_back(ConversionResult(it->ItemId(), it->Placement(),new OpenCascadeShape(a), (!!styles[1] ? styles[1] : it->StylePtr())));
 			} else {
@@ -334,7 +334,7 @@ bool IfcGeom::util::apply_layerset(const ConversionResults& items, const std::ve
 
 		for (ConversionResults::const_iterator it = items.begin(); it != items.end(); ++it) {
 
-			const TopoDS_Shape& s = ((OpenCascadeShape*)it->Shape())->shape();
+			const TopoDS_Shape& s = std::static_pointer_cast<OpenCascadeShape>(it->Shape())->shape();
 			TopoDS_Solid sld;
 			ensure_fit_for_subtraction(s, sld, tol);
 
diff --git a/src/ifcgeom/mapping/IfcObjectPlacement.cpp b/src/ifcgeom/mapping/IfcObjectPlacement.cpp
index 5c9539bfda..ae817bb424 100644
--- a/src/ifcgeom/mapping/IfcObjectPlacement.cpp
+++ b/src/ifcgeom/mapping/IfcObjectPlacement.cpp
@@ -67,13 +67,13 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcObjectPlacement* inst) {
 
 	taxonomy::ptr result;
 	if (!parent_placement_ignored && relative_to) {
-		// The parent placement of the current is a placement for a type that is
-		// being ignored (Site or Building) or it is the host element of an opening.
         result = taxonomy::make<taxonomy::matrix4>(
 			taxonomy::cast<taxonomy::matrix4>(map(relative_to))->ccomponents() *
 			taxonomy::cast<taxonomy::matrix4>(map(transform))->ccomponents()
 		);
 	} else {
+		// The parent placement of the current is a placement for a type that is
+		// being ignored (Site or Building) or it is the host element of an opening.
 		result = map(transform);
 	}
 
diff --git a/src/ifcwrap/IfcGeomWrapper.i b/src/ifcwrap/IfcGeomWrapper.i
index e538ac1f70..466a03ba54 100644
--- a/src/ifcwrap/IfcGeomWrapper.i
+++ b/src/ifcwrap/IfcGeomWrapper.i
@@ -60,6 +60,7 @@
 
 %include "../ifcgeom/ifc_geom_api.h"
 %include "../ifcgeom/Converter.h"
+%include "../ifcgeom/ConversionResult.h"
 %include "../ifcgeom/IteratorSettings.h"
 %include "../ifcgeom/IfcGeomElement.h"
 %include "../ifcgeom/IfcGeomRepresentation.h"
@@ -602,6 +603,9 @@ struct ShapeRTTI : public boost::static_visitor<PyObject*>
 %template(svg_loop) std::vector<std::array<double, 2>>;
 %template(svg_loops) std::vector<std::vector<std::array<double, 2>>>;
 
+%template(OpaqueCoordinate_3) IfcGeom::OpaqueCoordinate<3>;
+%template(OpaqueCoordinate_4) IfcGeom::OpaqueCoordinate<4>;
+
 %naturalvar svgfill::polygon_2::boundary;
 %naturalvar svgfill::polygon_2::inner_boundaries;
 %naturalvar svgfill::polygon_2::point_inside;
diff --git a/src/ifcwrap/IfcPython.i b/src/ifcwrap/IfcPython.i
index f7f233e4ba..0ab31e71c0 100644
--- a/src/ifcwrap/IfcPython.i
+++ b/src/ifcwrap/IfcPython.i
@@ -52,6 +52,10 @@
 %include "std_vector.i"
 %include "std_string.i"
 %include "exception.i"
+%include "std_shared_ptr.i"
+
+%shared_ptr(IfcGeom::OpaqueNumber);
+%ignore IfcGeom::NumberNativeDouble;
 
 // General python-specific rename rules for comparison operators.
 // Mostly to silence warnings, but might be of use some time.
@@ -130,6 +134,8 @@
 	#include "../ifcparse/IfcSchema.h"
 	#include "../ifcparse/utils.h"
 
+	#include "../ifcgeom/ConversionResult.h"
+
 	#include "../svgfill/src/svgfill.h"
 %}
 
@@ -152,6 +158,7 @@
 	#include <BRepTools_ShapeSet.hxx>
 #endif
 	#include "../ifcgeom/Iterator.h"
+	#include "../ifcgeom/ConversionResult.h"
 
 	#include "../serializers/SvgSerializer.h"
 	#include "../serializers/WavefrontObjSerializer.h"
diff --git a/src/ifcwrap/utils/type_conversion.i b/src/ifcwrap/utils/type_conversion.i
index 73013b2b1e..3ecd6c6660 100644
--- a/src/ifcwrap/utils/type_conversion.i
+++ b/src/ifcwrap/utils/type_conversion.i
@@ -138,6 +138,7 @@
 	PyObject* pythonize(const IfcParse::inverse_attribute* t) { return SWIG_NewPointerObj(SWIG_as_voidptr(t), SWIGTYPE_p_IfcParse__inverse_attribute, 0); }
 	PyObject* pythonize(const IfcParse::entity* t)      { return SWIG_NewPointerObj(SWIG_as_voidptr(t), SWIGTYPE_p_IfcParse__entity, 0);             }
 	PyObject* pythonize(const IfcParse::declaration* t) { return SWIG_NewPointerObj(SWIG_as_voidptr(t), declaration_type_to_swig(t), 0);             }
+	PyObject* pythonize(const IfcGeom::ConversionResultShape* t) { return SWIG_NewPointerObj(SWIG_as_voidptr(t), SWIGTYPE_p_IfcGeom__ConversionResultShape, 0); }
 	// NB: This cannot be temporary as a Python object is constructed from a pointer to the address of this object
 	// PyObject* pythonize(const IfcGeom::Material& t)     { return SWIG_NewPointerObj(SWIG_as_voidptr(&t), SWIGTYPE_p_IfcGeom__Material, 0);           }
 	
diff --git a/src/ifcwrap/utils/typemaps_out.i b/src/ifcwrap/utils/typemaps_out.i
index 3574e49ed8..e8ea58f2b0 100644
--- a/src/ifcwrap/utils/typemaps_out.i
+++ b/src/ifcwrap/utils/typemaps_out.i
@@ -146,3 +146,4 @@ CREATE_VECTOR_TYPEMAP_OUT(IfcParse::attribute const *)
 CREATE_VECTOR_TYPEMAP_OUT(IfcParse::inverse_attribute const *)
 CREATE_VECTOR_TYPEMAP_OUT(IfcParse::entity const *)
 CREATE_VECTOR_TYPEMAP_OUT(IfcParse::declaration const *)
+CREATE_VECTOR_TYPEMAP_OUT(IfcGeom::ConversionResultShape *)
diff --git a/src/serializers/HdfSerializer.cpp b/src/serializers/HdfSerializer.cpp
index f2235c5a36..e0b40ed4bf 100644
--- a/src/serializers/HdfSerializer.cpp
+++ b/src/serializers/HdfSerializer.cpp
@@ -593,7 +593,7 @@ void HdfSerializer::write(const IfcGeom::BRepElement* o) {
 		}
 
 		brep_strings.emplace_back();
-		write_shape(((ifcopenshell::geometry::OpenCascadeShape*)it->Shape())->shape(), brep_strings.back());
+		write_shape(std::static_pointer_cast<ifcopenshell::geometry::OpenCascadeShape>(it->Shape())->shape(), brep_strings.back());
 		
 		parts[i].surface_style = { "", "", 0, {nan,nan,nan}, {nan,nan,nan}, nan, nan };
 		if (it->hasStyle()) {