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 all 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
28 changes: 9 additions & 19 deletions kratos/geometries/line_2d_2.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1086,28 +1086,18 @@ 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

const double length = Length();

const double length_1 = std::sqrt( std::pow(rPoint[0] - r_first_point[0], 2)
+ std::pow(rPoint[1] - r_first_point[1], 2));
// 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_2 = std::sqrt( std::pow(rPoint[0] - r_second_point[0], 2)
+ std::pow(rPoint[1] - r_second_point[1], 2));
// 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;

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 @@ -1429,4 +1419,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