Quantities in convert and geomserver. Faceset helper for creating edge pairs.

src/ifcconvert/IfcConvert.cpp

@@ -69,7 +69,7 @@ void print_usage(bool suggest_help = true)
 {
     std::cout << "Usage: IfcConvert [options] <input.ifc> [<output>]\n"
         << "\n"
-        << "Converts the geometry in an IFC file into one of the following formats:\n"
+        << "Converts (the geometry in) an IFC file into one of the following formats:\n"
         << "  .obj   WaveFront OBJ  (a .mtl file is also created)\n"
 #ifdef WITH_OPENCOLLADA
         << "  .dae   Collada        Digital Assets Exchange\n"
@@ -78,7 +78,8 @@ void print_usage(bool suggest_help = true)
         << "  .igs   IGES           Initial Graphics Exchange Specification\n"
         << "  .xml   XML            Property definitions and decomposition tree\n"
         << "  .svg   SVG            Scalable Vector Graphics (2D floor plan)\n"
-        << "\n"
+		<< "  .ifc   IFC-SPF        Industry Foundation Classes\n"
+		<< "\n"
         << "If no output filename given, <input>." + DEFAULT_EXTENSION + " will be used as the output file.\n";
     if (suggest_help) {
         std::cout << "\nRun 'IfcConvert --help' for more information.";
@@ -120,6 +121,7 @@ bool rename_file(const std::string& old_filename, const std::string& new_filenam
 
 static std::stringstream log_stream;
 void write_log(bool);
+void fix_quantities(IfcParse::IfcFile&, bool, bool, bool);
 
 /// @todo make the filters non-global
 IfcGeom::entity_filter entity_filter; // Entity filter is used always by default.
@@ -181,7 +183,12 @@ int main(int argc, char** argv)
     std::string filter_filename;
     std::string default_material_filename;
 	
-    po::options_description geom_options("Geometry options");
+	po::options_description ifc_options("IFC options");
+	ifc_options.add_options()
+		("calculate-quantities", "Calculate or fix the physical quantity definitions "
+			"based on an interpretation of the geometry when exporting IFC");
+    
+	po::options_description geom_options("Geometry options");
 	geom_options.add_options()
 		("plan",
 			"Specifies whether to include curves in the output result. Typically "
@@ -466,23 +473,44 @@ int main(int argc, char** argv)
 
 	IfcParse::IfcFile* ifc_file = 0;
 
-    if (output_extension == ".xml") {
+	// @todo clean up serializer selection
-        int exit_code = EXIT_FAILURE;
+	// @todo detect program options that conflict with the chosen serializer
-        try {
+	if (output_extension == ".xml") {
-            if (init_input_file(input_filename, ifc_file, no_progress || quiet, mmap)) {
+		int exit_code = EXIT_FAILURE;
-                XmlSerializer s(ifc_file, output_temp_filename);
+		try {
-                Logger::Status("Writing XML output...");
+			if (init_input_file(input_filename, ifc_file, no_progress || quiet, mmap)) {
-                s.finalize();
+				XmlSerializer s(ifc_file, output_temp_filename);
-                Logger::Status("Done!");
+				Logger::Status("Writing XML output...");
-                rename_file(output_temp_filename, output_filename);
+				s.finalize();
-                exit_code = EXIT_SUCCESS;
+				Logger::Status("Done!");
-            }
+				rename_file(output_temp_filename, output_filename);
-        } catch (const std::exception& e) {
+				exit_code = EXIT_SUCCESS;
+			}
+		} catch (const std::exception& e) {
 			Logger::Error(e);
 		}
-        write_log(!quiet);
+		write_log(!quiet);
-        return exit_code;
+		return exit_code;
-    }
+	} else if (output_extension == ".ifc") {
+		int exit_code = EXIT_FAILURE;
+		try {
+			if (init_input_file(input_filename, ifc_file, no_progress || quiet, mmap)) {
+				std::ofstream fs(output_filename.c_str());
+				if (fs.is_open()) {
+					if (vmap.count("calculate-quantities")) {
+						fix_quantities(*ifc_file, no_progress, quiet, stderr_progress);
+					}
+					fs << *ifc_file;
+				} else {
+					Logger::Error("Unable to open output file for writing");
+				}
+			}
+		} catch (const std::exception& e) {
+			Logger::Error(e);
+		}
+		write_log(!quiet);
+		return exit_code;
+	}
 
     if (!filter_filename.empty()) {
         size_t num_filters = read_filters_from_file(filter_filename, include_filter, include_traverse_filter, exclude_filter, exclude_traverse_filter);
@@ -984,3 +1012,205 @@ std::vector<IfcGeom::filter_t> setup_filters(const std::vector<geom_filter>& fil
 
     return filter_funcs;
 }
+
+namespace latebound_access {
+
+	template <typename T>
+	void set(IfcUtil::IfcBaseClass* inst, const std::string& attr, T t);
+
+	template <typename T>
+	void set_enumeration(IfcUtil::IfcBaseClass*, const std::string&, const IfcParse::enumeration_type*, T) {}
+
+	template <>
+	void set_enumeration(IfcUtil::IfcBaseClass* inst, const std::string& attr, const IfcParse::enumeration_type* enum_type, std::string t) {
+		std::vector<std::string>::const_iterator it = std::find(
+			enum_type->enumeration_items().begin(),
+			enum_type->enumeration_items().end(),
+			t);
+
+		return set(inst, attr, IfcWrite::IfcWriteArgument::EnumerationReference(it - enum_type->enumeration_items().begin(), it->c_str()));
+	}
+
+	template <typename T>
+	void set(IfcUtil::IfcBaseClass* inst, const std::string& attr, T t) {
+		auto decl = inst->declaration().as_entity();
+		auto i = decl->attribute_index(attr);
+
+		auto attr_type = decl->attribute_by_index(i)->type_of_attribute();
+		if (attr_type->as_named_type() && attr_type->as_named_type()->declared_type()->as_enumeration_type()) {
+			set_enumeration(inst, attr, attr_type->as_named_type()->declared_type()->as_enumeration_type(), t);
+		}
+
+		IfcWrite::IfcWriteArgument* a = new IfcWrite::IfcWriteArgument;
+		a->set(t);
+		inst->data().attributes()[i] = a;
+	}
+
+	IfcUtil::IfcBaseClass* create(IfcParse::IfcFile& f, const std::string& entity) {
+		auto decl = f.schema()->declaration_by_name(entity);
+		auto data = new IfcEntityInstanceData(decl);
+		auto inst = f.schema()->instantiate(data);
+		if (decl->is("IfcRoot")) {
+			IfcParse::IfcGlobalId guid;
+			latebound_access::set(inst, "GlobalId", (std::string) guid);
+		}
+		return f.addEntity(inst);
+	}
+}
+
+void fix_quantities(IfcParse::IfcFile& f, bool no_progress, bool quiet, bool stderr_progress) {
+	{
+		auto delete_reversed = [&f](const IfcEntityList::ptr& insts) {
+			if (!insts) {
+				return;
+			}
+			// Lists are traversed back to front as the list may be mutated when
+			// instances are removed from the grouping by type.
+			for (auto it = insts->end() - 1; it >= insts->begin(); --it) {
+				IfcUtil::IfcBaseClass* const inst = *it;
+				f.removeEntity(inst);
+			}
+		};
+
+		// Delete quantities
+		auto quantities = f.instances_by_type("IfcPhysicalQuantity");
+		if (quantities) {
+			quantities = quantities->filtered({ f.schema()->declaration_by_name("IfcPhysicalComplexQuantity") });
+			delete_reversed(quantities);
+		}
+
+		// Delete complexes
+		delete_reversed(f.instances_by_type("IfcPhysicalComplexQuantity"));
+
+		auto element_quantities = f.instances_by_type("IfcElementQuantity");
+
+		// Capture relationship nodes
+		std::vector<IfcUtil::IfcBaseClass*> relationships;
+		auto IfcRelDefinesByProperties = f.schema()->declaration_by_name("IfcRelDefinesByProperties");
+		for (auto& eq : *element_quantities) {
+			auto rels = eq->data().getInverse(IfcRelDefinesByProperties, -1);
+			for (auto& rel : *rels) {
+				relationships.push_back(rel);
+			}
+		}
+
+		// Delete element quantities
+		delete_reversed(element_quantities);
+
+		// Delete relationship nodes
+		for (auto& rel : relationships) {
+			f.removeEntity(rel);
+		}
+	}
+
+	IfcGeom::IteratorSettings settings;
+	settings.set(IfcGeom::IteratorSettings::USE_WORLD_COORDS, false);
+	settings.set(IfcGeom::IteratorSettings::WELD_VERTICES, false);
+	settings.set(IfcGeom::IteratorSettings::SEW_SHELLS, true);
+	settings.set(IfcGeom::IteratorSettings::CONVERT_BACK_UNITS, true);
+	settings.set(IfcGeom::IteratorSettings::DISABLE_TRIANGULATION, true);
+
+	IfcGeom::Iterator<double> context_iterator(settings, &f);
+
+	if (!context_iterator.initialize()) {
+		return;
+	}
+
+	size_t num_created = 0;
+	int old_progress = quiet ? 0 : -1;
+
+	auto person = latebound_access::create(f, "IfcPerson");
+	latebound_access::set(person, "FamilyName", std::string("IfcOpenShell"));
+	latebound_access::set(person, "GivenName", std::string("IfcOpenShell"));
+	
+	auto org = latebound_access::create(f, "IfcOrganization");
+	latebound_access::set(org, "Name", std::string("IfcOpenShell"));
+	
+	auto pando = latebound_access::create(f, "IfcPersonAndOrganization");
+	latebound_access::set(pando, "ThePerson", person);
+	latebound_access::set(pando, "TheOrganization", org);
+	
+	auto application = latebound_access::create(f, "IfcApplication");
+	latebound_access::set(application, "ApplicationDeveloper", org);
+	latebound_access::set(application, "Version", std::string(IFCOPENSHELL_VERSION));
+	latebound_access::set(application, "ApplicationFullName", std::string("IfcConvert"));
+	latebound_access::set(application, "ApplicationIdentifier", std::string("IfcConvert" IFCOPENSHELL_VERSION));
+	
+	auto ownerhist = latebound_access::create(f, "IfcOwnerHistory");
+	latebound_access::set(ownerhist, "OwningUser", pando);
+	latebound_access::set(ownerhist, "OwningApplication", application);
+	latebound_access::set(ownerhist, "ChangeAction", std::string("MODIFIED"));
+	latebound_access::set(ownerhist, "CreationDate", (int)time(0));
+
+	IfcUtil::IfcBaseClass* quantity = nullptr;
+	IfcEntityList::ptr objects;
+	boost::shared_ptr<IfcGeom::Representation::BRep> previous_geometry_pointer;
+
+	do {
+		IfcGeom::BRepElement<double>* geom_object = context_iterator.get_native();
+
+		if (geom_object->geometry_pointer() == previous_geometry_pointer) {
+			objects->push(geom_object->product());
+		} else {
+			if (quantity) {
+				auto rel = latebound_access::create(f, "IfcRelDefinesByProperties");
+				latebound_access::set(rel, "OwnerHistory", ownerhist);
+				latebound_access::set(rel, "RelatedObjects", objects);
+				latebound_access::set(rel, "RelatingPropertyDefinition", quantity);
+			}
+
+			IfcEntityList::ptr quantities(new IfcEntityList);
+
+			double a, b, c;
+			if (geom_object->geometry().calculate_surface_area(a)) {
+				auto quantity_area = latebound_access::create(f, "IfcQuantityArea");
+				latebound_access::set(quantity_area, "Name", std::string("Total Surface Area"));
+				latebound_access::set(quantity_area, "AreaValue", a);
+				quantities->push(quantity_area);
+			}
+			
+			if (geom_object->geometry().calculate_volume(a)) {
+				auto quantity_volume = latebound_access::create(f, "IfcQuantityVolume");
+				latebound_access::set(quantity_volume, "Name", std::string("Volume"));
+				latebound_access::set(quantity_volume, "VolumeValue", a);
+				quantities->push(quantity_volume);
+			}
+
+			if (geom_object->calculate_projected_surface_area(a, b, c)) {
+				auto quantity_area = latebound_access::create(f, "IfcQuantityArea");
+				latebound_access::set(quantity_area, "Name", std::string("Footprint Area"));
+				latebound_access::set(quantity_area, "AreaValue", c);
+				quantities->push(quantity_area);
+			}
+
+			if (quantities->size()) {
+				quantity = latebound_access::create(f, "IfcElementQuantity");
+				latebound_access::set(quantity, "OwnerHistory", ownerhist);
+				latebound_access::set(quantity, "Quantities", quantities);
+			}
+
+			objects.reset(new IfcEntityList);
+			objects->push(geom_object->product());
+		}
+
+		previous_geometry_pointer = geom_object->geometry_pointer();
+
+		if (!no_progress) {
+			if (quiet) {
+				const int progress = context_iterator.progress();
+				for (; old_progress < progress; ++old_progress) {
+					std::cout << ".";
+					if (stderr_progress)
+						std::cerr << ".";
+				}
+				std::cout << std::flush;
+				if (stderr_progress)
+					std::cerr << std::flush;
+			} else {
+				const int progress = context_iterator.progress() / 2;
+				if (old_progress != progress) Logger::ProgressBar(progress);
+				old_progress = progress;
+			}
+		}
+	} while (++num_created, context_iterator.next());
+}

src/ifcgeom/IfcGeom.h

@@ -45,6 +45,8 @@ inline static bool ALMOST_THE_SAME(const T& a, const T& b, double tolerance=ALMO
 #include <TColgp_SequenceOfPnt.hxx>
 #include <TopTools_ListOfShape.hxx>
 #include <BOPAlgo_Operation.hxx>
+#include <BRep_Builder.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
 
 #include "../ifcparse/macros.h"
 #include "../ifcparse/IfcParse.h"
@@ -107,6 +109,82 @@ public:
 class IFC_GEOM_API MAKE_TYPE_NAME(Kernel) : public IfcGeom::Kernel {
 private:
 
+	/*
+	faceset_helper traverses the forward instance references of IfcConnectedFaceSet and then provides a mapping
+	M of (IfcCartesianPoint, IfcCartesianPoint) -> TopoDS_Edge, where M(a, b) is a partner of M(b, a), ie share
+	the same underlying edge but with orientation reversed. This then later speeds op the process of creating a
+	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).
+	*/
+	class faceset_helper {
+	private:
+		MAKE_TYPE_NAME(Kernel)* kernel_;
+		std::map<int, int> vertex_mapping_;
+		std::map<std::pair<int, int>, TopoDS_Edge> edges_;
+
+		template <typename Fn>
+		void loop_(IfcSchema::IfcCartesianPoint::list::ptr& ps, const Fn& callback) {
+			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 C = vertex_mapping_[A], D = vertex_mapping_[B];
+				bool fwd = C < D;
+				if (!fwd) {
+					std::swap(C, D);
+				}
+				if (C != D) {
+					callback(C, D, fwd);
+					A = B;
+				}
+			}
+		}
+	public:
+		faceset_helper(MAKE_TYPE_NAME(Kernel)* kernel, const IfcSchema::IfcConnectedFaceSet* l);
+
+		~faceset_helper();
+
+		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()];
+			if (A == B) {
+				return false;
+			}
+
+			return edge(A, B, e);
+		}
+
+		bool edge(int A, int B, TopoDS_Edge& e) {
+			e = edges_[{A, B}];
+			return true;
+		}
+
+		bool wire(const IfcSchema::IfcPolyLoop* loop, TopoDS_Wire& wire) {
+			BRep_Builder builder;
+			builder.MakeWire(wire);
+			bool valid;
+			auto ps = loop->Polygon();
+			loop_(ps, [this, &builder, &wire, &valid](int A, int B, bool fwd) {
+				TopoDS_Edge e;
+				if (edge(A, B, e)) {
+					if (!fwd) {
+						e.Reverse();
+					}
+					builder.Add(wire, e);
+					valid = true;
+				}
+			});
+			if (valid) {
+				wire.Closed(true);
+			}
+			return valid;
+		}
+	};
+
 	double deflection_tolerance;
 	double wire_creation_tolerance;
 	double point_equality_tolerance;
@@ -115,6 +193,7 @@ private:
 	double ifc_planeangle_unit;
 	double modelling_precision;
 	double dimensionality;
+	faceset_helper* faceset_helper_;
 
 #ifndef NO_CACHE
 	MAKE_TYPE_NAME(Cache) cache;
@@ -139,6 +218,7 @@ public:
 		, modelling_precision(0.00001)
 		, dimensionality(1.)
 		, placement_rel_to(0)
+		, faceset_helper_(nullptr)
 	{}
 
 	MAKE_TYPE_NAME(Kernel)(const MAKE_TYPE_NAME(Kernel)& other) : IfcGeom::Kernel(0) {

src/ifcgeom/IfcGeomElement.h

@@ -171,6 +171,13 @@ namespace IfcGeom {
 			: Element<P, PP>(geometry->settings() ,id, parent_id, name, type, guid, context, trsf, product)
 			, _geometry(geometry)
 		{}
+
+		bool calculate_projected_surface_area(double& along_x, double& along_y, double& along_z) const {
+			const auto& trsf = this->transformation().data();
+			const gp_Mat& mat = trsf.HVectorialPart();
+			gp_Ax3 ax(trsf.TranslationPart(), mat.Column(3), mat.Column(1));
+			return geometry().calculate_projected_surface_area(ax, along_x, along_y, along_z);
+		}
 	private:
 		BRepElement(const BRepElement& other);
 		BRepElement& operator=(const BRepElement& other);		

src/ifcgeom/IfcGeomFaces.cpp

@@ -178,35 +178,45 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcFace* l, TopoDS_Shape& face) {
 			if (is_interior == !process_interior) continue;
 		
 			TopoDS_Wire wire;
-			if (!convert_wire(loop, wire)) {
+			if (faceset_helper_ && loop->as<IfcSchema::IfcPolyLoop>()) {
+				faceset_helper_->wire(loop->as<IfcSchema::IfcPolyLoop>(), wire);
+			} else if (!convert_wire(loop, wire)) {
 				Logger::Message(Logger::LOG_ERROR, "Failed to process face boundary loop", loop);
 				delete mf;
 				return false;
 			}
 			
-			/*
+			// The approach below does not result in a significant speed-up
-			The approach below does not result in a significant speed-up
+			if (loop->as<IfcSchema::IfcPolyLoop>() && processed == 0 && face_surface.IsNull()) {
-			if (loop->declaration().is(IfcSchema::IfcPolyLoop::Class()) && processed == 0 && face_surface.IsNull()) {
+				TopExp_Explorer exp(wire, TopAbs_EDGE);
-				IfcSchema::IfcPolyLoop* polyloop = (IfcSchema::IfcPolyLoop*) loop;
+				int count = 0;
-				IfcSchema::IfcCartesianPoint::list::ptr points = polyloop->Polygon();
+				TopoDS_Edge edges[2];
-
+				for (; exp.More(); exp.Next(), count++) {
-				if (points->size() == 3) {
+					if (count < 2) {
-					// Help Open Cascade by finding the plane more efficiently
+						edges[count] = TopoDS::Edge(exp.Current());
-					IfcSchema::IfcCartesianPoint::list::it point_iterator = points->begin();
-					gp_Pnt a, b, c;
-					convert(*point_iterator++, a);
-					convert(*point_iterator++, b);
-					convert(*point_iterator++, c);
-					const gp_XYZ ab = (b.XYZ() - a.XYZ());
-					const gp_XYZ ac = (c.XYZ() - a.XYZ());
-					const gp_Vec cross = ab.Crossed(ac);
-					if (cross.SquareMagnitude() > ALMOST_ZERO) {
-						const gp_Dir n = cross;
-						face_surface = new Geom_Plane(a, n);
 					}
 				}
+				
+				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;
+						face_surface = new Geom_Plane(c1->Position().Location(), n);
+					}
+				} else {
+					gp_Pln pln;
+					approximate_plane_through_wire(wire, pln);
+					face_surface = new Geom_Plane(pln);
+				}
 			}
-			*/
 		
 			if (!same_sense) {
 				wire.Reverse();
@@ -297,16 +307,16 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcFace* l, TopoDS_Shape& face) {
 						TopTools_ListOfShape face_list;
 						triangulate_wire(wire, face_list);
 
-						TopoDS_Compound compound;
+						TopoDS_Compound triangulation_compound;
-						BRep_Builder builder;
+						BRep_Builder triangulation_builder;
-						builder.MakeCompound(compound);
+						triangulation_builder.MakeCompound(triangulation_compound);
 
 						TopTools_ListIteratorOfListOfShape face_iterator;
 						for (face_iterator.Initialize(face_list); face_iterator.More(); face_iterator.Next()) {
-							builder.Add(compound, face_iterator.Value());
+							triangulation_builder.Add(triangulation_compound, face_iterator.Value());
 						}
 
-						face = compound;
+						face = triangulation_compound;
 
 						return true;
 					}

src/ifcgeom/IfcGeomFunctions.cpp

@@ -297,30 +297,6 @@ namespace {
 		return M;
 	}
 
-	bool is_manifold(const TopoDS_Shape& a) {
-		TopTools_IndexedDataMapOfShapeListOfShape map;
-		TopExp::MapShapesAndAncestors(a, TopAbs_EDGE, TopAbs_FACE, map);
-
-		for (int i = 1; i <= map.Extent(); ++i) {
-			if (map.FindFromIndex(i).Extent() != 2) {
-				return false;
-			}
-		}
-
-		return true;
-	}
-
-	bool is_manifold(const TopTools_ListOfShape& l) {
-		TopTools_ListOfShape r;
-		TopTools_ListIteratorOfListOfShape it(l);
-		for (; it.More(); it.Next()) {
-			if (!is_manifold(it.Value())) {
-				return false;
-			}
-		}
-		return true;
-	}
-
 	void bounding_box_overlap(double p, const TopoDS_Shape& a, const TopTools_ListOfShape& b, TopTools_ListOfShape& c) {
 		Bnd_Box A;
 		BRepBndLib::Add(a, A);
@@ -3081,8 +3057,6 @@ bool IfcGeom::Kernel::wire_intersections(const TopoDS_Wire& wire, TopTools_ListO
 			// Only check non-consecutive edges
 			if (i == n - 1 && j == 0) continue;
 			
-			bool unbounded_intersects;
-
 			double u11, u12, u21, u22, U1, U2;
 			GeomAPI_ExtremaCurveCurve ecc(
 				BRep_Tool::Curve(wd->Edge(i + 1), u11, u12),
@@ -3090,87 +3064,89 @@ bool IfcGeom::Kernel::wire_intersections(const TopoDS_Wire& wire, TopTools_ListO
 			);
 
 			// @todo: extend this to work in case of multiple extrema and curved segments.
-			if ((unbounded_intersects = (ecc.NbExtrema() == 1 && ecc.Distance(1) < eps))) {
+			const bool unbounded_intersects = (ecc.NbExtrema() == 1 && ecc.Distance(1) < eps);
+			if (unbounded_intersects) {
 				ecc.Parameters(1, U1, U2);
-			}
 
-			if (u11 > u12) {
+				if (u11 > u12) {
-				std::swap(u11, u12);
+					std::swap(u11, u12);
-			}
+				}
-			if (u21 > u22) {
+				if (u21 > u22) {
-				std::swap(u21, u22);
+					std::swap(u21, u22);
-			}
+				}
-
+
-			/// @todo: tfk: probably need different thresholds on non-linear curves
+				/// @todo: tfk: probably need different thresholds on non-linear curves
-			u11 -= eps;
+				u11 -= eps;
-			u12 += eps;
+				u12 += eps;
-			u21 -= eps;
+				u21 -= eps;
-			u22 += eps;
+				u22 += eps;
-
+
-			// tfk: code below is for ShapeAnalysis_Wire::CheckIntersectingEdges()
+				// tfk: code below is for ShapeAnalysis_Wire::CheckIntersectingEdges()
-			// IntRes2d_SequenceOfIntersectionPoint points2d;
+				// IntRes2d_SequenceOfIntersectionPoint points2d;
-			// TColgp_SequenceOfPnt points3d;
+				// TColgp_SequenceOfPnt points3d;
-			// TColStd_SequenceOfReal errors;
+				// TColStd_SequenceOfReal errors;
-			// if (saw.CheckIntersectingEdges(i + 1, j + 1, points2d, points3d, errors)) {
+				// if (saw.CheckIntersectingEdges(i + 1, j + 1, points2d, points3d, errors)) {
-
+
-			if (unbounded_intersects &&	u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) {
+				if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) {
-
+
-				intersected = true;
+					intersected = true;
-
+
-				// Explore a forward and backward cycle from the intersection point
+					// Explore a forward and backward cycle from the intersection point
-				for (int fb = 0; fb <= 1; ++fb) {
+					for (int fb = 0; fb <= 1; ++fb) {
-					const bool forward = fb == 0;
+						const bool forward = fb == 0;
-
+
-					BRepBuilderAPI_MakeWire mw;
+						BRepBuilderAPI_MakeWire mw;
-					bool first = true;
+						bool first = true;
-
+
-					for (bounded_int k(j, n);;) {
+						for (bounded_int k(j, n);;) {
-						bool intersecting = k == j || k == i;
+							bool intersecting = k == j || k == i;
-						if (intersecting) {
+							if (intersecting) {
-							TopoDS_Edge e = wd->Edge(k + 1);
+								TopoDS_Edge e = wd->Edge(k + 1);
-							
+
-							TopoDS_Vertex v1, v2;
+								TopoDS_Vertex v1, v2;
-							TopExp::Vertices(e, v1, v2);
+								TopExp::Vertices(e, v1, v2);
-							const TopoDS_Vertex* v = first == forward ? &v2 : &v1;
+								const TopoDS_Vertex* v = first == forward ? &v2 : &v1;
-
+
-							// gp_Pnt p2 = points3d.Value(1);
+								// gp_Pnt p2 = points3d.Value(1);
-							
+
-							gp_Pnt p1 = BRep_Tool::Pnt(*v);
+								gp_Pnt p1 = BRep_Tool::Pnt(*v);
-							gp_Pnt pp1, pp2;
+								gp_Pnt pp1, pp2;
-							ecc.Points(1, pp1, pp2);
+								ecc.Points(1, pp1, pp2);
-							const gp_Pnt& p2 = k == i ? pp1 : pp2;
+								const gp_Pnt& p2 = k == i ? pp1 : pp2;
-							
+
-							// Substitute with a new edge from/to the intersection point
+								// Substitute with a new edge from/to the intersection point
-							if (p1.Distance(p2) > getValue(GV_PRECISION) * 2) {
+								if (p1.Distance(p2) > getValue(GV_PRECISION) * 2) {
-								double _, __;
+									double _, __;
-								Handle_Geom_Curve crv = BRep_Tool::Curve(e, _, __);
+									Handle_Geom_Curve crv = BRep_Tool::Curve(e, _, __);
-								BRepBuilderAPI_MakeEdge me(crv, p1, p2);
+									BRepBuilderAPI_MakeEdge me(crv, p1, p2);
-								TopoDS_Edge ed = me.Edge();
+									TopoDS_Edge ed = me.Edge();
-								mw.Add(ed);
+									mw.Add(ed);
-							}
+								}
-
+
-							first = false;
+								first = false;
-						} else {
+							} else {
-							// Re-use original edge
+								// Re-use original edge
-							mw.Add(wd->Edge(k+1));
+								mw.Add(wd->Edge(k + 1));
-						}
+							}
-
+
-						if (k == i) {
+							if (k == i) {
-							break;
+								break;
-						}
+							}
-
+
-						if (forward) {
+							if (forward) {
-							++k;
+								++k;
-						} else {
+							} else {
-							--k;
+								--k;
-						}
+							}
-					}
+						}
-
+
-					// Recursively process both cuts
+						// Recursively process both cuts
-					wire_intersections(mw.Wire(), wires);
+						wire_intersections(mw.Wire(), wires);
+					}
+
+					return true;
 				}
 
-				return true;
 			}
 		}
 	}
@@ -3473,3 +3449,78 @@ bool IfcGeom::Kernel::boolean_operation(const TopoDS_Shape& a, const TopoDS_Shap
 	return boolean_operation(a, bs, op, result, fuzziness);
 }
 #endif
+
+namespace {
+	void find_neighbours(IfcGeom::impl::tree<int>& tree, std::vector<gp_Pnt>& pnts, std::set<int>& visited, int p, double eps) {
+		visited.insert(p);
+
+		Bnd_Box b;
+		b.Set(pnts[p]);
+		b.Enlarge(eps);
+
+		std::vector<int> js = tree.select_box(b, false);
+		for (int j : js) {
+			if (visited.find(j) == visited.end()) {
+				find_neighbours(tree, pnts, visited, j, eps);
+			}
+		}
+	}
+}
+
+IfcGeom::Kernel::faceset_helper::~faceset_helper() {
+	kernel_->faceset_helper_ = nullptr;
+}
+
+IfcGeom::Kernel::faceset_helper::faceset_helper(Kernel* kernel, const IfcSchema::IfcConnectedFaceSet* l)
+	: kernel_(kernel) 
+{
+	kernel->faceset_helper_ = this;
+
+	IfcSchema::IfcCartesianPoint::list::ptr points = IfcParse::traverse((IfcUtil::IfcBaseClass*) l)->as<IfcSchema::IfcCartesianPoint>();
+	std::vector<gp_Pnt> pnts(std::distance(points->begin(), points->end()));
+	std::vector<TopoDS_Vertex> vertices(pnts.size());
+
+	BRep_Builder B;
+
+	const double eps = kernel->getValue(GV_PRECISION);
+	IfcGeom::impl::tree<int> tree;
+	{
+		int i = 0;
+		for (auto& pt : *points) {
+			if (kernel->convert(pt, pnts[i])) {
+				B.MakeVertex(vertices[i], pnts[i], Precision::Confusion());
+				tree.add(i, vertices[i]);
+				i++;
+			}
+		}
+	}
+
+	std::map<std::pair<int, int>, int> edge_use;
+
+	for (int i = 0; i < pnts.size(); ++i) {
+		std::set<int> vs;
+		find_neighbours(tree, pnts, vs, i, eps);
+
+		for (int v : vs) {
+			if (v <= i) {
+				auto pt = *(points->begin() + v);
+				vertex_mapping_.insert({pt->data().id(), i});
+			}
+		}
+	}
+
+	IfcSchema::IfcPolyLoop::list::ptr loops = IfcParse::traverse((IfcUtil::IfcBaseClass*)l)->as<IfcSchema::IfcPolyLoop>();
+
+	for (auto& loop : *loops) {
+		auto ps = loop->Polygon();
+		loop_(ps, [&edge_use](int C, int D, bool) {
+			edge_use[{C, D}] ++;
+		});
+	}
+
+	for (auto& p : edge_use) {
+		int a, b;
+		std::tie(a, b) = p.first;
+		edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]);
+	}
+}

src/ifcgeom/IfcGeomRepresentation.cpp

@@ -22,6 +22,9 @@
 #include <BRep_Builder.hxx>
 
 #include <TopoDS_Compound.hxx>
+#include <Geom_Plane.hxx>
+#include <GProp_GProps.hxx>
+#include <BRepGProp.hxx>
 
 #include "../ifcgeom/IfcGeom.h"
 
@@ -98,3 +101,135 @@ TopoDS_Compound IfcGeom::Representation::BRep::as_compound() const {
 	}
 	return compound;
 }
+
+namespace {
+	void accumulate(const gp_Ax3& ax, const gp_Dir& normal, double area, double& along_x, double& along_y, double& along_z) {
+		along_x += area * ax.XDirection().Dot(normal);
+		along_y += area * ax.YDirection().Dot(normal);
+		along_z += area * ax.Direction().Dot(normal);
+	}
+
+	void surface_area_along_direction(double tol, const TopoDS_Shape& s, const gp_Ax3& ax, double& along_x, double& along_y, double& along_z) {
+		along_x = along_y = along_z = 0.;
+
+		bool meshed = false;
+
+		// todo check whether manifold and divide by 2
+
+		TopExp_Explorer exp(s, TopAbs_FACE);
+		for (; exp.More(); exp.Next()) {
+			const TopoDS_Face& face = TopoDS::Face(exp.Current());
+			Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
+			Handle(Geom_Plane) plane = Handle(Geom_Plane)::DownCast(surf);
+
+			if (surf->DynamicType() == STANDARD_TYPE(Geom_Plane)) {
+				GProp_GProps prop_area;
+				BRepGProp::SurfaceProperties(face, prop_area);
+				const double area = prop_area.Mass();
+
+				accumulate(ax, plane->Position().Direction(), area, along_x, along_y, along_z);
+			} else {
+
+				if (!meshed) {
+					try {
+						BRepMesh_IncrementalMesh(s, tol);
+					} catch (...) {
+						Logger::Message(Logger::LOG_ERROR, "Failed to triangulate shape");
+						return;
+					}
+					meshed = true;
+				}
+
+				TopLoc_Location loc;
+				Handle(Poly_Triangulation) tri = BRep_Tool::Triangulation(face, loc);
+				if (!tri.IsNull()) {
+					const TColgp_Array1OfPnt& nodes = tri->Nodes();
+					std::vector<gp_XYZ> coords;
+					coords.reserve(nodes.Length());
+
+					for (int i = 1; i <= nodes.Length(); ++i) {
+						coords.push_back(nodes(i).Transformed(loc).XYZ());
+					}
+
+					const Poly_Array1OfTriangle& triangles = tri->Triangles();
+					for (int i = 1; i <= triangles.Length(); ++i) {
+						int n1, n2, n3;
+
+						if (face.Orientation() == TopAbs_REVERSED) {
+							triangles(i).Get(n3, n2, n1);
+						} else {
+							triangles(i).Get(n1, n2, n3);
+						}
+
+						const gp_XYZ& pt1 = coords[n1 - 1];
+						const gp_XYZ& pt2 = coords[n2 - 1];
+						const gp_XYZ& pt3 = coords[n3 - 1];
+						const gp_Vec v1 = pt2 - pt1;
+						const gp_Vec v2 = pt3 - pt2;
+						const gp_Vec v3 = pt1 - pt3;
+						gp_Dir normal = gp_Dir(v1^v2);
+
+						double edge_lengths[3] = { v1.Magnitude(), v2.Magnitude(), v3.Magnitude() };
+						std::sort(&edge_lengths[0], &edge_lengths[2]);
+
+						const double& a = edge_lengths[0];
+						const double& b = edge_lengths[1];
+						const double& c = edge_lengths[2];
+
+						const double area = 0.25 * sqrt((a + (b + c))*(c - (a - b))*(c + (a - b))*(a + (b - c)));
+						accumulate(ax, normal, area, along_x, along_y, along_z);
+					}
+				}
+			}
+		}
+	}
+}
+
+bool IfcGeom::Representation::BRep::calculate_surface_area(double& area) const {
+	area = 0.;
+	
+	for (IfcGeom::IfcRepresentationShapeItems::const_iterator it = begin(); it != end(); ++it) {
+		GProp_GProps prop;
+		BRepGProp::SurfaceProperties(it->Shape(), prop);
+		area += prop.Mass();
+	}
+
+	return true;
+}
+
+bool IfcGeom::Representation::BRep::calculate_volume(double& volume) const {
+	volume = 0.;
+
+	for (IfcGeom::IfcRepresentationShapeItems::const_iterator it = begin(); it != end(); ++it) {
+		if (Kernel::is_manifold(it->Shape())) {
+			GProp_GProps prop;
+			BRepGProp::VolumeProperties(it->Shape(), prop);
+			volume += prop.Mass();
+		} else {
+			return false;
+		}
+	}
+	
+	return true;
+}
+
+bool IfcGeom::Representation::BRep::calculate_projected_surface_area(const gp_Ax3 & ax, double & along_x, double & along_y, double & along_z) const {
+	along_x = along_y = along_z = 0.;
+
+	for (IfcGeom::IfcRepresentationShapeItems::const_iterator it = begin(); it != end(); ++it) {
+		double x, y, z;
+		surface_area_along_direction(settings().deflection_tolerance(), it->Shape(), ax, x, y, z);
+
+		if (Kernel::is_manifold(it->Shape())) {
+			x /= 2.;
+			y /= 2.;
+			z /= 2.;
+		}
+
+		along_x += x;
+		along_y += y;
+		along_z += z;
+	}
+
+	return true;	
+}

src/ifcgeom/IfcGeomRepresentation.h

@@ -76,6 +76,10 @@ namespace IfcGeom {
 			const IfcGeom::IfcRepresentationShapeItems& shapes() const { return shapes_; }
 			const std::string& id() const { return id_; }
 			TopoDS_Compound as_compound() const;
+
+			bool calculate_volume(double&) const;
+			bool calculate_surface_area(double&) const;
+			bool calculate_projected_surface_area(const gp_Ax3& ax, double& along_x, double& along_y, double& along_z) const;
 		};
 
 		class IFC_GEOM_API Serialization : public Representation  {

src/ifcgeom/IfcGeomShapes.cpp

@@ -103,6 +103,8 @@
 
 #include "../ifcgeom/IfcGeom.h"
 
+#include <memory>
+
 #define Kernel MAKE_TYPE_NAME(Kernel)
 
 bool IfcGeom::Kernel::convert(const IfcSchema::IfcExtrudedAreaSolid* l, TopoDS_Shape& shape) {
@@ -592,6 +594,12 @@ 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;
+
+	if (getValue(GV_MAX_FACES_TO_SEW) != -1) {
+		helper_scope.reset(new faceset_helper(this, l));
+	}
+
 	IfcSchema::IfcFace::list::ptr faces = l->CfsFaces();
 
 	TopTools_ListOfShape face_list;

src/ifcgeom_schema_agnostic/Kernel.cpp

@@ -1,5 +1,7 @@
 #include "Kernel.h"
 
+#include <TopExp.hxx>
+
 IfcGeom::Kernel::Kernel(IfcParse::IfcFile* file) {
 	if (file != 0) {
 		if (file->schema() == 0) {
@@ -172,3 +174,15 @@ std::map<std::string, IfcUtil::IfcBaseEntity*> IfcGeom::Kernel::get_layers(IfcUt
 	}
 }
 
+bool IfcGeom::Kernel::is_manifold(const TopoDS_Shape& a) {
+	TopTools_IndexedDataMapOfShapeListOfShape map;
+	TopExp::MapShapesAndAncestors(a, TopAbs_EDGE, TopAbs_FACE, map);
+
+	for (int i = 1; i <= map.Extent(); ++i) {
+		if (map.FindFromIndex(i).Extent() != 2) {
+			return false;
+		}
+	}
+
+	return true;
+}

src/ifcgeom_schema_agnostic/Kernel.h

@@ -85,6 +85,7 @@ namespace IfcGeom {
 
 		static int count(const TopoDS_Shape&, TopAbs_ShapeEnum);
 
+		static bool is_manifold(const TopoDS_Shape& a);
 		static IfcUtil::IfcBaseEntity* get_decomposing_entity(IfcUtil::IfcBaseEntity*);
 		static std::map<std::string, IfcUtil::IfcBaseEntity*> get_layers(IfcUtil::IfcBaseEntity*);
 	};

src/ifcgeomserver/IfcGeomServer.cpp

@@ -51,7 +51,11 @@
 #include <BRepGProp.hxx>
 #include <Geom_Plane.hxx>
 
+#include <boost/property_tree/ptree.hpp>
+#include <boost/property_tree/json_parser.hpp>
+
 using namespace boost;
+using boost::property_tree::ptree;
 
 template <typename T>
 union data_field {
@@ -358,15 +362,16 @@ protected:
 		swrite(s, value_);
 	}
 public:
-	Setting(uint32_t k = 0, uint32_t v = 0) : Command(DEFLECTION), id_(k), value_(v) {};
+	Setting(uint32_t k = 0, uint32_t v = 0) : Command(SETTING), id_(k), value_(v) {};
 	uint32_t id() const { return id_; }
 	uint32_t value() const { return value_; }
 };
 
 static const std::string TOTAL_SURFACE_AREA = "TOTAL_SURFACE_AREA";
 static const std::string TOTAL_SHAPE_VOLUME = "TOTAL_SHAPE_VOLUME";
-static const std::string WALKABLE_SURFACE_AREA = "WALKABLE_SURFACE_AREA";
+static const std::string SURFACE_AREA_ALONG_X = "SURFACE_AREA_ALONG_X";
-static const double MAX_WALKABLE_SURFACE_ANGLE_DEGREES = 15.;
+static const std::string SURFACE_AREA_ALONG_Y = "SURFACE_AREA_ALONG_Y";
+static const std::string SURFACE_AREA_ALONG_Z = "SURFACE_AREA_ALONG_Z";
 
 class QuantityWriter : public EntityExtension {
 private:
@@ -376,73 +381,25 @@ public:
 		elem_(elem)
 	{}
 	void write_contents(std::ostream& s) {
+		ptree pt;
+		double a, b, c;
 
-		double total_surface_area = 0.;
+		if (elem_->geometry().calculate_surface_area(a)) {
-		double total_shape_volume = 0.;
+			pt.put(TOTAL_SURFACE_AREA, a);
-		double walkable_surface_area = 0.;
-
-		TopoDS_Shape moved_shape = elem_->geometry().as_compound();
-			
-		{
-			GProp_GProps prop_area;
-			BRepGProp::SurfaceProperties(moved_shape, prop_area);
-			total_surface_area += prop_area.Mass();
 		}
 
-		{
+		if (elem_->geometry().calculate_volume(a)) {
-			GProp_GProps prop_volume;
+			pt.put(TOTAL_SHAPE_VOLUME, a);
-			BRepGProp::VolumeProperties(moved_shape, prop_volume);
-			total_shape_volume += prop_volume.Mass();
 		}
 
-		if (elem_->type() == "IfcSpace") {
+		if (elem_->calculate_projected_surface_area(a, b, c)) {
-			TopExp_Explorer exp(moved_shape, TopAbs_FACE);
+			pt.put(SURFACE_AREA_ALONG_X, a);
-			for (; exp.More(); exp.Next()) {
+			pt.put(SURFACE_AREA_ALONG_X, b);
-				const TopoDS_Face& face = TopoDS::Face(exp.Current());
+			pt.put(SURFACE_AREA_ALONG_X, c);
-				Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
-
-				// Assume we can only walk on planar surfaces
-				if (surf->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
-					continue;
-				}
-
-				BRepGProp_Face prop(face);
-				double u0, u1, v0, v1;
-				BRepTools::UVBounds(face, u0, u1, v0, v1);
-				gp_Pnt p;
-				gp_Vec normal_direction;
-				prop.Normal((u0 + u1) / 2., (v0 + v1) / 2., p, normal_direction);
-
-				gp_Vec normal(0., 0., 0.);
-				if (normal_direction.Magnitude() > 1.e-5) {
-					normal = gp_Dir(normal_direction.XYZ());
-				}
-
-				if (normal.Angle(gp::DZ()) < (MAX_WALKABLE_SURFACE_ANGLE_DEGREES * M_PI / 180.0)) {
-					GProp_GProps prop_face;
-					BRepGProp::SurfaceProperties(face, prop_face);
-					walkable_surface_area += prop_face.Mass();
-				}
-			}
 		}
 
-		// TODO: Manual JSON formatting is always a bad idea
 		std::ostringstream ss;
-		ss.write("{", 1);
+		boost::property_tree::write_json(ss, pt, false);
-		ss << format_json(TOTAL_SURFACE_AREA);
-		ss.write(":", 1);
-		ss << format_json(total_surface_area);
-		ss.write(",", 1);
-		ss << format_json(TOTAL_SHAPE_VOLUME);
-		ss.write(":", 1);
-		ss << format_json(total_shape_volume);
-		if (elem_->type() == "IfcSpace") {
-			ss.write(",", 1);
-			ss << format_json(WALKABLE_SURFACE_AREA);
-			ss.write(":", 1);
-			ss << format_json(walkable_surface_area);
-		}
-		ss.write("}", 1);
 
 		// We do a 4-byte manual alignment
 		std::string payload = ss.str();
@@ -461,6 +418,8 @@ int main () {
 	stdout_orig = std::cout.rdbuf();
 	std::cout.rdbuf(stdout_redir);
 
+	bool emit_quantities = false;
+
 #ifdef SET_BINARY_STREAMS
 	_setmode(_fileno(stdout), _O_BINARY);
 	std::cout.setf(std::ios_base::binary);
@@ -496,6 +455,11 @@ int main () {
 			std::vector< std::pair<uint32_t, uint32_t> >::const_iterator it = setting_pairs.begin();
 			for (; it != setting_pairs.end(); ++it) {
 				settings.set(it->first, it->second != 0);
+				if (it->first == IfcGeom::IteratorSettings::SEW_SHELLS && it->second) {
+					// Quantities (especially volume) can be emitted if there are proper
+					// topologically valid geometries being created.
+					emit_quantities = true;
+				}
 			}
 
 			settings.set_deflection_tolerance(deflection);
@@ -514,8 +478,11 @@ int main () {
 				break;
 			}
 			const IfcGeom::TriangulationElement<float, double>* geom = static_cast<const IfcGeom::TriangulationElement<float, double>*>(iterator->get());
-			QuantityWriter eext(iterator->get_native());
+			std::unique_ptr<QuantityWriter> eext;
-			Entity(geom, &eext).write(std::cout);
+			if (emit_quantities) {
+				eext = std::make_unique<QuantityWriter>(iterator->get_native());
+			}
+			Entity(geom, eext.get()).write(std::cout);
 			continue;
 		}
 		case NEXT: {