From 3e81c09e214e1c18e00a030343b03c4bb8819b81 Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Sat, 22 Jul 2023 17:08:36 +0200 Subject: [PATCH 1/7] revert to thrust 1.16 --- cmake/thrust.cmake | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cmake/thrust.cmake b/cmake/thrust.cmake index 4c6ebb8..a7ef8bf 100644 --- a/cmake/thrust.cmake +++ b/cmake/thrust.cmake @@ -11,7 +11,7 @@ set(THRUST_ENABLE_EXAMPLES "OFF") # Set standard CPP Dialect to 17 (default of thrust would be 14) set(THRUST_CPP_DIALECT 17) -find_package(Thrust 2.1.0 QUIET) +find_package(Thrust 1.16.0 QUIET) if (${Thrust_FOUND}) @@ -19,10 +19,10 @@ if (${Thrust_FOUND}) else() message(STATUS "Using thrust from git repository") - # Fetches the version 2.1.0 of the official NVIDIA Thrust repository + # Fetches the version 1.16.0 of the official NVIDIA Thrust repository FetchContent_Declare(thrust GIT_REPOSITORY https://github.com/NVIDIA/thrust.git - GIT_TAG 2.1.0 + GIT_TAG 1.16.0 ) FetchContent_MakeAvailable(thrust) endif() From e92431d21797bea80bf5a11a778e146020b81bf7 Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Sat, 22 Jul 2023 18:39:28 +0200 Subject: [PATCH 2/7] add missing includes --- src/polyhedralGravityPython/PolyhedralGravityPython.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp index 7160587..9f5c05c 100644 --- a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp +++ b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp @@ -1,5 +1,8 @@ #include #include +#include +#include +#include #include "pybind11/pybind11.h" #include "pybind11/stl.h" From 33d14984b094bd2dffdcc706847d8d0633738343 Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Sat, 22 Jul 2023 19:19:52 +0200 Subject: [PATCH 3/7] relax xsimd --- cmake/xsimd.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/xsimd.cmake b/cmake/xsimd.cmake index dd3e820..7fa01a1 100644 --- a/cmake/xsimd.cmake +++ b/cmake/xsimd.cmake @@ -3,7 +3,7 @@ include(FetchContent) message(STATUS "Setting up xsimd via CMake") -find_package(xsimd 11.1.0 QUIET) +find_package(xsimd 11.1 QUIET) if (${xsimd_FOUND}) From eb95898250e730878ad31db3eb72b9f20520769f Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Mon, 24 Jul 2023 20:26:25 +0200 Subject: [PATCH 4/7] option _LIBCPP_DISABLE_AVAILABILITY --- CMakeLists.txt | 5 +++++ setup.py | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4148467..06e8a93 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,6 +46,11 @@ option(BUILD_POLYHEDRAL_GRAVITY_TESTS "Set to on if the tests should be built (D message(STATUS "BUILD_POLYHEDRAL_GRAVITY_TESTS = ${BUILD_POLYHEDRAL_GRAVITY_TESTS}") +IF(_LIBCPP_DISABLE_AVAILABILITY) + message(STATUS "Disabling availability macros for libc++") + add_definitions(-D_LIBCPP_DISABLE_AVAILABILITY) +endif () + ####################################################### # Including dependencies needed across multiple targets ####################################################### diff --git a/setup.py b/setup.py index 980b820..36ddaf5 100644 --- a/setup.py +++ b/setup.py @@ -87,6 +87,12 @@ def build_extension(self, ext): f"-D{option_name}={final_value}" ] + # Disable availability of standard libc++ on macOS if requested + if os.environ.get("_LIBCPP_DISABLE_AVAILABILITY"): + cmake_args += [ + "-D_LIBCPP_DISABLE_AVAILABILITY=ON" + ] + # Sets the CMake Generator if specified (this is separate from the other variables since it is given to # CMake vie the -G prefix final_generator = os.environ.get("CMAKE_GENERATOR", CMAKE_GENERATOR) From 00d20ab09f52c444f2f0a2e3f52dedcf31d76236 Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Mon, 31 Jul 2023 20:00:03 +0200 Subject: [PATCH 5/7] resolve cyclic dependency --- .../calculation/GravityEvaluable.h | 4 +- .../calculation/GravityModel.cpp | 394 ----------------- .../calculation/GravityModel.h | 200 +-------- .../calculation/GravityModelDetail.cpp | 395 ++++++++++++++++++ .../calculation/GravityModelDetail.h | 207 +++++++++ .../calculation/MeshChecking.h | 2 +- src/polyhedralGravity/model/Polyhedron.h | 2 +- 7 files changed, 608 insertions(+), 596 deletions(-) create mode 100644 src/polyhedralGravity/calculation/GravityModelDetail.cpp create mode 100644 src/polyhedralGravity/calculation/GravityModelDetail.h diff --git a/src/polyhedralGravity/calculation/GravityEvaluable.h b/src/polyhedralGravity/calculation/GravityEvaluable.h index 9ebd704..1b053ca 100644 --- a/src/polyhedralGravity/calculation/GravityEvaluable.h +++ b/src/polyhedralGravity/calculation/GravityEvaluable.h @@ -6,12 +6,12 @@ #include #include "thrust/transform.h" -#include "polyhedralGravity/calculation/GravityModel.h" +#include "thrust/execution_policy.h" +#include "polyhedralGravity/calculation/GravityModelDetail.h" #include "polyhedralGravity/input/TetgenAdapter.h" #include "polyhedralGravity/model/GravityModelData.h" #include "polyhedralGravity/model/Polyhedron.h" -#include "thrust/execution_policy.h" namespace polyhedralGravity { diff --git a/src/polyhedralGravity/calculation/GravityModel.cpp b/src/polyhedralGravity/calculation/GravityModel.cpp index a88fcdf..e95735e 100644 --- a/src/polyhedralGravity/calculation/GravityModel.cpp +++ b/src/polyhedralGravity/calculation/GravityModel.cpp @@ -21,398 +21,4 @@ namespace polyhedralGravity::GravityModel { }, polyhedron)); } - - Array3Triplet detail::buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2) { - using util::operator-; - //Calculate G_ij - return {vertex1 - vertex0, vertex2 - vertex1, vertex0 - vertex2}; - } - - Array3 detail::buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2) { - //Calculate N_i as (G_i1 * G_i2) / |G_i1 * G_i2| with * being the cross product - return util::normal(segmentVector1, segmentVector2); - } - - Array3Triplet - detail::buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal) { - Array3Triplet segmentUnitNormal{}; - //Calculate n_ij as (G_ij * N_i) / |G_ig * N_i| with * being the cross product - std::transform(segmentVectors.cbegin(), segmentVectors.end(), segmentUnitNormal.begin(), - [&planeUnitNormal](const Array3 &segmentVector) -> Array3 { - return util::normal(segmentVector, planeUnitNormal); - }); - return segmentUnitNormal; - } - - double detail::computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0) { - using namespace util; - //Calculate N_i * -G_i1 where * is the dot product and then use the inverted sgn - //We abstain on the double multiplication with -1 in the line above and beyond since two - //times multiplying with -1 equals no change - return sgn(dot(planeUnitNormal, vertex0), util::EPSILON); - } - - HessianPlane detail::computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r) { - using namespace util; - constexpr Array3 origin{0.0, 0.0, 0.0}; - const auto crossProduct = cross(p - q, p - r); - const auto res = (origin - p) * crossProduct; - const double d = res[0] + res[1] + res[2]; - - return {crossProduct[0], crossProduct[1], crossProduct[2], d}; - } - - double detail::distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane) { - //Compute h_p as D/sqrt(A^2 + B^2 + C^2) - return std::abs(hessianPlane.d / std::sqrt( - hessianPlane.a * hessianPlane.a + hessianPlane.b * hessianPlane.b + hessianPlane.c * hessianPlane.c)); - } - - Array3 detail::projectPointOrthogonallyOntoPlane(const Array3 &planeUnitNormal, double planeDistance, - const HessianPlane &hessianPlane) { - using namespace util; - //Calculate the projection point by (22) P'_ = N_i / norm(N_i) * h_i - // norm(N_i) is always 1 since N_i is a "normed" vector --> we do not need this division - Array3 orthogonalProjectionPoint = planeUnitNormal * planeDistance; - - //Calculate alpha, beta and gamma as D/A, D/B and D/C (Notice that we "forget" the minus before those - // divisions. In consequence, the conditions for signs are reversed below!!!) - // These values represent the intersections of each polygonal plane with the axes - // Comparison x == 0.0 is ok, since we only want to avoid nan values - Array3 intersections = {hessianPlane.a == 0.0 ? 0.0 : hessianPlane.d / hessianPlane.a, - hessianPlane.b == 0.0 ? 0.0 : hessianPlane.d / hessianPlane.b, - hessianPlane.c == 0.0 ? 0.0 : hessianPlane.d / hessianPlane.c}; - - //Determine the signs of the coordinates of P' according to the intersection values alpha, beta, gamma - // denoted as __ below, i.e. -alpha, -beta, -gamma denoted -__ - for (unsigned int index = 0; index < 3; ++index) { - if (intersections[index] < 0) { - //If -__ >= 0 --> __ < 0 then coordinates are positive, we calculate abs(orthogonalProjectionPoint[..]) - orthogonalProjectionPoint[index] = std::abs(orthogonalProjectionPoint[index]); - } else { - //The coordinates need to be controlled - if (planeUnitNormal[index] > 0) { - //If -__ < 0 --> __ >= 0 then the coordinate is negative -orthogonalProjectionPoint[..] - orthogonalProjectionPoint[index] = -1.0 * orthogonalProjectionPoint[index]; - } else { - //Else the coordinate is positive orthogonalProjectionPoint[..] - orthogonalProjectionPoint[index] = orthogonalProjectionPoint[index]; - } - } - } - return orthogonalProjectionPoint; - } - - Array3 - detail::computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices, const Array3 &projectionPointOnPlane, - const Array3Triplet &segmentUnitNormalsForPlane) { - using namespace util; - std::array segmentNormalOrientations{}; - //Equation (23) - //Calculate x_P' - x_ij^1 (x_P' is the projectionPoint and x_ij^1 is the first vertices of one segment, - //i.e. the coordinates of the training-planes' nodes --> projectionPointOnPlane - vertex - //Calculate n_ij * x_ij with * being the dot product and use the inverted sgn to determine the value of sigma_pq - //--> sgn((dot(segmentNormalOrientation, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0 - std::transform(segmentUnitNormalsForPlane.cbegin(), segmentUnitNormalsForPlane.cend(), vertices.cbegin(), - segmentNormalOrientations.begin(), - [&projectionPointOnPlane](const Array3 &segmentUnitNormal, const Array3 &vertex) { - using namespace util; - return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0; - }); - return segmentNormalOrientations; - } - - Array3Triplet detail::projectPointOrthogonallyOntoSegments(const Array3 &projectionPointOnPlane, - const Array3 &segmentNormalOrientations, - const Array3Triplet &face) { - auto counterJ = thrust::counting_iterator(0); - Array3Triplet orthogonalProjectionPointOnSegmentPerPlane{}; - - //Running over the segments of this plane - thrust::transform(segmentNormalOrientations.begin(), segmentNormalOrientations.end(), counterJ, - orthogonalProjectionPointOnSegmentPerPlane.begin(), - [&](const double &segmentNormalOrientation, const unsigned int j) { - //We actually only accept +0.0 or -0.0, so the equal comparison is ok - if (segmentNormalOrientation == 0.0) { - //Geometrically trivial case, in neither of the half space --> already on segment - return projectionPointOnPlane; - } else { - //In one of the half space, evaluate the projection point P'' for the segment - //with the endpoints v1 and v2 - const auto &vertex1 = face[j]; - const auto &vertex2 = face[(j + 1) % 3]; - return projectPointOrthogonallyOntoSegment(vertex1, vertex2, projectionPointOnPlane); - } - }); - return orthogonalProjectionPointOnSegmentPerPlane; - } - - Array3 detail::projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2, - const Array3 &orthogonalProjectionPointOnPlane) { - using namespace util; - //Preparing our the planes/ equations in matrix form - const Array3 matrixRow1 = vertex2 - vertex1; - const Array3 matrixRow2 = cross(vertex1 - orthogonalProjectionPointOnPlane, matrixRow1); - const Array3 matrixRow3 = cross(matrixRow2, matrixRow1); - const Array3 d = {dot(matrixRow1, orthogonalProjectionPointOnPlane), - dot(matrixRow2, orthogonalProjectionPointOnPlane), dot(matrixRow3, vertex1)}; - Matrix columnMatrix = transpose(Matrix{matrixRow1, matrixRow2, matrixRow3}); - //Calculation and solving the equations of above - const double determinant = det(columnMatrix); - return Array3{det(Matrix{d, columnMatrix[1], columnMatrix[2]}), - det(Matrix{columnMatrix[0], d, columnMatrix[2]}), - det(Matrix{columnMatrix[0], columnMatrix[1], d})} / determinant; - } - - Array3 detail::distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane, - const Array3Triplet &orthogonalProjectionPointOnSegments) { - std::array segmentDistances{}; - //The inner loop with the running j --> iterating over the segments - //Using the values P'_i and P''_ij for the calculation of the distance - std::transform(orthogonalProjectionPointOnSegments.cbegin(), orthogonalProjectionPointOnSegments.cend(), - segmentDistances.begin(), [&](const Array3 &orthogonalProjectionPointOnSegment) { - using namespace util; - return euclideanNorm(orthogonalProjectionPointOnSegment - orthogonalProjectionPointOnPlane); - }); - return segmentDistances; - } - - std::array detail::distancesToSegmentEndpoints(const Array3Triplet &segmentVectorsForPlane, - const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane, - const Array3Triplet &face) { - std::array distancesForPlane{}; - auto counter = thrust::counting_iterator(0); - auto zip = util::zipPair(segmentVectorsForPlane, orthogonalProjectionPointsOnSegmentForPlane); - - thrust::transform(zip.first, zip.second, counter, distancesForPlane.begin(), - [&face](const auto &tuple, unsigned int j) { - using namespace util; - Distance distance{}; - //segment vector G_pq - const Array3 &segmentVector = thrust::get<0>(tuple); - //orthogonal projection point on segment P'' for plane p and segment q - const Array3 &orthogonalProjectionPointsOnSegment = thrust::get<1>(tuple); - - //Calculate the 3D distances between P (0, 0, 0) and - // the segment endpoints face[j] and face[(j + 1) % 3]) - distance.l1 = euclideanNorm(face[j]); - distance.l2 = euclideanNorm(face[(j + 1) % 3]); - //Calculate the 1D distances between P'' (every segment has its own) and - // the segment endpoints face[j] and face[(j + 1) % 3]) - distance.s1 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[j]); - distance.s2 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[(j + 1) % 3]); - - /* - * Additional remark: - * Details on these conditions are in the second paper referenced in the README.md (Tsoulis, 2021) - * The numbering of these conditions is equal to the numbering scheme of the paper - * Assign a sign to those magnitudes depending on the relative position of P'' to the two - * segment endpoints - */ - - //4. Option: |s1 - l1| == 0 && |s2 - l2| == 0 Computation point P is located from the beginning on - // the direction of a specific segment (P coincides with P' and P'') - if (std::abs(distance.s1 - distance.l1) < EPSILON && - std::abs(distance.s2 - distance.l2) < EPSILON) { - //4. Option - Case 2: P is located on the segment from its right side - // s1 = -|s1|, s2 = -|s2|, l1 = -|l1|, l2 = -|l2| - if (distance.s2 < distance.s1) { - distance.s1 *= -1.0; - distance.s2 *= -1.0; - distance.l1 *= -1.0; - distance.l2 *= -1.0; - return distance; - } else if (std::abs(distance.s2 - distance.s1) < EPSILON) { - //4. Option - Case 1: P is located inside the segment (s2 == s1) - // s1 = -|s1|, s2 = |s2|, l1 = -|l1|, l2 = |l2| - distance.s1 *= -1.0; - distance.l1 *= -1.0; - return distance; - } - //4. Option - Case 3: P is located on the segment from its left side - // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do! - } else { - const double norm = euclideanNorm(segmentVector); - if (distance.s1 < norm && distance.s2 < norm) { - //1. Option: |s1| < |G_ij| && |s2| < |G_ij| Point P'' is situated inside the segment - // s1 = -|s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| - distance.s1 *= -1.0; - return distance; - } else if (distance.s2 < distance.s1) { - //2. Option: |s2| < |s1| Point P'' is on the right side of the segment - // s1 = -|s1|, s2 = -|s2|, l1 = |l1|, l2 = |l2| - distance.s1 *= -1.0; - distance.s2 *= -1.0; - return distance; - } - //3. Option: |s1| < |s2| Point P'' is on the left side of the segment - // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do! - } - return distance; - }); - return distancesForPlane; - } - - std::array - detail::computeTranscendentalExpressions(const std::array &distancesForPlane, double planeDistance, - const Array3 &segmentDistancesForPlane, - const Array3 &segmentNormalOrientationsForPlane, - const Array3 &projectionPointVertexNorms) { - std::array transcendentalExpressionsForPlane{}; - - //Zip iterator consisting of 3D and 1D distances l1/l2 and s1/2 for this plane | h_pq | sigma_pq for this plane - auto zip = util::zipPair(distancesForPlane, segmentDistancesForPlane, segmentNormalOrientationsForPlane); - auto counter = thrust::counting_iterator(0); - - thrust::transform(zip.first, zip.second, counter, transcendentalExpressionsForPlane.begin(), - [&](const auto &tuple, const unsigned int j) { - using namespace util; - //distances l1, l2, s1, s1 for this segment q of plane p - const Distance &distance = thrust::get<0>(tuple); - //segment distance h_pq for this segment q of plane p - const double segmentDistance = thrust::get<1>(tuple); - //segment normal orientation sigma_pq for this segment q of plane p - const double segmentNormalOrientation = thrust::get<2>(tuple); - - //Result for this segment - TranscendentalExpression transcendentalExpressionPerSegment{}; - - //Computation of the norm of P' and segment endpoints - // If the one of the norms == 0 then P' lies on the corresponding vertex and coincides with P'' - const double r1Norm = projectionPointVertexNorms[(j + 1) % 3]; - const double r2Norm = projectionPointVertexNorms[j]; - - //Compute LN_pq according to (14) - // If sigma_pq == 0 && either of the distances of P' to the two segment endpoints == 0 OR - // the 1D and 3D distances are smaller than some EPSILON - // then LN_pq can be set to zero - if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON || r2Norm < EPSILON)) || - (std::abs(distance.s1 + distance.s2) < EPSILON && - std::abs(distance.l1 + distance.l2) < EPSILON)) { - transcendentalExpressionPerSegment.ln = 0.0; - } else { - //Implementation of - // log((s2_pq + l2_pq) / (s1_pq + l1_pq)) - transcendentalExpressionPerSegment.ln = std::log( - (distance.s2 + distance.l2) / (distance.s1 + distance.l1)); - } - - //Compute AN_pq according to (15) - // If h_p == 0 or h_pq == 0 then AN_pq is zero, too (distances are always positive!) - if (planeDistance < EPSILON || segmentDistance < EPSILON) { - transcendentalExpressionPerSegment.an = 0.0; - } else { - //Implementation of: - // atan(h_p * s2_pq / h_pq * l2_pq) - atan(h_p * s1_pq / h_pq * l1_pq) - // in a vectorized manner to reduce the number of atan(..) and divisions(..) - - //The last subtraction is implemented via -distance.l1 (the minus/ inversion - // is sustained through all operations, even the atan(..) - xsimd::batch reg1(distance.s2, distance.s1); - xsimd::batch reg2(distance.l2, -distance.l1); - - reg1 = xsimd::mul(reg1, planeDistance); - reg2 = xsimd::mul(reg2, segmentDistance); - - reg1 = xsimd::div(reg1, reg2); - reg1 = xsimd::atan(reg1); - - // "Subtraction" masked as addition (l1 was previously inverted with a minus) - transcendentalExpressionPerSegment.an = xsimd::reduce_add(reg1); - } - - return transcendentalExpressionPerSegment; - }); - return transcendentalExpressionsForPlane; - } - - std::pair detail::computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane, - const Array3 &segmentNormalOrientationForPlane, - const Array3 &projectionPointVertexNorms, - const Array3 &planeUnitNormal, double planeDistance, - double planeNormalOrientation) { - //1. Case: If all sigma_pq for a given plane p are 1.0 then P' lies inside the plane S_p - if (std::all_of(segmentNormalOrientationForPlane.cbegin(), segmentNormalOrientationForPlane.cend(), - [](const double sigma) { return sigma == 1.0; })) { - using namespace util; - return std::make_pair( - -1.0 * util::PI2 * planeDistance, //sing alpha = -2pi*h_p - planeUnitNormal * (-1.0 * util::PI2 * planeNormalOrientation)); //sing beta = -2pi*sigma_p*N_p - } - //2. Case: If sigma_pq == 0 AND norm(P' - v1) < norm(G_ij) && norm(P' - v2) < norm(G_ij) with G_ij - // as the vector of v1 and v2 - // then P' is located on one line segment G_p of plane p, but not on any of its vertices - auto counter2 = thrust::counting_iterator(0); - auto secondCaseBegin = util::zip(segmentVectorsForPlane.begin(), segmentNormalOrientationForPlane.begin(), - counter2); - auto secondCaseEnd = util::zip(segmentVectorsForPlane.end(), segmentNormalOrientationForPlane.end(), - counter2 + 3); - if (std::any_of(secondCaseBegin, secondCaseEnd, [&](const auto &tuple) { - using namespace util; - const Array3 &segmentVector = thrust::get<0>(tuple); - const double segmentNormalOrientation = thrust::get<1>(tuple); - const unsigned int j = thrust::get<2>(tuple); - - //segmentNormalOrientation != 0.0 - if (std::abs(segmentNormalOrientation) > EPSILON) { - return false; - } - - const double segmentVectorNorm = euclideanNorm(segmentVector); - return projectionPointVertexNorms[(j + 1) % 3] < segmentVectorNorm && - projectionPointVertexNorms[j] < segmentVectorNorm; - })) { - using namespace util; - return std::make_pair(-1.0 * util::PI * planeDistance, //sing alpha = -pi*h_p - planeUnitNormal * - (-1.0 * util::PI * planeNormalOrientation)); //sing beta = -pi*sigma_p*N_p - } - //3. Case If sigma_pq == 0 AND norm(P' - v1) < 0 || norm(P' - v2) < 0 - // then P' is located at one of G_p's vertices - auto counter3 = thrust::counting_iterator(0); - auto thirdCaseBegin = util::zip(segmentNormalOrientationForPlane.begin(), counter3); - auto thirdCaseEnd = util::zip(segmentNormalOrientationForPlane.end(), counter3 + 3); - double r1Norm; - double r2Norm; - unsigned int j; - if (std::any_of(thirdCaseBegin, thirdCaseEnd, [&](const auto &tuple) { - using namespace util; - const double segmentNormalOrientation = thrust::get<0>(tuple); - j = thrust::get<1>(tuple); - - //segmentNormalOrientation != 0.0 - if (std::abs(segmentNormalOrientation) > EPSILON) { - return false; - } - - r1Norm = projectionPointVertexNorms[(j + 1) % 3]; - r2Norm = projectionPointVertexNorms[j]; - //r1Norm == 0.0 || r2Norm == 0.0 - return r1Norm < EPSILON || r2Norm < EPSILON; - })) { - using namespace util; - //Two segment vectors G_1 and G_2 of this plane - const Array3 &g1 = r1Norm == 0.0 ? segmentVectorsForPlane[j] : segmentVectorsForPlane[(j - 1 + 3) % 3]; - const Array3 &g2 = r1Norm == 0.0 ? segmentVectorsForPlane[(j + 1) % 3] : segmentVectorsForPlane[j]; - // theta = arcos((G_2 * -G_1) / (|G_2| * |G_1|)) - const double gdot = dot(g1 * -1.0, g2); - const double theta = gdot == 0.0 ? util::PI_2 : std::acos(gdot / (euclideanNorm(g1) * euclideanNorm(g2))); - return std::make_pair(-1.0 * theta * planeDistance, //sing alpha = -theta*h_p - planeUnitNormal * - (-1.0 * theta * planeNormalOrientation)); //sing beta = -theta*sigma_p*N_p - } - //4. Case Otherwise P' is located outside the plane S_p and then the singularity equals zero - return std::make_pair(0.0, //sing alpha = 0 - Array3{0.0, 0.0, 0.0}); //sing beta = 0 - } - - Array3 detail::computeNormsOfProjectionPointAndVertices(const Array3 &orthogonalProjectionPointOnPlane, - const Array3Triplet &face) { - using namespace util; - return {euclideanNorm(orthogonalProjectionPointOnPlane - face[0]), - euclideanNorm(orthogonalProjectionPointOnPlane - face[1]), - euclideanNorm(orthogonalProjectionPointOnPlane - face[2])}; - } - - } diff --git a/src/polyhedralGravity/calculation/GravityModel.h b/src/polyhedralGravity/calculation/GravityModel.h index 021d71d..adc8b01 100644 --- a/src/polyhedralGravity/calculation/GravityModel.h +++ b/src/polyhedralGravity/calculation/GravityModel.h @@ -6,23 +6,9 @@ #include #include -#include "polyhedralGravity/input/TetgenAdapter.h" -#include "polyhedralGravity/model/Polyhedron.h" -#include "polyhedralGravity/model/GravityModelData.h" #include "polyhedralGravity/calculation/GravityEvaluable.h" -#include "polyhedralGravity/calculation/PolyhedronTransform.h" -#include "polyhedralGravity/util/UtilityConstants.h" -#include "polyhedralGravity/util/UtilityContainer.h" -#include "spdlog/spdlog.h" -#include "thrust/iterator/zip_iterator.h" -#include "thrust/iterator/transform_iterator.h" -#include "thrust/iterator/counting_iterator.h" -#include "thrust/transform.h" -#include "thrust/transform_reduce.h" -#include "thrust/execution_policy.h" -#include "polyhedralGravity/util/UtilityThrust.h" -#include "polyhedralGravity/output/Logging.h" -#include "xsimd/xsimd.hpp" +#include "polyhedralGravity/model/GravityModelData.h" +#include "polyhedralGravity/model/Polyhedron.h" /** * Namespace containing the methods used to evaluate the polyhedrale Gravity Model @@ -57,186 +43,4 @@ namespace polyhedralGravity::GravityModel { std::vector evaluate(const PolyhedronOrSource &polyhedron, double density, const std::vector &computationPoints, bool parallel = true); - namespace detail { - - /** - * Computes the segment vectors G_ij for one plane of the polyhedron according to Tsoulis (18). - * The segment vectors G_ij represent the vector from one vertex of the face to the neighboring vertex and - * depict every line segment of the triangular face (A-B-C) - * @param vertex0 - the first vertex A - * @param vertex1 - the second vertex B - * @param vertex2 - the third vertex C - * @return the segment vectors for a plane - */ - Array3Triplet buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2); - - /** - * Computes the plane unit normal N_p for one plane p of the polyhedron according to Tsoulis (19). - * The plane unit normal is the outward pointing normal of the face from the polyhedron. - * @param segmentVector1 - first edge - * @param segmentVector2 - second edge - * @return plane unit normal - */ - Array3 buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2); - - /** - * Computes the segment unit normals n_pq for one plane p of the polyhedron according to Tsoulis (20). - * The segment unit normal n_pq represent the normal of one line segment of a polyhedrale face. - * @param segmentVectors - the segment vectors of the face G_p(0-2) - * @param planeUnitNormal - the plane unit normal N_p - * @return segment unit normals n_pq for plane p with q = {0, 1, 2} - */ - Array3Triplet buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal); - - /** - * Computes the plane unit normal orientation/ direction sigma_p for one plane p of the polyhedron - * according to Tsoulis (21). - * The plane unit normal orientation values represents the relative position of computation point P - * with respect to the pointing direction of N_p. E. g. if N_p points to the half-space containing P, the - * inner product of N_p and -G_i1 will be positive, leading to a negative sigma_p. - * If sigma_p is zero than P and P' lie geometrically in the same plane --> P == P'. - * @param planeUnitNormal - the plane unit normal N_p - * @param vertex0 - the first vertex of the plane - * @return plane normal orientation - */ - double computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0); - - /** - * Calculates the Hessian Plane form spanned by three given points p, q, and r. - * @param p - first point on the plane - * @param q - second point on the plane - * @param r - third point on the plane - * @return HessianPlane - * @related Cross-Product method https://tutorial.math.lamar.edu/classes/calciii/eqnsofplanes.aspx - */ - HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r); - - /** - * Calculates the (plane) distances h_p of computation point P to the plane S_p given in Hessian Form - * according to the following equation: - * h_p = D / sqrt(A^2+B^2+C^2) - * @param hessianPlane - Hessian Plane Form of S_p - * @return plane distance h_p - */ - double distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane); - - /** - * Computes P' for a given plane p according to equation (22) of Tsoulis paper. - * P' is the orthogonal projection of the computation point P onto the plane S_p. - * @param planeUnitNormal - the plane unit normal N_p - * @param planeDistance - the distance from P to the plane h_p - * @param hessianPlane - the Hessian Plane Form - * @return P' for this plane - */ - Array3 projectPointOrthogonallyOntoPlane(const Array3 &planeUnitNormal, double planeDistance, - const HessianPlane &hessianPlane); - - /** - * Computes the segment normal orientations/ directions sigma_pq for a given plane p. - * If sigma_pq is negative, this denotes that n_pq points to the half-plane containing P'. Nn case - * sigma_pq is positive, P' resides in the other half-plane and if sigma_pq is zero, then P' lies directly - * on the segment pq. - * @param vertices - the vertices of this plane - * @param projectionPointOnPlane - the projection point P' for this plane - * @param segmentUnitNormalsForPlane - the segment unit normals sigma_pq for this plane - * @return the segment normal orientations for the plane p - */ - Array3 - computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices, const Array3 &projectionPointOnPlane, - const Array3Triplet &segmentUnitNormalsForPlane); - - /** - * Computes the orthogonal projection Points P'' foreach segment q of a given plane p. - * @param projectionPointOnPlane - the projection Point P' - * @param segmentNormalOrientations - the segment normal orientations sigma_pq for this plane p - * @param face - the vertices of the plane p - * @return the orthogonal projection points of P on the segment P'' foreach segment q of p - */ - Array3Triplet projectPointOrthogonallyOntoSegments(const Array3 &projectionPointOnPlane, - const Array3 &segmentNormalOrientations, - const Array3Triplet &face); - - /** - * Calculates the point P'' for a given Segment consisting of vertices v1 and v2 and the orthogonal projection - * point P' for the plane consisting of those vertices. Solves the three equations given in (24), (25) and (26). - * @param vertex1 - first endpoint of segment - * @param vertex2 - second endpoint of segment - * @param orthogonalProjectionPointOnPlane - the orthogonal projection P' of P on this plane - * @return P'' for this segment - * @note If sigma_pq is zero then P'' == P', this is not checked by this method, but has to be assured first - */ - Array3 projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2, - const Array3 &orthogonalProjectionPointOnPlane); - - /** - * Computes the (segment) distances h_pq between P' for a given plane p and P'' for a given segment q of plane p. - * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for p - * @param orthogonalProjectionPointOnSegments - the orthogonal projection points P'' for each segment q of p - * @return distances h_pq for plane p - */ - Array3 distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane, - const Array3Triplet &orthogonalProjectionPointOnSegments); - - /** - * Computes the 3D distances l1_pq and l2_pq between the computation point P and the line - * segment endpoints of each polyhedral segment for one plane. - * Computes the 1D distances s1_pq and s2_pq between orthogonal projection of P on the line - * segment P''_pq and the line segment endpoints for each polyhedral segment for one plane - * @param segmentVectorsForPlane - the segment vectors G_pq for plane p - * @param orthogonalProjectionPointsOnSegmentForPlane - the orthogonal projection Points P'' for plane p - * @param face - the vertices of plane p - * @return distances l1_pq and l2_pq and s1_pq and s2_pq foreach segment q of plane p - */ - std::array distancesToSegmentEndpoints(const Array3Triplet &segmentVectorsForPlane, - const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane, - const Array3Triplet &face); - - /** - * Calculates the Transcendental Expressions LN_pq and AN_pq for every line segment of the polyhedron for - * a given plane p. - * LN_pq is calculated according to (14) using the natural logarithm and AN_pq is calculated according - * to (15) using the arctan. - * @param distancesForPlane - the distances l1, l2, s1, s2 foreach segment q of plane p - * @param planeDistance - the plane distance h_p for plane p - * @param segmentDistancesForPlane - the segment distance h_pq for segment q of plane p - * @param segmentNormalOrientationsForPlane - the segment normal orientations n_pq for a plane p - * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for plane p - * @param face - the vertices of plane p - * @return LN_pq and AN_pq foreach segment q of plane p - */ - std::array - computeTranscendentalExpressions(const std::array &distancesForPlane, double planeDistance, - const Array3 &segmentDistancesForPlane, - const Array3 &segmentNormalOrientationsForPlane, - const Array3 &projectionPointVertexNorms); - - /** - * Calculates the singularities (correction) terms according to the Flow text for a given plane p. - * @param segmentVectorsForPlane - the segment vectors for a given plane - * @param segmentNormalOrientationForPlane - the segment orientation sigma_pq - * @param projectionPointVertexNorms - the projection point P' - * @param planeUnitNormal - the plane unit normal N_p - * @param planeDistance - the plane distance h_p - * @param planeNormalOrientation - the plane normal orientation sigma_p - * @param face - the vertices of plane p - * @return the singularities for a plane p - */ - std::pair computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane, - const Array3 &segmentNormalOrientationForPlane, - const Array3 &projectionPointVertexNorms, - const Array3 &planeUnitNormal, double planeDistance, - double planeNormalOrientation); - - - /** - * Computes the L2 norms of the orthogonal projection point P' on a plane p with each vertex of that plane p. - * The values are later used to determine if P' is situated at a vertex. - * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' - * @param face - the vertices of plane p - * @return the norms of p and each vertex - */ - Array3 computeNormsOfProjectionPointAndVertices(const Array3 &orthogonalProjectionPointOnPlane, - const Array3Triplet &face); - - } } diff --git a/src/polyhedralGravity/calculation/GravityModelDetail.cpp b/src/polyhedralGravity/calculation/GravityModelDetail.cpp new file mode 100644 index 0000000..6c0038b --- /dev/null +++ b/src/polyhedralGravity/calculation/GravityModelDetail.cpp @@ -0,0 +1,395 @@ +#include "GravityModelDetail.h" + +namespace polyhedralGravity::GravityModel::detail { + + Array3Triplet buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2) { + using util::operator-; + //Calculate G_ij + return {vertex1 - vertex0, vertex2 - vertex1, vertex0 - vertex2}; + } + + Array3 buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2) { + //Calculate N_i as (G_i1 * G_i2) / |G_i1 * G_i2| with * being the cross product + return util::normal(segmentVector1, segmentVector2); + } + + Array3Triplet buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal) { + Array3Triplet segmentUnitNormal{}; + //Calculate n_ij as (G_ij * N_i) / |G_ig * N_i| with * being the cross product + std::transform(segmentVectors.cbegin(), segmentVectors.end(), segmentUnitNormal.begin(), + [&planeUnitNormal](const Array3 &segmentVector) -> Array3 { + return util::normal(segmentVector, planeUnitNormal); + }); + return segmentUnitNormal; + } + + double computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0) { + using namespace util; + //Calculate N_i * -G_i1 where * is the dot product and then use the inverted sgn + //We abstain on the double multiplication with -1 in the line above and beyond since two + //times multiplying with -1 equals no change + return sgn(dot(planeUnitNormal, vertex0), util::EPSILON); + } + + HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r) { + using namespace util; + constexpr Array3 origin{0.0, 0.0, 0.0}; + const auto crossProduct = cross(p - q, p - r); + const auto res = (origin - p) * crossProduct; + const double d = res[0] + res[1] + res[2]; + + return {crossProduct[0], crossProduct[1], crossProduct[2], d}; + } + + double distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane) { + //Compute h_p as D/sqrt(A^2 + B^2 + C^2) + return std::abs(hessianPlane.d / std::sqrt( + hessianPlane.a * hessianPlane.a + hessianPlane.b * hessianPlane.b + hessianPlane.c * hessianPlane.c)); + } + + Array3 projectPointOrthogonallyOntoPlane(const Array3 &planeUnitNormal, double planeDistance, + const HessianPlane &hessianPlane) { + using namespace util; + //Calculate the projection point by (22) P'_ = N_i / norm(N_i) * h_i + // norm(N_i) is always 1 since N_i is a "normed" vector --> we do not need this division + Array3 orthogonalProjectionPoint = planeUnitNormal * planeDistance; + + //Calculate alpha, beta and gamma as D/A, D/B and D/C (Notice that we "forget" the minus before those + // divisions. In consequence, the conditions for signs are reversed below!!!) + // These values represent the intersections of each polygonal plane with the axes + // Comparison x == 0.0 is ok, since we only want to avoid nan values + Array3 intersections = {hessianPlane.a == 0.0 ? 0.0 : hessianPlane.d / hessianPlane.a, + hessianPlane.b == 0.0 ? 0.0 : hessianPlane.d / hessianPlane.b, + hessianPlane.c == 0.0 ? 0.0 : hessianPlane.d / hessianPlane.c}; + + //Determine the signs of the coordinates of P' according to the intersection values alpha, beta, gamma + // denoted as __ below, i.e. -alpha, -beta, -gamma denoted -__ + for (unsigned int index = 0; index < 3; ++index) { + if (intersections[index] < 0) { + //If -__ >= 0 --> __ < 0 then coordinates are positive, we calculate abs(orthogonalProjectionPoint[..]) + orthogonalProjectionPoint[index] = std::abs(orthogonalProjectionPoint[index]); + } else { + //The coordinates need to be controlled + if (planeUnitNormal[index] > 0) { + //If -__ < 0 --> __ >= 0 then the coordinate is negative -orthogonalProjectionPoint[..] + orthogonalProjectionPoint[index] = -1.0 * orthogonalProjectionPoint[index]; + } else { + //Else the coordinate is positive orthogonalProjectionPoint[..] + orthogonalProjectionPoint[index] = orthogonalProjectionPoint[index]; + } + } + } + return orthogonalProjectionPoint; + } + + Array3 computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices, const Array3 &projectionPointOnPlane, + const Array3Triplet &segmentUnitNormalsForPlane) { + using namespace util; + std::array segmentNormalOrientations{}; + //Equation (23) + //Calculate x_P' - x_ij^1 (x_P' is the projectionPoint and x_ij^1 is the first vertices of one segment, + //i.e. the coordinates of the training-planes' nodes --> projectionPointOnPlane - vertex + //Calculate n_ij * x_ij with * being the dot product and use the inverted sgn to determine the value of sigma_pq + //--> sgn((dot(segmentNormalOrientation, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0 + std::transform(segmentUnitNormalsForPlane.cbegin(), segmentUnitNormalsForPlane.cend(), vertices.cbegin(), + segmentNormalOrientations.begin(), + [&projectionPointOnPlane](const Array3 &segmentUnitNormal, const Array3 &vertex) { + using namespace util; + return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0; + }); + return segmentNormalOrientations; + } + + Array3Triplet + projectPointOrthogonallyOntoSegments(const Array3 &projectionPointOnPlane, const Array3 &segmentNormalOrientations, + const Array3Triplet &face) { + auto counterJ = thrust::counting_iterator(0); + Array3Triplet orthogonalProjectionPointOnSegmentPerPlane{}; + + //Running over the segments of this plane + thrust::transform(segmentNormalOrientations.begin(), segmentNormalOrientations.end(), counterJ, + orthogonalProjectionPointOnSegmentPerPlane.begin(), + [&](const double &segmentNormalOrientation, const unsigned int j) { + //We actually only accept +0.0 or -0.0, so the equal comparison is ok + if (segmentNormalOrientation == 0.0) { + //Geometrically trivial case, in neither of the half space --> already on segment + return projectionPointOnPlane; + } else { + //In one of the half space, evaluate the projection point P'' for the segment + //with the endpoints v1 and v2 + const auto &vertex1 = face[j]; + const auto &vertex2 = face[(j + 1) % 3]; + return projectPointOrthogonallyOntoSegment(vertex1, vertex2, projectionPointOnPlane); + } + }); + return orthogonalProjectionPointOnSegmentPerPlane; + } + + Array3 projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2, + const Array3 &orthogonalProjectionPointOnPlane) { + using namespace util; + //Preparing our the planes/ equations in matrix form + const Array3 matrixRow1 = vertex2 - vertex1; + const Array3 matrixRow2 = cross(vertex1 - orthogonalProjectionPointOnPlane, matrixRow1); + const Array3 matrixRow3 = cross(matrixRow2, matrixRow1); + const Array3 d = {dot(matrixRow1, orthogonalProjectionPointOnPlane), + dot(matrixRow2, orthogonalProjectionPointOnPlane), dot(matrixRow3, vertex1)}; + Matrix columnMatrix = transpose(Matrix < double, 3, 3 > {matrixRow1, matrixRow2, matrixRow3}); + //Calculation and solving the equations of above + const double determinant = det(columnMatrix); + return Array3{det(Matrix < double, 3, 3 > {d, columnMatrix[1], columnMatrix[2]}), + det(Matrix < double, 3, 3 > {columnMatrix[0], d, columnMatrix[2]}), + det(Matrix < double, 3, 3 > {columnMatrix[0], columnMatrix[1], d})} / determinant; + } + + Array3 distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane, + const Array3Triplet &orthogonalProjectionPointOnSegments) { + std::array segmentDistances{}; + //The inner loop with the running j --> iterating over the segments + //Using the values P'_i and P''_ij for the calculation of the distance + std::transform(orthogonalProjectionPointOnSegments.cbegin(), orthogonalProjectionPointOnSegments.cend(), + segmentDistances.begin(), [&](const Array3 &orthogonalProjectionPointOnSegment) { + using namespace util; + return euclideanNorm(orthogonalProjectionPointOnSegment - orthogonalProjectionPointOnPlane); + }); + return segmentDistances; + } + + std::array distancesToSegmentEndpoints(const Array3Triplet &segmentVectorsForPlane, + const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane, + const Array3Triplet &face) { + std::array distancesForPlane{}; + auto counter = thrust::counting_iterator(0); + auto zip = util::zipPair(segmentVectorsForPlane, orthogonalProjectionPointsOnSegmentForPlane); + + thrust::transform(zip.first, zip.second, counter, distancesForPlane.begin(), + [&face](const auto &tuple, unsigned int j) { + using namespace util; + Distance distance{}; + //segment vector G_pq + const Array3 &segmentVector = thrust::get<0>(tuple); + //orthogonal projection point on segment P'' for plane p and segment q + const Array3 &orthogonalProjectionPointsOnSegment = thrust::get<1>(tuple); + + //Calculate the 3D distances between P (0, 0, 0) and + // the segment endpoints face[j] and face[(j + 1) % 3]) + distance.l1 = euclideanNorm(face[j]); + distance.l2 = euclideanNorm(face[(j + 1) % 3]); + //Calculate the 1D distances between P'' (every segment has its own) and + // the segment endpoints face[j] and face[(j + 1) % 3]) + distance.s1 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[j]); + distance.s2 = euclideanNorm(orthogonalProjectionPointsOnSegment - face[(j + 1) % 3]); + + /* + * Additional remark: + * Details on these conditions are in the second paper referenced in the README.md (Tsoulis, 2021) + * The numbering of these conditions is equal to the numbering scheme of the paper + * Assign a sign to those magnitudes depending on the relative position of P'' to the two + * segment endpoints + */ + + //4. Option: |s1 - l1| == 0 && |s2 - l2| == 0 Computation point P is located from the beginning on + // the direction of a specific segment (P coincides with P' and P'') + if (std::abs(distance.s1 - distance.l1) < EPSILON && + std::abs(distance.s2 - distance.l2) < EPSILON) { + //4. Option - Case 2: P is located on the segment from its right side + // s1 = -|s1|, s2 = -|s2|, l1 = -|l1|, l2 = -|l2| + if (distance.s2 < distance.s1) { + distance.s1 *= -1.0; + distance.s2 *= -1.0; + distance.l1 *= -1.0; + distance.l2 *= -1.0; + return distance; + } else if (std::abs(distance.s2 - distance.s1) < EPSILON) { + //4. Option - Case 1: P is located inside the segment (s2 == s1) + // s1 = -|s1|, s2 = |s2|, l1 = -|l1|, l2 = |l2| + distance.s1 *= -1.0; + distance.l1 *= -1.0; + return distance; + } + //4. Option - Case 3: P is located on the segment from its left side + // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do! + } else { + const double norm = euclideanNorm(segmentVector); + if (distance.s1 < norm && distance.s2 < norm) { + //1. Option: |s1| < |G_ij| && |s2| < |G_ij| Point P'' is situated inside the segment + // s1 = -|s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| + distance.s1 *= -1.0; + return distance; + } else if (distance.s2 < distance.s1) { + //2. Option: |s2| < |s1| Point P'' is on the right side of the segment + // s1 = -|s1|, s2 = -|s2|, l1 = |l1|, l2 = |l2| + distance.s1 *= -1.0; + distance.s2 *= -1.0; + return distance; + } + //3. Option: |s1| < |s2| Point P'' is on the left side of the segment + // s1 = |s1|, s2 = |s2|, l1 = |l1|, l2 = |l2| --> Nothing to do! + } + return distance; + }); + return distancesForPlane; + } + + std::array + computeTranscendentalExpressions(const std::array &distancesForPlane, double planeDistance, + const Array3 &segmentDistancesForPlane, + const Array3 &segmentNormalOrientationsForPlane, + const Array3 &projectionPointVertexNorms) { + std::array transcendentalExpressionsForPlane{}; + + //Zip iterator consisting of 3D and 1D distances l1/l2 and s1/2 for this plane | h_pq | sigma_pq for this plane + auto zip = util::zipPair(distancesForPlane, segmentDistancesForPlane, segmentNormalOrientationsForPlane); + auto counter = thrust::counting_iterator(0); + + thrust::transform(zip.first, zip.second, counter, transcendentalExpressionsForPlane.begin(), + [&](const auto &tuple, const unsigned int j) { + using namespace util; + //distances l1, l2, s1, s1 for this segment q of plane p + const Distance &distance = thrust::get<0>(tuple); + //segment distance h_pq for this segment q of plane p + const double segmentDistance = thrust::get<1>(tuple); + //segment normal orientation sigma_pq for this segment q of plane p + const double segmentNormalOrientation = thrust::get<2>(tuple); + + //Result for this segment + TranscendentalExpression transcendentalExpressionPerSegment{}; + + //Computation of the norm of P' and segment endpoints + // If the one of the norms == 0 then P' lies on the corresponding vertex and coincides with P'' + const double r1Norm = projectionPointVertexNorms[(j + 1) % 3]; + const double r2Norm = projectionPointVertexNorms[j]; + + //Compute LN_pq according to (14) + // If sigma_pq == 0 && either of the distances of P' to the two segment endpoints == 0 OR + // the 1D and 3D distances are smaller than some EPSILON + // then LN_pq can be set to zero + if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON || r2Norm < EPSILON)) || + (std::abs(distance.s1 + distance.s2) < EPSILON && + std::abs(distance.l1 + distance.l2) < EPSILON)) { + transcendentalExpressionPerSegment.ln = 0.0; + } else { + //Implementation of + // log((s2_pq + l2_pq) / (s1_pq + l1_pq)) + transcendentalExpressionPerSegment.ln = std::log( + (distance.s2 + distance.l2) / (distance.s1 + distance.l1)); + } + + //Compute AN_pq according to (15) + // If h_p == 0 or h_pq == 0 then AN_pq is zero, too (distances are always positive!) + if (planeDistance < EPSILON || segmentDistance < EPSILON) { + transcendentalExpressionPerSegment.an = 0.0; + } else { + //Implementation of: + // atan(h_p * s2_pq / h_pq * l2_pq) - atan(h_p * s1_pq / h_pq * l1_pq) + // in a vectorized manner to reduce the number of atan(..) and divisions(..) + + //The last subtraction is implemented via -distance.l1 (the minus/ inversion + // is sustained through all operations, even the atan(..) + xsimd::batch reg1(distance.s2, distance.s1); + xsimd::batch reg2(distance.l2, -distance.l1); + + reg1 = xsimd::mul(reg1, planeDistance); + reg2 = xsimd::mul(reg2, segmentDistance); + + reg1 = xsimd::div(reg1, reg2); + reg1 = xsimd::atan(reg1); + + // "Subtraction" masked as addition (l1 was previously inverted with a minus) + transcendentalExpressionPerSegment.an = xsimd::reduce_add(reg1); + } + + return transcendentalExpressionPerSegment; + }); + return transcendentalExpressionsForPlane; + } + + std::pair + computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane, const Array3 &segmentNormalOrientationForPlane, + const Array3 &projectionPointVertexNorms, const Array3 &planeUnitNormal, + double planeDistance, double planeNormalOrientation) { + //1. Case: If all sigma_pq for a given plane p are 1.0 then P' lies inside the plane S_p + if (std::all_of(segmentNormalOrientationForPlane.cbegin(), segmentNormalOrientationForPlane.cend(), + [](const double sigma) { return sigma == 1.0; })) { + using namespace util; + return std::make_pair( + -1.0 * util::PI2 * planeDistance, //sing alpha = -2pi*h_p + planeUnitNormal * (-1.0 * util::PI2 * planeNormalOrientation)); //sing beta = -2pi*sigma_p*N_p + } + //2. Case: If sigma_pq == 0 AND norm(P' - v1) < norm(G_ij) && norm(P' - v2) < norm(G_ij) with G_ij + // as the vector of v1 and v2 + // then P' is located on one line segment G_p of plane p, but not on any of its vertices + auto counter2 = thrust::counting_iterator(0); + auto secondCaseBegin = util::zip(segmentVectorsForPlane.begin(), segmentNormalOrientationForPlane.begin(), + counter2); + auto secondCaseEnd = util::zip(segmentVectorsForPlane.end(), segmentNormalOrientationForPlane.end(), + counter2 + 3); + if (std::any_of(secondCaseBegin, secondCaseEnd, [&](const auto &tuple) { + using namespace util; + const Array3 &segmentVector = thrust::get<0>(tuple); + const double segmentNormalOrientation = thrust::get<1>(tuple); + const unsigned int j = thrust::get<2>(tuple); + + //segmentNormalOrientation != 0.0 + if (std::abs(segmentNormalOrientation) > EPSILON) { + return false; + } + + const double segmentVectorNorm = euclideanNorm(segmentVector); + return projectionPointVertexNorms[(j + 1) % 3] < segmentVectorNorm && + projectionPointVertexNorms[j] < segmentVectorNorm; + })) { + using namespace util; + return std::make_pair(-1.0 * util::PI * planeDistance, //sing alpha = -pi*h_p + planeUnitNormal * + (-1.0 * util::PI * planeNormalOrientation)); //sing beta = -pi*sigma_p*N_p + } + //3. Case If sigma_pq == 0 AND norm(P' - v1) < 0 || norm(P' - v2) < 0 + // then P' is located at one of G_p's vertices + auto counter3 = thrust::counting_iterator(0); + auto thirdCaseBegin = util::zip(segmentNormalOrientationForPlane.begin(), counter3); + auto thirdCaseEnd = util::zip(segmentNormalOrientationForPlane.end(), counter3 + 3); + double r1Norm; + double r2Norm; + unsigned int j; + if (std::any_of(thirdCaseBegin, thirdCaseEnd, [&](const auto &tuple) { + using namespace util; + const double segmentNormalOrientation = thrust::get<0>(tuple); + j = thrust::get<1>(tuple); + + //segmentNormalOrientation != 0.0 + if (std::abs(segmentNormalOrientation) > EPSILON) { + return false; + } + + r1Norm = projectionPointVertexNorms[(j + 1) % 3]; + r2Norm = projectionPointVertexNorms[j]; + //r1Norm == 0.0 || r2Norm == 0.0 + return r1Norm < EPSILON || r2Norm < EPSILON; + })) { + using namespace util; + //Two segment vectors G_1 and G_2 of this plane + const Array3 &g1 = r1Norm == 0.0 ? segmentVectorsForPlane[j] : segmentVectorsForPlane[(j - 1 + 3) % 3]; + const Array3 &g2 = r1Norm == 0.0 ? segmentVectorsForPlane[(j + 1) % 3] : segmentVectorsForPlane[j]; + // theta = arcos((G_2 * -G_1) / (|G_2| * |G_1|)) + const double gdot = dot(g1 * -1.0, g2); + const double theta = gdot == 0.0 ? util::PI_2 : std::acos(gdot / (euclideanNorm(g1) * euclideanNorm(g2))); + return std::make_pair(-1.0 * theta * planeDistance, //sing alpha = -theta*h_p + planeUnitNormal * + (-1.0 * theta * planeNormalOrientation)); //sing beta = -theta*sigma_p*N_p + } + //4. Case Otherwise P' is located outside the plane S_p and then the singularity equals zero + return std::make_pair(0.0, //sing alpha = 0 + Array3{0.0, 0.0, 0.0}); //sing beta = 0 + } + + Array3 computeNormsOfProjectionPointAndVertices(const Array3 &orthogonalProjectionPointOnPlane, + const Array3Triplet &face) { + using namespace util; + return {euclideanNorm(orthogonalProjectionPointOnPlane - face[0]), + euclideanNorm(orthogonalProjectionPointOnPlane - face[1]), + euclideanNorm(orthogonalProjectionPointOnPlane - face[2])}; + } + + +} // namespace polyhedralGravity::GravityModel::detail diff --git a/src/polyhedralGravity/calculation/GravityModelDetail.h b/src/polyhedralGravity/calculation/GravityModelDetail.h new file mode 100644 index 0000000..b7da2e6 --- /dev/null +++ b/src/polyhedralGravity/calculation/GravityModelDetail.h @@ -0,0 +1,207 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "spdlog/spdlog.h" +#include "thrust/iterator/zip_iterator.h" +#include "thrust/iterator/transform_iterator.h" +#include "thrust/iterator/counting_iterator.h" +#include "thrust/transform.h" +#include "thrust/transform_reduce.h" +#include "thrust/execution_policy.h" +#include "xsimd/xsimd.hpp" + +#include "polyhedralGravity/calculation/PolyhedronTransform.h" +#include "polyhedralGravity/input/TetgenAdapter.h" +#include "polyhedralGravity/model/Polyhedron.h" +#include "polyhedralGravity/model/GravityModelData.h" +#include "polyhedralGravity/util/UtilityConstants.h" +#include "polyhedralGravity/util/UtilityContainer.h" +#include "polyhedralGravity/util/UtilityThrust.h" +#include "polyhedralGravity/output/Logging.h" + +namespace polyhedralGravity::GravityModel::detail { + + /** + * Computes the segment vectors G_ij for one plane of the polyhedron according to Tsoulis (18). + * The segment vectors G_ij represent the vector from one vertex of the face to the neighboring vertex and + * depict every line segment of the triangular face (A-B-C) + * @param vertex0 - the first vertex A + * @param vertex1 - the second vertex B + * @param vertex2 - the third vertex C + * @return the segment vectors for a plane + */ + Array3Triplet buildVectorsOfSegments(const Array3 &vertex0, const Array3 &vertex1, const Array3 &vertex2); + + /** + * Computes the plane unit normal N_p for one plane p of the polyhedron according to Tsoulis (19). + * The plane unit normal is the outward pointing normal of the face from the polyhedron. + * @param segmentVector1 - first edge + * @param segmentVector2 - second edge + * @return plane unit normal + */ + Array3 buildUnitNormalOfPlane(const Array3 &segmentVector1, const Array3 &segmentVector2); + + /** + * Computes the segment unit normals n_pq for one plane p of the polyhedron according to Tsoulis (20). + * The segment unit normal n_pq represent the normal of one line segment of a polyhedrale face. + * @param segmentVectors - the segment vectors of the face G_p(0-2) + * @param planeUnitNormal - the plane unit normal N_p + * @return segment unit normals n_pq for plane p with q = {0, 1, 2} + */ + Array3Triplet buildUnitNormalOfSegments(const Array3Triplet &segmentVectors, const Array3 &planeUnitNormal); + + /** + * Computes the plane unit normal orientation/ direction sigma_p for one plane p of the polyhedron + * according to Tsoulis (21). + * The plane unit normal orientation values represents the relative position of computation point P + * with respect to the pointing direction of N_p. E. g. if N_p points to the half-space containing P, the + * inner product of N_p and -G_i1 will be positive, leading to a negative sigma_p. + * If sigma_p is zero than P and P' lie geometrically in the same plane --> P == P'. + * @param planeUnitNormal - the plane unit normal N_p + * @param vertex0 - the first vertex of the plane + * @return plane normal orientation + */ + double computeUnitNormalOfPlaneDirection(const Array3 &planeUnitNormal, const Array3 &vertex0); + + /** + * Calculates the Hessian Plane form spanned by three given points p, q, and r. + * @param p - first point on the plane + * @param q - second point on the plane + * @param r - third point on the plane + * @return HessianPlane + * @related Cross-Product method https://tutorial.math.lamar.edu/classes/calciii/eqnsofplanes.aspx + */ + HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r); + + /** + * Calculates the (plane) distances h_p of computation point P to the plane S_p given in Hessian Form + * according to the following equation: + * h_p = D / sqrt(A^2+B^2+C^2) + * @param hessianPlane - Hessian Plane Form of S_p + * @return plane distance h_p + */ + double distanceBetweenOriginAndPlane(const HessianPlane &hessianPlane); + + /** + * Computes P' for a given plane p according to equation (22) of Tsoulis paper. + * P' is the orthogonal projection of the computation point P onto the plane S_p. + * @param planeUnitNormal - the plane unit normal N_p + * @param planeDistance - the distance from P to the plane h_p + * @param hessianPlane - the Hessian Plane Form + * @return P' for this plane + */ + Array3 projectPointOrthogonallyOntoPlane(const Array3 &planeUnitNormal, double planeDistance, + const HessianPlane &hessianPlane); + + /** + * Computes the segment normal orientations/ directions sigma_pq for a given plane p. + * If sigma_pq is negative, this denotes that n_pq points to the half-plane containing P'. Nn case + * sigma_pq is positive, P' resides in the other half-plane and if sigma_pq is zero, then P' lies directly + * on the segment pq. + * @param vertices - the vertices of this plane + * @param projectionPointOnPlane - the projection point P' for this plane + * @param segmentUnitNormalsForPlane - the segment unit normals sigma_pq for this plane + * @return the segment normal orientations for the plane p + */ + Array3 computeUnitNormalOfSegmentsDirections(const Array3Triplet &vertices, const Array3 &projectionPointOnPlane, + const Array3Triplet &segmentUnitNormalsForPlane); + + /** + * Computes the orthogonal projection Points P'' foreach segment q of a given plane p. + * @param projectionPointOnPlane - the projection Point P' + * @param segmentNormalOrientations - the segment normal orientations sigma_pq for this plane p + * @param face - the vertices of the plane p + * @return the orthogonal projection points of P on the segment P'' foreach segment q of p + */ + Array3Triplet + projectPointOrthogonallyOntoSegments(const Array3 &projectionPointOnPlane, const Array3 &segmentNormalOrientations, + const Array3Triplet &face); + + /** + * Calculates the point P'' for a given Segment consisting of vertices v1 and v2 and the orthogonal projection + * point P' for the plane consisting of those vertices. Solves the three equations given in (24), (25) and (26). + * @param vertex1 - first endpoint of segment + * @param vertex2 - second endpoint of segment + * @param orthogonalProjectionPointOnPlane - the orthogonal projection P' of P on this plane + * @return P'' for this segment + * @note If sigma_pq is zero then P'' == P', this is not checked by this method, but has to be assured first + */ + Array3 projectPointOrthogonallyOntoSegment(const Array3 &vertex1, const Array3 &vertex2, + const Array3 &orthogonalProjectionPointOnPlane); + + /** + * Computes the (segment) distances h_pq between P' for a given plane p and P'' for a given segment q of plane p. + * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for p + * @param orthogonalProjectionPointOnSegments - the orthogonal projection points P'' for each segment q of p + * @return distances h_pq for plane p + */ + Array3 distancesBetweenProjectionPoints(const Array3 &orthogonalProjectionPointOnPlane, + const Array3Triplet &orthogonalProjectionPointOnSegments); + + /** + * Computes the 3D distances l1_pq and l2_pq between the computation point P and the line + * segment endpoints of each polyhedral segment for one plane. + * Computes the 1D distances s1_pq and s2_pq between orthogonal projection of P on the line + * segment P''_pq and the line segment endpoints for each polyhedral segment for one plane + * @param segmentVectorsForPlane - the segment vectors G_pq for plane p + * @param orthogonalProjectionPointsOnSegmentForPlane - the orthogonal projection Points P'' for plane p + * @param face - the vertices of plane p + * @return distances l1_pq and l2_pq and s1_pq and s2_pq foreach segment q of plane p + */ + std::array distancesToSegmentEndpoints(const Array3Triplet &segmentVectorsForPlane, + const Array3Triplet &orthogonalProjectionPointsOnSegmentForPlane, + const Array3Triplet &face); + + /** + * Calculates the Transcendental Expressions LN_pq and AN_pq for every line segment of the polyhedron for + * a given plane p. + * LN_pq is calculated according to (14) using the natural logarithm and AN_pq is calculated according + * to (15) using the arctan. + * @param distancesForPlane - the distances l1, l2, s1, s2 foreach segment q of plane p + * @param planeDistance - the plane distance h_p for plane p + * @param segmentDistancesForPlane - the segment distance h_pq for segment q of plane p + * @param segmentNormalOrientationsForPlane - the segment normal orientations n_pq for a plane p + * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' for plane p + * @param face - the vertices of plane p + * @return LN_pq and AN_pq foreach segment q of plane p + */ + std::array + computeTranscendentalExpressions(const std::array &distancesForPlane, double planeDistance, + const Array3 &segmentDistancesForPlane, + const Array3 &segmentNormalOrientationsForPlane, + const Array3 &projectionPointVertexNorms); + + /** + * Calculates the singularities (correction) terms according to the Flow text for a given plane p. + * @param segmentVectorsForPlane - the segment vectors for a given plane + * @param segmentNormalOrientationForPlane - the segment orientation sigma_pq + * @param projectionPointVertexNorms - the projection point P' + * @param planeUnitNormal - the plane unit normal N_p + * @param planeDistance - the plane distance h_p + * @param planeNormalOrientation - the plane normal orientation sigma_p + * @param face - the vertices of plane p + * @return the singularities for a plane p + */ + std::pair + computeSingularityTerms(const Array3Triplet &segmentVectorsForPlane, const Array3 &segmentNormalOrientationForPlane, + const Array3 &projectionPointVertexNorms, const Array3 &planeUnitNormal, + double planeDistance, double planeNormalOrientation); + + + /** + * Computes the L2 norms of the orthogonal projection point P' on a plane p with each vertex of that plane p. + * The values are later used to determine if P' is situated at a vertex. + * @param orthogonalProjectionPointOnPlane - the orthogonal projection point P' + * @param face - the vertices of plane p + * @return the norms of p and each vertex + */ + Array3 + computeNormsOfProjectionPointAndVertices(const Array3 &orthogonalProjectionPointOnPlane, const Array3Triplet &face); + + +} // namespace polyhedralGravity::GravityModel::detail diff --git a/src/polyhedralGravity/calculation/MeshChecking.h b/src/polyhedralGravity/calculation/MeshChecking.h index 59702d0..c1b71b0 100644 --- a/src/polyhedralGravity/calculation/MeshChecking.h +++ b/src/polyhedralGravity/calculation/MeshChecking.h @@ -8,7 +8,7 @@ #include "polyhedralGravity/model/GravityModelData.h" #include "polyhedralGravity/util/UtilityContainer.h" #include "polyhedralGravity/util/UtilityConstants.h" -#include "polyhedralGravity/calculation/GravityModel.h" +#include "polyhedralGravity/calculation/PolyhedronTransform.h" namespace polyhedralGravity::MeshChecking { diff --git a/src/polyhedralGravity/model/Polyhedron.h b/src/polyhedralGravity/model/Polyhedron.h index 97df911..2b35094 100644 --- a/src/polyhedralGravity/model/Polyhedron.h +++ b/src/polyhedralGravity/model/Polyhedron.h @@ -9,7 +9,7 @@ #include #include #include -#include "GravityModelData.h" +#include "polyhedralGravity/model/GravityModelData.h" namespace polyhedralGravity { From 3a1ab7c2130d89c7e51a68cc9cef8ed91d884254 Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Mon, 31 Jul 2023 20:17:38 +0200 Subject: [PATCH 6/7] using namespace in lambda --- src/polyhedralGravityPython/PolyhedralGravityPython.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp index 9f5c05c..4cd4e5a 100644 --- a/src/polyhedralGravityPython/PolyhedralGravityPython.cpp +++ b/src/polyhedralGravityPython/PolyhedralGravityPython.cpp @@ -24,12 +24,15 @@ PYBIND11_MODULE(polyhedral_gravity, m) { m.def("evaluate", [](const PolyhedralSource &polyhedralSource, double density, const std::variant> &computationPoints, bool parallel) -> std::variant> { + using namespace polyhedralGravity; if (std::holds_alternative(computationPoints)) { return std::visit([&density, &computationPoints, parallel](const auto &source) { + using namespace polyhedralGravity; return GravityModel::evaluate(source, density, std::get(computationPoints), parallel); }, polyhedralSource); } else { return std::visit([&density, &computationPoints, parallel](const auto &source) { + using namespace polyhedralGravity; return GravityModel::evaluate(source, density, std::get>(computationPoints), parallel); }, polyhedralSource); From f5d9a285d55a478489e071c5b5ef681d809640e4 Mon Sep 17 00:00:00 2001 From: Jonas Schuhmacher Date: Mon, 31 Jul 2023 21:17:19 +0200 Subject: [PATCH 7/7] update README.md --- README.md | 42 +++++++++++++++++++++++++++++++++--------- docs/build.rst | 43 +++++++++++++++++++++++++++++++------------ setup.py | 1 + 3 files changed, 65 insertions(+), 21 deletions(-) diff --git a/README.md b/README.md index 15d8633..e2e63b6 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,16 @@ ![Build and Test](https://github.com/schuhmaj/polyhedral-gravity-model-cpp/actions/workflows/ctest.yml/badge.svg) [![Documentation Status](https://readthedocs.org/projects/polyhedral-gravity-model-cpp/badge/?version=latest)](https://polyhedral-gravity-model-cpp.readthedocs.io/en/latest/?badge=latest) +![GitHub](https://img.shields.io/github/license/esa/polyhedral-gravity-model) + +![PyPI](https://img.shields.io/pypi/v/polyhedral-gravity) +![Static Badge](https://img.shields.io/badge/platform-linux--64_%7C_windows--64_%7C_osx--64_%7C_linux--arm64_%7C_osx--arm64-lightgrey) +![PyPI - Downloads](https://img.shields.io/pypi/dm/polyhedral-gravity) + +![Conda](https://img.shields.io/conda/v/conda-forge/polyhedral-gravity-model) +![Conda](https://img.shields.io/conda/pn/conda-forge/polyhedral-gravity-model) +![Conda](https://img.shields.io/conda/dn/conda-forge/polyhedral-gravity-model) + This code is a validated implementation in C++17 of the Polyhedral Gravity Model by Tsoulis et al.. It was created in a collaborative project between @@ -55,21 +65,35 @@ The python interface can be easily installed with conda install -c conda-forge polyhedral-gravity-model -This is currently only supported for `x86-64` systems since -one of the dependencies is not available on conda for `aarch64`. -However, building from source with `pip` can also be done -on `aarch64` as shown below. - ### pip -Use pip to install the python interface in your local python runtime. -The module will be build using CMake and the using the above -requirements. Just execute in repository root: +As a second option, you can also install the python interface with pip. + + pip install polyhedral-gravity + +Binaries for the most common platforms are available on PyPI including +Windows, Linux and macOS. For macOS and Linux, binaries for +`x86_64` and `aarch64` are provided. +In case `pip` uses the source distribution, please make sure that +you have a C++17 capable compiler, CMake and ninja-build installed. + +### From source + +The module will be build using a C++17 capable compiler, +CMake and ninja-build. Just execute the following command in +the repository root folder: pip install . To modify the build options (like parallelization) have a look -at the `setupy.py` and the [next paragraph](#build-c). +at the `setupy.py` and the [next paragraph](#build-c). The options +are modified by setting the environment variables before executing +the `pip install .` command, e.g.: + + export POLYHEDRAL_GRAVITY_PARALLELIZATION="TBB" + pip install . + + (Optional: For a faster build you can install all dependencies available for your system in your local python environment. That way, they won't be fetched from GitHub.) diff --git a/docs/build.rst b/docs/build.rst index ead5b0c..d150e4b 100644 --- a/docs/build.rst +++ b/docs/build.rst @@ -23,7 +23,7 @@ The requirements (see below) are set-up automatically during the build process. Use the instructions below to build the project, from the repository's root directory: -.. code-block:: +.. code-block:: bash mkdir build cd build @@ -43,8 +43,34 @@ BUILD_POLYHEDRAL_GRAVITY_TESTS (:code:`ON`) Build the Tests BUILD_POLYHEDRAL_PYTHON_INTERFACE (:code:`ON`) Build the Python interface ================================================ =================================================================================================================================== -Build & Installation with pip ------------------------------ +Installation with conda +----------------------- + +The python interface can be easily installed with `conda `__: + +.. code-block:: bash + + conda install -c conda-forge polyhedral-gravity-model + +Installation with pip +--------------------- + +As a second option, you can also install the python interface with pip from +`PyPi `__: + +.. code-block:: bash + + pip install polyhedral-gravity + +Binaries for the most common platforms are available on PyPI including +Windows, Linux and macOS. For macOS and Linux, binaries for +:code:`x86_64` and :code:`aarch64` are provided. +In case :code:`pip` uses the source distribution, please make sure that +you have a C++17 capable compiler, CMake and ninja-build installed. + + +Build & Installation from source +-------------------------------- Use pip to install the python interface in your local python runtime. The module will be build using CMake. Just execute in repository root: @@ -54,19 +80,12 @@ The module will be build using CMake. Just execute in repository root: pip install . To modify the build options (like parallelization) have a look -at the :code:`setupy.py`. As simple example, to modify the parallelization, -just set the environment variable like below: +at the :code:`setupy.py`. The options are modified by setting the +environment variables before executing the :code:`pip install .` command, e.g.: .. code-block:: export POLYHEDRAL_GRAVITY_PARALLELIZATION="TBB" -Installation with conda ------------------------ -The python interface can be easily installed with `conda `__: - -.. code-block:: - - conda install -c conda-forge polyhedral-gravity-model diff --git a/setup.py b/setup.py index 36ddaf5..d5858a1 100644 --- a/setup.py +++ b/setup.py @@ -158,6 +158,7 @@ def build_extension(self, ext): """, ext_modules=[CMakeExtension("polyhedral_gravity")], cmdclass={"build_ext": CMakeBuild}, + license="GPLv3", license_file="LICENSE", zip_safe=False, python_requires=">=3.6",