diff --git a/src/ifcgeom/infra_sweep_helper.cpp b/src/ifcgeom/infra_sweep_helper.cpp index 3c65dd3d60..b6c0dc109b 100644 --- a/src/ifcgeom/infra_sweep_helper.cpp +++ b/src/ifcgeom/infra_sweep_helper.cpp @@ -2,7 +2,9 @@ #include "infra_sweep_helper.h" #include "function_item_evaluator.h" +#include #include +#include using namespace ifcopenshell::geom; @@ -72,6 +74,53 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett longitudes.push_back(x.dist_along); } longitudes.push_back(std::numeric_limits::infinity()); + + // Directrix frame (col0 = tangent, col1 = lateral, col2 = up, col3 = position) at + // every cross section station. Only needed to reconcile a cross section's own + // Axis / RefDirection against the curve; skip the work when no placement has one. + std::vector section_directrix_frames; + if (std::any_of(cross_sections.begin(), cross_sections.end(), + [](const cross_section& cs) { return cs.rotation.has_value(); })) { + section_directrix_frames.reserve(cross_sections.size()); + for (const auto& cs : cross_sections) { + section_directrix_frames.push_back(evaluator.evaluate(std::min(std::max(cs.dist_along, start), end))); + } + } + + // The frame a cross section is placed in, given the directrix frame at its station + // and its own IfcAxis2PlacementLinear: profile X, profile Y = Axis, profile normal + // = RefDirection -- or, when RefDirection was not authored, the curve tangent, so + // the section stays perpendicular to the path (buildingSMART IFC4.x-IF #147). When + // the placement carries no direction vectors the section just follows the curve + // (lateral, up, tangent). + const auto profile_basis = + [](const std::optional& rotation, + const std::optional& ref_direction, + const Eigen::Matrix4d& directrix_frame) -> Eigen::Matrix3d { + const Eigen::Vector3d tangent = directrix_frame.col(0).head<3>().normalized(); + const Eigen::Vector3d lateral = directrix_frame.col(1).head<3>().normalized(); + const Eigen::Vector3d up = directrix_frame.col(2).head<3>().normalized(); + Eigen::Matrix3d B; + if (!rotation) { + B.col(0) = lateral; + B.col(1) = up; + B.col(2) = tangent; + return B; + } + const Eigen::Vector3d axis = rotation->col(2).normalized(); + const Eigen::Vector3d normal = ref_direction ? ref_direction->normalized() : tangent; + Eigen::Vector3d x = axis.cross(normal); + if (x.norm() < 1.e-9) { + // Axis parallel to the normal: fall back to the curve's own lateral. + x = lateral - lateral.dot(axis) * axis; + } + x.normalize(); + B.col(0) = x; + B.col(1) = axis; + B.col(2) = x.cross(axis); + return B; + }; + auto profile_index = longitudes.begin(); for (size_t i = 0; i <= num_steps; ++i) { auto dist_along = start + delta_step * i; @@ -88,6 +137,7 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett const auto& profile_a = cross_sections[std::distance(longitudes.begin(), profile_index)].section_geometry; const auto& offset_a = cross_sections[std::distance(longitudes.begin(), profile_index)].offset; const auto& rotation_a = cross_sections[std::distance(longitudes.begin(), profile_index)].rotation; + const auto& ref_direction_a = cross_sections[std::distance(longitudes.begin(), profile_index)].ref_direction; taxonomy::geom_item::ptr interpolated = nullptr; @@ -98,26 +148,37 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett (profile_index + 1 < longitudes.end()) && (relative_dist_along >= 1.e-9 || offset_a.cwiseAbs().maxCoeff() > 0. || rotation_a); + // When both bracketing placements ask for the same orientation, drive the sweep + // frame from it directly (relative to the curve). When they disagree, keep the + // sweep frame on the curve and fold each section's own orientation into its + // profile points via section_basis_a / section_basis_b so the ends still land + // exactly as authored without twisting the body between them. std::optional interpolated_rotation; + std::optional interpolated_ref_direction; + Eigen::Matrix3d section_basis_a = Eigen::Matrix3d::Identity(); + Eigen::Matrix3d section_basis_b = Eigen::Matrix3d::Identity(); if (should_interpolate) { taxonomy::geom_item::ptr profile_b; Eigen::Vector3d offset_b; std::optional rotation_b; + std::optional ref_direction_b; if ((profile_index + 1 < longitudes.end())) { profile_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1].section_geometry; offset_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1].offset; rotation_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1].rotation; + ref_direction_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1].ref_direction; } else { profile_b = profile_a; offset_b = offset_a; rotation_b = rotation_a; + ref_direction_b = ref_direction_a; } // Only interpolate if the profiles are different or either of the offsets is non-zero bool should_interpolate2 = (profile_a->instance != profile_b->instance) || - (offset_a.cwiseAbs().maxCoeff() > 0. || offset_b.cwiseAbs().maxCoeff() > 0. || rotation_b); + (offset_a.cwiseAbs().maxCoeff() > 0. || offset_b.cwiseAbs().maxCoeff() > 0. || rotation_a || rotation_b); if (should_interpolate2) { @@ -160,12 +221,33 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett } auto interpolated_offset = lerp(offset_a, offset_b, relative_dist_along); - if (rotation_a == rotation_b && rotation_a) { - // @todo we don't support an overridden rotation on only one of the placements - // in which case we would need to lerp with the rotation component below in m4b. - interpolated_rotation = lerp(*rotation_a, *rotation_b, relative_dist_along); - } else if (rotation_a != rotation_b) { - logger.error("GEO", 42, "Direction vectors on cross section placements only supported when used consistently"); + if (rotation_a == rotation_b) { + // Same orientation on both placements (including both absent): the + // sweep frame carries it, built against the curve just below. + interpolated_rotation = rotation_a; + interpolated_ref_direction = ref_direction_a; + } else if (rotation_a || rotation_b) { + // The two placements disagree -- in practice they share an Axis but + // only one carries a RefDirection (a raked end against a square run). + // Drive the sweep frame from the shared Axis with the profile normal + // on the curve tangent -- identical to the neighbouring consistent + // segments, so m4b stays continuous across the boundary and the body + // never flips -- then fold each section's *own* authored orientation + // into its profile points through section_basis_a / section_basis_b, + // so each end cap still lands exactly as authored. + const auto ia = static_cast(std::distance(longitudes.begin(), profile_index)); + assert(ia + 1 < section_directrix_frames.size()); + const std::optional no_ref; + const auto base_a = profile_basis(rotation_a, no_ref, section_directrix_frames[ia]); + const auto base_b = profile_basis(rotation_b, no_ref, section_directrix_frames[ia + 1]); + section_basis_a = + base_a.transpose() * + profile_basis(rotation_a, ref_direction_a, section_directrix_frames[ia]); + section_basis_b = + base_b.transpose() * + profile_basis(rotation_b, ref_direction_b, section_directrix_frames[ia + 1]); + interpolated_rotation = rotation_a ? rotation_a : rotation_b; + interpolated_ref_direction = no_ref; } taxonomy::loop::ptr w1, w2; @@ -337,7 +419,9 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett const auto& tagged_point_on_w1 = tag_to_point_on_w1[t]; const auto& tagged_point_on_w2 = tag_to_point_on_w2[t]; - auto p3 = (lerp(tagged_point_on_w1->ccomponents(), tagged_point_on_w2->ccomponents(), relative_dist_along) + interpolated_offset).eval(); + const Eigen::Vector3d rebased_w1 = section_basis_a * tagged_point_on_w1->ccomponents(); + const Eigen::Vector3d rebased_w2 = section_basis_b * tagged_point_on_w2->ccomponents(); + auto p3 = (lerp(rebased_w1, rebased_w2, relative_dist_along) + interpolated_offset).eval(); std::set tags_for_this_point_on_subsequent_profile = {t}; @@ -357,7 +441,9 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett } else { for (auto tmp__ : boost::combine(w1_points, w2_points)) { boost::tie(p1, p2) = tmp__; - auto p3 = (lerp(p1->ccomponents(), p2->ccomponents(), relative_dist_along) + interpolated_offset).eval(); + const Eigen::Vector3d rebased_1 = section_basis_a * p1->ccomponents(); + const Eigen::Vector3d rebased_2 = section_basis_b * p2->ccomponents(); + auto p3 = (lerp(rebased_1, rebased_2, relative_dist_along) + interpolated_offset).eval(); points.push_back(taxonomy::make(p3)); } } @@ -399,17 +485,13 @@ taxonomy::loft::ptr ifcopenshell::geom::make_loft(const ifcopenshell::geom::sett std::wcout << "#" << pwf->instance.data().id() << " " << dist_along << ": " << m4.col(3).row(2).value() << std::endl; }*/ + // Sweep frame at this station: the profile orientation asked for by the + // (consistent) placements, built against the curve here so it follows the + // directrix. Falls back to the plain curve frame (lateral, up, tangent) when + // no placement carries direction vectors. Inconsistent placements keep this + // on the curve and are reconciled through section_basis_a / section_basis_b. Eigen::Matrix4d m4b = Eigen::Matrix4d::Identity(); - if (interpolated_rotation) { - // direction vectors on the linear placement overwrite the placement otherwise inferred from the tangent - m4b.col(0).head<3>() = interpolated_rotation->col(1); - m4b.col(1).head<3>() = interpolated_rotation->col(2); - m4b.col(2).head<3>() = interpolated_rotation->col(0); - } else { - m4b.col(0).head<3>() = m4.col(1).head<3>().normalized(); - m4b.col(1).head<3>() = m4.col(2).head<3>().normalized(); - m4b.col(2).head<3>() = m4.col(0).head<3>().normalized(); - } + m4b.block<3, 3>(0, 0) = profile_basis(interpolated_rotation, interpolated_ref_direction, m4); m4b.col(3).head<3>() = m4.col(3).head<3>(); if (interpolated) { diff --git a/src/ifcgeom/infra_sweep_helper.h b/src/ifcgeom/infra_sweep_helper.h index cb0c06ac39..709292a2ca 100644 --- a/src/ifcgeom/infra_sweep_helper.h +++ b/src/ifcgeom/infra_sweep_helper.h @@ -14,7 +14,13 @@ namespace ifcopenshell { double dist_along; taxonomy::geom_item::ptr section_geometry; Eigen::Vector3d offset; + // rotation: the IfcAxis2PlacementLinear basis [X | Y | Axis], or nullopt when + // the placement carries neither Axis nor RefDirection. std::optional rotation; + // ref_direction: the raw RefDirection, only when it was actually authored on the + // placement (rotation carries a default otherwise). When absent the profile + // normal comes from the directrix tangent -- see buildingSMART IFC4.x-IF #147. + std::optional ref_direction; bool operator <(const cross_section& other) const { return dist_along < other.dist_along; diff --git a/src/ifcgeom/mapping/IfcSectionedSolidHorizontal.cpp b/src/ifcgeom/mapping/IfcSectionedSolidHorizontal.cpp index 546c8712d6..e790b36a35 100644 --- a/src/ifcgeom/mapping/IfcSectionedSolidHorizontal.cpp +++ b/src/ifcgeom/mapping/IfcSectionedSolidHorizontal.cpp @@ -49,7 +49,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal& in // The longitudes determine the range of the sweep and the offsets are interpolated in between // sweep segments. std::vector profile_offsets; - std::vector> profile_rotations; + std::vector> profile_axes; + std::vector> profile_ref_directions; std::vector longitudes; for (auto& cs : css) { @@ -72,23 +73,27 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal& in profile_offsets.push_back(po); std::optional rot; + std::optional ref_direction; if (csp.Axis() && csp.RefDirection()) { + ref_direction = taxonomy::cast(map(csp.RefDirection()))->ccomponents(); rot = taxonomy::matrix4( Eigen::Vector3d(0, 0, 0), taxonomy::cast(map(csp.Axis()))->ccomponents(), - taxonomy::cast(map(csp.RefDirection()))->ccomponents()).ccomponents().block<3,3>(0,0); + *ref_direction).ccomponents().block<3,3>(0,0); } else if (csp.Axis()) { rot = taxonomy::matrix4( Eigen::Vector3d(0, 0, 0), taxonomy::cast(map(csp.Axis()))->ccomponents()).ccomponents().block<3, 3>(0, 0); - } else if (csp.RefDirection()) { - rot = taxonomy::matrix4( - Eigen::Vector3d(0, 0, 0), - Eigen::Vector3d(0, 0, 1), - taxonomy::cast(map(csp.RefDirection()))->ccomponents() + } else if (csp.RefDirection()) { + ref_direction = taxonomy::cast(map(csp.RefDirection()))->ccomponents(); + rot = taxonomy::matrix4( + Eigen::Vector3d(0, 0, 0), + Eigen::Vector3d(0, 0, 1), + *ref_direction ).ccomponents().block<3, 3>(0, 0); - } - profile_rotations.push_back(rot); + } + profile_axes.push_back(rot); + profile_ref_directions.push_back(ref_direction); } if (faces.size() != profile_offsets.size()) { logger_.warning("GEO", 286, "Expected CrossSections and CrossSectionPositions to be equal length, but got " + std::to_string(faces.size()) + " and " + std::to_string(profile_offsets.size()) + " respectively", inst); @@ -100,7 +105,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal& in } for (size_t i = 0; i < faces.size(); ++i) { - cross_sections.push_back({ longitudes[i], faces[i], profile_offsets[i], profile_rotations[i]}); + cross_sections.push_back({longitudes[i], faces[i], profile_offsets[i], profile_axes[i], profile_ref_directions[i]}); } #else return nullptr; diff --git a/src/ifcgeom/mapping/IfcSectionedSurface.cpp b/src/ifcgeom/mapping/IfcSectionedSurface.cpp index a0ef5f727b..32423d9f7e 100644 --- a/src/ifcgeom/mapping/IfcSectionedSurface.cpp +++ b/src/ifcgeom/mapping/IfcSectionedSurface.cpp @@ -50,7 +50,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSurface& inst) { // The longitudes determine the range of the sweep and the offsets are interpolated in between // sweep segments. std::vector profile_offsets; - std::vector> profile_rotations; + std::vector> profile_axes; + std::vector> profile_ref_directions; std::vector longitudes; for (auto& cs : css) { @@ -74,24 +75,28 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSurface& inst) { profile_offsets.push_back(po); std::optional rot; + std::optional ref_direction; if (csp.Axis() && csp.RefDirection()) { + ref_direction = taxonomy::cast(map(csp.RefDirection()))->ccomponents(); rot = taxonomy::matrix4( Eigen::Vector3d(0, 0, 0), taxonomy::cast(map(csp.Axis()))->ccomponents(), - taxonomy::cast(map(csp.RefDirection()))->ccomponents()).ccomponents().block<3, 3>(0, 0); + *ref_direction).ccomponents().block<3, 3>(0, 0); } else if (csp.Axis()) { rot = taxonomy::matrix4( Eigen::Vector3d(0, 0, 0), taxonomy::cast(map(csp.Axis()))->ccomponents()).ccomponents().block<3, 3>(0, 0); - } else if (csp.RefDirection()) { - rot = taxonomy::matrix4( - Eigen::Vector3d(0, 0, 0), - Eigen::Vector3d(0, 0, 1), - taxonomy::cast(map(csp.RefDirection()))->ccomponents()) - .ccomponents() - .block<3, 3>(0, 0); - } - profile_rotations.push_back(rot); + } else if (csp.RefDirection()) { + ref_direction = taxonomy::cast(map(csp.RefDirection()))->ccomponents(); + rot = taxonomy::matrix4( + Eigen::Vector3d(0, 0, 0), + Eigen::Vector3d(0, 0, 1), + *ref_direction) + .ccomponents() + .block<3, 3>(0, 0); + } + profile_axes.push_back(rot); + profile_ref_directions.push_back(ref_direction); } #else return nullptr; @@ -106,8 +111,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSurface& inst) { } for (size_t i = 0; i < faces.size(); ++i) { - cross_sections.push_back({ longitudes[i], faces[i], profile_offsets[i], profile_rotations[i] }); - } + cross_sections.push_back({longitudes[i], faces[i], profile_offsets[i], profile_axes[i], profile_ref_directions[i]}); + } } return make_loft(settings_, inst, fn, cross_sections); diff --git a/src/ifcgeom/tests/test_ifcopenshell_geometry.cpp b/src/ifcgeom/tests/test_ifcopenshell_geometry.cpp index 63f20538dc..7933bcb8cf 100644 --- a/src/ifcgeom/tests/test_ifcopenshell_geometry.cpp +++ b/src/ifcgeom/tests/test_ifcopenshell_geometry.cpp @@ -1,7 +1,9 @@ #include +#include #include #include +#include #include #include "ifcgeom/converter.h" @@ -148,8 +150,152 @@ const ifcopenshell::geom::geometry_conversion_task* task_for_product( return nullptr; } +// A minimal IfcSectionedSolidHorizontal (IFC4X3_ADD2) whose two +// IfcAxis2PlacementLinear cross section positions use direction vectors +// inconsistently: the near position carries a raked RefDirection and a +// 1/cos(theta) wider profile, the far position carries neither. Before the +// make_loft() fix this logged GEO 42 and dropped the rotation, lofting a wedge. +constexpr const char* RAKED_SECTIONED_SOLID_SPF = R"IFC(ISO-10303-21; +HEADER; +FILE_DESCRIPTION((''),'2;1'); +FILE_NAME('','',(''),(''),'','',''); +FILE_SCHEMA(('IFC4X3_ADD2')); +ENDSEC; +DATA; +#1=IFCPROJECT('0RYK8PV8D0ee9DDm77xcTZ',$,'T',$,$,$,$,(#6),$); +#2=IFCCARTESIANPOINT((0.,0.,0.)); +#3=IFCDIRECTION((0.,0.,1.)); +#4=IFCDIRECTION((1.,0.,0.)); +#5=IFCAXIS2PLACEMENT3D(#2,#3,#4); +#6=IFCGEOMETRICREPRESENTATIONCONTEXT($,'Model',3,1.E-05,#5,$); +#7=IFCGEOMETRICREPRESENTATIONSUBCONTEXT('Body','Model',*,*,*,*,#6,$,.MODEL_VIEW.,$); +#8=IFCCARTESIANPOINT((0.,0.,0.)); +#9=IFCCARTESIANPOINT((40.,0.,0.)); +#10=IFCPOLYLINE((#8,#9)); +#11=IFCDIRECTION((0.,0.,1.)); +#12=IFCDIRECTION((0.9034641832977311,-0.42866358545853134,0.)); +#13=IFCPOINTBYDISTANCEEXPRESSION(IFCLENGTHMEASURE(0.),$,$,$,#10); +#14=IFCAXIS2PLACEMENTLINEAR(#13,#11,#12); +#15=IFCPOINTBYDISTANCEEXPRESSION(IFCLENGTHMEASURE(40.),$,$,$,#10); +#16=IFCAXIS2PLACEMENTLINEAR(#15,#11,$); +#17=IFCCARTESIANPOINTLIST2D(((-0.9223756168081689,0.),(0.9223756168081689,0.),(0.9223756168081689,6.),(-0.9223756168081689,6.),(-0.9223756168081689,0.)),$); +#18=IFCINDEXEDPOLYCURVE(#17,$,.F.); +#19=IFCARBITRARYCLOSEDPROFILEDEF(.AREA.,$,#18); +#20=IFCCARTESIANPOINTLIST2D(((-0.8333333333333334,0.),(0.8333333333333334,0.),(0.8333333333333334,6.),(-0.8333333333333334,6.),(-0.8333333333333334,0.)),$); +#21=IFCINDEXEDPOLYCURVE(#20,$,.F.); +#22=IFCARBITRARYCLOSEDPROFILEDEF(.AREA.,$,#21); +#23=IFCSECTIONEDSOLIDHORIZONTAL(#10,(#19,#22),(#14,#16)); +#24=IFCBUILDINGELEMENTPROXY('3cWcr4$892GAeKwryhDILR',$,'wingwall',$,$,$,#26,$,$); +#25=IFCSHAPEREPRESENTATION(#7,'Body','AdvancedSweptSolid',(#23)); +#26=IFCPRODUCTDEFINITIONSHAPE($,$,(#25)); +ENDSEC; +END-ISO-10303-21; +)IFC"; + +struct sectioned_solid_result { + std::size_t geo42_count = 0; + bool produced_brep = false; + double projected_area_x = 0.0; + double projected_area_y = 0.0; + double projected_area_z = 0.0; +}; + +sectioned_solid_result convert_sectioned_solid(const std::string& spf) { + std::istringstream stream(spf); + ifcopenshell::logger log; + log.output_format(ifcopenshell::logger::FMT_INMEMORY); + ifcopenshell::file file(stream, static_cast(spf.size()), log); + REQUIRE(file.good()); + + ifcopenshell::geom::settings settings; + ifcopenshell::geom::converter converter( + ifcopenshell::geom::kernels::construct(&file, "opencascade", settings, log), &file, settings, log); + + std::vector tasks; + std::vector filters; + converter.mapping()->get_representations(tasks, filters); + REQUIRE(!tasks.empty()); + + sectioned_solid_result result; + for (const auto& task : tasks) { + REQUIRE(!task.products.empty()); + auto* elem = converter.create_brep_for_representation_and_product(task.representation, task.products.front()); + if (elem) { + result.produced_brep = true; + elem->calculate_projected_surface_area( + result.projected_area_x, result.projected_area_y, result.projected_area_z); + } + delete elem; + } + result.geo42_count = log.count("GEO42"); + return result; +} + +std::size_t count_geo42_converting(const std::string& spf) { + return convert_sectioned_solid(spf).geo42_count; +} + +// A minimal IfcSectionedSolidHorizontal (IFC4X3_ADD2) whose cross section +// placements carry an explicit Axis = (0,0,1) but no RefDirection, on a +// directrix that runs along +Y (not the global +X). Per buildingSMART +// IFC4.x-IF #147 the profile normal follows the directrix tangent, so the +// 12 x 0.5 rectangle sweeps 60 along +Y: a plan (Z) projected area of ~720. +// Before the fix the profile was placed with a fixed axis permutation that +// ignored the directrix and the solid collapsed. +constexpr const char* AXIS_ALIGNED_SECTIONED_SOLID_SPF = R"IFC(ISO-10303-21; +HEADER; +FILE_DESCRIPTION((''),'2;1'); +FILE_NAME('','',(''),(''),'','',''); +FILE_SCHEMA(('IFC4X3_ADD2')); +ENDSEC; +DATA; +#1=IFCPROJECT('0RYK8PV8D0ee9DDm77xcTZ',$,'T',$,$,$,$,(#6),$); +#2=IFCCARTESIANPOINT((0.,0.,0.)); +#3=IFCDIRECTION((0.,0.,1.)); +#4=IFCDIRECTION((1.,0.,0.)); +#5=IFCAXIS2PLACEMENT3D(#2,#3,#4); +#6=IFCGEOMETRICREPRESENTATIONCONTEXT($,'Model',3,1.E-05,#5,$); +#7=IFCGEOMETRICREPRESENTATIONSUBCONTEXT('Body','Model',*,*,*,*,#6,$,.MODEL_VIEW.,$); +#8=IFCCARTESIANPOINT((0.,0.,0.)); +#9=IFCCARTESIANPOINT((0.,60.,0.)); +#10=IFCPOLYLINE((#8,#9)); +#11=IFCDIRECTION((0.,0.,1.)); +#13=IFCPOINTBYDISTANCEEXPRESSION(IFCLENGTHMEASURE(0.),$,$,$,#10); +#14=IFCAXIS2PLACEMENTLINEAR(#13,#11,$); +#15=IFCPOINTBYDISTANCEEXPRESSION(IFCLENGTHMEASURE(60.),$,$,$,#10); +#16=IFCAXIS2PLACEMENTLINEAR(#15,#11,$); +#17=IFCCARTESIANPOINTLIST2D(((-6.,0.),(6.,0.),(6.,0.5),(-6.,0.5),(-6.,0.)),$); +#18=IFCINDEXEDPOLYCURVE(#17,$,.F.); +#19=IFCARBITRARYCLOSEDPROFILEDEF(.AREA.,$,#18); +#23=IFCSECTIONEDSOLIDHORIZONTAL(#10,(#19,#19),(#14,#16)); +#24=IFCBUILDINGELEMENTPROXY('3cWcr4$892GAeKwryhDILR',$,'pavement',$,$,$,#26,$,$); +#25=IFCSHAPEREPRESENTATION(#7,'Body','AdvancedSweptSolid',(#23)); +#26=IFCPRODUCTDEFINITIONSHAPE($,$,(#25)); +ENDSEC; +END-ISO-10303-21; +)IFC"; + } // namespace +TEST_CASE("IfcSectionedSolidHorizontal raked end cut does not log GEO 42", "[ifcgeom][infra-sweep]") { + if (std::string(STRINGIFY(IfcSchema)) != "Ifc4x3_add2") { + SKIP("fixture is authored for IFC4X3_ADD2"); + } + CHECK(count_geo42_converting(RAKED_SECTIONED_SOLID_SPF) == 0); +} + +TEST_CASE("IfcSectionedSolidHorizontal follows a non-axis-aligned directrix", "[ifcgeom][infra-sweep]") { + if (std::string(STRINGIFY(IfcSchema)) != "Ifc4x3_add2") { + SKIP("fixture is authored for IFC4X3_ADD2"); + } + const auto result = convert_sectioned_solid(AXIS_ALIGNED_SECTIONED_SOLID_SPF); + CHECK(result.geo42_count == 0); + REQUIRE(result.produced_brep); + // Plan projection is the top plus the bottom of the slab, 2 x width x length; + // a collapsed sweep (the pre-fix behaviour) is nowhere near this. + CHECK(result.projected_area_z == Catch::Approx(2.0 * 12.0 * 60.0).margin(2.0)); +} + TEST_CASE("IfcGeom C++ fixture creates walls below and above the void limit", "[ifcgeom][voids]") { hierarchy_helper below_limit_file; const auto below_limit_wall = create_wall_with_voids(below_limit_file, MAX_VOIDS - 1); diff --git a/src/ifcopenshell-python/test/test_create_shape.py b/src/ifcopenshell-python/test/test_create_shape.py index 4118685dde..f9f9fd1af5 100644 --- a/src/ifcopenshell-python/test/test_create_shape.py +++ b/src/ifcopenshell-python/test/test_create_shape.py @@ -1,10 +1,12 @@ import functools import itertools +import math import multiprocessing import operator import os from typing import get_args +import numpy as np import pytest import ifcopenshell @@ -270,6 +272,276 @@ def test_logging(): ] +class TestSectionedSolidHorizontalRakedEndCut(test.bootstrap.IFC4X3): + """An IfcSectionedSolidHorizontal whose two IfcAxis2PlacementLinear cross + section positions use direction vectors inconsistently -- one raked + RefDirection + width scale, the other plain -- must still loft a uniform + prism (raked at one end, square at the other), not a wedge, and must not + log GEO 42.""" + + def _build_wingwall(self, theta): + f = self.file + ifcopenshell.api.root.create_entity(f, ifc_class="IfcProject", name="Test") + ctx = ifcopenshell.api.context.add_context(f, context_type="Model") + body = ifcopenshell.api.context.add_context( + f, context_type="Model", context_identifier="Body", target_view="MODEL_VIEW", parent=ctx + ) + + self.L, self.w, self.h = 40.0, 20.0 / 12.0, 6.0 + width_scale = 1.0 / math.cos(theta) + + directrix = f.createIfcPolyline( + Points=[f.createIfcCartesianPoint((0.0, 0.0, 0.0)), f.createIfcCartesianPoint((self.L, 0.0, 0.0))] + ) + + def rect(half_width): + coords = ( + (-half_width, 0.0), + (half_width, 0.0), + (half_width, self.h), + (-half_width, self.h), + (-half_width, 0.0), + ) + return f.createIfcArbitraryClosedProfileDef( + ProfileType="AREA", + OuterCurve=f.createIfcIndexedPolyCurve( + Points=f.createIfcCartesianPointList2D(coords), Segments=None, SelfIntersect=False + ), + ) + + axis = f.createIfcDirection((0.0, 0.0, 1.0)) + # Rakes the near end cap about the vertical axis by theta while keeping + # the section's width axis (Axis x RefDirection) perpendicular extent at w. + ref_direction = f.createIfcDirection((math.cos(theta), -math.sin(theta), 0.0)) + + def location(distance_along): + return f.createIfcPointByDistanceExpression( + DistanceAlong=f.createIfcLengthMeasure(distance_along), BasisCurve=directrix + ) + + near = f.createIfcAxis2PlacementLinear(Location=location(0.0), Axis=axis, RefDirection=ref_direction) + far = f.createIfcAxis2PlacementLinear(Location=location(self.L), Axis=axis) + + solid = f.createIfcSectionedSolidHorizontal( + Directrix=directrix, + CrossSections=[rect(0.5 * self.w * width_scale), rect(0.5 * self.w)], + CrossSectionPositions=[near, far], + ) + return f.createIfcShapeRepresentation( + ContextOfItems=body, + RepresentationIdentifier="Body", + RepresentationType="AdvancedSweptSolid", + Items=[solid], + ) + + @pytest.mark.skipif( + not ifcopenshell.geom.has_geometry_library("opencascade"), reason="requires the OpenCASCADE kernel" + ) + def test_raked_end_keeps_uniform_perpendicular_thickness(self): + theta = math.radians(25.0 + 22.0 / 60.0 + 58.0 / 3600.0) + representation = self._build_wingwall(theta) + + logger = ifcopenshell.logger() + logger.output_format(logger.FMT_INMEMORY) + settings = ifcopenshell.geom.settings() + settings.set("use-world-coords", True) + geometry = ifcopenshell.geom.create_shape(settings, representation, logger=logger) + + assert "GEO42" not in [msg.code for msg in logger] + + v = ifcopenshell.util.shape.get_vertices(geometry) + xs = v[:, 0] + + # Uniform thickness perpendicular to the (X aligned) directrix. The whole + # solid spans exactly w in Y; a wedge (the pre fix behaviour) does not. + assert np.ptp(v[:, 1]) == pytest.approx(self.w, abs=1e-4) + for x0 in np.linspace(3.0, self.L - 3.0, 12): + slab = v[np.abs(xs - x0) < 1.0] + if len(slab) < 4: + continue + assert np.ptp(slab[:, 1]) == pytest.approx(self.w, abs=1e-4), f"thickness at x={x0:.1f}" + + # Far end square to the directrix (all vertices at x == L), near end raked + # about the vertical axis by exactly theta (its extreme vertex sits at + # -w/2 * tan(theta) along the directrix). + assert xs.max() == pytest.approx(self.L, abs=1e-4) + assert np.ptp(v[xs > xs.max() - 1e-4][:, 0]) < 1e-4 + assert xs.min() == pytest.approx(-0.5 * self.w * math.tan(theta), abs=1e-4) + + +class TestSectionedSolidHorizontalHonoursAxis(test.bootstrap.IFC4X3): + """An IfcSectionedSolidHorizontal whose cross section placements carry an + explicit Axis = (0,0,1) but no RefDirection must sweep the profile with its + Y axis on that Axis and its normal on the directrix tangent (buildingSMART + IFC4.x-IF #147) -- i.e. following the curve. Regression test for the case + where the directrix does not run along the global X axis: the profile used + to be placed with a fixed [e_y | e_z | e_x] permutation that ignored the + directrix, collapsing the swept solid (a road pavement running north-south + would come out a sliver a few centimetres wide).""" + + def _build(self, theta=0.0): + f = self.file + ifcopenshell.api.root.create_entity(f, ifc_class="IfcProject", name="Test") + ctx = ifcopenshell.api.context.add_context(f, context_type="Model") + body = ifcopenshell.api.context.add_context( + f, context_type="Model", context_identifier="Body", target_view="MODEL_VIEW", parent=ctx + ) + + self.L, self.width, self.height = 60.0, 12.0, 0.5 + + # Directrix runs along +Y, not +X. + directrix = f.createIfcPolyline( + Points=[f.createIfcCartesianPoint((0.0, 0.0, 0.0)), f.createIfcCartesianPoint((0.0, self.L, 0.0))] + ) + + def rect(half_width): + coords = ( + (-half_width, 0.0), + (half_width, 0.0), + (half_width, self.height), + (-half_width, self.height), + (-half_width, 0.0), + ) + return f.createIfcArbitraryClosedProfileDef( + ProfileType="AREA", + OuterCurve=f.createIfcIndexedPolyCurve( + Points=f.createIfcCartesianPointList2D(coords), Segments=None, SelfIntersect=False + ), + ) + + axis = f.createIfcDirection((0.0, 0.0, 1.0)) + + def location(distance_along): + return f.createIfcPointByDistanceExpression( + DistanceAlong=f.createIfcLengthMeasure(distance_along), BasisCurve=directrix + ) + + near = f.createIfcAxis2PlacementLinear(Location=location(0.0), Axis=axis) + if theta: + # Directrix tangent is +Y; rake the far cap about the vertical Axis. + ref_direction = f.createIfcDirection((math.sin(theta), math.cos(theta), 0.0)) + far = f.createIfcAxis2PlacementLinear(Location=location(self.L), Axis=axis, RefDirection=ref_direction) + far_half = 0.5 * self.width / math.cos(theta) + else: + far = f.createIfcAxis2PlacementLinear(Location=location(self.L), Axis=axis) + far_half = 0.5 * self.width + + solid = f.createIfcSectionedSolidHorizontal( + Directrix=directrix, + CrossSections=[rect(0.5 * self.width), rect(far_half)], + CrossSectionPositions=[near, far], + ) + return f.createIfcShapeRepresentation( + ContextOfItems=body, + RepresentationIdentifier="Body", + RepresentationType="AdvancedSweptSolid", + Items=[solid], + ) + + def _shape(self, representation): + logger = ifcopenshell.logger() + logger.output_format(logger.FMT_INMEMORY) + settings = ifcopenshell.geom.settings() + settings.set("use-world-coords", True) + geometry = ifcopenshell.geom.create_shape(settings, representation, logger=logger) + assert "GEO42" not in [msg.code for msg in logger] + return ifcopenshell.util.shape.get_vertices(geometry) + + @pytest.mark.skipif( + not ifcopenshell.geom.has_geometry_library("opencascade"), reason="requires the OpenCASCADE kernel" + ) + def test_square_run_follows_the_directrix(self): + v = self._shape(self._build()) + # width across the road -> world X, crown -> world Z, length -> world Y. + assert np.ptp(v[:, 0]) == pytest.approx(self.width, abs=1e-4) + assert np.ptp(v[:, 1]) == pytest.approx(self.L, abs=1e-4) + assert np.ptp(v[:, 2]) == pytest.approx(self.height, abs=1e-4) + + @pytest.mark.skipif( + not ifcopenshell.geom.has_geometry_library("opencascade"), reason="requires the OpenCASCADE kernel" + ) + def test_raked_far_end_on_a_non_axis_aligned_directrix(self): + theta = math.radians(25.0 + 22.0 / 60.0 + 58.0 / 3600.0) + v = self._shape(self._build(theta)) + ys = v[:, 1] + # Uniform width perpendicular to the directrix, square near end, raked far end. + assert np.ptp(v[:, 0]) == pytest.approx(self.width, abs=1e-4) + for y0 in np.linspace(5.0, self.L - 5.0, 10): + slab = v[np.abs(ys - y0) < 1.0] + if len(slab) >= 4: + assert np.ptp(slab[:, 0]) == pytest.approx(self.width, abs=1e-4), f"width at y={y0:.1f}" + assert ys.min() == pytest.approx(0.0, abs=1e-4) + assert np.ptp(v[ys < ys.min() + 1e-4][:, 1]) < 1e-4 + assert ys.max() == pytest.approx(self.L + 0.5 * self.width * math.tan(theta), abs=1e-4) + + +class TestSectionedSolidHorizontalOffsetUnits(test.bootstrap.IFC4X3): + """IfcPointByDistanceExpression.OffsetLateral / OffsetVertical are + IfcLengthMeasure and must be scaled by the model length unit, exactly like + DistanceAlong. Regression test: in a millimetre model a 3000 mm lateral + offset must move the swept solid 3 m, not 3000 m.""" + + @pytest.mark.skipif( + not ifcopenshell.geom.has_geometry_library("opencascade"), reason="requires the OpenCASCADE kernel" + ) + def test_lateral_offset_is_scaled_by_the_length_unit(self): + f = self.file + ifcopenshell.api.root.create_entity(f, ifc_class="IfcProject", name="Test") + unit = ifcopenshell.api.unit.add_si_unit(f, unit_type="LENGTHUNIT", prefix="MILLI") + ifcopenshell.api.unit.assign_unit(f, units=[unit]) + ctx = ifcopenshell.api.context.add_context(f, context_type="Model") + body = ifcopenshell.api.context.add_context( + f, context_type="Model", context_identifier="Body", target_view="MODEL_VIEW", parent=ctx + ) + + length_mm, width_mm, offset_mm = 20000.0, 4000.0, 3000.0 + directrix = f.createIfcPolyline( + Points=[f.createIfcCartesianPoint((0.0, 0.0, 0.0)), f.createIfcCartesianPoint((length_mm, 0.0, 0.0))] + ) + hw = 0.5 * width_mm + profile = f.createIfcArbitraryClosedProfileDef( + ProfileType="AREA", + OuterCurve=f.createIfcIndexedPolyCurve( + Points=f.createIfcCartesianPointList2D(((-hw, 0.0), (hw, 0.0), (hw, 500.0), (-hw, 500.0), (-hw, 0.0))), + Segments=None, + SelfIntersect=False, + ), + ) + axis = f.createIfcDirection((0.0, 0.0, 1.0)) + + def position(distance_along): + return f.createIfcAxis2PlacementLinear( + Location=f.createIfcPointByDistanceExpression( + DistanceAlong=f.createIfcLengthMeasure(distance_along), + OffsetLateral=offset_mm, + BasisCurve=directrix, + ), + Axis=axis, + ) + + solid = f.createIfcSectionedSolidHorizontal( + Directrix=directrix, + CrossSections=[profile, profile], + CrossSectionPositions=[position(0.0), position(length_mm)], + ) + representation = f.createIfcShapeRepresentation( + ContextOfItems=body, + RepresentationIdentifier="Body", + RepresentationType="AdvancedSweptSolid", + Items=[solid], + ) + + settings = ifcopenshell.geom.settings() + settings.set("use-world-coords", True) + geometry = ifcopenshell.geom.create_shape(settings, representation) + v = ifcopenshell.util.shape.get_vertices(geometry) # metres + + # Directrix +X, Axis +Z -> profile local x is world +Y; +OffsetLateral shifts there. + assert np.ptp(v[:, 0]) == pytest.approx(length_mm / 1000.0, abs=1e-4) + assert np.ptp(v[:, 1]) == pytest.approx(width_mm / 1000.0, abs=1e-4) + assert v[:, 1].mean() == pytest.approx(offset_mm / 1000.0, abs=1e-3) + + if __name__ == "__main__": import pytest