Skip to content
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

[CORE] Fixed issue with projection onto 2D2 line when input does not lie exactly on the line. #12637

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
+51 −19
Open

[CORE] Fixed issue with projection onto 2D2 line when input does not lie exactly on the line. #12637

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
d049fd6
Fixed issue with projection onto 2D2 line and added two unit tests to…
rfaasse Aug 20, 2024
230a6f6
Slight modification to the wording of a comment
rfaasse Aug 20, 2024
7de6ae2
Fixed some alignment issues.
rfaasse Aug 20, 2024
977132d
Removed redundant comment
rfaasse Aug 20, 2024
e9d9224
Recover the tolerance, as it is necessary for certain functionality.
rfaasse Aug 22, 2024
64ee075
Explicitly use array_1d<double, 3> instead of auto, to fix unexpected…
rfaasse Aug 26, 2024
ddd6f86
Reduce recomputation
loumalouomega Aug 27, 2024
f317523
Removed a comment since code is self-explanatory
rfaasse Nov 21, 2024
c276994
Merge remote-tracking branch 'origin/master' into geo/fix-2-noded-lin…
rfaasse Nov 21, 2024
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
30 changes: 9 additions & 21 deletions kratos/geometries/line_2d_2.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1066,28 +1066,16 @@ class Line2D2 : public Geometry<TPointType>
const TPointType& r_first_point = BaseType::GetPoint(0);
const TPointType& r_second_point = BaseType::GetPoint(1);

// Project point
const double tolerance = 1e-14; // Tolerance
// Project the point on the line in global space
const array_1d<double, 3> vector_from_first_point_to_input = rPoint - r_first_point;
const array_1d<double, 3> unity_line_direction = (r_second_point - r_first_point) / Length();
const auto projection_on_line = inner_prod(vector_from_first_point_to_input, unity_line_direction);

const double length = Length();
// Conversion to local space
constexpr double tolerance = 1e-14;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i understand that this tolerance was hardcoded before, but it is definitely now too nice ( for example it may fail for very small geometries).

Not telling that we should consider this as a blocker, but this is not a very robust solution.

Our point here is: why do we need the tolerance at all? anyhow you are assuming lenght to be different from zero, otherwise it would fail already in 1093

Copy link
Contributor Author

@rfaasse rfaasse Dec 9, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I completely agree that we shouldn't have the tolerance here at all. I tried before to remove it, but the IsInside function seems to expect the results from PointLocalCoordinates to be as they are (as @loumalouomega mentioned), so tests started failing.

To keep the scope of this PR limited, I only fixed the projection itself, but kept the conversion of the projection to a $\xi$ value as is.

As proposed by @loumalouomega and @matekelemen, the IsInside function should be changed (i.e. the tolerance should be applied there instead of in PointLocalCoordinates), so I would propose to investigate and change this in a separate small-scoped PR.

rResult[0] = 2.0 * projection_on_line/(Length() + tolerance) - 1.0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to compute Length() twice.

Besides, the job of that tolerance is probably to prevent division by zero errors but there's a naked division on line 1071. I'd either get rid of tolerance completely (preferred - degenerate lines should be caught at construction) or apply it there as well to remain consistent.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right about avoiding computing Length twice, but the tolerance is not because is degenerated, it is because the tolerance on the checks to avoid false negatives.

Copy link
Contributor

@matekelemen matekelemen Aug 27, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

False negatives?

As I understand this function doesn't do point membership tests but transforms a point from global to local space.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, but that is used in the IsInside function, maybe the proper way should to check tolerance there (where Toelrance is an input)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yup, I like that idea better.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By the way IsInside should probably be changed after this PR because an infinitely large cylinder would project to this line.

Copy link
Contributor Author

@rfaasse rfaasse Aug 27, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, at first I removed the tolerance completely, but it failed some test cases in different applications, so it would require quite some effort to fix (simply adding the tolerance when calling IsInside didn't work unfortunately).

I can have a look again after I fixed the issue where using auto for an array1d gave issues, maybe the tolerance can now be removed safely.

However, if that's not the case and we need to change quite a bit of code, I would suggest to do that in a separate PR and keep this one focused on just fixing the projection and keep the tolerance that was already there.


const double length_1 = std::sqrt( std::pow(rPoint[0] - r_first_point[0], 2)
+ std::pow(rPoint[1] - r_first_point[1], 2));

const double length_2 = std::sqrt( std::pow(rPoint[0] - r_second_point[0], 2)
+ std::pow(rPoint[1] - r_second_point[1], 2));

if (length_1 <= (length + tolerance) && length_2 <= (length + tolerance)) {
rResult[0] = 2.0 * length_1/(length + tolerance) - 1.0;
} else {
if (length_1 > length_2) {
rResult[0] = 2.0 * length_1/(length + tolerance) - 1.0;
} else {
rResult[0] = -2.0 * length_1/(length + tolerance) - 1.0;
}
}

return rResult ;
return rResult;
}

///@}
Expand Down Expand Up @@ -1409,4 +1397,4 @@ const GeometryData Line2D2<TPointType>::msGeometryData(
template<class TPointType>
const GeometryDimension Line2D2<TPointType>::msGeometryDimension(2, 1);

} // namespace Kratos.
} // namespace Kratos.
42 changes: 42 additions & 0 deletions kratos/tests/cpp_tests/geometries/test_line_2d_2.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,6 +270,48 @@ namespace Testing {
KRATOS_EXPECT_NEAR(test_point_local_coords(2), 0.0, TOLERANCE);
}

/** Checks the point local coordinates for a given point with respect to the line.
* A point 'above' the line is selected (meaning it moved a bit in the direction normal to the line).
*/
KRATOS_TEST_CASE_IN_SUITE(Line2D2PointLocalCoordinates_NotAligned, KratosCoreGeometriesFastSuite)
{
// Create the test line
auto geom = GeneratePointsDiagonalLine2D2();

// Set the point to be checked, starting at the center and move it in the normal direction
Point test_point(0.5 - 0.1,0.5 + 0.1,0.0);

// Compute the centre local coordinates
array_1d<double, 3> test_point_local_coords;
geom->PointLocalCoordinates(test_point_local_coords, test_point);

// Since we started at the center and moved in the normal direction, we expect the local
// coordinate to be (0, 0, 0)
KRATOS_EXPECT_NEAR(test_point_local_coords(0), 0.0, TOLERANCE);
KRATOS_EXPECT_NEAR(test_point_local_coords(1), 0.0, TOLERANCE);
KRATOS_EXPECT_NEAR(test_point_local_coords(2), 0.0, TOLERANCE);
}

/** Checks the point local coordinates for a given point with respect to the line.
* A point outside and not aligned with the line is selected. In this case the distance from first node is larger than the length
*/
KRATOS_TEST_CASE_IN_SUITE(Line2D2PointLocalCoordinatesOutsidePointNotAligned, KratosCoreGeometriesFastSuite)
{
// Create the test line
auto geom = GeneratePointsDiagonalLine2D2();

// Set the point to be checked (moved in the normal direction)
Point test_point(-1.5 - 0.1,-1.5 + 0.1,0.0);

// Compute the centre local coordinates
array_1d<double, 3> test_point_local_coords;
geom->PointLocalCoordinates(test_point_local_coords, test_point);

KRATOS_EXPECT_NEAR(test_point_local_coords(0), -4.0, TOLERANCE);
KRATOS_EXPECT_NEAR(test_point_local_coords(1), 0.0, TOLERANCE);
KRATOS_EXPECT_NEAR(test_point_local_coords(2), 0.0, TOLERANCE);
}

/** Tests the Jacobian determinants using 'GI_GAUSS_1' integration method.
* Tests the Jacobian determinants using 'GI_GAUSS_1' integration method.
*/
Expand Down
Loading