Skip to content

fix: correct vertex singularity detection in computeSingularityTerms #52

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
darioizzo merged 2 commits into from
Jun 17, 2026
Merged

fix: correct vertex singularity detection in computeSingularityTerms #52

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ac4a74e
fix(GravityModelDetail): exclude near-vertex points from edge singula…
sasso-effe May 27, 2026
4af8df1
fix(GravityModelDetail): close gap between Case 2 and Case 3 singular…
sasso-effe May 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions src/polyhedralGravity/model/GravityModelDetail.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -337,7 +337,9 @@ namespace polyhedralGravity::GravityModel::detail {

const double segmentVectorNorm = euclideanNorm(segmentVector);
return projectionPointVertexNorms[(j + 1) % 3] < segmentVectorNorm &&
projectionPointVertexNorms[j] < segmentVectorNorm;
projectionPointVertexNorms[j] < segmentVectorNorm &&
projectionPointVertexNorms[(j + 1) % 3] >= EPSILON_ZERO_OFFSET &&
projectionPointVertexNorms[j] >= EPSILON_ZERO_OFFSET;
})) {
using namespace util;
return std::make_pair(-1.0 * util::PI * planeDistance, //sing alpha = -pi*h_p
Expand Down Expand Up @@ -369,8 +371,8 @@ namespace polyhedralGravity::GravityModel::detail {
})) {
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];
const Array3 &g1 = r1Norm < EPSILON_ZERO_OFFSET ? segmentVectorsForPlane[j] : segmentVectorsForPlane[(j - 1 + 3) % 3];
const Array3 &g2 = r1Norm < EPSILON_ZERO_OFFSET ? 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)));
Expand Down
Loading