#include <map>

#include <TopoDS.hxx>
#include <TopExp.hxx>
#include <BRepGProp.hxx>
#include <GProp_GProps.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_Plane.hxx>
#include <BRepTools_WireExplorer.hxx>
#include <TopoDS_Compound.hxx>
#include <BRep_Builder.hxx>

#include "opencascade_conversion_result.h"

#include "../../../ifcparse/logger.h"
#include "../../../ifcgeom/representation.h"
#include "base_utils.h"
#include "boolean_utils.h"

#include <Standard_Version.hxx>

#include <iostream>
#include <vector>
#include <unordered_map>
#include <tuple>
#include <algorithm>

#if OCC_VERSION_HEX >= 0x70600
#include <TopTools_FormatVersion.hxx>
#endif

using ifcopenshell::geom::opaque_number;
using ifcopenshell::geom::opaque_coordinate;
using ifcopenshell::geom::conversion_result_shape;

namespace {
	// We bypass the conversion to gp_GTrsf, because it does not work
	void taxonomy_transform(const Eigen::Matrix4d* m, gp_XYZ& xyz) {
		if (m) {
			Eigen::Vector4d v(xyz.X(), xyz.Y(), xyz.Z(), 1.0);
			auto v2 = (*m * v).eval();
			xyz.ChangeData()[0] = v2(0);
			xyz.ChangeData()[1] = v2(1);
			xyz.ChangeData()[2] = v2(2);
		}
	}
}

ifcopenshell::geom::open_cascade_shape::open_cascade_shape(const TopoDS_Shape& shape)
	: shape_(shape) {}

ifcopenshell::geom::open_cascade_shape::open_cascade_shape(TopoDS_Shape&& shape)
	: shape_(std::move(shape)) {}

const TopoDS_Shape& ifcopenshell::geom::open_cascade_shape::shape() const {
	return shape_;
}

ifcopenshell::geom::open_cascade_shape::operator const TopoDS_Shape& () {
	return shape_;
}

std::string_view ifcopenshell::geom::open_cascade_shape::backend_id() const {
	return "opencascade";
}

ifcopenshell::geom::conversion_result_shape* ifcopenshell::geom::open_cascade_shape::clone() const {
	return new open_cascade_shape(shape_);
}

void ifcopenshell::geom::open_cascade_shape::triangulate(ifcopenshell::geom::settings settings, const ifcopenshell::geom::taxonomy::matrix4& place, ifcopenshell::geom::triangulation* t, int item_id, int surface_style_id, ifcopenshell::logger& logger) const {

	// @todo remove duplication with open_cascade_kernel::convert(const taxonomy::matrix4::ptr matrix, gp_GTrsf& trsf);
	// above can be static?

	// A 3x3 matrix to rotate the vertex normals
	std::optional<gp_Mat> rotation_matrix;

	if (place.components_) {
		const auto& m = *place.components_;
		rotation_matrix.emplace(
			m(0, 0), m(0, 1), m(0, 2),
			m(1, 0), m(1, 1), m(1, 2),
			m(2, 0), m(2, 1), m(2, 2)
		);
	}

	// When welding vertices, vertex coords will be shared among faces so we need to per-shape set
	// to keep track of which edges were already emitted.
	std::set<std::pair<int, int>> emitted_edges;

	// Do our own check if there are triangulations. Any will do. This is faster than the OCCT incremental check which compares the deflection tolerances and initialized a bunch of state
	bool has_triangulation = false;
	{
		TopExp_Explorer exp;
		for (exp.Init(shape_, TopAbs_FACE); exp.More(); exp.Next()) {
			TopLoc_Location loc;
			const Handle(Poly_Triangulation)& tri =
				BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), loc);
			if (tri) {
				has_triangulation = true;
				break;
			}
		}
	}
	if (!has_triangulation) {
		// triangulate the shape
		try {
			BRepMesh_IncrementalMesh(shape_, settings.get<settings::MesherLinearDeflection>().get(), false, settings.get<settings::MesherAngularDeflection>().get());
		} catch (...) {
			ifcopenshell::logger::root().message(ifcopenshell::logger::LOG_ERROR, "GEO", 183, "Failed to triangulate shape");
			return;
		}
	}

	// Iterates over the faces of the shape
	int num_faces = 0;
	TopExp_Explorer exp;
	for (exp.Init(shape_, TopAbs_FACE); exp.More(); exp.Next(), ++num_faces) {
		TopoDS_Face face = TopoDS::Face(exp.Current());

		size_t num_bounds = 0;
		for (TopoDS_Iterator it(face); it.More(); it.Next(), ++num_bounds) {}

		const bool is_planar = BRep_Tool::Surface(face) && BRep_Tool::Surface(face)->DynamicType() == STANDARD_TYPE(Geom_Plane);
		const bool has_inner_bounds = num_bounds > 1;

		const bool polyhedral_output_with_holes = settings.get<settings::TriangulationType>().get() == settings::POLYHEDRON_WITH_HOLES && is_planar;
		const bool polyhedral_output_without_holes = settings.get<settings::TriangulationType>().get() == settings::POLYHEDRON_WITHOUT_HOLES && is_planar && !has_inner_bounds;

		std::vector<std::tuple<int, int, int>> triangle_indices;

		TopLoc_Location loc;
		opencascade::handle<Poly_Triangulation> tri = BRep_Tool::Triangulation(face, loc);

		if (tri.IsNull()) {
			ifcopenshell::logger::root().message(ifcopenshell::logger::LOG_ERROR, "GEO", 184, "Triangulation missing for face");
		} else {
			// Keep track of the number of times an edge is used
			// Manifold edges (i.e. edges used twice) are deemed invisible
			std::map<std::pair<int, int>, int> edgecount;

			std::vector<gp_XYZ> coords;
			BRepGProp_Face prop(face);
			std::map<int, int> dict;

			// Vertex normals are only calculated if vertices are not welded and calculation is not disable explicitly.
			const bool calculate_normals = !settings.get<settings::WeldVertices>().get() &&
				!settings.get<settings::DontEmitNormals>().get();

			for (int i = 1; i <= tri->NbNodes(); ++i) {
				coords.push_back(tri->Node(i).Transformed(loc).XYZ());
				taxonomy_transform(place.components_, *coords.rbegin());
				const gp_XYZ& last = *coords.rbegin();
				dict[i] = t->addVertex(item_id, surface_style_id, last.X(), last.Y(), last.Z());

				if (calculate_normals) {
					const gp_Pnt2d& uv = tri->UVNode(i);
					gp_Pnt p;
					gp_Vec normal_direction;
					prop.Normal(uv.X(), uv.Y(), p, normal_direction);
					gp_Vec normal(0., 0., 0.);
					if (normal_direction.Magnitude() > 1.e-9) {
						if (rotation_matrix) {
							normal = gp_Dir(normal_direction.XYZ() * *rotation_matrix);
						} else {
							normal = normal_direction;
						}
					} else {
						opencascade::handle<Geom_Surface> surf = BRep_Tool::Surface(face);
						// Special case the normal at the poles of a spherical surface
						if (surf->DynamicType() == STANDARD_TYPE(Geom_SphericalSurface)) {
							if (fabs(fabs(uv.Y()) - M_PI / 2.) < 1.e-9) {
								const bool is_top = uv.Y() > 0;
								const bool is_forward = face.Orientation() == TopAbs_FORWARD;
								const double z = (is_top == is_forward) ? 1. : -1.;
								if (rotation_matrix) {
									normal = gp_Dir(gp_XYZ(0, 0, z) * *rotation_matrix);
								} else {
									normal = gp_Dir(gp_XYZ(0, 0, z));
								}
							}
						}
						// TODO: Do the same for conical surfaces, but they are rare in IFC.
					}
					t->addNormal(normal.X(), normal.Y(), normal.Z());
				}
			}

			for (int i = 1; i <= tri->NbTriangles(); ++i) {
				int n1, n2, n3;
				if (face.Orientation() == TopAbs_REVERSED)
					tri->Triangle(i).Get(n3, n2, n1);
				else tri->Triangle(i).Get(n1, n2, n3);

				if (dict[n1] == dict[n2] || dict[n2] == dict[n3] || dict[n3] == dict[n1]) {
					logger.warning("GEO", 185, "Mesher generated a degenerate triangle, ignoring");
					continue;
				}

				/* An alternative would be to calculate normals based
				* on the coordinates of the mesh vertices */
				/*
				const gp_XYZ pt1 = coords[n1-1];
				const gp_XYZ pt2 = coords[n2-1];
				const gp_XYZ pt3 = coords[n3-1];
				const gp_XYZ v1 = pt2-pt1;
				const gp_XYZ v2 = pt3-pt2;
				gp_Dir normal = gp_Dir(v1^v2);
				_normals.push_back((float)normal.X());
				_normals.push_back((float)normal.Y());
				_normals.push_back((float)normal.Z());
				*/

				if (polyhedral_output_without_holes || polyhedral_output_with_holes) {
					triangle_indices.push_back({ dict[n1], dict[n2], dict[n3] });
				} else {
					if (settings.get<settings::TriangulationType>().get() == settings::POLYHEDRON_WITHOUT_HOLES) {
						t->addFace(item_id, surface_style_id, std::vector<int>{ dict[n1], dict[n2], dict[n3] });
					} else if (settings.get<settings::TriangulationType>().get() == settings::POLYHEDRON_WITH_HOLES) {
						t->addFace(item_id, surface_style_id, std::vector<std::vector<int>>{{ dict[n1], dict[n2], dict[n3] }});
					} else {
						t->addFace(item_id, surface_style_id, dict[n1], dict[n2], dict[n3]);

						t->registerEdgeCount(dict[n1], dict[n2], edgecount);
						t->registerEdgeCount(dict[n2], dict[n3], edgecount);
						t->registerEdgeCount(dict[n3], dict[n1], edgecount);
					}
				}
			}
			for (auto& p : edgecount) {
				// @todo should be != 2?
				if (p.second == 1 && emitted_edges.find(p.first) == emitted_edges.end()) {
					// non manifold edge, face boundary
					t->registerEdge(item_id, p.first.first, p.first.second);
					if (settings.get<settings::WeldVertices>().get()) {
						// only relevant while welding, because otherwise vertices are not shared among distinct faces
						emitted_edges.insert(p.first);
					}
				}
			}
		}

		if (polyhedral_output_without_holes || polyhedral_output_with_holes) {
			auto loops = ifcopenshell::geom::util::find_boundary_loops(t->verts(), triangle_indices);
			if (polyhedral_output_without_holes) {
				if (!loops.empty() && !loops[0].empty()) {
					t->addFace(item_id, surface_style_id, loops[0]);
				}
			} else {
				if (!loops.empty()) {
					t->addFace(item_id, surface_style_id, loops);
				}
			}
		}
	}

	if (!t->normals().empty() && settings.get<settings::GenerateUvs>().get()) {
		t->uvs_ref() = ifcopenshell::geom::triangulation::box_project_uvs(t->verts(), t->normals());
	}

	if (num_faces == 0) {
		// Edges are only emitted if there are no faces. A mixed representation of faces
		// and loose edges is discouraged by the standard. An alternative would be to use
		// TopExp_Explorer texp(s, TopAbs_EDGE, TopAbs_FACE) to find edges that do not
		// belong to any face.

		NCollection_List<TopoDS_Shape> edges;
		// First collect edges part of wire in order
		for (TopExp_Explorer texp(shape_, TopAbs_WIRE); texp.More(); texp.Next()) {
			BRepTools_WireExplorer wexp(TopoDS::Wire(texp.Current()));
			for (; wexp.More(); wexp.Next()) {
				edges.Append(wexp.Current());
			}
		}

		// Then collect edges not part of wire
		for (TopExp_Explorer texp(shape_, TopAbs_EDGE, TopAbs_WIRE); texp.More(); texp.Next()) {
			edges.Append(texp.Current());
		}

		for (NCollection_List<TopoDS_Shape>::Iterator texp(edges); texp.More(); texp.Next()) {
			BRepAdaptor_Curve crv(TopoDS::Edge(texp.Value()));
			GCPnts_QuasiUniformDeflection tessellater(crv, settings.get<settings::MesherLinearDeflection>().get());
			int n = tessellater.NbPoints();
			int previous = -1;
			const bool reversed = texp.Value().Orientation() == TopAbs_REVERSED;
			bool first = true;

			gp_Pnt p0, p1;
			double u0 = std::numeric_limits<double>::quiet_NaN(), u1 = std::numeric_limits<double>::quiet_NaN();
			if (auto crv = BRep_Tool::Curve(TopoDS::Edge(texp.Value()), u0, u1)) {
				TopoDS_Vertex v0, v1;
				TopExp::Vertices(TopoDS::Edge(texp.Value()), v0, v1, false);
				if (!v0.IsNull() && !v1.IsNull()) {
					p0 = BRep_Tool::Pnt(v0);
					p1 = BRep_Tool::Pnt(v1);
				} else {
					u0 = u1 = std::numeric_limits<double>::quiet_NaN();
				}
			}

			for (int i = (reversed ? n : 1); reversed ? (i >= 1) : (i <= n); i += reversed ? -1 : 1) {
				gp_XYZ p;
				if (std::fabs(tessellater.Parameter(i) - u0) < 1.e-7) {
					// Use the exact points from the topology when parameter is close to the begin or end of the parametric range
					// This guarantees points are properly welded, because the GCPnts_QuasiUniformDeflection could otherwise introduce
					// minor differences between the approximated points from shared vertices.
					// @todo Using GCPnts_QuasiUniformDeflection on linear edges is pure lazyness
					p = p0.XYZ();
				} else if (std::fabs(tessellater.Parameter(i) - u1) < 1.e-7) {
					p = p1.XYZ();
				} else {
					p = tessellater.Value(i).XYZ();
				}

				auto p_local = p;
				taxonomy_transform(place.components_, p);

				int current = t->addVertex(item_id, surface_style_id, p.X(), p.Y(), p.Z());

				std::vector<std::pair<int, int>> segments;
				if (!first) {
					segments.push_back(std::make_pair(previous, current));
				}
				first = false;

				if (settings.get<settings::EdgeArrows>().get()) {
					// In case you want direction arrows on your edges
					double u = tessellater.Parameter(i);
					gp_XYZ p2, p3;
					gp_Pnt tmp;
					gp_Vec tmp2;
					crv.D1(u, tmp, tmp2);
					gp_Dir d1, d2, d3, d4;
					d1 = tmp2;
					if (reversed) {
						d1 = -d1;
					}
					if (fabs(d1.Z()) < 0.5) {
						d2 = d1.Crossed(gp::DZ());
					} else {
						d2 = d1.Crossed(gp::DY());
					}
					d3 = d1.XYZ() + d2.XYZ();
					d4 = d1.XYZ() - d2.XYZ();
					p2 = p_local - d3.XYZ() / 10.;
					p3 = p_local - d4.XYZ() / 10.;

					taxonomy_transform(place.components_, p2);
					taxonomy_transform(place.components_, p3);

					int left = t->addVertex(item_id, surface_style_id, p2.X(), p2.Y(), p2.Z());
					int right = t->addVertex(item_id, surface_style_id, p3.X(), p3.Y(), p3.Z());

					segments.push_back(std::make_pair(left, current));
					segments.push_back(std::make_pair(right, current));
				}

				for (auto& sgmt : segments) {
					t->addEdge(item_id, surface_style_id, sgmt.first, sgmt.second);
				}

				previous = current;
			}
		}
	}

	if (!settings.get<settings::OcctNoCleanTriangulation>().get()) {
		BRepTools::Clean(shape_);
	}
}

void ifcopenshell::geom::open_cascade_shape::serialize(const ifcopenshell::geom::taxonomy::matrix4& place, std::string& r) const {
	auto s = ifcopenshell::geom::util::apply_transformation(shape_, place);
	std::stringstream sstream;
#if OCC_VERSION_HEX >= 0x70600
	BRepTools::Write(s, sstream, false, false, TopTools_FormatVersion_VERSION_2);
#else
	BRepTools::Write(s, sstream);
#endif
	r = sstream.str();
}

int ifcopenshell::geom::open_cascade_shape::surface_genus() const {
	return ifcopenshell::geom::util::surface_genus(shape_);
}

bool ifcopenshell::geom::open_cascade_shape::is_manifold() const {
	return ifcopenshell::geom::util::is_manifold(shape_);
}

int ifcopenshell::geom::open_cascade_shape::num_vertices() const
{
	return ifcopenshell::geom::util::count(shape_, TopAbs_VERTEX);
}

int ifcopenshell::geom::open_cascade_shape::num_edges() const
{
	return ifcopenshell::geom::util::count(shape_, TopAbs_EDGE);
}

int ifcopenshell::geom::open_cascade_shape::num_faces() const
{
	return ifcopenshell::geom::util::count(shape_, TopAbs_FACE);
}

opaque_number ifcopenshell::geom::open_cascade_shape::open_cascade_shape::length()
{
	GProp_GProps prop;
	BRepGProp::LinearProperties(shape_, prop);
	double l = prop.Mass();
	return opaque_number(l);
}

opaque_number ifcopenshell::geom::open_cascade_shape::area()
{
	GProp_GProps prop;
	BRepGProp::SurfaceProperties(shape_, prop);
	double l = prop.Mass();
	return opaque_number(l);
}

opaque_number ifcopenshell::geom::open_cascade_shape::volume()
{
	GProp_GProps prop;
	BRepGProp::VolumeProperties(shape_, prop);
	double l = prop.Mass();
	return opaque_number(l);
}

#include <Geom_Plane.hxx>

opaque_coordinate<3> ifcopenshell::geom::open_cascade_shape::position()
{
	if (shape_.ShapeType() == TopAbs_FACE) {
		auto surf = BRep_Tool::Surface(TopoDS::Face(shape_));
		auto plane = Handle(Geom_Plane)::DownCast(surf);
		if (plane) {
			auto loc = plane->Location();
			return opaque_coordinate<3>(
				opaque_number(loc.X()),
				opaque_number(loc.Y()),
				opaque_number(loc.Z())
			);
		}
	}
	throw std::runtime_error("Invalid shape type");
}

opaque_coordinate<3> ifcopenshell::geom::open_cascade_shape::axis()
{
	if (shape_.ShapeType() == TopAbs_FACE) {
		auto surf = BRep_Tool::Surface(TopoDS::Face(shape_));
		auto plane = Handle(Geom_Plane)::DownCast(surf);
		if (plane) {
			auto dir = plane->Axis().Direction();
			return opaque_coordinate<3>(
				opaque_number(dir.X()),
				opaque_number(dir.Y()),
				opaque_number(dir.Z())
			);
		}
	}
	throw std::runtime_error("Invalid shape type");
}

opaque_coordinate<4> ifcopenshell::geom::open_cascade_shape::plane_equation()
{
	if (shape_.ShapeType() == TopAbs_FACE) {
		auto surf = BRep_Tool::Surface(TopoDS::Face(shape_));
		auto plane = Handle(Geom_Plane)::DownCast(surf);
		if (plane) {
			double a, b, c, d;
			plane->Pln().Coefficients(a, b, c, d);
			return opaque_coordinate<4>(
				opaque_number(a),
				opaque_number(b),
				opaque_number(c),
				opaque_number(d)
			);
		}
	}
	throw std::runtime_error("Invalid shape type");
}

std::vector<conversion_result_shape*> ifcopenshell::geom::open_cascade_shape::convex_decomposition()
{
	throw std::runtime_error("Not implemented");
}

conversion_result_shape * ifcopenshell::geom::open_cascade_shape::halfspaces()
{
	throw std::runtime_error("Not implemented");
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::solid()
{
	throw std::runtime_error("Not implemented");
}

conversion_result_shape * ifcopenshell::geom::open_cascade_shape::box()
{
	throw std::runtime_error("Not implemented");
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::wrap_in_compound()
{
	TopoDS_Compound compound;
	BRep_Builder builder;
	builder.MakeCompound(compound);
	builder.Add(compound, shape_);
	return new open_cascade_shape(std::move(compound));
}

std::vector<conversion_result_shape*> ifcopenshell::geom::open_cascade_shape::vertices()
{
    NCollection_IndexedMap<TopoDS_Shape, TopTools_ShapeMapHasher> map;
	TopExp::MapShapes(shape_, TopAbs_VERTEX, map);
	std::vector<conversion_result_shape*> vec;
	for (int i = 1; i <= map.Extent(); ++i) {
		vec.push_back(new open_cascade_shape(map.FindKey(i)));
	}
	return vec;
}

std::vector<conversion_result_shape*> ifcopenshell::geom::open_cascade_shape::edges()
{
    NCollection_IndexedMap<TopoDS_Shape, TopTools_ShapeMapHasher> map;
	TopExp::MapShapes(shape_, TopAbs_EDGE, map);
	std::vector<conversion_result_shape*> vec;
	for (int i = 1; i <= map.Extent(); ++i) {
		vec.push_back(new open_cascade_shape(map.FindKey(i)));
	}
	return vec;
}

std::vector<conversion_result_shape*> ifcopenshell::geom::open_cascade_shape::facets()
{
    NCollection_IndexedMap<TopoDS_Shape, TopTools_ShapeMapHasher> map;
	TopExp::MapShapes(shape_, TopAbs_FACE, map);
	std::vector<conversion_result_shape*> vec;
	for (int i = 1; i <= map.Extent(); ++i) {
		vec.push_back(new open_cascade_shape(map.FindKey(i)));
	}
	return vec;
}

namespace {
	conversion_result_shape* boolean_op(BOPAlgo_Operation op, const TopoDS_Shape& shape_, const TopoDS_Shape& other_shape) {
		ifcopenshell::geom::util::boolean_settings st;
		st.attempt_2d = true;
		st.debug = false;
		st.precision = 1.e-5;

		TopoDS_Shape result;
		if (ifcopenshell::geom::util::boolean_operation(st, shape_, other_shape, op, result)) {
			return new ifcopenshell::geom::open_cascade_shape(result);
		} else {
			throw std::runtime_error("Failed to process boolean operation");
		}
	}
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::add(conversion_result_shape* other)
{
	return boolean_op(BOPAlgo_FUSE, shape_, ((ifcopenshell::geom::open_cascade_shape*)other)->shape_);
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::subtract(conversion_result_shape* other)
{
	return boolean_op(BOPAlgo_CUT, shape_, ((ifcopenshell::geom::open_cascade_shape*)other)->shape_);
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::intersect(conversion_result_shape* other)
{
	return boolean_op(BOPAlgo_COMMON, shape_, ((ifcopenshell::geom::open_cascade_shape*)other)->shape_);
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::concat(conversion_result_shape* other)
{
	TopoDS_Compound compound;
	BRep_Builder builder;

	auto& left = shape_;
	auto& right = ((ifcopenshell::geom::open_cascade_shape*)other)->shape_;

	// This reads a bit strange, but we want to specifically avoid compounds of faces that are
	// the result of shell instances that are not sewn into a shell (yet).
	if (left.ShapeType() == TopAbs_COMPOUND && !ifcopenshell::geom::util::is_compound_of_faces(left)) {
		compound = TopoDS::Compound(left);
	} else {
		builder.MakeCompound(compound);
		builder.Add(compound, left);
	}

	builder.Add(compound, right);

	return new open_cascade_shape(std::move(compound));
}

std::pair<opaque_coordinate<3>, opaque_coordinate<3>> ifcopenshell::geom::open_cascade_shape::bounding_box() const
{
	throw std::runtime_error("Not implemented");
}

conversion_result_shape* ifcopenshell::geom::open_cascade_shape::moved(ifcopenshell::geom::taxonomy::matrix4::ptr t) const
{
	return new open_cascade_shape(ifcopenshell::geom::util::apply_transformation(shape_, *t));
}

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 * fabs(ax.XDirection().Dot(normal));
		along_y += area * fabs(ax.YDirection().Dot(normal));
		along_z += area * fabs(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;

		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 (...) {
						ifcopenshell::logger::root().message(ifcopenshell::logger::LOG_ERROR, "GEO", 186, "Failed to triangulate shape");
						return;
					}
					meshed = true;
				}

				TopLoc_Location loc;
				Handle(Poly_Triangulation) tri = BRep_Tool::Triangulation(face, loc);
				if (!tri.IsNull()) {
					std::vector<gp_XYZ> coords;
					coords.reserve(tri->NbNodes());

					for (int i = 1; i <= tri->NbNodes(); ++i) {
						coords.push_back(tri->Node(i).Transformed(loc).XYZ());
					}

					for (int i = 1; i <= tri->NbTriangles(); ++i) {
						int n1, n2, n3;

						if (face.Orientation() == TopAbs_REVERSED) {
							tri->Triangle(i).Get(n3, n2, n1);
						} else {
							tri->Triangle(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;
						const gp_Vec normal_vector = v1 ^ v2;
						if (normal_vector.Magnitude() > 1.e-7) {
							gp_Dir normal = gp_Dir();

							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 ifcopenshell::geom::open_cascade_shape::surface_area_along_direction(double tol, const ifcopenshell::geom::taxonomy::matrix4::ptr& place, double& along_x, double& along_y, double& along_z) const
{
	gp_GTrsf trsf;

	if (place->components_) {
		gp_Trsf tr;
		const auto& m = place->ccomponents();
		tr.SetValues(
			m(0, 0), m(0, 1), m(0, 2), m(0, 3),
			m(1, 0), m(1, 1), m(1, 2), m(1, 3),
			m(2, 0), m(2, 1), m(2, 2), m(2, 3)
		);
		trsf = tr;
	}

	gp_Mat mat = trsf.Trsf().HVectorialPart();
	gp_Ax3 ax(trsf.TranslationPart(), mat.Column(3), mat.Column(1));

	surface_area_along_direction_(tol, shape_, ax, along_x, along_y, along_z);

	return true;
}

std::size_t ifcopenshell::geom::open_cascade_shape::map(opaque_coordinate<4>&, opaque_coordinate<4>&) {
	throw std::runtime_error("Not implemented");
}

std::size_t ifcopenshell::geom::open_cascade_shape::map(const std::vector<opaque_coordinate<4>>&, const std::vector<opaque_coordinate<4>>&) {
	throw std::runtime_error("Not implemented");
}