|
32 | 32 | #include "MeshKernel/Constants.hpp" |
33 | 33 | #include "MeshKernel/Exceptions.hpp" |
34 | 34 | #include "MeshKernel/LandBoundary.hpp" |
| 35 | +#include "MeshKernel/MeshFaceCenters.hpp" |
35 | 36 | #include "MeshKernel/Operations.hpp" |
36 | 37 | #include "MeshKernel/Polygon.hpp" |
37 | 38 | #include "MeshKernel/Vector.hpp" |
@@ -751,147 +752,8 @@ std::tuple<double, meshkernel::Point, meshkernel::TraversalDirection> meshkernel |
751 | 752 | double area = 0.0; |
752 | 753 |
|
753 | 754 | const double minArea = 1e-8; |
754 | | - const auto numberOfPointsOpenedPolygon = static_cast<UInt>(polygon.size()) - 1; |
755 | 755 |
|
756 | | - const double updateStepSize = 0.1; |
757 | | - |
758 | | - Point centreOfMass(0.0, 0.0); |
759 | | - |
760 | | - for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) |
761 | | - { |
762 | | - centreOfMass += polygon[n]; |
763 | | - } |
764 | | - |
765 | | - centreOfMass *= 1.0 / static_cast<double>(numberOfPointsOpenedPolygon); |
766 | | - // Will be non-unity for spherical coordinates only |
767 | | - const double xTransformation = projection == Projection::cartesian ? 1.0 : 1.0 / std::cos(centreOfMass.y * constants::conversion::degToRad); |
768 | | - const double circumcentreTolerance = constants::geometric::circumcentreTolerance * (projection == Projection::cartesian ? 1.0 : 1.0 / (constants::geometric::earth_radius * constants::conversion::degToRad)); |
769 | | - |
770 | | - if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) |
771 | | - { |
772 | | - |
773 | | - Point midPoint1 = 0.5 * (polygon[0] + polygon[1]); |
774 | | - Point midPoint2 = 0.5 * (polygon[1] + polygon[2]); |
775 | | - Point midPoint3 = 0.5 * (polygon[2] + polygon[0]); |
776 | | - |
777 | | - Vector edgeVector1 = static_cast<Vector>(NormalVector(polygon[0], polygon[1], midPoint1, projection)); |
778 | | - Vector edgeVector2 = static_cast<Vector>(NormalVector(polygon[1], polygon[2], midPoint2, projection)); |
779 | | - Vector edgeVector3 = static_cast<Vector>(NormalVector(polygon[2], polygon[0], midPoint3, projection)); |
780 | | - |
781 | | - edgeVector1.normalise(); |
782 | | - edgeVector2.normalise(); |
783 | | - edgeVector3.normalise(); |
784 | | - |
785 | | - Vector edgeVectorSum = edgeVector1 + edgeVector2 + edgeVector3; |
786 | | - |
787 | | - double edgeVectorSumLength = edgeVectorSum.length(); |
788 | | - |
789 | | - for (UInt i = 1; i <= MaximumNumberOfCircumcentreIterations; ++i) |
790 | | - { |
791 | | - Vector delta1 = GetDelta(midPoint1, centreOfMass, projection); |
792 | | - Vector delta2 = GetDelta(midPoint2, centreOfMass, projection); |
793 | | - Vector delta3 = GetDelta(midPoint3, centreOfMass, projection); |
794 | | - |
795 | | - double ds = dot(delta1, edgeVector1) + dot(delta2, edgeVector2) + dot(delta3, edgeVector3); |
796 | | - |
797 | | - if (projection != Projection::cartesian) |
798 | | - { |
799 | | - ds *= constants::conversion::radToDeg * constants::geometric::inverse_earth_radius; |
800 | | - } |
801 | | - |
802 | | - centreOfMass.x -= updateStepSize * ds * edgeVectorSum.x() * xTransformation; |
803 | | - centreOfMass.y -= updateStepSize * ds * edgeVectorSum.y(); |
804 | | - |
805 | | - if (ds * edgeVectorSumLength < circumcentreTolerance || i == MaximumNumberOfCircumcentreIterations) |
806 | | - { |
807 | | - break; |
808 | | - } |
809 | | - } |
810 | | - } |
811 | | - else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) |
812 | | - { |
813 | | - |
814 | | - Point midPoint1 = 0.5 * (polygon[0] + polygon[1]); |
815 | | - Point midPoint2 = 0.5 * (polygon[1] + polygon[2]); |
816 | | - Point midPoint3 = 0.5 * (polygon[2] + polygon[3]); |
817 | | - Point midPoint4 = 0.5 * (polygon[3] + polygon[0]); |
818 | | - |
819 | | - Vector edgeVector1 = static_cast<Vector>(NormalVector(polygon[0], polygon[1], midPoint1, projection)); |
820 | | - Vector edgeVector2 = static_cast<Vector>(NormalVector(polygon[1], polygon[2], midPoint2, projection)); |
821 | | - Vector edgeVector3 = static_cast<Vector>(NormalVector(polygon[2], polygon[3], midPoint3, projection)); |
822 | | - Vector edgeVector4 = static_cast<Vector>(NormalVector(polygon[3], polygon[0], midPoint4, projection)); |
823 | | - |
824 | | - edgeVector1.normalise(); |
825 | | - edgeVector2.normalise(); |
826 | | - edgeVector3.normalise(); |
827 | | - edgeVector4.normalise(); |
828 | | - |
829 | | - Vector edgeVectorSum = edgeVector1 + edgeVector2 + edgeVector3 + edgeVector4; |
830 | | - |
831 | | - if (projection != Projection::cartesian) |
832 | | - { |
833 | | - edgeVectorSum.x() *= constants::conversion::radToDeg * constants::geometric::inverse_earth_radius; |
834 | | - } |
835 | | - |
836 | | - double edgeVectorSumLength = edgeVectorSum.length(); |
837 | | - |
838 | | - for (UInt i = 1; i <= MaximumNumberOfCircumcentreIterations; ++i) |
839 | | - { |
840 | | - Vector delta1 = GetDelta(midPoint1, centreOfMass, projection); |
841 | | - Vector delta2 = GetDelta(midPoint2, centreOfMass, projection); |
842 | | - Vector delta3 = GetDelta(midPoint3, centreOfMass, projection); |
843 | | - Vector delta4 = GetDelta(midPoint4, centreOfMass, projection); |
844 | | - |
845 | | - double ds = dot(delta1, edgeVector1) + dot(delta2, edgeVector2) + dot(delta3, edgeVector3) + dot(delta4, edgeVector4); |
846 | | - |
847 | | - if (projection != Projection::cartesian) |
848 | | - { |
849 | | - ds *= constants::conversion::radToDeg * constants::geometric::inverse_earth_radius; |
850 | | - } |
851 | | - |
852 | | - centreOfMass.x -= updateStepSize * ds * edgeVectorSum.x() * xTransformation; |
853 | | - centreOfMass.y -= updateStepSize * ds * edgeVectorSum.y(); |
854 | | - |
855 | | - if (ds * edgeVectorSumLength < circumcentreTolerance) |
856 | | - { |
857 | | - break; |
858 | | - } |
859 | | - } |
860 | | - } |
861 | | - else |
862 | | - { |
863 | | - for (UInt j = 1; j <= MaximumNumberOfCircumcentreIterations; ++j) |
864 | | - { |
865 | | - Vector edgeVectorSum(0.0, 0.0); |
866 | | - double ds = 0.0; |
867 | | - |
868 | | - for (UInt i = 0; i < numberOfPointsOpenedPolygon; ++i) |
869 | | - { |
870 | | - const auto nextNode = NextCircularForwardIndex(i, numberOfPointsOpenedPolygon); |
871 | | - |
872 | | - Point midPoint = 0.5 * (polygon[i] + polygon[nextNode]); |
873 | | - Vector edgeVector = static_cast<Vector>(NormalVector(polygon[i], polygon[nextNode], midPoint, projection)); |
874 | | - Vector delta = GetDelta(midPoint, centreOfMass, projection); |
875 | | - |
876 | | - edgeVector.normalise(); |
877 | | - ds += dot(delta, edgeVector); |
878 | | - edgeVectorSum += edgeVector; |
879 | | - } |
880 | | - |
881 | | - if (projection != Projection::cartesian) |
882 | | - { |
883 | | - ds *= constants::conversion::radToDeg * constants::geometric::inverse_earth_radius; |
884 | | - } |
885 | | - |
886 | | - centreOfMass.x -= updateStepSize * ds * edgeVectorSum.x() * xTransformation; |
887 | | - centreOfMass.y -= updateStepSize * ds * edgeVectorSum.y(); |
888 | | - |
889 | | - if (ds * edgeVectorSum.length() < circumcentreTolerance) |
890 | | - { |
891 | | - break; |
892 | | - } |
893 | | - } |
894 | | - } |
| 756 | + Point centreOfMass(algo::ComputeMassCentre(polygon)); |
895 | 757 |
|
896 | 758 | area = ComputeArea(polygon, projection); |
897 | 759 | TraversalDirection direction = area > 0.0 ? TraversalDirection::AntiClockwise : TraversalDirection::Clockwise; |
|
0 commit comments