wip

src/ifcgeom/kernels/opencascade/OpenCascadeKernel.h

@@ -136,6 +136,10 @@ public:
 
 	typedef boost::variant<boost::blank, Handle(Geom_Curve), TopoDS_Wire> curve_creation_visitor_result_type;
 	curve_creation_visitor_result_type convert_curve(const ifcopenshell::geometry::taxonomy::ptr);
+
+	std::map<uint32_t, curve_creation_visitor_result_type> curve_cache_, edge_curve_cache_;
+    std::map<std::tuple<double, double, double>, TopoDS_Vertex> vertex_cache_;
+
 	Handle(Geom_Surface) convert_surface(const ifcopenshell::geometry::taxonomy::ptr);
 
 	template <typename T, typename U>

src/ifcgeom/kernels/opencascade/face.cpp

@@ -47,6 +47,9 @@
 #include <Geom_SurfaceOfLinearExtrusion.hxx>
 #include <Geom_SurfaceOfRevolution.hxx>
 #include <BRepPrimAPI_MakeRevol.hxx>
+#include <GeomProjLib.hxx>
+#include <BRepTools_WireExplorer.hxx>
+#include <Geom2d_TrimmedCurve.hxx>
 
 #if OCC_VERSION_HEX < 0x70600
 #include <BRepAdaptor_HCompCurve.hxx>
@@ -431,6 +434,155 @@ bool OpenCascadeKernel::convert(const taxonomy::face::ptr face, TopoDS_Shape& re
 		if (reversed_surface) {
 			surf = surf->UReversed();
 		}
+
+		if (!(surf->DynamicType() == STANDARD_TYPE(Geom_Plane))) {
+            for (auto& w : fd.wires()) {
+                bool has_seam = false;
+                double u0, u1;
+                std::set<Geom_Curve*> curves;
+                for (TopoDS_Iterator it(w); it.More(); it.Next()) {
+                    const auto& ed = TopoDS::Edge(it.Value());
+                    auto crv = BRep_Tool::Curve(ed, u0, u1);
+                    if (curves.find(crv.get()) != curves.end()) {
+                        has_seam = true;
+                    } else {
+                        curves.insert(crv.get());
+                    }
+                }
+
+				if (has_seam) {
+                    BRep_Builder BB;
+
+					std::vector<TopoDS_Edge> edges;
+                    for (BRepTools_WireExplorer it(w); it.More(); it.Next()) {
+                        const auto& ed = it.Current();
+                        edges.push_back(ed);
+                    }
+
+					for (auto it = edges.begin(); it != edges.end(); ++it) {
+                        const auto& ed = *it;
+
+                        TopLoc_Location L3d;
+                        Standard_Real f3d = 0.0, l3d = 0.0;
+                        Handle(Geom_Curve) c3d = BRep_Tool::Curve(ed, L3d, f3d, l3d); // may be null
+                        
+						std::ostringstream oss;
+
+						oss << "Ed orien: " << ed.Orientation() << std::endl;
+
+                        oss << "3d: " << f3d << " - " << l3d << " ";
+                        c3d->DumpJson(oss);
+
+						for (Standard_Integer idx = 1;; ++idx) {
+                            Handle(Geom2d_Curve) c2d;
+                            Handle(Geom_Surface) s;
+                            TopLoc_Location L;
+                            Standard_Real f = 0.0, l = 0.0;
+
+                            BRep_Tool::CurveOnSurface(ed, c2d, s, L, f, l, idx);
+
+							// "C and S are null if the index is out of range"
+                            if (c2d.IsNull() || s.IsNull()) {
+                                break;
+                            }
+							
+							oss << idx << " ";
+                            c2d->DumpJson(oss);
+                        }
+
+                        // Build 2D curve for the segment [f3d,l3d] on surface S
+                        Handle(Geom2d_Curve) base2d = GeomProjLib::Curve2d(c3d, f3d, l3d, fd.surface());
+                        if (base2d.IsNull()) {
+                            throw std::runtime_error("Failed to project 3D curve to 2D on surface with seam");
+                        }
+
+						Handle(Geom2d_TrimmedCurve) trim2d = new Geom2d_TrimmedCurve(base2d, f3d, l3d);
+
+						auto Pf = trim2d->Value(trim2d->FirstParameter());
+                        auto Pl = trim2d->Value(trim2d->LastParameter());
+
+												oss << std::endl
+                            << Pf.X() << "," << Pf.Y() << " to " << Pl.X() << "," << Pl.Y() << std::endl;
+
+
+						BB.UpdateEdge(ed, trim2d, fd.surface(), L3d, 1.e-7, Pf, Pl);
+
+						oss << "2d: ";
+                        trim2d->DumpJson(oss);
+
+						auto s = oss.str();
+                        std::wcout << s.c_str() << "\n\n\n\n" << std::endl;
+                    }
+
+					for (auto it = edges.begin(); it != edges.end(); ++it) {
+                        const auto& ed = *it;
+                        const auto& prev = (it == edges.begin()) ? *(edges.end() - 1) : *(it - 1);
+                        const auto& next = (it == (edges.end() - 1)) ? *(edges.begin()) : *(it + 1);
+
+                        std::ostringstream oss;
+
+						TopLoc_Location L3d;
+                        Standard_Real f3d = 0.0, l3d = 0.0;
+                        Handle(Geom_Curve) c3d = BRep_Tool::Curve(ed, L3d, f3d, l3d); // may be null
+
+                        auto idx = 1;
+                        Handle(Geom2d_Curve) c2d, prev_c2d;
+                        Handle(Geom_Surface) surf;
+                        TopLoc_Location L;
+                        Standard_Real f = 0.0, l = 0.0;
+
+                        BRep_Tool::CurveOnSurface(ed, c2d, surf, L, f, l, idx);
+
+						if (c2d.IsNull()) {
+							throw std::runtime_error("Failed to retrieve 2D curve on surface with seam");
+                        }
+
+						if (surf != fd.surface()) {
+							throw std::runtime_error("Unexpected surface on edge with seam");
+                        }
+
+						auto Pf = c2d->Value(c2d->FirstParameter());
+                        auto Pl = c2d->Value(c2d->LastParameter());
+
+						if (ed.Orientation() == TopAbs_REVERSED) {
+                            std::swap(Pf, Pl);
+                        }
+
+												oss << std::endl
+                            << Pf.X() << "," << Pf.Y() << " to " << Pl.X() << "," << Pl.Y() << std::endl;
+
+
+						BRep_Tool::CurveOnSurface(prev, prev_c2d, surf, L, f, l, idx);
+                                                auto PPf = prev_c2d->Value(prev_c2d->FirstParameter());
+						auto PPl = prev_c2d->Value(prev_c2d->LastParameter());
+
+						if (prev.Orientation() == TopAbs_REVERSED) {
+							std::swap(PPf, PPl);
+                        }
+
+						if (Pl.Distance(PPl) < 1.e-7) {
+							continue;
+                        }
+
+						gp_Trsf2d tr;
+                        tr.SetTranslation(PPl, Pf);
+                        auto updated = Handle(Geom2d_Curve)::DownCast(c2d->Transformed(tr));
+
+                        updated->DumpJson(oss);
+                        BB.UpdateEdge(ed, updated, fd.surface(), L3d, 1.e-7, Pf, Pl);
+						
+						oss << std::endl << Pf.X() << "," << Pf.Y() << " to " << Pl.X() << "," << Pl.Y() << std::endl;
+
+						auto s = oss.str();
+                        std::wcout << s.c_str() << "\n\n\n\n"
+                                   << std::endl;
+
+                    }
+                }
+            }
+
+		}
+
 		BRepBuilderAPI_MakeFace mf(surf, fd.outer_wire());
 
 		if (mf.IsDone()) {
@@ -501,6 +653,8 @@ bool OpenCascadeKernel::convert(const taxonomy::face::ptr face, TopoDS_Shape& re
 			// For planar faces, Open Cascade generates p-curves on the fly.
 
 			for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) {
+                continue;
+
 				ShapeFix_Shape sfs(it.Value());
 
 				Handle(ShapeExtend_MsgRegistrator) msg;
@@ -534,6 +688,9 @@ bool OpenCascadeKernel::convert(const taxonomy::face::ptr face, TopoDS_Shape& re
 		}
 
 		for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) {
+            continue;
+
+
 			const TopoDS_Face& occ_face = TopoDS::Face(it.Value());
 
 			ShapeFix_Face sfs(TopoDS::Face(occ_face));
@@ -559,6 +716,8 @@ bool OpenCascadeKernel::convert(const taxonomy::face::ptr face, TopoDS_Shape& re
 		}
 
 		for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) {
+            continue;
+
 			TopoDS_Face& occ_face = TopoDS::Face(it.Value());
 
 			bool all_reversed = true;

src/ifcgeom/kernels/opencascade/loop.cpp

@@ -19,6 +19,10 @@
 #include <TopTools_ListOfShape.hxx>
 #include <BRepTools_WireExplorer.hxx>
 #include <BRepBuilderAPI_Transform.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <ShapeFix_Shape.hxx>
+
 
 #include <Standard_Version.hxx>
 #if OCC_VERSION_HEX < 0x70600
@@ -35,6 +39,18 @@ namespace {
 		OpenCascadeKernel* kernel;
 		OpenCascadeKernel::curve_creation_visitor_result_type result;
 
+		const TopoDS_Vertex create_cached_vertex(const taxonomy::point3::ptr& pnt) {
+            auto cs = pnt->components();
+            auto it = kernel->vertex_cache_.find({cs(0), cs(1), cs(2)});
+            if (it == kernel->vertex_cache_.end()) {
+				TopoDS_Vertex v = BRepBuilderAPI_MakeVertex(OpenCascadeKernel::convert_xyz<gp_Pnt>(*pnt)).Vertex();
+				kernel->vertex_cache_[{cs(0), cs(1), cs(2)}] = v;
+				return v;
+			} else {
+				return it->second;
+            }
+        }
+
 		OpenCascadeKernel::curve_creation_visitor_result_type operator()(const taxonomy::bspline_curve::ptr& bc) {
 
 			const bool is_rational = !!bc->weights;
@@ -155,13 +171,44 @@ namespace {
 				if (e_start.which() == 0) {
 					E = BRepBuilderAPI_MakeEdge(curve).Edge();
 				} else if (e_start.which() == 1) {
-					auto p1 = OpenCascadeKernel::convert_xyz<gp_Pnt>(*boost::get<taxonomy::point3::ptr>(e_start));
-					auto p2 = OpenCascadeKernel::convert_xyz<gp_Pnt>(*boost::get<taxonomy::point3::ptr>(e_end));
+                    auto p1 = create_cached_vertex(boost::get<taxonomy::point3::ptr>(e_start));
+                    auto p2 = create_cached_vertex(boost::get<taxonomy::point3::ptr>(e_end));
 
-					if (curve->IsClosed() && p1.Distance(p2) <= kernel->settings().get<settings::Precision>().get()) {
-						E = BRepBuilderAPI_MakeEdge(curve).Edge();
+					// Closed or Periodic, what do in case of Bspline?
+					// p1.Distance(p2) <= kernel->settings().get<settings::Precision>().get()
+					if (curve->IsClosed() && p1.IsSame(p2)) {
+                        if (curve->IsPeriodic()) {
+                            auto pp1 = BRep_Tool::Pnt(p1);
+							GeomAPI_ProjectPointOnCurve proj(pp1, curve);
+							auto u0 = proj.Parameter(1);
+							auto u1 = u0 + curve->Period();
+
+							E = BRepBuilderAPI_MakeEdge(curve, p1, p1, u0, u1);
+                        } else {
+                            E = BRepBuilderAPI_MakeEdge(curve);
+						}
 					} else {
-						E = BRepBuilderAPI_MakeEdge(curve, p1, p2).Edge();
+                        double u0 = 0., u1 = 0.;
+                        auto pp1 = BRep_Tool::Pnt(p1);
+                        auto pp2 = BRep_Tool::Pnt(p2);
+
+                            try {
+                                GeomAPI_ProjectPointOnCurve proj(pp1, curve);
+								u0 = proj.Parameter(1);
+                            }
+                            catch (const Standard_Failure&) {
+                            }
+                            try {
+                                GeomAPI_ProjectPointOnCurve proj(pp2, curve);
+                                u1 = proj.Parameter(1);
+                            } catch (const Standard_Failure&) {
+                            }
+
+						if (u1 > u0) {                            
+							E = BRepBuilderAPI_MakeEdge(curve, p1, p2).Edge();
+						} else {
+                            E = TopoDS::Edge(BRepBuilderAPI_MakeEdge(curve, p2, p1).Edge().Reversed());
+                        }
 					}
 				} else if (e_start.which() == 2) {
 					auto v1 = boost::get<double>(e_start);
@@ -222,8 +269,33 @@ namespace {
 }
 
 OpenCascadeKernel::curve_creation_visitor_result_type OpenCascadeKernel::convert_curve(const taxonomy::ptr curve) {
+    auto it = curve_cache_.find(curve->identity());
+    if (it != curve_cache_.end()) {
+		return it->second;
+    }
+	// hack hack
+    auto has_edge_geom = curve->kind() == taxonomy::EDGE && std::static_pointer_cast<taxonomy::edge>(curve)->basis;
+    if (has_edge_geom) {
+		auto edge_curve = std::static_pointer_cast<taxonomy::edge>(curve)->basis;
+        auto jt = edge_curve_cache_.find(edge_curve->identity());
+        if (jt != edge_curve_cache_.end()) {
+            auto& w = boost::get<TopoDS_Wire>(jt->second);
+            TopoDS_Iterator it(w);
+            auto& e = it.Value();
+            auto partner = e.Reversed();
+            BRep_Builder B;
+            TopoDS_Wire W;
+            B.MakeWire(W);
+            B.Add(W, partner);
+            return W;
+        }
+    }
 	curve_creation_visitor v{ this };
 	if (dispatch_curve_creation<curve_creation_visitor, 0>::dispatch(curve, v)) {
+        curve_cache_[curve->identity()] = v.result;
+        if (curve->kind() == taxonomy::EDGE && std::static_pointer_cast<taxonomy::edge>(curve)->basis) {
+            edge_curve_cache_[std::static_pointer_cast<taxonomy::edge>(curve)->basis->identity()] = v.result;
+        }
 		return v.result;
 	} else {
 		throw std::runtime_error("No curve created");
@@ -265,6 +337,19 @@ bool OpenCascadeKernel::convert(const taxonomy::loop::ptr loop, TopoDS_Wire& wir
 		converted_segments.Append(segment_wire);
 	}
 
+	{
+        for (auto& seg : converted_segments) {
+            TopoDS_Iterator it_seg(seg);
+            auto& e = TopoDS::Edge(it_seg.Value());
+            TopoDS_Vertex v0, v1;
+            TopExp::Vertices(e, v0, v1, true);
+            gp_Pnt p0 = BRep_Tool::Pnt(v0);
+            gp_Pnt p1 = BRep_Tool::Pnt(v1);
+            std::wcout << "or: " << e.Orientation() << " p0 " << p0.X() << " " << p0.Y() << " " << p0.Z() << " p1 " << p1.X() << " " << p1.Y() << " " << p1.Z() << std::endl;
+		}
+	}
+
+	/*
 	if (converted_segments.Extent() == 0) {
 		Logger::Message(Logger::LOG_ERROR, "No segment successfully converted:", loop->instance);
 		return false;
@@ -362,6 +447,71 @@ bool OpenCascadeKernel::convert(const taxonomy::loop::ptr loop, TopoDS_Wire& wir
 
 		wire = mw.Wire();
 	}
+	*/
+
+	{
+		// Rebuild to use partners in case of seam curves
+		// This swaps the edge, but fails to maintain partners in vertices
+		// Therefore see hack hack.
+		
+		std::map<Geom_Curve*, TopoDS_Edge> curve_edge_map;
+        
+        BRep_Builder BB;
+        TopoDS_Wire new_wire;
+        BB.MakeWire(new_wire);
+
+        for (auto& seg : converted_segments) {
+            TopoDS_Iterator it_seg(seg);
+            auto& edge = TopoDS::Edge(it_seg.Value());
+
+            double u0, u1;
+            if (auto crv = BRep_Tool::Curve(TopoDS::Edge(edge), u0, u1)) {
+                auto crvp = crv.get();
+                auto it = curve_edge_map.find(crvp);
+
+				TopoDS_Vertex v0, v1;
+                TopExp::Vertices(edge, v0, v1, true);
+                gp_Pnt p0 = BRep_Tool::Pnt(v0);
+                gp_Pnt p1 = BRep_Tool::Pnt(v1);
+
+				std::wcout << "or: " << edge.Orientation() << " p0 " << p0.X() << " " << p0.Y() << " " << p0.Z() << " p1 " << p1.X() << " " << p1.Y() << " " << p1.Z() << std::endl;
+
+                if (it == curve_edge_map.end()) {
+                    curve_edge_map[crvp] = edge;
+                    BB.Add(new_wire, edge);
+                } else {
+                    auto partner = TopoDS::Edge(it->second.Reversed());
+                    TopoDS_Vertex v0, v1;
+                    TopExp::Vertices(partner, v0, v1, true);
+                    gp_Pnt p0 = BRep_Tool::Pnt(v0);
+                    gp_Pnt p1 = BRep_Tool::Pnt(v1);
+
+                    std::wcout << ">>> " << partner.Orientation() << " p0 " << p0.X() << " " << p0.Y() << " " << p0.Z() << " p1 " << p1.X() << " " << p1.Y() << " " << p1.Z() << std::endl;
+
+
+                    BB.Add(new_wire, TopoDS::Edge(it->second.Reversed()));
+                }
+            }
+        }
+
+		new_wire.Closed(true);
+        ShapeFix_Shape sw(new_wire);
+        sw.Perform();
+        wire = TopoDS::Wire(sw.Shape());
+
+		/*
+		new_wire.Checked(wire.Checked());
+        new_wire.Closed(wire.Closed());
+        new_wire.Convex(wire.Convex());
+        new_wire.Free(wire.Free());
+        new_wire.Infinite(wire.Infinite());
+        new_wire.Locked(wire.Locked());
+        new_wire.Modified(wire.Modified());
+        new_wire.Orientable(wire.Orientable());
+		*/
+
+        wire = new_wire;
+	}
 
 	if (loop->matrix && !loop->matrix->is_identity()) {
 		const auto& m = loop->matrix->ccomponents();

src/ifcgeom/kernels/opencascade/solid.cpp

@@ -25,6 +25,7 @@
 #include <BRepAlgoAPI_Common.hxx>
 #include <ShapeFix_Solid.hxx>
 #include <BRepPrimAPI_MakeSphere.hxx>
+#include <BRepPrimAPI_MakeCylinder.hxx>
 
 using namespace ifcopenshell::geometry;
 using namespace ifcopenshell::geometry::kernels;
@@ -98,6 +99,19 @@ bool OpenCascadeKernel::convert(const taxonomy::solid::ptr solid, TopoDS_Shape&
 	if (!S.IsNull()) {
 		result = S;
 	}
+
+	{
+        std::ofstream f("C:\\temp\\my_solid.brep");
+        BRepTools::Dump(result, f);
+	}
+
+	{
+        auto test = BRepPrimAPI_MakeCylinder(1., 10.).Solid();
+        std::ofstream f("C:\\temp\\ref_solid.brep");
+        BRepTools::Dump(test, f);
+	}
+	
+
 	return !result.IsNull();
 }
 

src/ifcgeom/mapping/IfcRightCircularCylinder.cpp

@@ -1,4 +1,4 @@
-/********************************************************************************
+/********************************************************************************
  *                                                                              *
  * This file is part of IfcOpenShell.                                           *
  *                                                                              *
@@ -21,21 +21,119 @@
 #define mapping POSTFIX_SCHEMA(mapping)
 using namespace ifcopenshell::geometry;
 
-taxonomy::ptr mapping::map_impl(const IfcSchema::IfcRightCircularCylinder* inst) {
-	// @todo
-	return nullptr;
-	/*
-
-	const double r = inst->Radius() * length_unit_;
-	const double h = inst->Height() * length_unit_;
-
-	BRepPrimAPI_MakeCylinder builder(r, h);
-	gp_Trsf trsf;
-	IfcGeom::Kernel::convert(inst->Position(),trsf);
-	
-	// IfcCsgPrimitive3D.Position has unit scale factor
-	shape = builder.Solid().Moved(trsf);
-
-	return true;
-	*/
+static taxonomy::point3::ptr mk_point(const Eigen::Vector3d& p) {
+    return taxonomy::make<taxonomy::point3>(p);
+}
+
+static taxonomy::circle::ptr mk_circle(const Eigen::Vector3d& C,
+                                       const Eigen::Vector3d& Z,
+                                       const Eigen::Vector3d& X,
+                                       double R) {
+    auto c = taxonomy::make<taxonomy::circle>();
+    c->radius = R;
+    c->matrix = taxonomy::make<taxonomy::matrix4>(C, Z, X); // u=0 aligned with X
+    return c;
+}
+
+static taxonomy::line::ptr mk_line(const Eigen::Vector3d& A, const Eigen::Vector3d& B) {
+    auto l = taxonomy::make<taxonomy::line>();
+    l->matrix = taxonomy::make<taxonomy::matrix4>(A, (B - A)); // direction inferred in matrix ctor
+    return l;
+}
+
+taxonomy::ptr mapping::map_impl(const IfcSchema::IfcRightCircularCylinder* inst) {
+    const double pi2 = 2.0 * boost::math::constants::pi<double>();
+
+    auto O = Eigen::Vector3d::Zero();
+    auto X = Eigen::Vector3d::UnitX();
+    auto Y = Eigen::Vector3d::UnitY();
+    auto Z = Eigen::Vector3d::UnitZ();
+
+    auto R = inst->Radius() * length_unit_;
+    auto H = inst->Height() * length_unit_;
+
+    // Seam vertices at u=0 (and u=2pi, coincident) on bottom/top
+    const Eigen::Vector3d P0 = O + X * R;         // (u=0, v=0)
+    const Eigen::Vector3d P1 = O + Z * H + X * R; // (u=0, v=H)
+
+    auto cyl = taxonomy::make<taxonomy::cylinder>();
+    cyl->radius = R;
+    cyl->matrix = taxonomy::make<taxonomy::matrix4>(O, Z, X);
+
+    auto pl_bot = taxonomy::make<taxonomy::plane>();
+    pl_bot->matrix = taxonomy::make<taxonomy::matrix4>(O, -Z, X); // outward normal
+
+    auto pl_top = taxonomy::make<taxonomy::plane>();
+    pl_top->matrix = taxonomy::make<taxonomy::matrix4>(O + Z * H, Z, X); // outward normal
+
+    // Boundary 3D curves
+    auto c_bot = mk_circle(O, Z, X, R);
+    auto c_top = mk_circle(O + Z * H, Z, X, R);
+    auto seam = mk_line(P0, P1); // coincident seam geometry for u=0 and u=2π
+
+    auto v00 = mk_point(P0);
+    auto v01 = mk_point(P1);
+
+    // Edges
+    auto e_bot = taxonomy::make<taxonomy::edge>();
+    e_bot->basis = c_bot;
+    e_bot->start = v00; // closed circle: start=end = seam vertex
+    e_bot->end = v00;
+    e_bot->curve_sense = true; // forward umin->umax
+
+    auto e_top = taxonomy::make<taxonomy::edge>();
+    e_top->basis = c_top;
+    e_top->start = v01;
+    e_top->end = v01;
+    e_top->curve_sense = false; // reversed umax->umin
+
+    auto e_seam_u0 = taxonomy::make<taxonomy::edge>();
+    e_seam_u0->basis = seam;
+    e_seam_u0->start = v00;
+    e_seam_u0->end = v01;
+
+    // Duplicate seam edge for periodic closure (u=2π). Same 3D geometry, distinct topological edge.
+    auto e_seam_u2pi = taxonomy::make<taxonomy::edge>();
+    e_seam_u2pi->basis = seam;
+    e_seam_u2pi->start = v01;
+    e_seam_u2pi->end = v00;
+
+    // Lateral loop (order matches the conceptual rectangle in (u,v))
+    auto loop_lat = taxonomy::make<taxonomy::loop>();
+    loop_lat->external = true;
+    loop_lat->closed = true;
+    loop_lat->children = {e_bot, e_seam_u0, e_top, e_seam_u2pi};
+
+    auto face_lat = taxonomy::make<taxonomy::face>();
+    face_lat->basis = cyl;
+    face_lat->children = {loop_lat};
+
+    // Bottom cap
+    auto loop_bot = taxonomy::make<taxonomy::loop>();
+    loop_bot->external = true;
+    loop_bot->closed = true;
+    loop_bot->children = {e_bot};
+
+    auto face_bot = taxonomy::make<taxonomy::face>();
+    face_bot->basis = pl_bot;
+    face_bot->children = {loop_bot};
+
+    // Top cap
+    auto loop_top = taxonomy::make<taxonomy::loop>();
+    loop_top->external = true;
+    loop_top->closed = true;
+    loop_top->children = {e_top};
+
+    auto face_top = taxonomy::make<taxonomy::face>();
+    face_top->basis = pl_top;
+    face_top->children = {loop_top};
+
+    // Shell + solid
+    auto sh = taxonomy::make<taxonomy::shell>();
+    sh->closed = true;
+    sh->children = {face_lat, face_bot, face_top};
+
+    auto so = taxonomy::make<taxonomy::solid>();
+    so->children = {sh};
+    return so;
 }