Skip to content

Commit

Permalink
[Core] adding PointLocalCoordinates and ShapeFunctionValues to triang…
Browse files Browse the repository at this point in the history 
…le 3d6 (#12961)

* adding PointLocalCoordinates and ShapeFunctionValues to triangle 3d6

* Fix edge definition in triangle 3D6

---------

Co-authored-by: Daniel Diez <ddiez@altair.com>
  • Loading branch information
@ddiezrod
ddiezrod and Daniel Diez authored Dec 23, 2024
1 parent bb98651 commit 3ddc3d6
Show file tree
Hide file tree
Showing 5 changed files with 226 additions and 49 deletions.
50 changes: 1 addition & 49 deletions kratos/geometries/triangle_3d_3.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -939,55 +939,7 @@ template<class TPointType> class Triangle3D3
const CoordinatesArrayType& rPoint
) const override
{
// Initialize
noalias(rResult) = ZeroVector(3);

// Tangent vectors
array_1d<double, 3> tangent_xi = this->GetPoint(1) - this->GetPoint(0);
tangent_xi /= norm_2(tangent_xi);
array_1d<double, 3> 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<double, 3, 3> 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<CoordinatesArrayType, 3> 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<double, 2, 2> 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);
}

///@}
Expand Down
73 changes: 73 additions & 0 deletions kratos/geometries/triangle_3d_6.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -915,6 +915,24 @@ template<class TPointType> 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
///@{
Expand Down Expand Up @@ -1111,6 +1129,61 @@ template<class TPointType> 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<std::array<size_t, 3>, 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<double>::Norm3(r_points[r_edge[0]] - r_points[r_edge[1]]);
const double b = MathUtils<double>::Norm3(r_points[r_edge[1]] - r_points[r_edge[2]]);
const double c = MathUtils<double>::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
///@{
Expand Down
69 changes: 69 additions & 0 deletions kratos/tests/cpp_tests/geometries/test_triangle_3d_6.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,21 @@ namespace {
TestAllShapeFunctionsLocalGradients(*geom);
}

KRATOS_TEST_CASE_IN_SUITE(Triangle3D6ShapeFunctionsValues, KratosCoreGeometriesFastSuite) {
auto geom = GenerateReferenceTriangle3D6();
array_1d<double, 3> 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();

Expand Down Expand Up @@ -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<Point> geom(
Kratos::make_shared<Point>(0.0, 0.0, 0.0),
Kratos::make_shared<Point>(1.0, 0.0, 0.0),
Kratos::make_shared<Point>(0.0, 1.0, 0.0),
Kratos::make_shared<Point>(0.5, 0.0, 0.0),
Kratos::make_shared<Point>(0.5, 0.5, 0.0),
Kratos::make_shared<Point>(0.0, 0.5, 0.0)
);

// Compute the global coordinates of the baricentre
array_1d<double, 3> baricentre;
baricentre[0] = 1.0/3.0; baricentre[1] = 1.0/3.0; baricentre[2] = 0.0;

// Compute the baricentre local coordinates
array_1d<double, 3> 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<Point> geom(
Kratos::make_shared<Point>(1.0, 0.0, 0.0),
Kratos::make_shared<Point>(0.0, 1.0, 0.0),
Kratos::make_shared<Point>(0.0, 0.0, 1.0),
Kratos::make_shared<Point>(0.6, 0.5, 0.0),
Kratos::make_shared<Point>(0.0, 0.5, 0.5),
Kratos::make_shared<Point>(0.5, 0.0, 0.5)
);

// Compute the global coordinates of the baricentre
array_1d<double, 3> 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<double, 3> 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.
22 changes: 22 additions & 0 deletions kratos/tests/cpp_tests/utilities/test_geometry_utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Point> geom(
Kratos::make_shared<Point>(0.0, 0.0, 0.0),
Kratos::make_shared<Point>(1.0, 0.0, 0.0),
Kratos::make_shared<Point>(0.0, 1.0, 0.0)
);

// Compute the global coordinates of the baricentre
array_1d<double, 3> baricentre;
baricentre[0] = 1.0/3.0; baricentre[1] = 1.0/3.0; baricentre[2] = 0.0;

// Compute the baricentre local coordinates
array_1d<double, 3> 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.
61 changes: 61 additions & 0 deletions kratos/utilities/geometry_utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<class TGeometryType>
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<double, 3> tangent_xi = rGeometry.GetPoint(1) - rGeometry.GetPoint(0);
tangent_xi /= norm_2(tangent_xi);
array_1d<double, 3> tangent_eta = rGeometry.GetPoint(2) - rGeometry.GetPoint(0);
tangent_eta /= norm_2(tangent_eta);

const auto center = rGeometry.Center();

BoundedMatrix<double, 3, 3> 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<typename TGeometryType::CoordinatesArrayType, 3> 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<double, 2, 2> 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
Expand Down

0 comments on commit 3ddc3d6

Please sign in to comment.