diff --git a/kratos/geometries/triangle_3d_3.h b/kratos/geometries/triangle_3d_3.h index 017f7d280ff0..e2913c3c7095 100644 --- a/kratos/geometries/triangle_3d_3.h +++ b/kratos/geometries/triangle_3d_3.h @@ -939,55 +939,7 @@ template class Triangle3D3 const CoordinatesArrayType& rPoint ) const override { - // Initialize - noalias(rResult) = ZeroVector(3); - - // Tangent vectors - array_1d tangent_xi = this->GetPoint(1) - this->GetPoint(0); - tangent_xi /= norm_2(tangent_xi); - array_1d tangent_eta = this->GetPoint(2) - this->GetPoint(0); - tangent_eta /= norm_2(tangent_eta); - - // The center of the geometry - const auto center = this->Center(); - - // Computation of the rotation matrix - BoundedMatrix rotation_matrix = ZeroMatrix(3, 3); - for (IndexType i = 0; i < 3; ++i) { - rotation_matrix(0, i) = tangent_xi[i]; - rotation_matrix(1, i) = tangent_eta[i]; - } - - // Destination point rotated - CoordinatesArrayType aux_point_to_rotate, destination_point_rotated; - noalias(aux_point_to_rotate) = rPoint - center.Coordinates(); - noalias(destination_point_rotated) = prod(rotation_matrix, aux_point_to_rotate) + center.Coordinates(); - - // Points of the geometry - array_1d points_rotated; - for (IndexType i = 0; i < 3; ++i) { - noalias(aux_point_to_rotate) = this->GetPoint(i).Coordinates() - center.Coordinates(); - noalias(points_rotated[i]) = prod(rotation_matrix, aux_point_to_rotate) + center.Coordinates(); - } - - // Compute the Jacobian matrix and its determinant - BoundedMatrix J; - J(0,0) = points_rotated[1][0] - points_rotated[0][0]; - J(0,1) = points_rotated[2][0] - points_rotated[0][0]; - J(1,0) = points_rotated[1][1] - points_rotated[0][1]; - J(1,1) = points_rotated[2][1] - points_rotated[0][1]; - const double det_J = J(0,0)*J(1,1) - J(0,1)*J(1,0); - - // Compute eta and xi - const double eta = (J(1,0)*(points_rotated[0][0] - destination_point_rotated[0]) + - J(0,0)*(destination_point_rotated[1] - points_rotated[0][1])) / det_J; - const double xi = (J(1,1)*(destination_point_rotated[0] - points_rotated[0][0]) + - J(0,1)*(points_rotated[0][1] - destination_point_rotated[1])) / det_J; - - rResult(0) = xi; - rResult(1) = eta; - - return rResult; + return GeometryUtils::PointLocalCoordinatesStraightEdgesTriangle(*this, rResult, rPoint); } ///@} diff --git a/kratos/geometries/triangle_3d_6.h b/kratos/geometries/triangle_3d_6.h index a5fe2f04427b..69656575e1cc 100644 --- a/kratos/geometries/triangle_3d_6.h +++ b/kratos/geometries/triangle_3d_6.h @@ -915,6 +915,24 @@ template class Triangle3D6 return rResult; } + /** + * @brief Returns the local coordinates of a given arbitrary point + * @param rResult The vector containing the local coordinates of the point + * @param rPoint The point in global coordinates + * @return The vector containing the local coordinates of the point + */ + CoordinatesArrayType& PointLocalCoordinates( + CoordinatesArrayType& rResult, + const CoordinatesArrayType& rPoint + ) const override + { + if (this->EdgesAreStraight()) { + return GeometryUtils::PointLocalCoordinatesStraightEdgesTriangle(*this, rResult, rPoint); + } else { + return BaseType::PointLocalCoordinates( rResult, rPoint ); + } + } + ///@} ///@name Friends ///@{ @@ -1111,6 +1129,61 @@ template class Triangle3D6 return shape_functions_local_gradients; } + /** This method gives all shape functions values + * evaluated at the rCoordinates provided + * @return Vector of values of shape functions \f$ F_{i} \f$ + * where i is the shape function index (for NURBS it is the index + * of the local enumeration in the element). + * + * @see ShapeFunctionValue + * @see ShapeFunctionsLocalGradients + * @see ShapeFunctionLocalGradient + */ + Vector& ShapeFunctionsValues ( + Vector &rResult, + const CoordinatesArrayType& rCoordinates) const override + { + if(rResult.size() != 6) { + rResult.resize(6, false); + } + + const double xi = rCoordinates[0]; + const double eta = rCoordinates[1]; + const double zeta = 1.0 - xi - eta; + + rResult[0] = zeta * (2.0 * zeta - 1.0); + rResult[1] = xi * (2.0 * xi - 1.0); + rResult[2] = eta * (2.0 * eta - 1.0); + rResult[3] = 4.0 * xi * zeta; + rResult[4] = 4.0 * xi * eta; + rResult[5] = 4.0 * eta * zeta; + + return rResult; + } + + + /** + * @brief Checks if edges are straight. We iterate though all edges + * and check that the sum of 0-2 and 2-1 segments is no bigger than 0-1. + * @return bool edges are straight or not + */ + bool EdgesAreStraight() const + { + constexpr double tol = 1e-6; + constexpr std::array, 3> edges{ + {{0, 1, 3}, {1, 2, 4}, {2, 0, 5}}}; + const auto& r_points = this->Points(); + for (const auto& r_edge : edges) { + const double a = MathUtils::Norm3(r_points[r_edge[0]] - r_points[r_edge[1]]); + const double b = MathUtils::Norm3(r_points[r_edge[1]] - r_points[r_edge[2]]); + const double c = MathUtils::Norm3(r_points[r_edge[2]] - r_points[r_edge[0]]); + if (b + c > a*(1.0+tol) ) { + return false; + } + } + return true; + } + ///@} ///@name Private Access ///@{ diff --git a/kratos/tests/cpp_tests/geometries/test_triangle_3d_6.cpp b/kratos/tests/cpp_tests/geometries/test_triangle_3d_6.cpp index 71c7e63fcc6d..1661307ab651 100644 --- a/kratos/tests/cpp_tests/geometries/test_triangle_3d_6.cpp +++ b/kratos/tests/cpp_tests/geometries/test_triangle_3d_6.cpp @@ -115,6 +115,21 @@ namespace { TestAllShapeFunctionsLocalGradients(*geom); } + KRATOS_TEST_CASE_IN_SUITE(Triangle3D6ShapeFunctionsValues, KratosCoreGeometriesFastSuite) { + auto geom = GenerateReferenceTriangle3D6(); + array_1d coord(3); + coord[0] = 0.5; + coord[1] = 0.125; + coord[2] = 0.0; + KRATOS_EXPECT_NEAR(geom->ShapeFunctionValue(0, coord), -0.09375, TOLERANCE); + KRATOS_EXPECT_NEAR(geom->ShapeFunctionValue(1, coord), 0.0, TOLERANCE); + KRATOS_EXPECT_NEAR(geom->ShapeFunctionValue(2, coord), -0.09375, TOLERANCE); + KRATOS_EXPECT_NEAR(geom->ShapeFunctionValue(3, coord), 0.75, TOLERANCE); + KRATOS_EXPECT_NEAR(geom->ShapeFunctionValue(4, coord), 0.25, TOLERANCE); + KRATOS_EXPECT_NEAR(geom->ShapeFunctionValue(5, coord), 0.1875, TOLERANCE); + CrossCheckShapeFunctionsValues(*geom); + } + KRATOS_TEST_CASE_IN_SUITE(Triangle3D6LumpingFactorsRegularShape, KratosCoreGeometriesFastSuite) { auto geom = GenerateReferenceTriangle3D6(); @@ -161,4 +176,58 @@ namespace { KRATOS_EXPECT_DOUBLE_EQ(geom->CalculateDistance(point2), 0.5); } + /* + * Computes point local coordinates from a given point. + * his triangle has straight edges so we can use the same + * implementation as the Triangle3D3 + */ + KRATOS_TEST_CASE_IN_SUITE(Triangle3D6PointLocalCoordinates, KratosCoreGeometriesFastSuite) { + Triangle3D6 geom( + Kratos::make_shared(0.0, 0.0, 0.0), + Kratos::make_shared(1.0, 0.0, 0.0), + Kratos::make_shared(0.0, 1.0, 0.0), + Kratos::make_shared(0.5, 0.0, 0.0), + Kratos::make_shared(0.5, 0.5, 0.0), + Kratos::make_shared(0.0, 0.5, 0.0) + ); + + // Compute the global coordinates of the baricentre + array_1d baricentre; + baricentre[0] = 1.0/3.0; baricentre[1] = 1.0/3.0; baricentre[2] = 0.0; + + // Compute the baricentre local coordinates + array_1d baricentre_local_coords; + geom.PointLocalCoordinates(baricentre_local_coords, baricentre); + + KRATOS_EXPECT_NEAR(baricentre_local_coords(0), 1.0/3.0, TOLERANCE); + KRATOS_EXPECT_NEAR(baricentre_local_coords(1), 1.0/3.0, TOLERANCE); + KRATOS_EXPECT_NEAR(baricentre_local_coords(2), 0.0, TOLERANCE); + } + + /* + * Computes point local coordinates from a given point. + * This triangle does not ghave straight edges so this will + * throw an exception + */ + KRATOS_TEST_CASE_IN_SUITE(Triangle3D6PointLocalCoordinatesError, KratosCoreGeometriesFastSuite) { + Triangle3D6 geom( + Kratos::make_shared(1.0, 0.0, 0.0), + Kratos::make_shared(0.0, 1.0, 0.0), + Kratos::make_shared(0.0, 0.0, 1.0), + Kratos::make_shared(0.6, 0.5, 0.0), + Kratos::make_shared(0.0, 0.5, 0.5), + Kratos::make_shared(0.5, 0.0, 0.5) + ); + + // Compute the global coordinates of the baricentre + array_1d baricentre; + baricentre[0] = 1.0/3.0; baricentre[1] = 1.0/3.0; baricentre[2] = 1.0/3.0; + + // Compute the baricentre local coordinates + array_1d baricentre_local_coords; + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + geom.PointLocalCoordinates(baricentre_local_coords, baricentre), + "ERROR:: Attention, the Point Local Coordinates must be specialized for the current geometry"); + } + } // namespace Kratos::Testing. diff --git a/kratos/tests/cpp_tests/utilities/test_geometry_utils.cpp b/kratos/tests/cpp_tests/utilities/test_geometry_utils.cpp index cb1342d88565..e1bf75691d7e 100644 --- a/kratos/tests/cpp_tests/utilities/test_geometry_utils.cpp +++ b/kratos/tests/cpp_tests/utilities/test_geometry_utils.cpp @@ -586,5 +586,27 @@ namespace Testing { KRATOS_EXPECT_EQ(GeometryUtils::ProjectedIsInside(triangle, *p_node_3, aux), true); KRATOS_EXPECT_EQ(GeometryUtils::ProjectedIsInside(triangle, *p_node_4, aux), false); } + + KRATOS_TEST_CASE_IN_SUITE(GeometryUtilsPointLocalCoordinatesStraightEdgesTriangle, KratosCoreGeometriesFastSuite) + { + const double tolerance = 1.0e-9; + Triangle3D3 geom( + Kratos::make_shared(0.0, 0.0, 0.0), + Kratos::make_shared(1.0, 0.0, 0.0), + Kratos::make_shared(0.0, 1.0, 0.0) + ); + + // Compute the global coordinates of the baricentre + array_1d baricentre; + baricentre[0] = 1.0/3.0; baricentre[1] = 1.0/3.0; baricentre[2] = 0.0; + + // Compute the baricentre local coordinates + array_1d baricentre_local_coords; + GeometryUtils::PointLocalCoordinatesStraightEdgesTriangle(geom, baricentre_local_coords, baricentre); + + KRATOS_EXPECT_NEAR(baricentre_local_coords(0), 1.0/3.0, tolerance); + KRATOS_EXPECT_NEAR(baricentre_local_coords(1), 1.0/3.0, tolerance); + KRATOS_EXPECT_NEAR(baricentre_local_coords(2), 0.0, tolerance); + } } // namespace Testing. } // namespace Kratos. diff --git a/kratos/utilities/geometry_utilities.h b/kratos/utilities/geometry_utilities.h index 0f9dc17664f3..ce450fcd0922 100644 --- a/kratos/utilities/geometry_utilities.h +++ b/kratos/utilities/geometry_utilities.h @@ -166,6 +166,67 @@ class KRATOS_API(KRATOS_CORE) GeometryUtils return rResult; } + /** + * @brief Returns the local coordinates of a given arbitrary point for a given linear triangle + * @param rGeometry The geometry to be considered + * @param rResult The vector containing the local coordinates of the point + * @param rPoint The point in global coordinates + * @return The vector containing the local coordinates of the point + */ + template + static inline typename TGeometryType::CoordinatesArrayType& PointLocalCoordinatesStraightEdgesTriangle( + const TGeometryType& rGeometry, + typename TGeometryType::CoordinatesArrayType& rResult, + const typename TGeometryType::CoordinatesArrayType& rPoint + ) + { + KRATOS_DEBUG_ERROR_IF_NOT(rGeometry.GetGeometryFamily() == GeometryData::KratosGeometryFamily::Kratos_Triangle) << + "Geometry should be a triangle in order to use PointLocalCoordinatesStraightEdgesTriangle" << std::endl; + + noalias(rResult) = ZeroVector(3); + + array_1d tangent_xi = rGeometry.GetPoint(1) - rGeometry.GetPoint(0); + tangent_xi /= norm_2(tangent_xi); + array_1d tangent_eta = rGeometry.GetPoint(2) - rGeometry.GetPoint(0); + tangent_eta /= norm_2(tangent_eta); + + const auto center = rGeometry.Center(); + + BoundedMatrix rotation_matrix = ZeroMatrix(3, 3); + for (IndexType i = 0; i < 3; ++i) { + rotation_matrix(0, i) = tangent_xi[i]; + rotation_matrix(1, i) = tangent_eta[i]; + } + + typename TGeometryType::CoordinatesArrayType aux_point_to_rotate, destination_point_rotated; + noalias(aux_point_to_rotate) = rPoint - center.Coordinates(); + noalias(destination_point_rotated) = prod(rotation_matrix, aux_point_to_rotate) + center.Coordinates(); + + array_1d points_rotated; + for (IndexType i = 0; i < 3; ++i) { + noalias(aux_point_to_rotate) = rGeometry.GetPoint(i).Coordinates() - center.Coordinates(); + noalias(points_rotated[i]) = prod(rotation_matrix, aux_point_to_rotate) + center.Coordinates(); + } + + // Compute the Jacobian matrix and its determinant + BoundedMatrix J; + J(0,0) = points_rotated[1][0] - points_rotated[0][0]; + J(0,1) = points_rotated[2][0] - points_rotated[0][0]; + J(1,0) = points_rotated[1][1] - points_rotated[0][1]; + J(1,1) = points_rotated[2][1] - points_rotated[0][1]; + const double det_J = J(0,0)*J(1,1) - J(0,1)*J(1,0); + + const double eta = (J(1,0)*(points_rotated[0][0] - destination_point_rotated[0]) + + J(0,0)*(destination_point_rotated[1] - points_rotated[0][1])) / det_J; + const double xi = (J(1,1)*(destination_point_rotated[0] - points_rotated[0][0]) + + J(0,1)*(points_rotated[0][1] - destination_point_rotated[1])) / det_J; + + rResult(0) = xi; + rResult(1) = eta; + + return rResult; + } + /** * @brief This function is designed to compute the shape function derivatives, shape functions and volume in 3D * @param rGeometry it is the array of nodes. It is expected to be a tetrahedra