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.

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 <algorithm>
+
+namespace
+{
+
+using namespace Kratos;
+
+void CheckThatVectorsHaveSameSizeAndAreNotEmpty(const std::vector<Matrix>& rBs,
+                                                const std::vector<Vector>& rStressVectors,
+                                                const std::vector<double>& 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<Matrix>& 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<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<Matrix>& rBs,
+                                              const std::vector<Vector>& rStressVectors,
+                                              const std::vector<double>& 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<Matrix>& rBs,
+                                                                  const std::vector<Vector>& rStressVectors,
+                                                                  const std::vector<double>& 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.*/

applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h

@@ -49,5 +49,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoEquationOfMotionUtilities
                                            const std::vector<Matrix>& rConstitutiveMatrices,
                                            const std::vector<double>& rIntegrationCoefficients);
 
+    static Vector CalculateInternalForceVector(const std::vector<Matrix>& rBs,
+                                               const std::vector<Vector>& rStressVectors,
+                                               const std::vector<double>& rIntegrationCoefficients);
+
 }; /* Class GeoTransportEquationUtilities*/
 } /* namespace Kratos.*/

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 <boost/numeric/ublas/assignment.hpp>
 
 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<Matrix>{b_matrix, b_matrix};
+    const auto stress_vector            = Vector{ScalarVector{2, 1.0}};
+    const auto stress_vectors           = std::vector<Vector>{stress_vector, stress_vector};
+    const auto integration_coefficients = std::vector<double>{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<Matrix>{b_matrix}; // Error: missing one matrix
+    const auto stress_vector  = Vector{ScalarVector{2, 1.0}};
+    const auto stress_vectors = std::vector<Vector>{stress_vector, stress_vector};
+    const auto integration_coefficients = std::vector<double>{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<Matrix>{};
+    const auto stress_vectors           = std::vector<Vector>{};
+    const auto integration_coefficients = std::vector<double>{};
+
+    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<Matrix>{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<Vector>{stress_vector, stress_vector};
+    const auto integration_coefficients = std::vector<double>{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<Matrix>{b_matrix, b_matrix};
+    const auto stress_vectors = std::vector<Vector>{ScalarVector{2, 1.0}, ScalarVector{3, 1.0}}; // Error: vectors have different sizes
+    const auto integration_coefficients = std::vector<double>{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<Matrix>{b_matrix, b_matrix};
+    const auto stress_vector            = Vector{ScalarVector{2, 1.0}};
+    const auto stress_vectors           = std::vector<Vector>{stress_vector, stress_vector};
+    const auto integration_coefficients = std::vector<double>{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