- added CalculateOnIntegrationPoints(Vector,.. for LinearTimoshenkoBe…

…amElement2D2N (#12948) - added TIMOSHENKO_BEAM_PRESTRESS_PK2 - added unit test for this function - moved InitializeConstitutiveLawValuesForStressCalculation to StructuralMechanicsElementUtilities
KratosMultiphysics · Dec 19, 2024 · cab78c8 · cab78c8
1 parent 62475a1
commit cab78c8
Show file tree

Hide file tree

Showing 11 changed files with 180 additions and 45 deletions.
diff --git a/...cturalMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp b/...cturalMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp
@@ -850,18 +850,11 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
         const auto &r_geometry = GetGeometry();
 
         ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
-        auto &r_cl_options = cl_values.GetOptions();
-        r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
-        r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);
-
         const double length = CalculateLength();
         const double Phi    = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);
-
-        // Let's initialize the cl values
         VectorType strain_vector(strain_size), stress_vector(strain_size);
-        strain_vector.clear();
-        cl_values.SetStrainVector(strain_vector);
-        cl_values.SetStressVector(stress_vector);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);
+
         VectorType nodal_values(mat_size);
         GetNodalValuesVector(nodal_values);
 
@@ -880,18 +873,10 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
         const auto &r_geometry = GetGeometry();
 
         ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
-        auto &r_cl_options = cl_values.GetOptions();
-        r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
-        r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);
-
         const double length = CalculateLength();
         const double Phi    = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);
-
-        // Let's initialize the cl values
         VectorType strain_vector(strain_size), stress_vector(strain_size);
-        strain_vector.clear();
-        cl_values.SetStrainVector(strain_vector);
-        cl_values.SetStressVector(stress_vector);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);
         VectorType nodal_values(mat_size);
         GetNodalValuesVector(nodal_values);
 
@@ -910,18 +895,10 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
         const auto &r_geometry = GetGeometry();
 
         ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
-        auto &r_cl_options = cl_values.GetOptions();
-        r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
-        r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);
-
         const double length = CalculateLength();
         const double Phi    = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);
-
-        // Let's initialize the cl values
         VectorType strain_vector(strain_size), stress_vector(strain_size);
-        strain_vector.clear();
-        cl_values.SetStrainVector(strain_vector);
-        cl_values.SetStressVector(stress_vector);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);
         VectorType nodal_values(mat_size);
         GetNodalValuesVector(nodal_values);
 
@@ -964,6 +941,44 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
     }
 }
 
+void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
+    const Variable<Vector>& rVariable,
+    std::vector<Vector>& rOutput,
+    const ProcessInfo& rProcessInfo
+    )
+{
+    const auto& integration_points = IntegrationPoints(GetIntegrationMethod());
+    const SizeType strain_size = mConstitutiveLawVector[0]->GetStrainSize();
+    const SizeType mat_size = GetDoFsPerNode() * GetGeometry().size();
+    rOutput.resize(integration_points.size());
+    const auto &r_props = GetProperties();
+
+    if (rVariable == PK2_STRESS_VECTOR) {
+
+        const auto &r_geometry = GetGeometry();
+
+        ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
+        const double length = CalculateLength();
+        const double Phi    = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);
+
+        VectorType strain_vector(strain_size), stress_vector(strain_size);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);
+
+        VectorType nodal_values(mat_size);
+        GetNodalValuesVector(nodal_values);
+
+        // Loop over the integration points
+        for (SizeType integration_point = 0; integration_point < integration_points.size(); ++integration_point) {
+            CalculateGeneralizedStrainsVector(strain_vector, length, Phi, integration_points[integration_point].X(), nodal_values);
+            mConstitutiveLawVector[integration_point]->CalculateMaterialResponsePK2(cl_values);
+            auto stress_vector = cl_values.GetStressVector();
+            if ( this->GetProperties().Has(TIMOSHENKO_BEAM_PRESTRESS_PK2)) {
+                stress_vector += this->GetProperties()[TIMOSHENKO_BEAM_PRESTRESS_PK2];
+            }
+            rOutput[integration_point] = stress_vector;
+        }
+    }
+}
 /***********************************************************************************/
 /***********************************************************************************/
 

diff --git a/...ructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h b/...ructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h
@@ -387,6 +387,18 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTimoshenkoBeamElement2D
         const ProcessInfo& rCurrentProcessInfo
         ) override;
 
+    /**
+     * @brief Calculate a double Variable on the Element Constitutive Law
+     * @param rVariable The variable we want to get
+     * @param rOutput The values obtained in the integration points
+     * @param rCurrentProcessInfo the current process info instance
+     */
+    void CalculateOnIntegrationPoints(
+        const Variable<Vector>& rVariable,
+        std::vector<Vector>& rOutput,
+        const ProcessInfo& rCurrentProcessInfo
+        ) override;
+
     /**
      * @brief Get on rVariable Constitutive Law from the element
      * @param rVariable The variable we want to get

diff --git a/...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp b/...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp
@@ -733,7 +733,7 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
         ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
         VectorType strain_vector(1), stress_vector(1);
         MatrixType C(1,1);
-        InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, C);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C);
 
         const double length = CalculateLength();
         SystemSizeBoundedArrayType nodal_values(SystemSize);
@@ -756,7 +756,7 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
         ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
         VectorType strain_vector(1), stress_vector(1);
         MatrixType C(1,1);
-        InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, C);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C);
 
         const double length = CalculateLength();
         SystemSizeBoundedArrayType nodal_values(SystemSize);
@@ -773,19 +773,6 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
     }
 }
 
-template <SizeType TDimension, SizeType TNNodes>
-void LinearTrussElement<TDimension, TNNodes>::InitializeConstitutiveLawValues(ConstitutiveLaw::Parameters& rValues,
-    VectorType& rStrainVector, VectorType& rStressVector, MatrixType rConstitutiveMatrix)
-{
-    rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
-    rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);
-
-    rStrainVector.clear();
-    rValues.SetStrainVector(rStrainVector);
-    rValues.SetStressVector(rStressVector);
-    rValues.SetConstitutiveMatrix(rConstitutiveMatrix);
-}
-
 /***********************************************************************************/
 /***********************************************************************************/
 
@@ -803,7 +790,7 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
         ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
         VectorType strain_vector(1), stress_vector(1);
         MatrixType constitutive_matrix(1,1);
-        InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, constitutive_matrix);
+        StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, constitutive_matrix);
 
         const double length = CalculateLength();
         SystemSizeBoundedArrayType nodal_values(SystemSize);

diff --git a/...ions/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h b/...ions/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h
@@ -477,8 +477,6 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement
     ///@}
     ///@name Private Operations
     ///@{
-    void InitializeConstitutiveLawValues(ConstitutiveLaw::Parameters& rValues,
-    VectorType& rStrainVector, VectorType& rStressVector, MatrixType rConstitutiveMatrix);
 
     ///@}
     ///@name Private  Access

diff --git a/.../StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp b/.../StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp
@@ -91,6 +91,7 @@ PYBIND11_MODULE(KratosStructuralMechanicsApplication,m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,I33)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,IZ)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,BEAM_INITIAL_STRAIN_VECTOR)
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,TIMOSHENKO_BEAM_PRESTRESS_PK2)
 
     // semi rigid beam variables
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROTATIONAL_STIFFNESS_AXIS_2)

diff --git a/...tructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp b/...tructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp
@@ -635,6 +635,24 @@ double CalculatePhi(const Properties& rProperties, const double L)
 /***********************************************************************************/
 /***********************************************************************************/
 
+void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
+    Vector& rStrainVector, Vector& rStressVector, Matrix& rConstitutiveMatrix)
+{
+    InitializeConstitutiveLawValuesForStressCalculation(rValues, rStrainVector, rStressVector);
+    rValues.SetConstitutiveMatrix(rConstitutiveMatrix);
+}
+
+void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
+    Vector& rStrainVector, Vector& rStressVector)
+{
+    rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
+    rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);
+
+    rStrainVector.clear();
+    rValues.SetStrainVector(rStrainVector);
+    rValues.SetStressVector(rStressVector);
+}
+
 } // namespace StructuralMechanicsElementUtilities.
 }  // namespace Kratos.
 

diff --git a/.../StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h b/.../StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h
@@ -350,5 +350,11 @@ double GetReferenceRotationAngle2D3NBeam(const GeometryType &rGeometry);
  */
 KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) double CalculatePhi(const Properties& rProperties, const double L);
 
+void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
+    Vector& rStrainVector, Vector& rStressVector, Matrix& rConstitutiveMatrix);
+
+void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
+    Vector& rStrainVector, Vector& rStressVector);
+
 } // namespace StructuralMechanicsElementUtilities.
 }  // namespace Kratos.
diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp
@@ -347,6 +347,8 @@ void KratosStructuralMechanicsApplication::Register() {
     KRATOS_REGISTER_VARIABLE(I33)
     KRATOS_REGISTER_VARIABLE(LUMPED_MASS_ROTATION_COEFFICIENT)
     KRATOS_REGISTER_VARIABLE(BEAM_INITIAL_STRAIN_VECTOR)
+    KRATOS_REGISTER_VARIABLE(TIMOSHENKO_BEAM_PRESTRESS_PK2)
+
 
     // semi rigid beam variables
     KRATOS_REGISTER_VARIABLE(ROTATIONAL_STIFFNESS_AXIS_2)

diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp
@@ -71,6 +71,8 @@ KRATOS_CREATE_VARIABLE(double, I22)
 KRATOS_CREATE_VARIABLE(double, I33)
 KRATOS_CREATE_VARIABLE(double, LUMPED_MASS_ROTATION_COEFFICIENT)
 KRATOS_CREATE_VARIABLE(Vector, BEAM_INITIAL_STRAIN_VECTOR)
+KRATOS_CREATE_VARIABLE(Vector, TIMOSHENKO_BEAM_PRESTRESS_PK2)
+
 
 // semi rigid beam variables
 KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_2)

diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h b/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h
@@ -92,6 +92,7 @@ namespace Kratos
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, I33)
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, LUMPED_MASS_ROTATION_COEFFICIENT)
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, BEAM_INITIAL_STRAIN_VECTOR)
+    KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, TIMOSHENKO_BEAM_PRESTRESS_PK2)
 
 
     // Semi rigid beam variables

diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp
@@ -0,0 +1,93 @@
+// KRATOS  ___|  |                   |                   |
+//       \___ \  __|  __| |   |  __| __| |   |  __| _` | |
+//             | |   |    |   | (    |   |   | |   (   | |
+//       _____/ \__|_|   \__,_|\___|\__|\__,_|_|  \__,_|_| MECHANICS
+//
+//  License:         BSD License
+//                   license: StructuralMechanicsApplication/license.txt
+//
+//  Main authors:    Gennady Markelov
+//
+
+// Project includes
+#include "containers/model.h"
+#include "structural_mechanics_fast_suite.h"
+#include "structural_mechanics_application_variables.h"
+
+#include <utility>
+#include <boost/numeric/ublas/assignment.hpp>
+
+namespace Kratos::Testing
+{
+
+void CreateBeamModel2N(std::string TimoshenkoBeamElementName)
+    {
+        Model current_model;
+        auto &r_model_part = current_model.CreateModelPart("ModelPart",1);
+        r_model_part.GetProcessInfo().SetValue(DOMAIN_SIZE, 3);
+        r_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);
+        r_model_part.AddNodalSolutionStepVariable(ROTATION_Z);
+
+        // Set the element properties
+        auto p_elem_prop = r_model_part.CreateNewProperties(0);
+        constexpr auto youngs_modulus = 2.0e+06;
+        p_elem_prop->SetValue(YOUNG_MODULUS, youngs_modulus);
+        constexpr auto cross_area = 1.0;
+        p_elem_prop->SetValue(CROSS_AREA, cross_area);
+        constexpr auto i33 = 1.0;
+        p_elem_prop->SetValue(I33, i33);
+        constexpr auto area_effective_y = 1.0;
+        p_elem_prop->SetValue(AREA_EFFECTIVE_Y, area_effective_y);
+
+        const auto &r_clone_cl = KratosComponents<ConstitutiveLaw>::Get("TimoshenkoBeamElasticConstitutiveLaw");
+        p_elem_prop->SetValue(CONSTITUTIVE_LAW, r_clone_cl.Clone());
+
+        // Create the test element
+        constexpr double directional_length = 2.0;
+        auto p_node_1 = r_model_part.CreateNewNode(1, 0.0, 0.0, 0.0);
+        auto p_node_2 = r_model_part.CreateNewNode(2, directional_length, directional_length, directional_length);
+
+        for (auto& r_node : r_model_part.Nodes()){
+            r_node.AddDof(DISPLACEMENT_X);
+            r_node.AddDof(DISPLACEMENT_Y);
+            r_node.AddDof(DISPLACEMENT_Z);
+            r_node.AddDof(ROTATION_Z);
+        }
+
+        std::vector<ModelPart::IndexType> element_nodes {1,2};
+        auto p_element = r_model_part.CreateNewElement(std::move(TimoshenkoBeamElementName), 1, element_nodes, p_elem_prop);
+        const auto& r_process_info = r_model_part.GetProcessInfo();
+        p_element->Initialize(r_process_info); // Initialize the element to initialize the constitutive law
+
+        constexpr auto induced_strain = 0.1;
+        p_element->GetGeometry()[1].FastGetSolutionStepValue(DISPLACEMENT) += ScalarVector(3, induced_strain * directional_length);
+        p_element->GetGeometry()[1].FastGetSolutionStepValue(ROTATION_Z) += 0.1;
+
+        std::vector<Vector> stress_vectors;
+        p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vectors, r_process_info);
+
+        constexpr auto expected_stress = induced_strain * youngs_modulus;
+        constexpr auto expected_shear_stress = -37500.0;
+
+        constexpr auto tolerance = 1.0e-5;
+        Vector expected_stress_vector(3);
+        expected_stress_vector <<= expected_stress, 29631.5,expected_shear_stress;
+        KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[0], tolerance);
+        expected_stress_vector <<= expected_stress, 70710.7,expected_shear_stress;
+        KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[1], tolerance);
+        expected_stress_vector <<= expected_stress, 111790,expected_shear_stress;
+        KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[2], tolerance);
+
+        Vector pre_stress(3);
+        pre_stress <<= 1.0e5, 1.0e4, 1.0e3;
+        p_element->GetProperties().SetValue(TIMOSHENKO_BEAM_PRESTRESS_PK2, pre_stress);
+        p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vectors, r_process_info);
+        expected_stress_vector += pre_stress;
+        KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[2], tolerance);
+    }
+
+    KRATOS_TEST_CASE_IN_SUITE(LinearTimoshenkoBeam2D2N_CalculatesPK2Stress, KratosStructuralMechanicsFastSuite)
+    {
+        CreateBeamModel2N("LinearTimoshenkoBeamElement2D2N");
+    }
+}