IfcPolynomialCurve

Core/src/IfcGeometryConverter/CurveConverter.h

@@ -1593,6 +1593,11 @@ namespace OpenInfraPlatform
 							point = getPointOnCurve(curveSegment->ParentCurve.as<typename IfcEntityTypesT::IfcCircle>(), runningLength);
 							direction = getDirectionOfCurve(curveSegment->ParentCurve.as<typename IfcEntityTypesT::IfcCircle>(), runningLength);
 						}
+						else if (curveSegment->ParentCurve.isOfType<typename IfcEntityTypesT::IfcPolynomialCurve>())
+						{
+							point = getPointOnCurve(curveSegment->ParentCurve.as<typename IfcEntityTypesT::IfcPolynomialCurve>(), runningLength);
+							direction = getDirectionOfCurve(curveSegment->ParentCurve.as<typename IfcEntityTypesT::IfcPolynomialCurve>(), runningLength);
+						}
 						segmentPoints.push_back(point);
 						segmentDirections.push_back(direction);
 						// determine next length
@@ -1618,6 +1623,7 @@ namespace OpenInfraPlatform
 
 					if (!segmentPoints.empty())
 					{
+						//TODO: implement different rotations around y and z axis
 						//get the local coordinate system
 						carve::geom::vector<3> tangent = segmentDirections[0].normalize();
 						tangent.y = -tangent.y;
@@ -2895,7 +2901,297 @@ namespace OpenInfraPlatform
 					return carve::geom::VECTOR(std::cos(angle), std::sin(angle), 0.);
 				}
 #endif
+#if defined(OIP_MODULE_EARLYBINDING_IFC4X3_RC4)
+				/*! \brief Calculates an angle of the polynomial curve.
+				* \param[in] polynomialCurve		A pointer to data from \c IfcPolynomialCurve.
+				* \param[in] parameter				The length.
+				* \return							The direction of the curve.
+				* \note
+				*/
+				template <>
+				carve::geom::vector<3> getDirectionOfCurve(const EXPRESSReference<typename IfcEntityTypesT::IfcPolynomialCurve>& polynomialCurve,
+					const typename IfcEntityTypesT::IfcParameterValue& parameter) const noexcept(false)
+				{
+					return getDirectionOfCurve(polynomialCurve, parameter * this->UnitConvert()->getLengthInMeterFactor());
+				}
+				template <>
+				carve::geom::vector<3> getDirectionOfCurve(const EXPRESSReference<typename IfcEntityTypesT::IfcPolynomialCurve>& polynomialCurve,
+					const typename IfcEntityTypesT::IfcNonNegativeLengthMeasure& parameter) const noexcept(false)
+				{
+					return getDirectionOfCurve(polynomialCurve, parameter * this->UnitConvert()->getLengthInMeterFactor());
+				}
+				template<>
+				carve::geom::vector<3> getDirectionOfCurve(const EXPRESSReference<typename IfcEntityTypesT::IfcPolynomialCurve>& polynomialCurve,
+					const double& parameter) const noexcept(false)
+				{				
+					//std::vector<typename IfcEntityTypesT::IfcReal> coefficientX, coefficientY, coefficientZ;
+					//std::vector<double> polynomialConstantX, polynomialConstantY, polynomialConstantZ;
+					int polynomialConstantCntX, polynomialConstantCntY, polynomialConstantCntZ;
+					double angleX, angleY, angleZ;
+					//get coefficients
+					if (polynomialCurve->CoefficientsX)
+					{
+						std::vector<double> polynomialConstantX;
+						std::vector<typename IfcEntityTypesT::IfcReal> coefficientX = polynomialCurve->CoefficientsX;
+						//convert to double
+						for (int i = 0; i < std::size(coefficientX); i++)
+						{
+							polynomialConstantX.push_back(double(coefficientX[i]));
+						}
+						polynomialConstantCntX = std::size(polynomialConstantX);
+						//calculate angle
+						angleX = SpiralUtils::AngleByAngleDeviationPolynomial(polynomialConstantX, polynomialConstantCntX, parameter);
+					}
+					if (polynomialCurve->CoefficientsY)
+					{
+						std::vector<double> polynomialConstantY;
+						std::vector<typename IfcEntityTypesT::IfcReal> coefficientY = polynomialCurve->CoefficientsY;
+						//convert to double
+						for (int i = 0; i < std::size(coefficientY); i++)
+						{
+							polynomialConstantY.push_back(double(coefficientY[i]));
+						}
+						polynomialConstantCntY = std::size(polynomialConstantY);
+						//calculate angle
+						angleY = SpiralUtils::AngleByAngleDeviationPolynomial(polynomialConstantY, polynomialConstantCntY, parameter);
+					}
+					if (polynomialCurve->CoefficientsZ)
+					{
+						std::vector<double> polynomialConstantZ;
+						std::vector<typename IfcEntityTypesT::IfcReal> coefficientZ = polynomialCurve->CoefficientsZ;
+						//convert to double
+						for (int i = 0; i < std::size(coefficientZ); i++)
+						{
+							polynomialConstantZ.push_back(double(coefficientZ[i]));
+						}
+						polynomialConstantCntZ = std::size(polynomialConstantZ);
+						//calculate angle
+						angleZ = SpiralUtils::AngleByAngleDeviationPolynomial(polynomialConstantZ, polynomialConstantCntZ, parameter);
+					}
+
+					// calculate angle between two polynomial curves
+					//double angle = std::atan2(angleY, angleX);
+
+					if (polynomialConstantCntX>0 && polynomialConstantCntY>0)
+					{
+						double angle = std::atan2(sin(angleY), cos(angleX));
+						//double angle = 0.;
+						return carve::geom::VECTOR(std::cos(angle), std::sin(angle), 0.);
+
+					}
+					// Only works if correctly rotation around y axis is in convertIfcCurveSegment 
+					else if (polynomialConstantCntX>0 && polynomialConstantCntZ>0)
+					{
+						//double angle = std::atan2(angleZ, angleX);
+						return carve::geom::VECTOR(std::cos(angleX), 0., -std::sin(angleZ));
+					}
+					// Only works if correctly rotation around x axis is in convertIfcCurveSegment 
+					else if (polynomialConstantCntY>0 && polynomialConstantCntZ>0)
+					{
+						//double angle = std::atan2(angleZ, angleY);
+						return carve::geom::VECTOR(0., std::cos(angleY), std::sin(angleZ));
+					}
+
+					//TODO: angle in 3D. At first implement additional code in convertIfcCurveSegment
+					/*else if (polynomialConstantCntX>0 && polynomialConstantCntY>0 && polynomialConstantCntZ>0)
+					{
+
+					}*/			
+				}
+
+#endif
+
+#if defined(OIP_MODULE_EARLYBINDING_IFC4X3_RC4)
+				/*! \brief Calculates a trimming point on the polynomial curve.
+				* \param[in] polynomialCurve	    A pointer to data from \c IfcPolynomialCurve.
+				* \param[in] parameter				A pointer to data from \c IfcCurveSegment.
+				* \return							The location of the trimming point.
+				* \note
+				*/
+				template <>
+				carve::geom::vector<3> getPointOnCurve(const EXPRESSReference<typename IfcEntityTypesT::IfcPolynomialCurve>& polynomialCurve,
+					const typename IfcEntityTypesT::IfcParameterValue& parameter) const noexcept(false)
+				{
+					return getPointOnCurve(polynomialCurve, parameter * this->UnitConvert()->getLengthInMeterFactor());
+				}
+				template <>
+				carve::geom::vector<3> getPointOnCurve(const EXPRESSReference<typename IfcEntityTypesT::IfcPolynomialCurve>& polynomialCurve,
+					const typename IfcEntityTypesT::IfcNonNegativeLengthMeasure& parameter) const noexcept(false)
+				{
+					return getPointOnCurve(polynomialCurve, parameter * this->UnitConvert()->getLengthInMeterFactor());
+				}
+
+				carve::geom::vector<3> getPointOnCurve(const EXPRESSReference<typename IfcEntityTypesT::IfcPolynomialCurve>& polynomialCurve,
+					const double& parameter) const noexcept(false)
+
+				{
+
+					std::vector<typename IfcEntityTypesT::IfcReal> coefficientX, coefficientY, coefficientZ;
+					//typename IfcEntityTypesT::IfcReal x, y, z;
+					double x, y, z;
+					int polynomialConstantCntX, polynomialConstantCntY, polynomialConstantCntZ;
+					std::vector<double> polynomialConstantX, polynomialConstantY, polynomialConstantZ;
+
+					if (polynomialCurve->CoefficientsX)
+					{
+						// Interpret coefficients for X
+						coefficientX = polynomialCurve->CoefficientsX;
+						polynomialConstantCntX = std::size(coefficientX);
+						//std::vector<double> polynomialConstantX;
+						//convert to double
+						for (int i = 0; i < std::size(coefficientX); i++)
+						{
+							polynomialConstantX.push_back(double(coefficientX[i]));
+						}
+
+						//x = SpiralUtils::XbyAngleDeviationPolynomial(polynomialConstantX, polynomialConstantCntX, parameter);
+						//x = SpiralUtils::XbyAngleDeviationPolynomialByTerms( 0., 0., 0., 0., 0., 0., 1., 0., parameter);
+
+					}
+					if (polynomialCurve->CoefficientsY)
+					{
+						// Interpret coefficients for Y
+						coefficientY = polynomialCurve->CoefficientsY;
+						// Implement polynomial term for y coordinate
+						polynomialConstantCntY = std::size(coefficientY);
+						//std::vector<double> polynomialConstantY;
+						//convert to double
+						for (int i = 0; i < std::size(coefficientY); i++)
+						{
+							polynomialConstantY.push_back(double(coefficientY[i]));
+						}
+						//y = SpiralUtils::YbyAngleDeviationPolynomial(polynomialConstantY, polynomialConstantCntY, parameter);
+						//y = SpiralUtils::XbyAngleDeviationPolynomialByTerms( 0., 0., 0., 0., 0., 1., 0., 0., parameter);
+
+					}
+					if (polynomialCurve->CoefficientsZ)
+					{
+						// Interpret coefficients for Z
+						coefficientZ = polynomialCurve->CoefficientsZ;
+						polynomialConstantCntZ = std::size(coefficientZ);
+						// Implement polynomial term for z coordinate/ We don't have any functions for supporting z
+						z = 0.;
+					}
+					// Impplement different cases for 2D and 3D
+					if (polynomialConstantCntX > 0 && polynomialConstantCntY > 0)
+					{
+						double t = integrateParameter(polynomialConstantX, polynomialConstantY, parameter);
+						x = calculatePolynomialCurve(polynomialConstantX, t);
+						y = calculatePolynomialCurve(polynomialConstantY, t);
+						return carve::geom::VECTOR(x, y, 0.);
+					}
+
+					else if (polynomialConstantCntX>0 && polynomialConstantCntZ>0) return carve::geom::VECTOR(x, 0., z);
+
+					else if (polynomialConstantCntY>0 && polynomialConstantCntZ>0) return carve::geom::VECTOR(0., y, z);
+
+					else if (polynomialConstantCntX>0 && polynomialConstantCntY>0 && polynomialConstantCntZ>0) return carve::geom::VECTOR(x, y, z);
+				}
+
+				/*! \brief Calculates a polinomial curve
+				* \param[in] coefficient			A vector with coeffitients corresponding ordinate from \c IfcPolynomialCurve.
+				* \param[in] parameter				A length during the curve from \c IfcCurveSegment.
+				* \return							The location of the ordinate point.
+				* \note
+				*/
+				typename IfcEntityTypesT::IfcReal calculatePolynomialCurve(std::vector<double> coefficient, const double parameter) const noexcept(false)
+				{
+					int n = std::size(coefficient);
+					typename IfcEntityTypesT::IfcReal sum = coefficient[0];
+					double factor = 1;
+
+					for (int i = 1; i < n; i++)
+					{
+						factor *= parameter;
+						sum += coefficient[i] * factor;
+					}
+					return sum;
+				}
+
+				/*! \brief Calculates derivative of the polynomial curve
+				* \param[in] polynomialConstant	    A vector with coeffitients corresponding ordinate from \c IfcPolynomialCurve.
+				* \param[in] b				        Parameter value from intagrateParameter
+				* \return							Derivative value.
+				* \note
+				*/
+				double calculatePolynomialDerivative(std::vector<double> polynomialConstant, double b)  const noexcept(false)
+				{
+					int polynomialConstantCnt = std::size(polynomialConstant);
+					double value = polynomialConstant[1];
+					int factor = 2;
+					// derivative of the polynomial
+					for (int i = 2; i < polynomialConstantCnt; i++)
+					{
+						value += polynomialConstant[i] * factor * b;
+						b *= b;
+						factor += 1;
+					}
+					return value;
+				}
+				/*! \brief Calculates parameter for parametric polynomial curve
+				* \param[in] polynomialConstantX	A vector with coeffitients axis X from \c IfcPolynomialCurve.
+				* \param[in] polynomialConstantY	A vector with coeffitients axis Y from \c IfcPolynomialCurve.
+				* \param[in] length				    Length of the curve cegment
+				* \return							Parameter value.
+				* \note
+				*/
+				double integrateParameter(std::vector<double> polynomialConstantX, std::vector<double> polynomialConstantY, const double length) const noexcept(false)
+				{
+					double a = 0.;
+					double b = 0.;
+					int n = 1;// small numbers return better values
+					double value = 0.;
+
+					if (length == 0) { return b; }
+					// wide steps 0.1
+					while (value <= length)
+					{
+						value = integral(polynomialConstantX, polynomialConstantY, a, b, n);
+						if (value == length) return b;
+						b += 0.1;
+					}
+					b -= 0.2;
+					value = integral(polynomialConstantX, polynomialConstantY, a, b, n);
+					// check the error and find the apropriate value
+					double error = length - value;
+					double border = 0.0001;
+					if (error > border)
+					{
+						while (value <= length)
+						{
+							value = integral(polynomialConstantX, polynomialConstantY, a, b, n);
+							if (value == length) return b;
+							b += border;
+						}
+					}
+					double parameter = b - 2 * border;
+					return parameter;
+				}
+
+				/*! \brief Calculates parameter for parametric polynomial curve
+				* \param[in] polynomialConstantX	A vector with coeffitients axis X from \c IfcPolynomialCurve.
+				* \param[in] polynomialConstantY	A vector with coeffitients axis Y from \c IfcPolynomialCurve.
+				* \param[in] a				        Low limit of the Integral
+				* \param[in] b				        Upper limit of the Integral
+				* \param[in] n				        Number of integral steps
+				* \return							Value of the integral.
+				* \note
+				*/
+				// source code https://helloacm.com/c-function-to-compute-numerical-integral-using-function-pointers/
 
+				double integral(std::vector<double> polynomialConstantX, std::vector<double> polynomialConstantY, double a, double b, int n) const noexcept(false)
+				{
+					double step = (b - a) / n;  // width of each small rectangle
+					double area = 0.0;  // signed area
+					for (int i = 0; i < n; i++) {
+						//source arcticle https://www.math.usm.edu/lambers/mat169/fall09/lecture31.pdf
+						area += sqrt(pow(calculatePolynomialDerivative(polynomialConstantY, (a + (i + 0.5) * step)),2)
+							+ pow(calculatePolynomialDerivative(polynomialConstantX, (a + (i + 0.5) * step)),2)) * step; // sum up each small rectangle
+					}
+					return area;
+				}
+
+#endif
 
 			protected:
 

UnitTests/Schemas/IFC4X3_RC4/CMakeLists.txt

@@ -0,0 +1,26 @@
+#
+#    Copyright (c) 2021 Technical University of Munich
+#    Chair of Computational Modeling and Simulation.
+#
+#    TUM Open Infra Platform is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License Version 3
+#    as published by the Free Software Foundation.
+#
+#    TUM Open Infra Platform is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+
+include(CreateUnitTests)
+
+# get all subdirectories (that is, unit tests)
+SUBDIRLIST(SUBDIRS ${CMAKE_CURRENT_LIST_DIR})
+
+# add to the solution
+FOREACH(subdir ${SUBDIRS})
+	ADD_SUBDIRECTORY(${subdir})
+ENDFOREACH()

UnitTests/Schemas/IFC4X3_RC4/bloss-curve_100.0_300_1000_1_Meter/CMakeLists.txt

@@ -0,0 +1,20 @@
+#
+#    Copyright (c) 2021 Technical University of Munich
+#    Chair of Computational Modeling and Simulation.
+#
+#    TUM Open Infra Platform is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License Version 3
+#    as published by the Free Software Foundation.
+#
+#    TUM Open Infra Platform is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+
+include(CreateUnitTests)
+
+CreateIfcFileVisualUnitTestForSchema(bloss-curve_100.0_300_1000_1_Meter IFC4X3_RC4)