Merge branch 'master' into scripts/docker/update-rocky-python-

loumalouomega · web-flow · commit 8e37b282e593 · 2025-11-28T15:19:19.000+01:00
diff --git a/applications/GeoMechanicsApplication/custom_elements/interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/interface_element.cpp
@@ -198,8 +198,11 @@ void InterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
         std::vector<std::optional<Vector>> interface_nodal_cauchy_stresses(interface_node_ids.size());
         auto&               r_neighbour_element = this->GetValue(NEIGHBOUR_ELEMENTS).front();
         std::vector<Vector> neighbour_cauchy_stresses;
+        // Note that the interface elements don't account for water pressures yet. Consequently,
+        // we need to consider the total stresses rather than the effective stresses to calculate
+        // the appropriate prestresses to be applied.
         r_neighbour_element.CalculateOnIntegrationPoints(
-            CAUCHY_STRESS_VECTOR, neighbour_cauchy_stresses, rCurrentProcessInfo);
+            TOTAL_STRESS_VECTOR, neighbour_cauchy_stresses, rCurrentProcessInfo);
         interface_nodal_cauchy_stresses = ExtrapolationUtilities::CalculateNodalVectors(
             interface_node_ids, r_neighbour_element.GetGeometry(), r_neighbour_element.GetIntegrationMethod(),
             neighbour_cauchy_stresses, r_neighbour_element.Id());
diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_interface_element.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_interface_element.cpp
@@ -16,6 +16,7 @@
 #include "custom_constitutive/interface_plane_strain.h"
 #include "custom_constitutive/interface_three_dimensional_surface.h"
 #include "custom_constitutive/plane_strain.h"
+#include "custom_constitutive/three_dimensional.h"
 #include "custom_elements/interface_element.h"
 #include "custom_elements/interface_stress_state.h"
 #include "custom_geometries/interface_geometry.h"
@@ -311,6 +312,75 @@ Matrix ExpectedLeftHandSideForTriangleElement()
     return expected_left_hand_side;
 }
 
+class MockElementWithTotalStressVectors : public Element
+{
+public:
+    MockElementWithTotalStressVectors(std::size_t                     ElementId,
+                                      std::shared_ptr<Geometry<Node>> pGeometry,
+                                      Properties::Pointer             pProperties);
+
+    void SetValuesOnIntegrationPoints(const Variable<Vector>&    rVariable,
+                                      const std::vector<Vector>& rValues,
+                                      const ProcessInfo&) override;
+    using Element::SetValuesOnIntegrationPoints;
+
+    void CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
+                                      std::vector<Vector>&    rOutput,
+                                      const ProcessInfo&) override;
+    using Element::CalculateOnIntegrationPoints;
+
+    IntegrationMethod GetIntegrationMethod() const override;
+    void              SetIntegrationMethod(IntegrationMethod CustomIntegrationMethod);
+
+private:
+    std::vector<Vector>              mTotalStressVectors;
+    std::optional<IntegrationMethod> mOptionalCustomIntegrationMethod;
+};
+
+MockElementWithTotalStressVectors::MockElementWithTotalStressVectors(std::size_t ElementId,
+                                                                     std::shared_ptr<Geometry<Node>> pGeometry,
+                                                                     Properties::Pointer pProperties)
+    : Element{ElementId, std::move(pGeometry), std::move(pProperties)}
+{
+}
+
+void MockElementWithTotalStressVectors::SetValuesOnIntegrationPoints(const Variable<Vector>& rVariable,
+                                                                     const std::vector<Vector>& rValues,
+                                                                     const ProcessInfo&)
+{
+    KRATOS_DEBUG_ERROR_IF_NOT(rVariable == TOTAL_STRESS_VECTOR)
+        << "This mock element can only set total stress vectors\n";
+
+    mTotalStressVectors = rValues;
+}
+
+void MockElementWithTotalStressVectors::CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
+                                                                     std::vector<Vector>& rOutput,
+                                                                     const ProcessInfo&)
+{
+    KRATOS_DEBUG_ERROR_IF_NOT(rVariable == TOTAL_STRESS_VECTOR)
+        << "This mock element can only calculate total stress vectors\n";
+
+    rOutput = mTotalStressVectors;
+}
+
+GeometryData::IntegrationMethod MockElementWithTotalStressVectors::GetIntegrationMethod() const
+{
+    return mOptionalCustomIntegrationMethod.value_or(GetGeometry().GetDefaultIntegrationMethod());
+}
+
+void MockElementWithTotalStressVectors::SetIntegrationMethod(IntegrationMethod CustomIntegrationMethod)
+{
+    mOptionalCustomIntegrationMethod = CustomIntegrationMethod;
+}
+
+template <typename DerivedElementPtrType>
+GlobalPointersVector<Element> MakeElementGlobalPtrContainerWith(const DerivedElementPtrType& rpElement)
+{
+    auto p_element = Element::Pointer{rpElement};
+    return {{GlobalPointer<Element>{p_element}}};
+}
+
 } // namespace
 
 namespace Kratos::Testing
@@ -1314,18 +1384,20 @@ KRATOS_TEST_CASE_IN_SUITE(LineInterfaceElement_InterpolatesNodalStresses, Kratos
     p_properties->SetValue(YOUNG_MODULUS, 1.000000e+07);
     p_properties->SetValue(POISSON_RATIO, 0.000000e+00);
     // create a triangle neighbour element
-    auto p_neighbour_element = ElementSetupUtilities::Create2D3NElement(nodes, p_properties, 2);
+    auto p_neighbour_element = make_intrusive<MockElementWithTotalStressVectors>(
+        2, std::make_shared<Triangle2D3<Node>>(nodes), p_properties);
+    p_neighbour_element->SetIntegrationMethod(GeometryData::IntegrationMethod::GI_GAUSS_2);
     ProcessInfo dummy_process_info;
     p_neighbour_element->Initialize(dummy_process_info);
-    std::vector<ConstitutiveLaw::StressVectorType> r_stress_vectors;
+    std::vector<ConstitutiveLaw::StressVectorType> total_stress_vectors;
     ConstitutiveLaw::StressVectorType              stress_vector(4);
     stress_vector <<= 3.0, 13.0 / 6.0, 4.0, 1.0;
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
     stress_vector <<= 3.0, 8.0 / 3.0, 4.0, 1.0;
-    r_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
 
-    p_neighbour_element->SetValuesOnIntegrationPoints(CAUCHY_STRESS_VECTOR, r_stress_vectors, dummy_process_info);
+    p_neighbour_element->SetValuesOnIntegrationPoints(TOTAL_STRESS_VECTOR, total_stress_vectors, dummy_process_info);
 
     nodes.clear();
     nodes.push_back(r_model_part.pGetNode(1));
@@ -1338,8 +1410,7 @@ KRATOS_TEST_CASE_IN_SUITE(LineInterfaceElement_InterpolatesNodalStresses, Kratos
     const auto     p_interface_properties =
         CreateElasticMaterialProperties<InterfacePlaneStrain>(normal_stiffness, shear_stiffness);
     auto interface_element = CreateInterfaceElementWithUDofs<Interface2D>(p_interface_properties, p_geometry);
-    GlobalPointersVector<Element> neighbours{p_neighbour_element};
-    interface_element.SetValue(NEIGHBOUR_ELEMENTS, neighbours);
+    interface_element.SetValue(NEIGHBOUR_ELEMENTS, MakeElementGlobalPtrContainerWith(p_neighbour_element));
 
     // Act
     interface_element.Initialize(dummy_process_info);
@@ -1359,20 +1430,20 @@ KRATOS_TEST_CASE_IN_SUITE(LineInterfaceElement_InterpolatesNodalStresses, Kratos
     nodes.push_back(r_model_part.pGetNode(4));
     nodes.push_back(r_model_part.pGetNode(5));
     nodes.push_back(r_model_part.CreateNewNode(6, 0.0, 2.0, 0.0));
-    auto p_other_neighbour_element = ElementSetupUtilities::Create2D3NElement(nodes, p_properties, 3);
+    auto p_other_neighbour_element = make_intrusive<MockElementWithTotalStressVectors>(
+        3, std::make_shared<Triangle2D3<Node>>(nodes), p_properties);
+    p_other_neighbour_element->SetIntegrationMethod(GeometryData::IntegrationMethod::GI_GAUSS_2);
     p_other_neighbour_element->Initialize(dummy_process_info);
-    r_stress_vectors.clear();
+    total_stress_vectors.clear();
     stress_vector <<= 7.0, 4.0, 11.0, 5.0 / 6.0;
-    r_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
     stress_vector <<= 7.0, 4.0, 11.0, 10.0 / 3.0;
-    r_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
     stress_vector <<= 7.0, 4.0, 11.0, 5.0 / 6.0;
-    r_stress_vectors.emplace_back(stress_vector);
-    p_other_neighbour_element->SetValuesOnIntegrationPoints(CAUCHY_STRESS_VECTOR, r_stress_vectors,
-                                                            dummy_process_info);
-
-    GlobalPointersVector<Element> other_neighbours{p_other_neighbour_element};
-    interface_element.SetValue(NEIGHBOUR_ELEMENTS, other_neighbours);
+    total_stress_vectors.emplace_back(stress_vector);
+    p_other_neighbour_element->SetValuesOnIntegrationPoints(
+        TOTAL_STRESS_VECTOR, total_stress_vectors, dummy_process_info);
+    interface_element.SetValue(NEIGHBOUR_ELEMENTS, MakeElementGlobalPtrContainerWith(p_other_neighbour_element));
 
     // Act
     interface_element.Initialize(dummy_process_info);
@@ -1404,27 +1475,29 @@ KRATOS_TEST_CASE_IN_SUITE(PlaneInterfaceElement_InterpolatesNodalStresses, Krato
     // properties for neighbour elements
     const auto p_properties = std::make_shared<Properties>();
     p_properties->SetValue(CONSTITUTIVE_LAW, std::make_shared<GeoIncrementalLinearElasticLaw>(
-                                                 std::make_unique<PlaneStrain>()));
+                                                 std::make_unique<ThreeDimensional>()));
     p_properties->SetValue(YOUNG_MODULUS, 1.000000e+07);
     p_properties->SetValue(POISSON_RATIO, 0.000000e+00);
     // create a hexagonal neighbour element
-    auto p_neighbour_element = ElementSetupUtilities::Create3D8NElement(nodes, p_properties, 2);
+    auto p_neighbour_element = make_intrusive<MockElementWithTotalStressVectors>(
+        2, std::make_shared<Hexahedra3D8<Node>>(nodes), p_properties);
+    p_neighbour_element->SetIntegrationMethod(GeometryData::IntegrationMethod::GI_GAUSS_2);
 
     ProcessInfo dummy_process_info;
     p_neighbour_element->Initialize(dummy_process_info);
-    std::vector<ConstitutiveLaw::StressVectorType> r_stress_vectors;
+    std::vector<ConstitutiveLaw::StressVectorType> total_stress_vectors;
     ConstitutiveLaw::StressVectorType              stress_vector(6);
     stress_vector <<= 3.0, (std::sqrt(3.0) - 1.0) / (2.0 * std::sqrt(3.0)), 4.0, 1.0, 2.0, 4.0;
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
     stress_vector <<= 3.0, (std::sqrt(3.0) + 1.0) / (2.0 * std::sqrt(3.0)), 4.0, 1.0, 2.0, 4.0;
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    p_neighbour_element->SetValuesOnIntegrationPoints(CAUCHY_STRESS_VECTOR, r_stress_vectors, dummy_process_info);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    p_neighbour_element->SetValuesOnIntegrationPoints(TOTAL_STRESS_VECTOR, total_stress_vectors, dummy_process_info);
 
     nodes.clear();
     nodes.push_back(r_model_part.pGetNode(2));
@@ -1441,8 +1514,7 @@ KRATOS_TEST_CASE_IN_SUITE(PlaneInterfaceElement_InterpolatesNodalStresses, Krato
     const auto     p_interface_properties =
         CreateElasticMaterialProperties<InterfaceThreeDimensionalSurface>(normal_stiffness, shear_stiffness);
     auto interface_element = CreateInterfaceElementWithUDofs<Interface3D>(p_interface_properties, p_geometry);
-    GlobalPointersVector<Element> neighbours{p_neighbour_element};
-    interface_element.SetValue(NEIGHBOUR_ELEMENTS, neighbours);
+    interface_element.SetValue(NEIGHBOUR_ELEMENTS, MakeElementGlobalPtrContainerWith(p_neighbour_element));
 
     // Act
     interface_element.Initialize(dummy_process_info);
@@ -1471,24 +1543,25 @@ KRATOS_TEST_CASE_IN_SUITE(PlaneInterfaceElement_InterpolatesNodalStresses, Krato
     nodes.push_back(r_model_part.CreateNewNode(16, 0.5, 1.5, 1.0));
     nodes.push_back(r_model_part.CreateNewNode(17, -0.5, 1.5, 1.0));
     nodes.push_back(r_model_part.pGetNode(18));
-    auto p_other_neighbour_element = ElementSetupUtilities::Create3D8NElement(nodes, p_properties, 3);
+    auto p_other_neighbour_element = make_intrusive<MockElementWithTotalStressVectors>(
+        3, std::make_shared<Hexahedra3D8<Node>>(nodes), p_properties);
+    p_neighbour_element->SetIntegrationMethod(GeometryData::IntegrationMethod::GI_GAUSS_2);
     p_other_neighbour_element->Initialize(dummy_process_info);
-    r_stress_vectors.clear();
+    total_stress_vectors.clear();
     stress_vector <<= 3.0, 4.0, 1.0, 2.0, (std::sqrt(3.0) - 1.0) / (2.0 * std::sqrt(3.0)), 4.0;
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
     stress_vector <<= 3.0, 4.0, 1.0, 2.0, (std::sqrt(3.0) + 1.0) / (2.0 * std::sqrt(3.0)), 4.0;
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    r_stress_vectors.emplace_back(stress_vector);
-    p_other_neighbour_element->SetValuesOnIntegrationPoints(CAUCHY_STRESS_VECTOR, r_stress_vectors,
-                                                            dummy_process_info);
-
-    GlobalPointersVector<Element> other_neighbours{p_other_neighbour_element};
-    interface_element.SetValue(NEIGHBOUR_ELEMENTS, other_neighbours);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    total_stress_vectors.emplace_back(stress_vector);
+    p_other_neighbour_element->SetValuesOnIntegrationPoints(
+        TOTAL_STRESS_VECTOR, total_stress_vectors, dummy_process_info);
+
+    interface_element.SetValue(NEIGHBOUR_ELEMENTS, MakeElementGlobalPtrContainerWith(p_other_neighbour_element));
 
     // Act
     interface_element.Initialize(dummy_process_info);
diff --git a/applications/GeoMechanicsApplication/tests/test_interface_prestress/MaterialParameters.json b/applications/GeoMechanicsApplication/tests/test_interface_prestress/MaterialParameters.json
@@ -8,7 +8,7 @@
             },
             "Variables": {
                 "IGNORE_UNDRAINED"         :  true,
-                "YOUNG_MODULUS"            :  120.0e6,
+                "YOUNG_MODULUS"            :  1.0e4,
                 "POISSON_RATIO"            :  0.0,
                 "DENSITY_SOLID"            :  0.0,
                 "DENSITY_WATER"            :  1.0e3,
@@ -29,8 +29,8 @@
                 "name"             :  "GeoIncrementalLinearElasticInterfaceLaw"
             },
             "Variables": {
-                "INTERFACE_NORMAL_STIFFNESS":  480.0e6,
-                "INTERFACE_SHEAR_STIFFNESS":   200.0e6
+                "INTERFACE_NORMAL_STIFFNESS":  2.0e3,
+                "INTERFACE_SHEAR_STIFFNESS":   1.0e3
             }
         }
     }]
diff --git a/applications/GeoMechanicsApplication/tests/test_interface_prestress/ProjectParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_interface_prestress/ProjectParameters_stage1.json
@@ -85,7 +85,7 @@
                     "variable_name":   "WATER_PRESSURE",
                     "active":          true,
                     "is_fixed":        true,
-                    "value":           0.0,
+                    "value":           -250.0,
                     "table":           0
                 }
             },
diff --git a/applications/GeoMechanicsApplication/tests/test_interface_prestress/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_interface_prestress/ProjectParameters_stage2.json
@@ -97,7 +97,7 @@
                     "variable_name": "WATER_PRESSURE",
                     "active": true,
                     "is_fixed": true,
-                    "value": 0.0,
+                    "value": -250.0,
                     "table": 0
                 }
             }
diff --git a/applications/GeoMechanicsApplication/tests/test_interface_prestress/README.md b/applications/GeoMechanicsApplication/tests/test_interface_prestress/README.md
@@ -46,13 +46,13 @@ o---o---o
 Next to the aforementioned elements and top load, the following constraints are applied:
 - The bottom nodes of the continuum elements are fixed in all directions.
 - The side nodes are fixed in the horizontal direction.
-- The water pressure is kept at zero in the entire domain.
+- The water pressure is kept at -250.0 Pa in the entire domain (to verify that the effect it has on the prestress calculation is properly taken into account).
 
 A linear elastic material is used for both the continuum and the interface elements.
 
 ## Assertions
 
 The following assertions are done to validate the prestress functionality:
 - In the second stage, all stage displacements of the nodes are checked to be zero (within a small tolerance). This would not be the case if the interfaces were not prestressed based on the first stage stresses in the continuum elements.
-- The stresses at the integration points of the interface elements are checked to be equal to the expected stresses at the integration points of the neighbouring continuum element. These are expected to be 1 kPa in the normal direction (following from the vertical equilibrium) and zero in the shear direction.
+- The traction components at the integration points of the interface elements are checked to be equal to the corresponding stress components at the integration points of the neighbouring continuum element. These are expected to be 1 kPa in the normal direction (following from the vertical equilibrium) and zero in the shear direction.
 - The relative (normal and shear) displacements of the interface elements are checked to be zero (within a small tolerance).
diff --git a/applications/OptimizationApplication/python_scripts/convergence_criteria/l2_conv_criterion.py b/applications/OptimizationApplication/python_scripts/convergence_criteria/l2_conv_criterion.py
@@ -1,6 +1,6 @@
 import typing
-import numpy
 import KratosMultiphysics as Kratos
+import KratosMultiphysics.OptimizationApplication as KratosOA
 from KratosMultiphysics.OptimizationApplication.utilities.optimization_problem import OptimizationProblem
 from KratosMultiphysics.OptimizationApplication.utilities.component_data_view import ComponentDataView
 from KratosMultiphysics.OptimizationApplication.utilities.optimization_problem_utilities import GetComponentHavingDataByFullName
@@ -57,7 +57,7 @@ def IsConverged(self) -> bool:
         if not hasattr(field, "Evaluate"):
             raise RuntimeError(f"The value represented by {self.__field_name} is not a field.")
 
-        self.__norm = numpy.linalg.norm(field.Evaluate().flatten())
+        self.__norm = KratosOA.ExpressionUtils.NormL2(field)
         self.__conv = self.__norm <= self.__tolerance
         self.__component_data_view.GetBufferedData().SetValue(self.__field_name.split(':')[0] + "_l2_norm", self.__norm)
         return self.__conv

Original file line number	Diff line number	Diff line change
`@@ -85,7 +85,7 @@`
`85`	`85`	`"variable_name": "WATER_PRESSURE",`
`86`	`86`	`"active": true,`
`87`	`87`	`"is_fixed": true,`
`88`		`- "value": 0.0,`
	`88`	`+ "value": -250.0,`
`89`	`89`	`"table": 0`
`90`	`90`	`}`
`91`	`91`	`},`
Original file line number	Diff line number	Diff line change
`@@ -97,7 +97,7 @@`
`97`	`97`	`"variable_name": "WATER_PRESSURE",`
`98`	`98`	`"active": true,`
`99`	`99`	`"is_fixed": true,`
`100`		`- "value": 0.0,`
	`100`	`+ "value": -250.0,`
`101`	`101`	`"table": 0`
`102`	`102`	`}`
`103`	`103`	`}`