[StructuralMechanicsApplication] Add PK2_STRESS_VECTOR in LinearTrussElement2D by markelov208 · Pull Request #12915 · KratosMultiphysics/Kratos

markelov208 · 2024-12-09T11:37:13Z

📝 Description
A brief description of the PR.

Added CalculateOnIntegrationPoints function with PK2_STRESS_VECTOR and TRUSS_PRESTRESS_PK2 for LinearTrussElement<3D/2D>
added two unit tests for this function (based on TrussElementLinear3D2N_CalculatesPK2Stress)

WPK4FEM

Hi Gennady,
Thank you for the implementation. On the implementation itself I have no comments, it follows what the element already is doing. Some nagging about following the style.
However, there is 1 thing to be aware of. To create also the 3D variants of this 2 and 3 node element Alejandro has a pull request pending ( #12894 ) that removes the file that you have done the implementation in. It may be wise to wait with your PR until Aljandro has merged and then put your implementation in the new file for both 2 and 3D. That would save doing it again.
Regards, Wijtze Pieter

WPK4FEM · 2024-12-10T13:27:15Z

...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp

+        SystemSizeBoundedArrayType B;
+
+        // Loop over the integration points
+        for (SizeType IP = 0; IP < integration_points.size(); ++IP) {


The Kratos style guide ( https://github.com/KratosMultiphysics/Kratos/wiki/Style-Guide ) says lower case with underscores for local variables. Here that would result in:, c, b, integration_point i.s.o. C, B, IP

Thank you for the comment. Fixed names and replaced C/B with longer names.

WPK4FEM · 2024-12-10T13:29:24Z

applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp

+        p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vector, r_process_info);
+        KRATOS_EXPECT_DOUBLE_EQ(expected_stress + pre_stress, stress_vector[0][0]);
+    }
+KRATOS_TEST_CASE_IN_SUITE(LienarTrussElement2D_CalculatesPK2Stress, KratosStructuralMechanicsFastSuite)


Lienar --> Linear

Many thanks. Fixed the typo.

markelov208 · 2024-12-10T21:06:48Z

Hi Gennady, Thank you for the implementation. On the implementation itself I have no comments, it follows what the element already is doing. Some nagging about following the style. However, there is 1 thing to be aware of. To create also the 3D variants of this 2 and 3 node element Alejandro has a pull request pending ( #12894 ) that removes the file that you have done the implementation in. It may be wise to wait with your PR until Aljandro has merged and then put your implementation in the new file for both 2 and 3D. That would save doing it again. Regards, Wijtze Pieter

Hi Wijtze-Pieter, thank you for the info. In any case, this PR would be merged only when Alejandro approves it.

rfaasse

Nice addition of PK2 stress to the new element with tests, I've got a few suggestions, but none of them blocking

rfaasse · 2024-12-10T12:39:51Z

...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp

+    if (rVariable == PK2_STRESS_VECTOR) {
+        ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
+        cl_values.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
+        cl_values.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true);


Since we only need the stress vector in the end, do you know if it's possible to set COMPUTE_CONSTITUTIVE_TENSOR to false?

yes, if the stiffness (E) is not required, you can set this to false

Actually, if I am not wrong, the stress could be retrieved by the CalculateOnIntegrationPoints(AXIAL_FORCE, --, --) and then divided by the GetPRoperties()[CROSS_AREA] right? I mean implementing the new CalculateOnIntegrationPoints(PK2_STRESS) and calling the CalculateOnIntegrationPoints(AXIAL_FORCE, --, --) inside

Richard, you are absolutely right.
Alejandro, thank you for the info on COMPUTE_CONSTITUTIVE_TENSOR. I also made it false for other CalculateOnIntegrationPoints functions and there are no any failed tests.
You are right it is possible to re-use CalculateOnIntegrationPoints(AXIAL_FORCE,); however, 1) I have added TRUSS_PRESTRESS_PK2, 2) I would like to keep theses functions separately because CalculateOnIntegrationPoints(AXIAL_FORCE,) could be changed in future, for example, PK2 could be replaced with other stress.
Next, Aljandro please could you advice how to get the stress vector for TImoshenko beams? As far as I can see there are two areas, CROSS_AREA and AREA_EFFECTIVE_Y.

The Timoshenko constitutive laws compute the axial force, moment and shear forces. To obtain the stress vector one should compose the bending and axial contributions as the normal stress and the shear stresses with respecto to the shear force and shear area.

The CROSS_AREA is the area of the cross section used in axial loading. The AREA_EFFECTIVE_Y is the reduced shear area for evaluating the shear forces.

applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp

rfaasse · 2024-12-11T08:12:38Z

applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp

+        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);
+
+        // 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);
+        const auto &r_clone_cl = KratosComponents<ConstitutiveLaw>::Get("TrussConstitutiveLaw");
+        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);
+
+        AddDisplacementDofsElement(r_model_part);
+
+        std::vector<ModelPart::IndexType> element_nodes {1,2};
+        auto p_element = r_model_part.CreateNewElement("LinearTrussElement2D2N", 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);


If I'm correct, the setup of this test is almost identical to the setup of the previous unit test. I'd propose to create a setup function which takes the name of the element as an input (so "LinearTrussElement2D2N" or "TrussLinearElement3D2N"), to avoid some duplication.

It is a very good proposal. Currently, one more test is added. Hopefully, it will be extended for Timoshenko beams.

rfaasse · 2024-12-11T08:14:50Z

...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp

+        ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
+        cl_values.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
+        cl_values.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true);
+
+        const double length = CalculateLength();
+
+        // Let's initialize the cl values
+        VectorType strain_vector(1), stress_vector(1);
+        MatrixType constitutive_matrix(1,1);
+        strain_vector.clear();
+        cl_values.SetStrainVector(strain_vector);
+        cl_values.SetStressVector(stress_vector);
+        cl_values.SetConstitutiveMatrix(constitutive_matrix);
+        SystemSizeBoundedArrayType nodal_values(SystemSize);
+        GetNodalValuesVector(nodal_values);
+
+        SystemSizeBoundedArrayType dN_dX;
+
+        // Loop over the integration points
+        for (SizeType integration_point = 0; integration_point < integration_points.size(); ++integration_point) {
+            GetFirstDerivativesShapeFunctionsValues(dN_dX, length, integration_points[integration_point].X());
+
+            strain_vector[0] = inner_prod(dN_dX, nodal_values);
+
+            mConstitutiveLawVector[integration_point]->CalculateMaterialResponsePK2(cl_values);
+            auto stress = cl_values.GetStressVector()[0];
+            if (GetProperties().Has(TRUSS_PRESTRESS_PK2)) {
+                stress += GetProperties()[TRUSS_PRESTRESS_PK2];
+            }
+
+            rOutput[integration_point] = ScalarVector(1, stress);


If this part of code is identical to another part in this class, we could extract a function for the common functionality

A few lines are extracted as a function.

AlejandroCornejo · 2024-12-11T08:34:55Z

HI! I just merged the linear_truss_element generic 2D/3D for 2N and 3N in the other PR. I'll review the changes in this PR soon

AlejandroCornejo · 2024-12-11T08:39:23Z

Actually, I believe that these additions should be done to the new linear_truss_element that I just added. Now the element in 2D has been generalized to 3D also. Please move the implementation to the new file.

AlejandroCornejo · 2024-12-11T08:40:39Z

...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp

 /***********************************************************************************/
 /***********************************************************************************/

+template<SizeType TNNodes>


now you'll need anothe template entry

AlejandroCornejo · 2024-12-11T08:41:16Z

...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp

+    if (rVariable == PK2_STRESS_VECTOR) {
+        ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
+        cl_values.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS             , true);
+        cl_values.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true);


yes, if the stiffness (E) is not required, you can set this to false

# Conflicts: # applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp

markelov208

Hi everybody, thank you very much for the suggestions and info.

rfaasse

Nicely reduced the duplication, thanks for the effort! from my side this is ready to go, but I leave the approval to @AlejandroCornejo / @loumalouomega 👍

markelov208 added 2 commits December 9, 2024 12:35

added CalculateOnIntegrationPoints for PK2_STRESS_VECTOR

cd818ae

added unit test for LinearTrussElement2D_PK2

30d3a84

markelov208 self-assigned this Dec 9, 2024

markelov208 requested a review from a team as a code owner December 9, 2024 11:37

used correct name LinearTrussElement2D2N

3195d1c

markelov208 changed the title ~~Geo/12868 pk2 linear truss element2 d~~ [StructuralMechanicsApplication] Add PK2_STRESS_VECTOR in LinearTrussElement2D Dec 9, 2024

markelov208 requested review from WPK4FEM and rfaasse December 9, 2024 16:13

WPK4FEM reviewed Dec 10, 2024

View reviewed changes

markelov208 added 2 commits December 10, 2024 19:16

used Kratos style

160c4c1

fixed typo

de1b040

rfaasse reviewed Dec 11, 2024

View reviewed changes

AlejandroCornejo reviewed Dec 11, 2024

View reviewed changes

markelov208 added 4 commits December 11, 2024 16:35

Merge branch 'master' into geo/12868-pk2_LinearTrussElement2D

d698313

# Conflicts: # applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp

added InitializeConstitutiveLawValues

4f7e8c2

added CreateTrussModel2N_and_CheckPK2Stress

e21c790

removed unused definition

4795532

markelov208 requested review from AlejandroCornejo, WPK4FEM and rfaasse December 11, 2024 20:42

markelov208 commented Dec 11, 2024

View reviewed changes

rfaasse reviewed Dec 12, 2024

View reviewed changes

AlejandroCornejo approved these changes Dec 13, 2024

View reviewed changes

markelov208 merged commit 5958552 into master Dec 13, 2024

markelov208 deleted the geo/12868-pk2_LinearTrussElement2D branch December 13, 2024 10:44

Conversation

markelov208 commented Dec 9, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

WPK4FEM left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

markelov208 commented Dec 10, 2024

Uh oh!

rfaasse left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

AlejandroCornejo commented Dec 11, 2024

Uh oh!

AlejandroCornejo commented Dec 11, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

markelov208 left a comment

Choose a reason for hiding this comment

Uh oh!

rfaasse left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

markelov208 commented Dec 9, 2024 •

edited

Loading

AlejandroCornejo commented Dec 11, 2024 •

edited

Loading