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

[FastPR][Core] Missing Tetrahedra3D10 second derivatives #12082

Merged
merged 3 commits into from
Feb 20, 2024
Merged

[FastPR][Core] Missing Tetrahedra3D10 second derivatives #12082

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
11a5789
Tetrahedra3D10 shape functions second derivatives
rubenzorrilla Feb 19, 2024
53e2f9f
Missing typename
rubenzorrilla Feb 19, 2024
79dfed8
Test
rubenzorrilla Feb 19, 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
73 changes: 73 additions & 0 deletions kratos/geometries/tetrahedra_3d_10.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,10 @@ class Tetrahedra3D10
/// type as its result.
using ShapeFunctionsGradientsType = typename BaseType::ShapeFunctionsGradientsType;

/// A third order tensor to hold shape functions' local second derivatives.
/// ShapeFunctionsSecondDerivatives function return this type as its result.
using ShapeFunctionsSecondDerivativesType = typename BaseType::ShapeFunctionsSecondDerivativesType;

/// Type of the normal vector used for normal to edges in geometry.
using NormalType = typename BaseType::NormalType;

Expand Down Expand Up @@ -742,6 +746,75 @@ class Tetrahedra3D10
return rResult;
}

ShapeFunctionsSecondDerivativesType& ShapeFunctionsSecondDerivatives(
ShapeFunctionsSecondDerivativesType& rResult,
const CoordinatesArrayType& rPoint) const override
{
// Check and resize results container
if (rResult.size() != this->PointsNumber()) {
rResult.resize(this->PointsNumber(), false);
}

for (IndexType i = 0; i < this->PointsNumber(); ++i) {
auto& r_DDN_i = rResult[i];
if (r_DDN_i.size1() != 3 || r_DDN_i.size2() != 3) {
r_DDN_i.resize(3,3, false);
}
}

// Node 0
rResult[0](0,0) = 4.0; rResult[0](0,1) = 4.0; rResult[0](0,2) = 4.0;
rResult[0](1,0) = 4.0; rResult[0](1,1) = 4.0; rResult[0](1,2) = 4.0;
rResult[0](2,0) = 4.0; rResult[0](2,1) = 4.0; rResult[0](2,2) = 4.0;

// Node 1
rResult[1](0,0) = 4.0; rResult[1](0,1) = 0.0; rResult[1](0,2) = 0.0;
rResult[1](1,0) = 0.0; rResult[1](1,1) = 0.0; rResult[1](1,2) = 0.0;
rResult[1](2,0) = 0.0; rResult[1](2,1) = 0.0; rResult[1](2,2) = 0.0;

// Node 2
rResult[2](0,0) = 0.0; rResult[2](0,1) = 0.0; rResult[2](0,2) = 0.0;
rResult[2](1,0) = 0.0; rResult[2](1,1) = 4.0; rResult[2](1,2) = 0.0;
rResult[2](2,0) = 0.0; rResult[2](2,1) = 0.0; rResult[2](2,2) = 0.0;

// Node 3
rResult[3](0,0) = 0.0; rResult[3](0,1) = 0.0; rResult[3](0,2) = 0.0;
rResult[3](1,0) = 0.0; rResult[3](1,1) = 0.0; rResult[3](1,2) = 0.0;
rResult[3](2,0) = 0.0; rResult[3](2,1) = 0.0; rResult[3](2,2) = 4.0;

// Node 4
rResult[4](0,0) = -8.0; rResult[4](0,1) = -4.0; rResult[4](0,2) = -4.0;
rResult[4](1,0) = -4.0; rResult[4](1,1) = 0.0; rResult[4](1,2) = 0.0;
rResult[4](2,0) = -4.0; rResult[4](2,1) = 0.0; rResult[4](2,2) = 0.0;

// Node 5
rResult[5](0,0) = 0.0; rResult[5](0,1) = 4.0; rResult[5](0,2) = 0.0;
rResult[5](1,0) = 4.0; rResult[5](1,1) = 0.0; rResult[5](1,2) = 0.0;
rResult[5](2,0) = 0.0; rResult[5](2,1) = 0.0; rResult[5](2,2) = 0.0;

// Node 6
rResult[6](0,0) = 0.0; rResult[6](0,1) = -4.0; rResult[6](0,2) = 0.0;
rResult[6](1,0) = -4.0; rResult[6](1,1) = -8.0; rResult[6](1,2) = -4.0;
rResult[6](2,0) = 0.0; rResult[6](2,1) = -4.0; rResult[6](2,2) = 0.0;

// Node 7
rResult[7](0,0) = 0.0; rResult[7](0,1) = 0.0; rResult[7](0,2) = -4.0;
rResult[7](1,0) = 0.0; rResult[7](1,1) = 0.0; rResult[7](1,2) = -4.0;
rResult[7](2,0) = -4.0; rResult[7](2,1) = -4.0; rResult[7](2,2) = -8.0;

// Node 8
rResult[8](0,0) = 0.0; rResult[8](0,1) = 0.0; rResult[8](0,2) = 4.0;
rResult[8](1,0) = 0.0; rResult[8](1,1) = 0.0; rResult[8](1,2) = 0.0;
rResult[8](2,0) = 4.0; rResult[8](2,1) = 0.0; rResult[8](2,2) = 0.0;

// Node 9
rResult[9](0,0) = 0.0; rResult[9](0,1) = 0.0; rResult[9](0,2) = 0.0;
rResult[9](1,0) = 0.0; rResult[9](1,1) = 0.0; rResult[9](1,2) = 4.0;
rResult[9](2,0) = 0.0; rResult[9](2,1) = 4.0; rResult[9](2,2) = 0.0;

return rResult;
}

/** Tests the intersection of the geometry with
* a 3D box defined by rLowPoint and rHighPoint.
* The method is only implemented for simple tets
Expand Down
44 changes: 43 additions & 1 deletion kratos/tests/cpp_tests/geometries/test_tetrahedra_3d_10.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,48 @@ namespace{
TestAllShapeFunctionsLocalGradients(*geom);
}

KRATOS_TEST_CASE_IN_SUITE(Tetrahedra3D10ShapeFunctionsSecondDerivatives, KratosCoreGeometriesFastSuite)
{
// Set a second order isoparametric tetrahedron so we avoid the Jacobian transformation
auto p_node_0 = Kratos::make_intrusive<NodeType>(1, 0.0, 0.0, 0.0);
auto p_node_1 = Kratos::make_intrusive<NodeType>(2, 1.0, 0.0, 0.0);
auto p_node_2 = Kratos::make_intrusive<NodeType>(3, 0.0, 1.0, 0.0);
auto p_node_3 = Kratos::make_intrusive<NodeType>(4, 0.0, 0.0, 1.0);
auto p_node_4 = Kratos::make_intrusive<NodeType>(5, 0.5, 0.0, 0.0);
auto p_node_5 = Kratos::make_intrusive<NodeType>(6, 0.5, 0.5, 0.0);
auto p_node_6 = Kratos::make_intrusive<NodeType>(7, 0.0, 0.5, 0.0);
auto p_node_7 = Kratos::make_intrusive<NodeType>(8, 0.0, 0.0, 0.5);
auto p_node_8 = Kratos::make_intrusive<NodeType>(9, 0.5, 0.0, 0.5);
auto p_node_9 = Kratos::make_intrusive<NodeType>(10, 0.0, 0.5, 0.5);
auto p_geom = Kratos::make_unique<Tetrahedra3D10<NodeType>>(
p_node_0, p_node_1, p_node_2, p_node_3, p_node_4, p_node_5, p_node_6, p_node_7, p_node_8, p_node_9);

// Compute the shape functions second derivatives in the barycencer
GeometryType::ShapeFunctionsSecondDerivativesType DDN_DDX;
array_1d<double,3> barycenter_coords({1.0/3.0, 1.0/3.0, 1.0/3.0});
p_geom->ShapeFunctionsSecondDerivatives(DDN_DDX, barycenter_coords);

// Evaluate the second derivatives of the auxiliary field xi^2 + eta^2 + zeta^2 at the barycenter
BoundedMatrix<double, 3, 3> expected_result = ZeroMatrix(3,3);
expected_result(0,0) = 2.0;
expected_result(1,1) = 2.0;
expected_result(2,2) = 2.0;

const auto field_func = [](array_1d<double,3>& rCoords){return std::pow(rCoords[0],2) + std::pow(rCoords[1],2) + std::pow(rCoords[2],2);};
BoundedMatrix<double, 3, 3> result = ZeroMatrix(3,3);
for (IndexType i = 0; i < 10; ++i) {
const auto& r_DDN_DDX_i = DDN_DDX[i];
const double field_i = field_func((*p_geom)[i]);
for (IndexType d1 = 0; d1 < 3; ++d1) {
for (IndexType d2 = 0; d2 < 3; ++d2) {
result(d1,d2) += r_DDN_DDX_i(d1,d2) * field_i;
}
}
}

KRATOS_CHECK_MATRIX_NEAR(expected_result, result, 1.0e-12)
}

KRATOS_TEST_CASE_IN_SUITE(Tetrahedra3D10AverageEdgeLength, KratosCoreGeometriesFastSuite) {
auto geom = GenerateReferenceTetrahedra3D10();
KRATOS_EXPECT_NEAR(geom->AverageEdgeLength(), 1.20710678119, 1e-7);
Expand All @@ -171,7 +213,7 @@ namespace{
KRATOS_TEST_CASE_IN_SUITE(Tetrahedra3D10HasIntersection, KratosCoreGeometriesFastSuite) {
const Point LowPoint(0.1,0.1,-0.1);
const Point HighPoint(0.1,0.1,1.1);

const Point OutLowPoint(1.1,0.1,-0.1);
const Point OutHighPoint(1.1,0.1,1.1);

Expand Down
Loading