From 01c86af747c3ed6401ed443cea49ae278510c204 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Wed, 11 Sep 2024 11:10:42 +0200 Subject: [PATCH] Add a utility function for calculating the internal force vector (#12669) Added a utility function for calculating the internal force vector. As input, it requires vectors of B-matrices, stress vectors, and integration coefficients (one for each integration point). When the input vectors have different sizes or when they are empty, an error is raised. --- .../equation_of_motion_utilities.cpp | 71 +++++++++++++++ .../equation_of_motion_utilities.h | 4 + .../cpp_tests/test_equation_of_motion.cpp | 87 +++++++++++++++++++ 3 files changed, 162 insertions(+) diff --git a/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp index 00aabe71950e..16bd8737a112 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp +++ b/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp @@ -17,6 +17,63 @@ #include "custom_utilities/element_utilities.hpp" #include "utilities/geometry_utilities.h" +#include + +namespace +{ + +using namespace Kratos; + +void CheckThatVectorsHaveSameSizeAndAreNotEmpty(const std::vector& rBs, + const std::vector& rStressVectors, + const std::vector& rIntegrationCoefficients) +{ + KRATOS_DEBUG_ERROR_IF((rBs.size() != rStressVectors.size()) || + (rBs.size() != rIntegrationCoefficients.size())) + << "Cannot calculate the internal force vector: input vectors have different sizes\n"; + KRATOS_DEBUG_ERROR_IF(rBs.empty()) + << "Cannot calculate the internal force vector: input vectors are empty\n"; +} + +void CheckThatBMatricesHaveSameSizes(const std::vector& rBs) +{ + auto has_inconsistent_sizes = [number_of_rows = rBs.front().size1(), + number_of_columns = rBs.front().size2()](const auto& rMatrix) { + return (rMatrix.size1() != number_of_rows) || (rMatrix.size2() != number_of_columns); + }; + KRATOS_DEBUG_ERROR_IF(std::any_of(rBs.begin() + 1, rBs.end(), has_inconsistent_sizes)) + << "Cannot calculate the internal force vector: B-matrices have different sizes\n"; +} + +void CheckThatStressVectorsHaveSameSize(const std::vector& rStressVectors) +{ + auto has_inconsistent_size = [size = rStressVectors.front().size()](const auto& rVector) { + return rVector.size() != size; + }; + KRATOS_DEBUG_ERROR_IF(std::any_of(rStressVectors.begin() + 1, rStressVectors.end(), has_inconsistent_size)) + << "Cannot calculate the internal force vector: stress vectors have different sizes\n"; +} + +void CheckThatTransposedMatrixVectorProductCanBeCalculated(const Matrix& rFirstB, const Vector& rFirstStressVector) +{ + // B-matrix will be transposed, so check against its number of rows rather than its number of columns! + KRATOS_DEBUG_ERROR_IF(rFirstB.size1() != rFirstStressVector.size()) + << "Cannot calculate the internal force vector: matrix-vector product cannot be calculated " + "due to size mismatch\n"; +} + +void CheckInputOfCalculateInternalForceVector(const std::vector& rBs, + const std::vector& rStressVectors, + const std::vector& rIntegrationCoefficients) +{ + CheckThatVectorsHaveSameSizeAndAreNotEmpty(rBs, rStressVectors, rIntegrationCoefficients); + CheckThatBMatricesHaveSameSizes(rBs); + CheckThatStressVectorsHaveSameSize(rStressVectors); + CheckThatTransposedMatrixVectorProductCanBeCalculated(rBs.front(), rStressVectors.front()); +} + +} // namespace + namespace Kratos { @@ -83,4 +140,18 @@ Matrix GeoEquationOfMotionUtilities::CalculateStiffnessMatrix(const std::vector< return result; } +Vector GeoEquationOfMotionUtilities::CalculateInternalForceVector(const std::vector& rBs, + const std::vector& rStressVectors, + const std::vector& rIntegrationCoefficients) +{ + CheckInputOfCalculateInternalForceVector(rBs, rStressVectors, rIntegrationCoefficients); + + auto result = Vector{ZeroVector{rBs.front().size2()}}; + for (auto i = std::size_t{0}; i < rBs.size(); ++i) { + result += prod(trans(rBs[i]), rStressVectors[i]) * rIntegrationCoefficients[i]; + } + + return result; +} + } /* namespace Kratos.*/ diff --git a/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h b/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h index f74a0df2ae56..a3ffce77cf98 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h +++ b/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h @@ -49,5 +49,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoEquationOfMotionUtilities const std::vector& rConstitutiveMatrices, const std::vector& rIntegrationCoefficients); + static Vector CalculateInternalForceVector(const std::vector& rBs, + const std::vector& rStressVectors, + const std::vector& rIntegrationCoefficients); + }; /* Class GeoTransportEquationUtilities*/ } /* namespace Kratos.*/ diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp index f40b7c997fe2..e74b2b0ee448 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp @@ -14,6 +14,7 @@ #include "geo_mechanics_application_variables.h" #include "geo_mechanics_fast_suite.h" #include "tests/cpp_tests/test_utilities/model_setup_utilities.h" + #include using namespace Kratos; @@ -183,4 +184,90 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateStiffnessMatrixGivesCorrectResults, KratosGeo KRATOS_CHECK_MATRIX_NEAR(stiffness_matrix, expected_stiffness_matrix, 1e-4) } +KRATOS_TEST_CASE_IN_SUITE(TheInternalForceVectorIsTheIntegralOfBTransposedTimesSigma, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto b_matrix = Matrix{ScalarMatrix{2, 8, 1.0}}; + const auto b_matrices = std::vector{b_matrix, b_matrix}; + const auto stress_vector = Vector{ScalarVector{2, 1.0}}; + const auto stress_vectors = std::vector{stress_vector, stress_vector}; + const auto integration_coefficients = std::vector{0.25, 0.4}; + + const auto expected_internal_force_vector = Vector{ScalarVector{8, 1.3}}; + constexpr auto relative_tolerance = 1.0e-6; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(GeoEquationOfMotionUtilities::CalculateInternalForceVector( + b_matrices, stress_vectors, integration_coefficients), + expected_internal_force_vector, relative_tolerance) +} + +// The following tests only raise errors when using debug builds +#ifdef KRATOS_DEBUG + +KRATOS_TEST_CASE_IN_SUITE(CalculatingTheInternalForceVectorFailsWhenTheInputVectorsHaveDifferentSizes, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto b_matrix = Matrix{ScalarMatrix{2, 8, 1.0}}; + const auto b_matrices = std::vector{b_matrix}; // Error: missing one matrix + const auto stress_vector = Vector{ScalarVector{2, 1.0}}; + const auto stress_vectors = std::vector{stress_vector, stress_vector}; + const auto integration_coefficients = std::vector{0.25, 0.4}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + GeoEquationOfMotionUtilities::CalculateInternalForceVector(b_matrices, stress_vectors, integration_coefficients), + "Cannot calculate the internal force vector: input vectors have different sizes") +} + +KRATOS_TEST_CASE_IN_SUITE(CalculatingTheInternalForceVectorFailsWhenAllInputVectorsAreEmpty, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto b_matrices = std::vector{}; + const auto stress_vectors = std::vector{}; + const auto integration_coefficients = std::vector{}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + GeoEquationOfMotionUtilities::CalculateInternalForceVector(b_matrices, stress_vectors, integration_coefficients), + "Cannot calculate the internal force vector: input vectors are empty") +} + +KRATOS_TEST_CASE_IN_SUITE(CalculatingTheInternalForceVectorFailsWhenBMatricesHaveDifferentSizes, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto b_matrices = std::vector{ScalarMatrix{2, 8, 1.0}, ScalarMatrix{1, 8, 1.0}}; // Error: matrices have different numbers of rows + const auto stress_vector = Vector{ScalarVector{2, 1.0}}; + const auto stress_vectors = std::vector{stress_vector, stress_vector}; + const auto integration_coefficients = std::vector{0.25, 0.4}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + GeoEquationOfMotionUtilities::CalculateInternalForceVector(b_matrices, stress_vectors, integration_coefficients), + "Cannot calculate the internal force vector: B-matrices have different sizes") +} + +KRATOS_TEST_CASE_IN_SUITE(CalculatingTheInternalForceVectorFailsWhenStressVectorsHaveDifferentSizes, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto b_matrix = Matrix{ScalarMatrix{2, 8, 1.0}}; + const auto b_matrices = std::vector{b_matrix, b_matrix}; + const auto stress_vectors = std::vector{ScalarVector{2, 1.0}, ScalarVector{3, 1.0}}; // Error: vectors have different sizes + const auto integration_coefficients = std::vector{0.25, 0.4}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + GeoEquationOfMotionUtilities::CalculateInternalForceVector(b_matrices, stress_vectors, integration_coefficients), + "Cannot calculate the internal force vector: stress vectors have different sizes") +} + +KRATOS_TEST_CASE_IN_SUITE(CalculatingTheInternalForceVectorFailsWhenTheMatrixVectorProductCantBeComputed, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + // Error: transpose of the B-matrix has more columns (3) than the number of stress components (2) + const auto b_matrix = Matrix{ScalarMatrix{3, 8, 1.0}}; + const auto b_matrices = std::vector{b_matrix, b_matrix}; + const auto stress_vector = Vector{ScalarVector{2, 1.0}}; + const auto stress_vectors = std::vector{stress_vector, stress_vector}; + const auto integration_coefficients = std::vector{0.25, 0.4}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + GeoEquationOfMotionUtilities::CalculateInternalForceVector(b_matrices, stress_vectors, integration_coefficients), "Cannot calculate the internal force vector: matrix-vector product cannot be calculated due to size mismatch") +} + +#endif + } // namespace Kratos::Testing \ No newline at end of file