diff --git a/src/ifcgeom/kernel_agnostic/AbstractKernel.h b/src/ifcgeom/kernel_agnostic/AbstractKernel.h index 3b3021b42f..45c3ed5bf5 100644 --- a/src/ifcgeom/kernel_agnostic/AbstractKernel.h +++ b/src/ifcgeom/kernel_agnostic/AbstractKernel.h @@ -45,7 +45,7 @@ namespace ifcopenshell { namespace geometry { namespace kernels { virtual bool convert_impl(const taxonomy::line*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::circle*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::ellipse*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } - virtual bool convert_impl(const taxonomy::bspline*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } + virtual bool convert_impl(const taxonomy::bspline_curve*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::edge*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::loop*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::shell*, ifcopenshell::geometry::ConversionResults&) { throw std::runtime_error("Not implemented"); } diff --git a/src/ifcgeom/kernels/opencascade/IfcGeomShapes.cpp b/src/ifcgeom/kernels/opencascade/IfcGeomShapes.cpp index 672789fc41..71c382b5c2 100644 --- a/src/ifcgeom/kernels/opencascade/IfcGeomShapes.cpp +++ b/src/ifcgeom/kernels/opencascade/IfcGeomShapes.cpp @@ -170,6 +170,646 @@ bool OpenCascadeKernel::convert(const taxonomy::extrusion* extrusion, TopoDS_Sha return !shape.IsNull(); } +namespace { + /* Returns whether wire conforms to a polyhedron, i.e. only edges with linear curves*/ + bool is_polyhedron(const TopoDS_Wire& wire) { + double a, b; + TopLoc_Location l; + + TopoDS_Iterator it(wire, false, false); + for (; it.More(); it.Next()) { + auto crv = BRep_Tool::Curve(TopoDS::Edge(it.Value()), l, a, b); + if (!crv || crv->DynamicType() != STANDARD_TYPE(Geom_Line)) { + return false; + } + } + + return true; + } + + /* Returns whether wire conforms to a polyhedron, i.e. only edges with linear curves*/ + bool is_polyhedron(const taxonomy::loop* wire) { + for (auto& edge : wire->children_as()) { + if (edge->basis) { + if (edge->basis->kind() != taxonomy::LINE) { + return false; + } + } + } + return true; + } + + /* A temporary structure to store the intermediate data for the face conversion */ + class face_definition { + private: + Handle(Geom_Surface) surface_; + std::vector wires_; + bool all_outer_; + public: + face_definition() : surface_(), all_outer_(false) {} + + typedef std::vector::const_iterator wire_it; + + bool& all_outer() { + return all_outer_; + } + + bool all_outer() const { + return all_outer_; + } + + Handle(Geom_Surface)& surface() { + return surface_; + } + + const Handle(Geom_Surface)& surface() const { + return surface_; + } + + std::vector& wires() { + return wires_; + } + + const TopoDS_Wire& outer_wire() const { + return wires_.front(); + } + + std::pair inner_wires() const { + return { wires_.begin() + 1, wires_.end() }; + } + }; +} + +#include +#include +#include +#include + +bool OpenCascadeKernel::convert(const taxonomy::face* face, TopoDS_Shape& result) { + std::vector bounds; + std::transform(face->children.begin(), face->children.end(), std::back_inserter(bounds), [](auto item){ + return static_cast(item); + }); + + face_definition fd; + + const bool is_face_surface = false; /* todo */ + + /* + if (is_face_surface) { + IfcSchema::IfcFaceSurface* fs = (IfcSchema::IfcFaceSurface*) l; + fs->FaceSurface(); + // FIXME: Surfaces are interpreted as a TopoDS_Shape + TopoDS_Shape surface_shape; + if (!convert_shape(fs->FaceSurface(), surface_shape)) return false; + + // FIXME: Assert this obtaines the only face + TopExp_Explorer exp(surface_shape, TopAbs_FACE); + if (!exp.More()) return false; + + TopoDS_Face surface = TopoDS::Face(exp.Current()); + fd.surface() = BRep_Tool::Surface(surface); + } + */ + + const int num_bounds = bounds.size(); + int num_outer_bounds = 0; + + for (auto& bound: bounds) { + if (bound->external.get_value_or(false)) { + num_outer_bounds++; + } + } + + // The number of outer bounds should be one according to the schema. Also Open Cascade + // expects this, but it is not strictly checked. Regardless, if the number is greater, + // the face will still be processed as long as there are no holes. A compound of faces + // is returned in that case. + if (num_bounds > 1 && num_outer_bounds > 1 && num_bounds != num_outer_bounds) { + Logger::Message(Logger::LOG_ERROR, "Invalid configuration of boundaries for:", face->instance); + return false; + } + + if (num_outer_bounds > 1) { + Logger::Message(Logger::LOG_WARNING, "Multiple outer boundaries for:", face->instance); + fd.all_outer() = true; + } + + TopTools_DataMapOfShapeInteger wire_senses; + + for (int process_interior = 0; process_interior <= 1; ++process_interior) { + for (auto& bound : bounds) { + bool same_sense = true; /* todo bound->Orientation(); */ + + const bool is_interior = + !bound->external.get_value_or(false) && + (num_bounds > 1) && + (num_outer_bounds < num_bounds); + + // The exterior face boundary is processed first + if (is_interior == !process_interior) continue; + + TopoDS_Wire wire; + if (faceset_helper_ && is_polyhedron(bound)) { + if (!faceset_helper_->wire(bound, wire)) { + Logger::Message(Logger::LOG_WARNING, "Face boundary loop not included", bound->instance); + continue; + } + } else if (!convert(bound, wire)) { + Logger::Message(Logger::LOG_ERROR, "Failed to process face boundary loop", bound->instance); + return false; + } + + if (!same_sense) { + wire.Reverse(); + } + + wire_senses.Bind(wire.Oriented(TopAbs_FORWARD), same_sense ? TopAbs_FORWARD : TopAbs_REVERSED); + + fd.wires().emplace_back(wire); + } + } + + if (fd.wires().empty()) { + Logger::Warning("Face with no boundaries", face->instance); + return false; + } + + if (fd.surface().IsNull()) { + // Use the first wire to find a plane manually for polygonal wires + const TopoDS_Wire& wire = fd.wires().front(); + if (is_polyhedron(wire)) { + TopExp_Explorer exp(wire, TopAbs_EDGE); + int count = 0; + TopoDS_Edge edges[2]; + for (; exp.More(); exp.Next(), count++) { + if (count < 2) { + edges[count] = TopoDS::Edge(exp.Current()); + } + } + + if (count == 3) { + // Help Open Cascade by finding the plane more efficiently + double _, __; + Handle(Geom_Line) c1 = Handle(Geom_Line)::DownCast(BRep_Tool::Curve(edges[0], _, __)); + Handle(Geom_Line) c2 = Handle(Geom_Line)::DownCast(BRep_Tool::Curve(edges[1], _, __)); + + const gp_Vec ab = c1->Position().Direction(); + const gp_Vec ac = c2->Position().Direction(); + const gp_Vec cross = ab.Crossed(ac); + + if (cross.SquareMagnitude() > ALMOST_ZERO) { + const gp_Dir n = cross; + fd.surface() = new Geom_Plane(c1->Position().Location(), n); + } + } else { + gp_Pln pln; + if (approximate_plane_through_wire(wire, pln)) { + fd.surface() = new Geom_Plane(pln); + } + } + } + } + + if (fd.surface().IsNull()) { + // BRepLib_FindSurface is used in case no surface is found or provided + + const TopoDS_Wire& wire = fd.wires().front(); + + BRepLib_FindSurface fs(wire, precision_, true, true); + if (fs.Found()) { + fd.surface() = fs.Surface(); + ShapeFix_ShapeTolerance ftol; + ftol.SetTolerance(wire, fs.ToleranceReached(), TopAbs_WIRE); + } + } + + TopTools_ListOfShape face_list; + + if (fd.surface().IsNull()) { + // The set of wires is triangulated in case no surface can be found + Logger::Message(Logger::LOG_WARNING, "Triangulating face boundaries for face", face->instance); + + if (fd.all_outer()) { + for (const auto& w : fd.wires()) { + TopTools_ListOfShape fl; + triangulate_wire({ w }, fl); + face_list.Append(fl); + } + } else { + triangulate_wire(fd.wires(), face_list); + } + } else if (!fd.all_outer()) { + BRepBuilderAPI_MakeFace mf(fd.surface(), fd.outer_wire()); + + if (mf.IsDone()) { + // Is this necessary + TopoDS_Face f = mf.Face(); + mf.Init(f); + + for (auto it = fd.inner_wires().first; it != fd.inner_wires().second; ++it) { + mf.Add(*it); + } + + face_list.Append(mf.Face()); + } + } else { + for (const auto& w : fd.wires()) { + BRepBuilderAPI_MakeFace mf(fd.surface(), w); + if (mf.IsDone()) { + face_list.Append(mf.Face()); + } + } + } + + if (!fd.surface().IsNull()) { + // Some fixes for orientation and p-curves. If we have no surface, it + // means the face has been triangulated in which case none of these + // fixes are necessary. + + if (fd.surface()->DynamicType() != STANDARD_TYPE(Geom_Plane)) { + // In case of (non-planar) face surface, p-curves need to be computed. + // For planar faces, Open Cascade generates p-curves on the fly. + + for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) { + // Small chance there are multiple faces + const TopoDS_Face& occ_face = TopoDS::Face(it.Value()); + for (TopExp_Explorer exp2(occ_face, TopAbs_EDGE); exp2.More(); exp2.Next()) { + const TopoDS_Edge& edge = TopoDS::Edge(exp2.Current()); + ShapeFix_Edge fix_edge; + fix_edge.FixAddPCurve(edge, occ_face, false, precision_); + } + } + } + + for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) { + const TopoDS_Face& occ_face = TopoDS::Face(it.Value()); + + ShapeFix_Face sfs(TopoDS::Face(occ_face)); + TopTools_DataMapOfShapeListOfShape wire_map; + sfs.FixOrientation(wire_map); + + TopoDS_Iterator jt(occ_face, false); + for (; jt.More(); jt.Next()) { + const TopoDS_Wire& w = TopoDS::Wire(jt.Value()); + // tfk: @todo if wire_map contains w, I would assume wire_senses also contains w, + // this is not the case in github issue #405. + if (wire_map.IsBound(w) && wire_senses.IsBound(w)) { + const TopTools_ListOfShape& shapes = wire_map.Find(w); + TopTools_ListIteratorOfListOfShape kt(shapes); + for (; kt.More(); kt.Next()) { + // Apparently the wire got reversed, so register it with opposite orientation in the map + wire_senses.Bind(kt.Value(), wire_senses.Find(w) == TopAbs_FORWARD ? TopAbs_REVERSED : TopAbs_FORWARD); + } + } + } + + it.Value() = sfs.Face(); + } + + for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) { + TopoDS_Face& occ_face = TopoDS::Face(it.Value()); + + bool all_reversed = true; + TopoDS_Iterator jt(occ_face, false); + for (; jt.More(); jt.Next()) { + const TopoDS_Wire& w = TopoDS::Wire(jt.Value()); + if (!wire_senses.IsBound(w.Oriented(TopAbs_FORWARD)) || (w.Orientation() == wire_senses.Find(w.Oriented(TopAbs_FORWARD)))) { + all_reversed = false; + } + } + + if (all_reversed) { + occ_face.Reverse(); + } + } + } + + if (face_list.Extent() > 1) { + TopoDS_Compound compound; + BRep_Builder builder; + builder.MakeCompound(compound); + for (TopTools_ListIteratorOfListOfShape it(face_list); it.More(); it.Next()) { + TopoDS_Face& occ_face = TopoDS::Face(it.Value()); + builder.Add(compound, occ_face); + } + result = compound; + } else { + result = face_list.First(); + } + + return true; +} + +#include +#include + +namespace { + /* A compile-time for loop over the curve kinds */ + template + struct dispatch_curve_creation { + static bool dispatch(const ifcopenshell::geometry::taxonomy::item* item, T visitor) { + // @todo it should be possible to eliminate this dynamic_cast when there is a static equivalent to kind() + const ifcopenshell::geometry::taxonomy::curves::type* v = dynamic_cast*>(item); + if (v) { + visitor(*v); + return true; + } else { + return dispatch_curve_creation::dispatch(item, visitor); + } + } + }; + + template + struct dispatch_curve_creation { + static bool dispatch(const ifcopenshell::geometry::taxonomy::item* item, T visitor) { + Logger::Error("No conversion for " + std::to_string(item->kind())); + return false; + } + }; + + template + T convert_xyz(const U& u) { + const double* vs = u.components.data(); + return T(vs[0], vs[1], vs[2]); + } + + struct curve_creation_visitor { + typedef boost::variant result_type; + result_type result; + + result_type operator()(const taxonomy::bspline_curve&) { + throw std::runtime_error("Not implemented"); + } + + result_type operator()(const taxonomy::line& l) { + return result = Handle(Geom_Curve)(new Geom_Line(convert_xyz(l.origin), convert_xyz(l.direction))); + } + + result_type operator()(const taxonomy::circle& c) { + return result = Handle(Geom_Curve)(new Geom_Circle(gp_Ax2(convert_xyz(c.origin), convert_xyz(c.z), convert_xyz(c.x)), c.radius)); + } + + result_type operator()(const taxonomy::ellipse& e) { + return result = Handle(Geom_Curve)(new Geom_Ellipse(gp_Ax2(convert_xyz(e.origin), convert_xyz(e.z), convert_xyz(e.x)), e.radius, e.radius2)); + } + }; + + curve_creation_visitor::result_type convert_curve(const taxonomy::item* curve) { + curve_creation_visitor v; + if (dispatch_curve_creation::dispatch(curve, v)) { + return v.result; + } else { + throw std::runtime_error("No curve created"); + } + } +} + +#include +#include + +namespace { + // Returns the other vertex of an edge + TopoDS_Vertex other(const TopoDS_Edge& e, const TopoDS_Vertex& v) { + TopoDS_Vertex a, b; + TopExp::Vertices(e, a, b); + return v.IsSame(b) ? a : b; + } + + TopoDS_Edge first_edge(const TopoDS_Wire& w) { + TopoDS_Vertex v1, v2; + TopExp::Vertices(w, v1, v2); + TopTools_IndexedDataMapOfShapeListOfShape wm; + TopExp::MapShapesAndAncestors(w, TopAbs_VERTEX, TopAbs_EDGE, wm); + return TopoDS::Edge(wm.FindFromKey(v1).First()); + } + + // Returns new wire with the edge replaced by a linear edge with the vertex v moved to p + TopoDS_Wire adjust(const TopoDS_Wire& w, const TopoDS_Vertex& v, const gp_Pnt& p) { + TopTools_IndexedDataMapOfShapeListOfShape map; + TopExp::MapShapesAndAncestors(w, TopAbs_VERTEX, TopAbs_EDGE, map); + + bool all_linear = true, single_circle = false, first = true; + + const TopTools_ListOfShape& edges = map.FindFromKey(v); + TopTools_ListIteratorOfListOfShape it(edges); + for (; it.More(); it.Next()) { + const TopoDS_Edge& e = TopoDS::Edge(it.Value()); + double _, __; + Handle(Geom_Curve) crv = BRep_Tool::Curve(e, _, __); + const bool is_line = crv->DynamicType() == STANDARD_TYPE(Geom_Line); + const bool is_circle = crv->DynamicType() == STANDARD_TYPE(Geom_Circle); + all_linear = all_linear && is_line; + single_circle = first && is_circle; + } + + if (all_linear) { + BRep_Builder b; + TopoDS_Vertex v2; + b.MakeVertex(v2, p, BRep_Tool::Tolerance(v)); + + ShapeBuild_ReShape reshape; + reshape.Replace(v.Oriented(TopAbs_FORWARD), v2); + + return TopoDS::Wire(reshape.Apply(w)); + } else if (single_circle) { + TopoDS_Vertex v1, v2; + TopExp::Vertices(w, v1, v2); + + gp_Pnt p1, p2, p3; + p1 = v.IsEqual(v1) ? p : BRep_Tool::Pnt(v1); + p3 = v.IsEqual(v2) ? p : BRep_Tool::Pnt(v2); + + double a, b; + Handle(Geom_Curve) crv = BRep_Tool::Curve(TopoDS::Edge(edges.First()), a, b); + crv->D0((a + b) / 2., p2); + + GC_MakeCircle mc(p1, p2, p3); + if (!mc.IsDone()) { + throw std::runtime_error("Failed to adjust circle"); + } + + TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(mc.Value(), p1, p3).Edge(); + BRepBuilderAPI_MakeWire builder; + builder.Add(edge); + return builder.Wire(); + } else { + throw std::runtime_error("Unexpected wire to adjust"); + } + } + + // A wrapper around BRepBuilderAPI_MakeWire that makes sure segments are connected either by moving end points or by adding intermediate segments + class wire_builder { + private: + BRepBuilderAPI_MakeWire mw_; + double p_; + bool override_next_; + gp_Pnt next_override_; + const IfcUtil::IfcBaseClass* inst_; + + public: + wire_builder(double p, const IfcUtil::IfcBaseClass* inst = 0) : p_(p), override_next_(false), inst_(inst) {} + + void operator()(const TopoDS_Shape& a) { + const TopoDS_Wire& w = TopoDS::Wire(a); + if (override_next_) { + override_next_ = false; + TopoDS_Edge e = first_edge(w); + mw_.Add(adjust(w, TopExp::FirstVertex(e, true), next_override_)); + } else { + mw_.Add(w); + } + } + + void operator()(const TopoDS_Shape& a, const TopoDS_Shape& b, bool last) { + TopoDS_Wire w1 = TopoDS::Wire(a); + const TopoDS_Wire& w2 = TopoDS::Wire(b); + + if (override_next_) { + override_next_ = false; + TopoDS_Edge e = first_edge(w1); + w1 = adjust(w1, TopExp::FirstVertex(e, true), next_override_); + } + + TopoDS_Vertex w11, w12, w21, w22; + TopExp::Vertices(w1, w11, w12); + TopExp::Vertices(w2, w21, w22); + + gp_Pnt p1 = BRep_Tool::Pnt(w12); + gp_Pnt p2 = BRep_Tool::Pnt(w21); + + double dist = p1.Distance(p2); + + // Distance is within tolerance, this is fine + if (dist < p_) { + mw_.Add(w1); + goto check; + } + + // Distance is too large for attempting to move end points, add intermediate edge + if (dist > 1000. * p_) { + mw_.Add(w1); + mw_.Add(BRepBuilderAPI_MakeEdge(p1, p2)); + Logger::Warning("Added additional segment to close gap with length " + boost::lexical_cast(dist) + " to:", inst_); + goto check; + } + + { + TopTools_IndexedDataMapOfShapeListOfShape wmap1, wmap2; + + // Find edges connected to end- and begin vertex + TopExp::MapShapesAndAncestors(w1, TopAbs_VERTEX, TopAbs_EDGE, wmap1); + TopExp::MapShapesAndAncestors(w2, TopAbs_VERTEX, TopAbs_EDGE, wmap2); + + const TopTools_ListOfShape& last_edges = wmap1.FindFromKey(w12); + const TopTools_ListOfShape& first_edges = wmap2.FindFromKey(w21); + + double _, __; + if (last_edges.Extent() == 1 && first_edges.Extent() == 1) { + Handle(Geom_Curve) c1 = BRep_Tool::Curve(TopoDS::Edge(last_edges.First()), _, __); + Handle(Geom_Curve) c2 = BRep_Tool::Curve(TopoDS::Edge(first_edges.First()), _, __); + + const bool is_line1 = c1->DynamicType() == STANDARD_TYPE(Geom_Line); + const bool is_line2 = c2->DynamicType() == STANDARD_TYPE(Geom_Line); + + const bool is_circle1 = c1->DynamicType() == STANDARD_TYPE(Geom_Circle); + const bool is_circle2 = c2->DynamicType() == STANDARD_TYPE(Geom_Circle); + + // Preferably adjust the segment that is linear + if (is_line1 || (is_circle1 && !is_line2)) { + mw_.Add(adjust(w1, w12, p2)); + Logger::Notice("Adjusted edge end-point with distance " + boost::lexical_cast(dist) + " on:", inst_); + } else if ((is_line2 || is_circle2) && !last) { + mw_.Add(w1); + override_next_ = true; + next_override_ = p1; + Logger::Notice("Adjusted edge end-point with distance " + boost::lexical_cast(dist) + " on:", inst_); + } else { + // In all other cases an edge is added + mw_.Add(w1); + mw_.Add(BRepBuilderAPI_MakeEdge(p1, p2)); + Logger::Warning("Added additional segment to close gap with length " + boost::lexical_cast(dist) + " to:", inst_); + } + } else { + Logger::Error("Internal error, inconsistent wire segments", inst_); + mw_.Add(w1); + } + } + + check: + if (mw_.Error() == BRepBuilderAPI_NonManifoldWire) { + Logger::Error("Non-manifold curve segments:", inst_); + } else if (mw_.Error() == BRepBuilderAPI_DisconnectedWire) { + Logger::Error("Failed to join curve segments:", inst_); + } + } + + const TopoDS_Wire& wire() { return mw_.Wire(); } + }; + + template + void shape_pair_enumerate(TopTools_ListIteratorOfListOfShape& it, Fn& fn, bool closed) { + bool is_first = true; + TopoDS_Shape first, previous, current; + for (; it.More(); it.Next(), is_first = false) { + current = it.Value(); + if (is_first) { + first = current; + } else { + fn(previous, current, false); + } + previous = current; + } + if (closed) { + fn(current, first, true); + } else { + fn(current); + } + } +} + +bool OpenCascadeKernel::convert(const taxonomy::loop* loop, TopoDS_Wire& wire) { + auto segments = loop->children_as(); + + TopTools_ListOfShape converted_segments; + + for (auto& segment : segments) { + TopoDS_Wire segment_wire = boost::get(convert_curve(segment)); + + if (!segment->orientation) { + segment_wire.Reverse(); + } + + ShapeFix_ShapeTolerance FTol; + FTol.SetTolerance(segment_wire, precision_, TopAbs_WIRE); + + converted_segments.Append(segment_wire); + } + + if (converted_segments.Extent() == 0) { + Logger::Message(Logger::LOG_ERROR, "No segment succesfully converted:", loop->instance); + return false; + } + + BRepBuilderAPI_MakeWire w; + TopoDS_Vertex wire_first_vertex, wire_last_vertex, edge_first_vertex, edge_last_vertex; + + TopTools_ListIteratorOfListOfShape it(converted_segments); + + /* + @todo + IfcEntityList::ptr profile = l->data().getInverse(&IfcSchema::IfcProfileDef::Class(), -1); + const bool force_close = profile && profile->size() > 0; + */ + const bool force_close = false; + + wire_builder bld(precision_, loop->instance); + shape_pair_enumerate(it, bld, force_close); + wire = bld.wire(); + + return true; +} + bool OpenCascadeKernel::convert_impl(const taxonomy::extrusion* extrusion, ifcopenshell::geometry::ConversionResults& results) { TopoDS_Shape shape; if (!convert(extrusion, shape)) { @@ -192,4 +832,276 @@ bool OpenCascadeKernel::convert(const taxonomy::matrix4* matrix, gp_GTrsf& trsf) } } return true; -} \ No newline at end of file +} + +#include + +bool OpenCascadeKernel::approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps) { + // Newell's Method is used for the normal calculation + // as a simple edge cross product can give opposite results + // for a concave face boundary. + // Reference: Graphics Gems III p. 231 + + const double eps_ = eps < 1. ? precision_ : eps; + const double eps2 = eps_ * eps_; + + double x = 0, y = 0, z = 0; + gp_Pnt current, previous, first; + gp_XYZ center; + int n = 0; + + BRepTools_WireExplorer exp(wire); + + for (;; exp.Next()) { + const bool has_more = exp.More() != 0; + if (has_more) { + const TopoDS_Vertex& v = exp.CurrentVertex(); + current = BRep_Tool::Pnt(v); + center += current.XYZ(); + } else { + current = first; + } + if (n) { + const double& xn = previous.X(); + const double& yn = previous.Y(); + const double& zn = previous.Z(); + const double& xn1 = current.X(); + const double& yn1 = current.Y(); + const double& zn1 = current.Z(); + x += (yn - yn1)*(zn + zn1); + y += (xn + xn1)*(zn - zn1); + z += (xn - xn1)*(yn + yn1); + } else { + first = current; + } + if (!has_more) { + break; + } + previous = current; + ++n; + } + + if (n < 3) { + return false; + } + + plane = gp_Pln(center / n, gp_Dir(x, y, z)); + + exp.Init(wire); + for (; exp.More(); exp.Next()) { + const TopoDS_Vertex& v = exp.CurrentVertex(); + current = BRep_Tool::Pnt(v); + if (plane.SquareDistance(current) > eps2) { + return false; + } + } + + return true; +} + + +bool OpenCascadeKernel::triangulate_wire(const std::vector& wires, TopTools_ListOfShape& faces) { + // This is a bit of a precarious approach, but seems to work for the + // versions of OCCT tested for. OCCT has a Delaunay triangulation function + // BRepMesh_Delaun, but it is notoriously hard to interpret the results + // (due to the Bowyer-Watson super triangle perhaps?). Therefore + // alternatively we use the regular OCCT incremental mesher on a new face + // created from the UV coordinates of the original wire. Pray to our gods + // that the vertex coordinates are unaffected by the meshing algorithm and + // map them back to 3d coordinates when iterating over the mesh triangles. + + // In addition, to maintain a manifold shell, we need to make sure that + // every edge from the input wire is used exactly once in the list of + // resulting faces. And that other internal edges are used twice. + + typedef std::pair uv_node; + + gp_Pln pln; + if (!approximate_plane_through_wire(wires.front(), pln, std::numeric_limits::infinity())) { + return false; + } + + const gp_XYZ& udir = pln.Position().XDirection().XYZ(); + const gp_XYZ& vdir = pln.Position().YDirection().XYZ(); + const gp_XYZ& pnt = pln.Position().Location().XYZ(); + + std::map mapping; + std::map, TopoDS_Edge> existing_edges, new_edges; + + std::unique_ptr mf; + + for (auto it = wires.begin(); it != wires.end(); ++it) { + const TopoDS_Wire& wire = *it; + BRepTools_WireExplorer exp(wire); + BRepBuilderAPI_MakePolygon mp; + + // Add UV coordinates to a newly created polygon + for (; exp.More(); exp.Next()) { + // Project onto plane + const TopoDS_Vertex& V = exp.CurrentVertex(); + gp_Pnt p = BRep_Tool::Pnt(V); + double u = (p.XYZ() - pnt).Dot(udir); + double v = (p.XYZ() - pnt).Dot(vdir); + mp.Add(gp_Pnt(u, v, 0.)); + + mapping.insert(std::make_pair(std::make_pair(u, v), V)); + + // Store existing edges in a map so that triangles can + // actually reference the preexisting edges. + const TopoDS_Edge& e = exp.Current(); + TopoDS_Vertex V0, V1; + TopExp::Vertices(e, V0, V1, true); + gp_Pnt p0 = BRep_Tool::Pnt(V0); + gp_Pnt p1 = BRep_Tool::Pnt(V1); + double u0 = (p0.XYZ() - pnt).Dot(udir); + double v0 = (p0.XYZ() - pnt).Dot(vdir); + double u1 = (p1.XYZ() - pnt).Dot(udir); + double v1 = (p1.XYZ() - pnt).Dot(vdir); + uv_node uv0 = std::make_pair(u0, v0); + uv_node uv1 = std::make_pair(u1, v1); + existing_edges.insert(std::make_pair(std::make_pair(uv0, uv1), e)); + existing_edges.insert(std::make_pair(std::make_pair(uv1, uv0), TopoDS::Edge(e.Reversed()))); + } + + // Not closed by default + mp.Close(); + + if (mf) { + if (it - 1 == wires.begin()) { + // @todo is this necessary? + TopoDS_Face f = mf->Face(); + mf->Init(f); + } + mf->Add(mp.Wire()); + } else { + mf.reset(new BRepBuilderAPI_MakeFace(mp.Wire())); + } + } + + const TopoDS_Face& face = mf->Face(); + + // Create a triangular mesh from the face + BRepMesh_IncrementalMesh(face, Precision::Confusion()); + + int n123[3]; + TopLoc_Location loc; + Handle_Poly_Triangulation tri = BRep_Tool::Triangulation(face, loc); + + if (!tri.IsNull()) { + const TColgp_Array1OfPnt& nodes = tri->Nodes(); + + const Poly_Array1OfTriangle& triangles = tri->Triangles(); + for (int i = 1; i <= triangles.Length(); ++i) { + if (face.Orientation() == TopAbs_REVERSED) + triangles(i).Get(n123[2], n123[1], n123[0]); + else triangles(i).Get(n123[0], n123[1], n123[2]); + + // Create polygons from the mesh vertices + BRepBuilderAPI_MakeWire mp2; + for (int j = 0; j < 3; ++j) { + + uv_node uvnodes[2]; + TopoDS_Vertex vs[2]; + + for (int k = 0; k < 2; ++k) { + const gp_Pnt& uv = nodes.Value(n123[(j + k) % 3]); + uvnodes[k] = std::make_pair(uv.X(), uv.Y()); + + auto it = mapping.find(uvnodes[k]); + if (it == mapping.end()) { + Logger::Error("Internal error: unable to unproject uv-mesh"); + return false; + } + + vs[k] = it->second; + } + + auto it = existing_edges.find(std::make_pair(uvnodes[0], uvnodes[1])); + if (it != existing_edges.end()) { + // This is a boundary edge, reuse existing edge from wire + mp2.Add(it->second); + } else { + auto jt = new_edges.find(std::make_pair(uvnodes[0], uvnodes[1])); + if (jt != new_edges.end()) { + // We have already added the reverse as part of another + // triangle, reuse this edge. + mp2.Add(TopoDS::Edge(jt->second)); + } else { + // This is a new internal edge. Register the reverse + // for reuse later. We need to be sure to reuse vertices + // for the edge construction because otherwise the wire + // builder will use geometrical proximity for vertex + // connections in which case the edge will be copied + // and no longer partner with other edges from the shell. + TopoDS_Edge ne = BRepBuilderAPI_MakeEdge(vs[0], vs[1]); + mp2.Add(ne); + // Store the reverse to be picked up later. + new_edges.insert(std::make_pair(std::make_pair(uvnodes[1], uvnodes[0]), TopoDS::Edge(ne.Reversed()))); + } + } + } + + BRepBuilderAPI_MakeFace mft(mp2.Wire()); + if (mft.IsDone()) { + TopoDS_Face triangle_face = mft.Face(); + TopoDS_Iterator jt(triangle_face, false); + for (; jt.More(); jt.Next()) { + const TopoDS_Wire& w = TopoDS::Wire(jt.Value()); + if (w.Orientation() != wires.front().Orientation()) { + triangle_face.Reverse(); + } + } + faces.Append(triangle_face); + } else { + Logger::Error("Internal error: missing face"); + return false; + } + } + } + + TopTools_IndexedDataMapOfShapeListOfShape mape, mapn; + for (auto& wire : wires) { + TopExp::MapShapesAndAncestors(wire, TopAbs_EDGE, TopAbs_WIRE, mape); + } + TopTools_ListIteratorOfListOfShape it(faces); + for (; it.More(); it.Next()) { + TopExp::MapShapesAndAncestors(it.Value(), TopAbs_EDGE, TopAbs_WIRE, mapn); + } + + // Validation + + for (int i = 1; i <= mape.Extent(); ++i) { +#if OCC_VERSION_HEX >= 0x70000 + TopTools_ListOfShape val; + if (!mapn.FindFromKey(mape.FindKey(i), val)) { +#else + bool contains = false; + try { + TopTools_ListOfShape val = mapn.FindFromKey(mape.FindKey(i)); + contains = true; + } catch (Standard_NoSuchObject&) {} + if (!contains) { +#endif + // All existing edges need to exist in the new faces + Logger::Error("Internal error, missing edge from triangulation"); + if (faceset_helper_ != nullptr) { + faceset_helper_->non_manifold() = true; + } + } + } + + for (int i = 1; i <= mapn.Extent(); ++i) { + const TopoDS_Shape& v = mapn.FindKey(i); + int n = mapn.FindFromIndex(i).Extent(); + // Existing edges are boundaries with use 1 + // New edges are internal with use 2 + if (n != (mape.Contains(v) ? 1 : 2)) { + Logger::Error("Internal error, non-manifold result from triangulation"); + if (faceset_helper_ != nullptr) { + faceset_helper_->non_manifold() = true; + } + } + } + + return true; + } \ No newline at end of file diff --git a/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.h b/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.h index ffce0fa956..101c6643c2 100644 --- a/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.h +++ b/src/ifcgeom/kernels/opencascade/OpenCascadeKernel.h @@ -235,8 +235,12 @@ namespace kernels { bool convert(const taxonomy::extrusion*, TopoDS_Shape&); bool convert(const taxonomy::face*, TopoDS_Shape&); + bool convert(const taxonomy::loop*, TopoDS_Wire&); bool convert(const taxonomy::matrix4*, gp_GTrsf&); + bool approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps = -1.); + bool triangulate_wire(const std::vector& wires, TopTools_ListOfShape& faces); + virtual bool convert_impl(const taxonomy::shell*, ifcopenshell::geometry::ConversionResults&); virtual bool convert_impl(const taxonomy::extrusion*, ifcopenshell::geometry::ConversionResults&); }; diff --git a/src/ifcgeom/taxonomy.h b/src/ifcgeom/taxonomy.h index 122f5a0fef..3be8ec2f6d 100644 --- a/src/ifcgeom/taxonomy.h +++ b/src/ifcgeom/taxonomy.h @@ -24,7 +24,7 @@ public: topology_error() : std::runtime_error("Generic topology error") {} }; -enum kinds { MATRIX4, POINT3, DIRECTION3, LINE, CIRCLE, ELLIPSE, BSPLINE, EDGE, LOOP, FACE, SHELL, EXTRUSION, NODE, COLLECTION, COLOUR, STYLE }; +enum kinds { MATRIX4, POINT3, DIRECTION3, LINE, CIRCLE, ELLIPSE, BSPLINE_CURVE, EDGE, LOOP, FACE, SHELL, EXTRUSION, NODE, COLLECTION, COLOUR, STYLE }; struct item { const IfcUtil::IfcBaseClass* instance; @@ -114,46 +114,74 @@ struct direction3 : public cartesian_base<3> { direction3(double x = 0., double y = 0., double z = 0.) : cartesian_base(x, y, z) {} }; -struct line : public geom_item { +struct curve : public geom_item {}; + +struct line : public curve { + point3 origin; + direction3 direction; + virtual item* clone() const { return new line(*this); } virtual kinds kind() const { return LINE; } }; -struct circle : public geom_item { +struct circle : public curve { + point3 origin; + direction3 x; + direction3 z; + double radius; + virtual item* clone() const { return new circle(*this); } virtual kinds kind() const { return CIRCLE; } }; -struct ellipse : public geom_item { +struct ellipse : public circle { + double radius2; + virtual item* clone() const { return new ellipse(*this); } virtual kinds kind() const { return ELLIPSE; } }; -struct bspline : public geom_item { - virtual item* clone() const { return new bspline(*this); } - virtual kinds kind() const { return BSPLINE; } +struct bspline_curve : public curve { + virtual item* clone() const { return new bspline_curve(*this); } + virtual kinds kind() const { return BSPLINE_CURVE; } }; -typedef boost::variant curve; - -struct edge : public geom_item { +struct trimmed_curve : public curve { boost::variant start, end; - boost::optional basis; + // @todo somehow account for the fact that curve in IFC can be trimmed curve, polyline and composite curve as well. + curve* basis; bool orientation; - edge() : orientation(true) {} - virtual item* clone() const { return new edge(*this); } - virtual kinds kind() const { return EDGE; } - + trimmed_curve() : basis(nullptr), orientation(true) {} + virtual void reverse() { std::swap(start, end); orientation = !orientation; } }; +struct edge : public trimmed_curve { + // @todo how to express similarity between trimmed_curve and edge? + virtual item* clone() const { return new edge(*this); } + virtual kinds kind() const { return EDGE; } +}; + struct collection : public geom_item { std::vector children; + template + std::vector children_as() const { + std::vector ts; + ts.reserve(children.size()); + std::for_each(children.begin(), children.end(), [&ts](item* i){ + auto v = dynamic_cast(i); + if (v) { + ts.push_back(v); + } + }); + return ts; + } + virtual item* clone() const { return new collection(*this); } virtual kinds kind() const { return COLLECTION; } virtual void reverse() { @@ -175,6 +203,8 @@ struct face : public collection { }; struct loop : public collection { + boost::optional external; + virtual item* clone() const { return new loop(*this); } virtual kinds kind() const { return LOOP; } }; @@ -205,7 +235,8 @@ struct node : public geom_item { }; namespace impl { - typedef std::tuple KindsTuple; + typedef std::tuple KindsTuple; + typedef std::tuple CurvesTuple; } struct type_by_kind { @@ -215,6 +246,14 @@ struct type_by_kind { static const size_t max = std::tuple_size< impl::KindsTuple>::value; }; +struct curves { + template + using type = typename std::tuple_element::type; + + static const size_t max = std::tuple_size< impl::CurvesTuple>::value; +}; + + }