new_helper

src/ifcgeom/IfcGeom.h

@@ -121,25 +121,29 @@ private:
 	manifold Shell / Solid from this set of faces. Only IfcPolyLoop instances are used. Points within the tolerance
 	threshiold are merged, so consider points a, b, c, distance(a, b) < eps then M(a, b) = Null, M(a, b) = M(a, c).
 	*/
+
+	template <typename CP=const IfcSchema::IfcCartesianPoint*, typename LP=const IfcSchema::IfcPolyLoop*>
 	class faceset_helper {
 	private:
 		MAKE_TYPE_NAME(Kernel)* kernel_;
-		std::set<const IfcSchema::IfcPolyLoop*> duplicates_;
-		std::map<int, int> vertex_mapping_;
+		std::set<LP> duplicates_;
+		std::map<const void*, int> vertex_mapping_;
 		std::map<std::pair<int, int>, TopoDS_Edge> edges_;
+		// not always in use
+		const std::vector<std::vector<double>>* points_ = nullptr;
 		double eps_;
 		bool non_manifold_;
 
 		template <typename Fn>
-		void loop_(IfcSchema::IfcCartesianPoint::list::ptr& ps, const Fn& callback) {
-			if (ps->size() < 3) {
+		void loop_(const LP& lp, const Fn& callback) {
+			auto ps = get_idxs(lp);
+
+			if (ps.size() < 3) {
 				return;
 			}
 
-			auto a = *(ps->end() - 1);
-			auto A = a->data().id();
-			for (auto& b : *ps) {
-				auto B = b->data().id();
+			auto A = ps.back();
+			for (auto& B : ps) {
 				auto C = vertex_mapping_[A], D = vertex_mapping_[B];
 				bool fwd = C < D;
 				if (!fwd) {
@@ -151,23 +155,66 @@ private:
 				}
 			}
 		}
+
+		bool construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l);
+		bool construct(const std::vector<double>& cp, gp_Pnt* l);
+
+		const void* get_idx(const IfcSchema::IfcCartesianPoint* cp) {
+			return cp;
+		}
+
+		const void* get_idx(const std::vector<double>& cp) {
+			return &cp;
+		}
+
+		std::vector<const void*> get_idxs(const IfcSchema::IfcPolyLoop* lp) {
+			auto poly = lp->Polygon();
+			std::vector<const void*> idxs;
+			std::transform(poly->begin(), poly->end(), std::back_inserter(idxs), [this](const IfcSchema::IfcCartesianPoint* p) {return get_idx(p); });
+			return idxs;
+		}
+
+		std::vector<const void*> get_idxs(const std::vector<int>& it) {
+			std::vector<const void*> idxs;
+			std::transform(it.begin(), it.end(), std::back_inserter(idxs), [this](int i) { return get_idx((*points_)[i - 1]); });
+			return idxs;
+		}
+
+		/*
+		std::vector<uintptr_t> get_idxs(std::vector<std::vector<int>>::const_iterator it) {
+			std::vector<uintptr_t> idxs;
+			std::transform(it->begin(), it->end(), std::back_inserter(idxs), [this](int i) {return get_idx(i); });
+			return idxs;
+		}
+		*/
+
 	public:
+		/*
 		faceset_helper(MAKE_TYPE_NAME(Kernel)* kernel, const IfcSchema::IfcConnectedFaceSet* l);
+		*/
+
+		faceset_helper(
+			MAKE_TYPE_NAME(Kernel)* kernel, 
+			const std::vector<CP>& points,
+			const std::vector<LP>& indices,
+			bool should_by_closed);
 
 		~faceset_helper();
 
 		bool non_manifold() const { return non_manifold_; }
 		bool& non_manifold() { return non_manifold_; }
 
-		bool edge(const IfcSchema::IfcCartesianPoint* a, const IfcSchema::IfcCartesianPoint* b, TopoDS_Edge& e) {
-			int A = vertex_mapping_[a->data().id()];
-			int B = vertex_mapping_[b->data().id()];
+		/*
+		bool edge(CP a, CP b, TopoDS_Edge& e) {
+			int A = vertex_mapping_[get_idx(a)];
+			int B = vertex_mapping_[get_idx(b)];
 			if (A == B) {
 				return false;
 			}
 
 			return edge(A, B, e);
 		}
+		*/
 
 		bool edge(int A, int B, TopoDS_Edge& e) {
 			auto it = edges_.find({A, B});
@@ -178,15 +225,14 @@ private:
 			return true;
 		}
 
-		bool wire(const IfcSchema::IfcPolyLoop* loop, TopoDS_Wire& wire) {
+		bool wire(LP loop, TopoDS_Wire& wire) {
 			if (duplicates_.find(loop) != duplicates_.end()) {
 				return false;
 			}
 			BRep_Builder builder;
 			builder.MakeWire(wire);
 			int count = 0;
-			auto ps = loop->Polygon();
-			loop_(ps, [this, &builder, &wire, &count](int A, int B, bool fwd) {
+			loop_(loop, [this, &builder, &wire, &count](int A, int B, bool fwd) {
 				TopoDS_Edge e;
 				if (edge(A, B, e)) {
 					if (!fwd) {
@@ -201,7 +247,7 @@ private:
 
 				TopTools_ListOfShape results;
 				if (kernel_->wire_intersections(wire, results)) {
-					Logger::Warning("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected", loop);
+					Logger::Warning("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected");
 					kernel_->select_largest(results, wire);
 					non_manifold_ = true;
 				}
@@ -231,7 +277,7 @@ private:
 	// For stopping PlacementRelTo recursion in convert(const IfcSchema::IfcObjectPlacement* l, gp_Trsf& trsf)
 	const IfcParse::declaration* placement_rel_to;
 
-	faceset_helper* faceset_helper_;
+	faceset_helper<>* faceset_helper_;
 	double disable_boolean_result;
 
 	gp_Vec offset = gp_Vec{0.0, 0.0, 0.0};

src/ifcgeom/IfcGeomFunctions.cpp

@@ -4426,11 +4426,32 @@ namespace {
 	}
 }
 
-IfcGeom::Kernel::faceset_helper::~faceset_helper() {
+template <typename CP, typename LP>
+IfcGeom::Kernel::faceset_helper<CP, LP>::~faceset_helper() {
+	// @todo this is super ugly, but how else can we be notified that the unique_ptr goes out of scope?
+	// Perhaps just supply a custom std::deleter?
 	kernel_->faceset_helper_ = nullptr;
 }
 
-IfcGeom::Kernel::faceset_helper::faceset_helper(Kernel* kernel, const IfcSchema::IfcConnectedFaceSet* l)
+
+template <typename CP, typename LP>
+bool IfcGeom::Kernel::faceset_helper<CP, LP>::construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l) {
+	return kernel_->convert(cp, *l);
+}
+
+template <typename CP, typename LP>
+bool IfcGeom::Kernel::faceset_helper<CP, LP>::construct(const std::vector<double>& cp, gp_Pnt* l) {
+	if (cp.size() != 3) {
+		return false;
+	}
+	auto LU = kernel_->getValue(GV_LENGTH_UNIT);
+	l->SetCoord(cp[0] * LU, cp[1] * LU, cp[2] * LU);
+}
+
+/*
+
+template <typename CP, typename LP>
+IfcGeom::Kernel::faceset_helper<CP, LP>::faceset_helper(Kernel* kernel, const IfcSchema::IfcConnectedFaceSet* l)
 	: kernel_(kernel) 
 	, non_manifold_(false)
 {
@@ -4519,7 +4540,7 @@ IfcGeom::Kernel::faceset_helper::faceset_helper(Kernel* kernel, const IfcSchema:
 				for (int v : vs) {
 					auto pt = *(points->begin() + v);
 					// NB: insert() ignores duplicate keys
-					vertex_mapping_.insert({ pt->data().id() , pnt_i });
+					vertex_mapping_.insert({ get_idx(pt), pnt_i });
 				}
 			}
 		}		
@@ -4576,3 +4597,165 @@ IfcGeom::Kernel::faceset_helper::faceset_helper(Kernel* kernel, const IfcSchema:
 		Logger::Warning(boost::lexical_cast<std::string>(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast<std::string>(loops_removed) + " loops removed and " + boost::lexical_cast<std::string>(non_manifold) + " non-manifold edges for:", l);
 	}
 }
+
+*/
+
+namespace {
+	const std::vector<std::vector<double>>* store_cache(const std::vector<std::vector<double>>& p) {
+		return &p;
+	}
+
+	const std::vector<std::vector<double>>* store_cache(const std::vector<const IfcSchema::IfcCartesianPoint*>& p) {
+		return nullptr;
+	}
+}
+
+template <typename CP, typename LP>
+IfcGeom::Kernel::faceset_helper<CP, LP>::faceset_helper(
+	Kernel* kernel,
+	const std::vector<CP>& points,
+	const std::vector<LP>& indices,
+	bool should_be_closed
+)
+	: kernel_(kernel)
+	, non_manifold_(false)
+	, points_(store_cache(points))
+{
+	std::vector<std::unique_ptr<gp_Pnt>> pnts(std::distance(points.begin(), points.end()));
+	std::vector<TopoDS_Vertex> vertices(pnts.size());
+
+	auto LU = kernel_->getValue(GV_LENGTH_UNIT);
+
+	IfcGeom::impl::tree<int> tree;
+
+	BRep_Builder B;
+
+	Bnd_Box box;
+	for (size_t i = 0; i < points.size(); ++i) {
+		gp_Pnt* p = new gp_Pnt;
+		if (construct(points[i], p)) {
+			pnts[i].reset(p);
+			B.MakeVertex(vertices[i], *p, Precision::Confusion());
+			tree.add(i, vertices[i]);
+			box.Add(*p);
+		} else {
+			delete p;
+		}
+	}
+
+	// Use the bbox diagonal to influence local epsilon
+	// double bdiff = std::sqrt(box.SquareExtent());
+
+	// @todo the bounding box diagonal is not used (see above)
+	// because we're explicitly interested in the miminal
+	// dimension of the element to limit the tolerance (for sheet-
+	// like elements for example). But the way below is very
+	// dependent on orientation due to the usage of the
+	// axis-aligned bounding box. Use PCA to find three non-aligned
+	// set of dimensions and use the one with the smallest eigenvalue.
+
+	// Find the minimal bounding box edge
+	double bmin[3], bmax[3];
+	box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]);
+	double bdiff = std::numeric_limits<double>::infinity();
+	for (size_t i = 0; i < 3; ++i) {
+		const double d = bmax[i] - bmin[i];
+		if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) {
+			bdiff = d;
+		}
+	}
+
+	eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff);
+
+	size_t loops_removed, non_manifold, duplicate_faces;
+
+	std::map<std::pair<int, int>, int> edge_use;
+
+	for (int i = 0; i < 3; ++i) {
+		// Some times files, have large tolerance values specified collapsing too many vertices.
+		// This case we detect below and re-run the loop with smaller epsilon. Normally
+		// the body of this loop would only be executed once.
+
+		loops_removed = 0;
+		non_manifold = 0;
+		duplicate_faces = 0;
+
+		vertex_mapping_.clear();
+		duplicates_.clear();
+
+		edge_use.clear();
+
+		if (eps_ < Precision::Confusion()) {
+			// occt uses some hard coded precision values, don't go smaller than that.
+			// @todo, can be reset though with BRepLib::Precision(double)
+			eps_ = Precision::Confusion();
+		}
+
+		for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) {
+			if (pnts[pnt_i]) {
+				std::set<int> vs;
+				find_neighbours(tree, pnts, vs, pnt_i, eps_);
+
+				for (int v : vs) {
+					// NB: insert() ignores duplicate keys
+					// v-1?
+					vertex_mapping_.insert({ get_idx(points[v]), pnt_i });
+				}
+			}
+		}
+
+		typedef std::array<int, 2> edge_t;
+		typedef std::set<edge_t> edge_set_t;
+		std::set<edge_set_t> edge_sets;
+
+		for (auto ps = indices.begin(); ps != indices.end(); ++ps) {
+			std::vector<std::pair<int, int> > segments;
+			edge_set_t segment_set;
+
+			loop_(*ps, [&segments, &segment_set](int C, int D, bool) {
+				segment_set.insert(edge_t{ C,D });
+				segments.push_back(std::make_pair(C, D));
+			});
+
+			if (edge_sets.find(segment_set) != edge_sets.end()) {
+				duplicate_faces++;
+				duplicates_.insert(*ps);
+				continue;
+			}
+			edge_sets.insert(segment_set);
+
+			if (segments.size() >= 3) {
+				for (auto& p : segments) {
+					edge_use[p] ++;
+				}
+			}
+			else {
+				loops_removed += 1;
+			}
+		}
+
+		if (edge_use.size() != 0) {
+			break;
+		}
+		else {
+			eps_ /= 10.;
+		}
+	}
+
+	for (auto& p : edge_use) {
+		int a, b;
+		std::tie(a, b) = p.first;
+		edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]);
+
+		if (p.second != 2) {
+			non_manifold += 1;
+		}
+	}
+
+	if (loops_removed || (non_manifold && should_be_closed)) {
+		Logger::Warning(boost::lexical_cast<std::string>(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast<std::string>(loops_removed) + " loops removed and " + boost::lexical_cast<std::string>(non_manifold) + " non-manifold edges");
+	}
+}
+
+template class IfcGeom::Kernel::faceset_helper<const IfcSchema::IfcCartesianPoint*, const IfcSchema::IfcPolyLoop*>;
+template class IfcGeom::Kernel::faceset_helper<std::vector<double>, std::vector<int>>;

src/ifcgeom/IfcGeomShapes.cpp

@@ -725,8 +725,22 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcBooleanResult* l, TopoDS_Shape
 }
 
 bool IfcGeom::Kernel::convert(const IfcSchema::IfcConnectedFaceSet* l, TopoDS_Shape& shape) {
-	std::unique_ptr<faceset_helper> helper_scope;
-	helper_scope.reset(new faceset_helper(this, l));
+	std::unique_ptr<faceset_helper<>> helper_scope;
+
+	IfcSchema::IfcCartesianPoint::list::ptr points = IfcParse::traverse((IfcUtil::IfcBaseClass*) l)->as<IfcSchema::IfcCartesianPoint>();
+	std::vector<const IfcSchema::IfcCartesianPoint*> points_(points->begin(), points->end());
+
+	IfcSchema::IfcPolyLoop::list::ptr loops = IfcParse::traverse((IfcUtil::IfcBaseClass*)l)->as<IfcSchema::IfcPolyLoop>();
+	std::vector<const IfcSchema::IfcPolyLoop*> loops_(loops->begin(), loops->end());
+
+	helper_scope.reset(new faceset_helper<>(
+		this,
+		points_,
+		loops_,
+		l->declaration().is(IfcSchema::IfcClosedShell::Class())
+	));
+
+	faceset_helper_ = helper_scope.get();
 
 	IfcSchema::IfcFace::list::ptr faces = l->CfsFaces();
 
@@ -1643,259 +1657,129 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcAdvancedBrep* l, TopoDS_Shape&
 
 bool IfcGeom::Kernel::convert(const IfcSchema::IfcTriangulatedFaceSet* l, TopoDS_Shape& shape) {
 	IfcSchema::IfcCartesianPointList3D* point_list = l->Coordinates();
-	const std::vector< std::vector<double> > coordinates = point_list->CoordList();
-	std::vector<gp_Pnt> points;
-	points.reserve(coordinates.size());
-	for (std::vector< std::vector<double> >::const_iterator it = coordinates.begin(); it != coordinates.end(); ++it) {
-		const std::vector<double>& coords = *it;
-		if (coords.size() != 3) {
-			Logger::Message(Logger::LOG_ERROR, "Invalid dimensions encountered on Coordinates", l);
-			return false;
-		}
-		points.push_back(gp_Pnt(coords[0] * getValue(GV_LENGTH_UNIT),
-		                        coords[1] * getValue(GV_LENGTH_UNIT),
-		                        coords[2] * getValue(GV_LENGTH_UNIT)));
-	}
+	auto coord_list = point_list->CoordList();
+	std::vector<std::vector<int>> indices = l->CoordIndex();
 
-	std::vector< std::vector<int> > indices = l->CoordIndex();
-	
-	std::vector<TopoDS_Face> faces;
-	faces.reserve(indices.size());
-
-	for(std::vector< std::vector<int> >::const_iterator it = indices.begin(); it != indices.end(); ++ it) {
-		const std::vector<int>& tri = *it;
-		if (tri.size() != 3) {
-			Logger::Message(Logger::LOG_ERROR, "Invalid dimensions encountered on CoordIndex", l);
-			return false;
-		}
-
-		const int min_index = *std::min_element(tri.begin(), tri.end());
-		const int max_index = *std::max_element(tri.begin(), tri.end());
-
-		if (min_index < 1 || max_index > (int) points.size()) {
-			Logger::Message(Logger::LOG_ERROR, "Contents of CoordIndex out of bounds", l);
-			return false;
-		}
-
-		if (tri[0] == tri[1] || tri[1] == tri[2] || tri[0] == tri[2]) {
-			auto tri_0 = boost::lexical_cast<std::string>(tri[0]);
-			auto tri_1 = boost::lexical_cast<std::string>(tri[1]);
-			auto tri_2 = boost::lexical_cast<std::string>(tri[2]);
-			Logger::Message(Logger::LOG_ERROR, "Degenerate triangle indices, skipping (" + tri_0 + "," + tri_1 + "," + tri_2 + ")", l);
-			continue;
-		}
-		
-		const gp_Pnt& a = points[tri[0] - 1]; // account for zero- vs
-		const gp_Pnt& b = points[tri[1] - 1]; // one-based indices in
-		const gp_Pnt& c = points[tri[2] - 1]; // c++ and express
-
-		BRepBuilderAPI_MakePolygon mp(a, b, c, true);
-
-		if (!mp.IsDone()) {
-			auto tri_0 = boost::lexical_cast<std::string>(tri[0]);
-			auto tri_1 = boost::lexical_cast<std::string>(tri[1]);
-			auto tri_2 = boost::lexical_cast<std::string>(tri[2]);
-			Logger::Message(Logger::LOG_ERROR, "Degenerate triangle, skipping (" + tri_0 + "," + tri_1 + "," + tri_2 + ")", l);
-			continue;
-		}
-
-		TopoDS_Wire wire = mp.Wire();
-		TopoDS_Face face = BRepBuilderAPI_MakeFace(wire).Face();
-
-		TopoDS_Iterator face_it(face, false);
-		const TopoDS_Wire& w = TopoDS::Wire(face_it.Value());
-		const bool reversed = w.Orientation() == TopAbs_REVERSED;
-		if (reversed) {
-			face.Reverse();
-		}
-
-		if (face_area(face) > getValue(GV_MINIMAL_FACE_AREA)) {
-			faces.push_back(face);
-		}
-	}
-
-	if (faces.empty()) return false;
-	
-	bool valid_shell = false;
-
-	// @todo Do this more efficiently by creating proper half-edge pairs.
-	BRepOffsetAPI_Sewing sewing_builder;
-	sewing_builder.SetTolerance(getValue(GV_POINT_EQUALITY_TOLERANCE));
-	sewing_builder.SetMaxTolerance(getValue(GV_POINT_EQUALITY_TOLERANCE));
-	sewing_builder.SetMinTolerance(getValue(GV_POINT_EQUALITY_TOLERANCE));
-		
-	for (std::vector<TopoDS_Face>::const_iterator it = faces.begin(); it != faces.end(); ++it) {
-		sewing_builder.Add(*it);
-	}
-
-	try {
-		sewing_builder.Perform();
-		shape = sewing_builder.SewedShape();
-		valid_shell = BRepCheck_Analyzer(shape).IsValid();
-	} catch(...) {}
-
-	if (valid_shell) {
-		try {
-			ShapeFix_Solid solid;
-			solid.LimitTolerance(getValue(GV_POINT_EQUALITY_TOLERANCE));
-			TopoDS_Solid solid_shape = solid.SolidFromShell(TopoDS::Shell(shape));
-			if (!solid_shape.IsNull()) {
-				try {
-					BRepClass3d_SolidClassifier classifier(solid_shape);
-					shape = solid_shape;
-				} catch (...) {}
-			}
-		} catch(...) {}
-	} else {
-		Logger::Message(Logger::LOG_WARNING, "Failed to sew faceset:", l);
-	}
-
-	if (!valid_shell) {
-		TopoDS_Compound compound;
-		BRep_Builder builder;
-		builder.MakeCompound(compound);
-		
-		for (std::vector<TopoDS_Face>::const_iterator it = faces.begin(); it != faces.end(); ++it) {
-			builder.Add(compound, *it);
-		}
-
-		shape = compound;
-	}
-
-	return true;
-}
-
-namespace {
-	bool make_indexed_polygon(IfcGeom::Kernel& k, const std::vector<gp_Pnt>& points, const std::vector<int>& indices, TopoDS_Wire& wire) {
-		TColgp_SequenceOfPnt polygon;
-		for (std::vector<int>::size_type j = 0; j != indices.size(); j++) {
-			const gp_Pnt& point = points[indices[j] - 1];
-			polygon.Append(point);
-		}
-		k.remove_duplicate_points_from_loop(polygon, true);
-
-		if (polygon.Length() < 3) {
-			return false;
-		}
-
-		BRepBuilderAPI_MakePolygon wire_builder;
-		for (int i = 1; i <= polygon.Length(); ++i) {
-			wire_builder.Add(polygon.Value(i));
-		}
-		wire_builder.Close();
-
-		wire = wire_builder.Wire();
-
-		TopoDS_Iterator it(wire);
-		for (; it.More(); it.Next()) {
-			BRepAdaptor_Curve ad(TopoDS::Edge(it.Value()));
-		}
-
-		ShapeFix_ShapeTolerance FTol;
-		FTol.SetTolerance(wire, k.getValue(IfcGeom::Kernel::GV_PRECISION), TopAbs_WIRE);
-
-		return true;
-	}
-}
-
-bool IfcGeom::Kernel::convert(const IfcSchema::IfcPolygonalFaceSet* pfs, TopoDS_Shape& shape) {
-    IfcSchema::IfcCartesianPointList3D* point_list = pfs->Coordinates();
-    const std::vector<std::vector<double> > coordinates = point_list->CoordList();
-    std::vector<gp_Pnt> points;
-    points.reserve(coordinates.size());
-
-    for (std::vector<std::vector<double> >::const_iterator it = coordinates.begin(); it != coordinates.end(); ++it) {
-        const std::vector<double>& coords = *it;
-        points.push_back(gp_Pnt(
-			coords[0] * getValue(GV_LENGTH_UNIT),
-			coords[1] * getValue(GV_LENGTH_UNIT),
-			coords[2] * getValue(GV_LENGTH_UNIT)));
-    }
-
-    auto polygonal_faces = pfs->Faces();
+	faceset_helper<
+		std::vector<double>,
+		std::vector<int>
+	> helper(this, coord_list, indices, l->hasClosed() ? l->Closed() : false);
 
 	TopTools_ListOfShape faces;
 
-    for (unsigned i = 0; i < polygonal_faces->size(); i++) {
-        IfcSchema::IfcIndexedPolygonalFace* la = (IfcSchema::IfcIndexedPolygonalFace*)*(polygonal_faces->begin() + i);
-        TopoDS_Face face;
-		TopoDS_Wire wire;
-		if (!make_indexed_polygon(*this, points, la->CoordIndex(), wire)) {
+	for (auto it = indices.begin(); it != indices.end(); ++it) {
+		TopoDS_Wire w;
+		if (helper.wire(*it, w)) {
+			BRepBuilderAPI_MakeFace mf(w);
+			if (mf.IsDone()) {
+				faces.Append(mf.Face());
+			}
+		}
+	}
+	
+	return create_solid_from_faces(faces, shape);
+}
+
+bool IfcGeom::Kernel::convert(const IfcSchema::IfcPolygonalFaceSet* pfs, TopoDS_Shape& shape) {
+	IfcSchema::IfcCartesianPointList3D* point_list = pfs->Coordinates();
+	auto coord_list = point_list->CoordList();
+	auto polygonal_faces = pfs->Faces();
+
+	std::vector<std::vector<int>> indices;
+	indices.reserve(polygonal_faces->size() * 2);
+
+	std::vector<std::vector<int>> loop_grouping;
+	loop_grouping.reserve(polygonal_faces->size());
+	
+	for (auto& f : *polygonal_faces) {
+		loop_grouping.emplace_back();
+		loop_grouping.back().push_back(indices.size());
+		indices.push_back(f->CoordIndex());
+		if (f->as<IfcSchema::IfcIndexedPolygonalFaceWithVoids>()) {
+			auto inner_coordinates = f->as<IfcSchema::IfcIndexedPolygonalFaceWithVoids>()->InnerCoordIndices();
+			for (auto& x : inner_coordinates) {
+				loop_grouping.back().push_back(indices.size());
+				indices.push_back(x);
+			}
+		}
+	}
+
+	faceset_helper<
+		std::vector<double>,
+		std::vector<int>
+	> helper(this, coord_list, indices, pfs->hasClosed() ? pfs->Closed() : false);
+
+	TopTools_ListOfShape faces;
+
+	for (auto& f : loop_grouping) {
+		bool not_planar = false;
+		
+		TopoDS_Wire w;
+		if (!helper.wire(indices[f[0]], w)) {
 			continue;
 		}
 
-        if (la->declaration().is(IfcSchema::IfcIndexedPolygonalFaceWithVoids::Class())) {
-            IfcSchema::IfcIndexedPolygonalFaceWithVoids* converted = (IfcSchema::IfcIndexedPolygonalFaceWithVoids*)la;
-            std::vector<std::vector<int> > innercoordinates = converted->InnerCoordIndices();
+		TopoDS_Face face;
+		std::vector<TopoDS_Wire> ws = { w };
 
-            BRepBuilderAPI_MakeFace facemaker = BRepBuilderAPI_MakeFace(wire);
-			std::vector<TopoDS_Wire> vectorofwires{ wire };
-            for (std::vector<std::vector<int> >::const_iterator it = innercoordinates.begin(); it != innercoordinates.end(); ++it) {
-				TopoDS_Wire inner_wire;
-				if (make_indexed_polygon(*this, points, *it, inner_wire)) {
-					vectorofwires.push_back(inner_wire);
-					facemaker.Add(inner_wire);
-				}
-            }
-
-			facemaker.Build();
-			if (facemaker.Error() == BRepBuilderAPI_FaceDone) {
-				face = facemaker.Face();
-			} else if (facemaker.Error() == BRepBuilderAPI_NotPlanar) {
-				TopTools_ListOfShape fs;
-				if (triangulate_wire(vectorofwires, fs)) {
-					Logger::Warning("Triangulated face boundary:", la);
-					TopTools_ListIteratorOfListOfShape it(fs);
-					for (; it.More(); it.Next()) {
-						const TopoDS_Face& tri = TopoDS::Face(it.Value());
-						if (face_area(tri) > getValue(GV_MINIMAL_FACE_AREA)) {
-							faces.Append(tri);
-						}
-					}
-					continue;
-				}
-			}
-        } else {
-			BRepBuilderAPI_MakeFace facemaker(wire);
-			facemaker.Build();
-			if (facemaker.Error() == BRepBuilderAPI_FaceDone) {
-				face = facemaker.Face();
-			} else if (facemaker.Error() == BRepBuilderAPI_NotPlanar) {
-				TopTools_ListOfShape fs;
-				if (triangulate_wire({ wire }, fs)) {
-					Logger::Warning("Triangulated face boundary:", la);
-					TopTools_ListIteratorOfListOfShape it(fs);
-					for (; it.More(); it.Next()) {
-						const TopoDS_Face& tri = TopoDS::Face(it.Value());
-						if (face_area(tri) > getValue(GV_MINIMAL_FACE_AREA)) {
-							faces.Append(tri);
-						}
-					}
-					continue;
-				}
-			}
-        }
-
-		if (face.IsNull()) {
-			Logger::Warning("Face creation failed:", la);
+		// @todo triangulate
+		BRepBuilderAPI_MakeFace mf(w);
+		if (mf.Error() == BRepBuilderAPI_NotPlanar) {
+			not_planar = true;
+		} else if (mf.IsDone()) {
+			face = mf.Face();
+		} else {
+			// todo log
 			continue;
 		}
 
-        TopoDS_Iterator face_it(face, false);
-        const TopoDS_Wire& w = TopoDS::Wire(face_it.Value());
-        const bool reversed = w.Orientation() == TopAbs_REVERSED;
-        if (reversed) {
-            face.Reverse();
-        }
+		if (f.size() > 1) {
 
-        if (face_area(face) > getValue(GV_MINIMAL_FACE_AREA)) {
-            faces.Append(face);
-        }
-    }
+			if (not_planar) {
+				for (auto it = f.begin() + 1; it != f.end(); ++it) {
+					TopoDS_Wire w2;
+					if (helper.wire(indices[*it], w2)) {
+						ws.push_back(w2);
+					}
+				}
+			} else {
+				BRepBuilderAPI_MakeFace mf2(face);
+				for (auto it = f.begin() + 1; it != f.end(); ++it) {
+					TopoDS_Wire w2;
+					if (helper.wire(indices[*it], w2)) {
+						mf2.Add(w2);
+						ws.push_back(w2);
+					}
 
-    if (faces.IsEmpty()) return false;
+				}
+				
+				if (mf2.Error() == BRepBuilderAPI_NotPlanar) {
+					not_planar = true;
+				} else if (mf2.IsDone()) {
+					face = mf2.Face();
+				}
+			}	
+			
+		}
 
-    return create_solid_from_faces(faces, shape);
+		if (not_planar) {
+			TopTools_ListOfShape fs;
+			if (triangulate_wire(ws, fs)) {
+				Logger::Warning("Triangulated face boundary:", pfs);
+				TopTools_ListIteratorOfListOfShape it(fs);
+				for (; it.More(); it.Next()) {
+					const TopoDS_Face& tri = TopoDS::Face(it.Value());
+					if (face_area(tri) > getValue(GV_MINIMAL_FACE_AREA)) {
+						faces.Append(tri);
+					}
+				}
+			}
+		} else {
+			faces.Append(face);
+		}
+	}
+
+	return create_solid_from_faces(faces, shape);
 }
 
 #endif