[GeoMechanicsApplication] Fix issues related to system integration fo…

…r interface element (#12691)

* Added a check for element being active + unit test
* Added test for multiple initializations of an element (emulating multi-stage)

applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp

@@ -193,6 +193,7 @@ void LineInterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
         GeoElementUtilities::EvaluateShapeFunctionsAtIntegrationPoints(
             mIntegrationScheme->GetIntegrationPoints(), GetGeometry());
 
+    mConstitutiveLaws.clear();
     for (const auto& r_shape_function_values : shape_function_values_at_integration_points) {
         mConstitutiveLaws.push_back(GetProperties()[CONSTITUTIVE_LAW]->Clone());
         mConstitutiveLaws.back()->InitializeMaterial(GetProperties(), GetGeometry(), r_shape_function_values);
@@ -204,13 +205,15 @@ int LineInterfaceElement::Check(const ProcessInfo& rCurrentProcessInfo) const
     int error = Element::Check(rCurrentProcessInfo);
     if (error != 0) return error;
 
-    KRATOS_ERROR_IF(mIntegrationScheme->GetNumberOfIntegrationPoints() != mConstitutiveLaws.size())
-        << "Number of integration points (" << mIntegrationScheme->GetNumberOfIntegrationPoints()
-        << ") and constitutive laws (" <<  mConstitutiveLaws.size() << ") do not match.\n";
+    if (this->IsActive()) {
+        KRATOS_ERROR_IF(mIntegrationScheme->GetNumberOfIntegrationPoints() != mConstitutiveLaws.size())
+            << "Number of integration points (" << mIntegrationScheme->GetNumberOfIntegrationPoints()
+            << ") and constitutive laws (" << mConstitutiveLaws.size() << ") do not match.\n";
 
-    for (const auto& r_constitutive_law : mConstitutiveLaws) {
-        error = r_constitutive_law->Check(GetProperties(), GetGeometry(), rCurrentProcessInfo);
-        if (error != 0) return error;
+        for (const auto& r_constitutive_law : mConstitutiveLaws) {
+            error = r_constitutive_law->Check(GetProperties(), GetGeometry(), rCurrentProcessInfo);
+            if (error != 0) return error;
+        }
     }
 
     return 0;

applications/GeoMechanicsApplication/tests/cpp_tests/custom_elements/test_line_interface_element.cpp

@@ -378,6 +378,29 @@ KRATOS_TEST_CASE_IN_SUITE(GetInitializedConstitutiveLawsAfterElementInitializati
     }
 }
 
+KRATOS_TEST_CASE_IN_SUITE(InterfaceElement_HasCorrectNumberOfConstitutiveLawsAfterMultipleInitializations,
+                          KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    // Arrange
+    const auto p_properties = std::make_shared<Properties>();
+    p_properties->GetValue(CONSTITUTIVE_LAW) = std::make_shared<GeoIncrementalLinearElasticInterfaceLaw>();
+
+    Model model;
+    auto element = CreateHorizontalUnitLength2Plus2NodedLineInterfaceElementWithUDofs(model, p_properties);
+
+    const auto dummy_process_info = ProcessInfo{};
+
+    // Multiple initializations emulate a multi-stage simulation
+    element.Initialize(dummy_process_info);
+    element.Initialize(dummy_process_info);
+    element.Initialize(dummy_process_info);
+
+    // Assert
+    auto constitutive_laws = std::vector<ConstitutiveLaw::Pointer>{};
+    element.CalculateOnIntegrationPoints(CONSTITUTIVE_LAW, constitutive_laws, dummy_process_info);
+    KRATOS_EXPECT_EQ(constitutive_laws.size(), 2);
+}
+
 KRATOS_TEST_CASE_IN_SUITE(LineInterfaceElement_CalculateLocalSystem_ReturnsExpectedLeftAndRightHandSide,
                           KratosGeoMechanicsFastSuiteWithoutKernel)
 {
@@ -539,7 +562,6 @@ KRATOS_TEST_CASE_IN_SUITE(3Plus3NodedLineInterfaceElement_CalculateLocalSystem_R
 
 KRATOS_TEST_CASE_IN_SUITE(InterfaceElement_CheckThrowsWhenElementIsNotInitialized, KratosGeoMechanicsFastSuiteWithoutKernel)
 {
-    // Arrange
     constexpr auto normal_stiffness = 20.0;
     constexpr auto shear_stiffness  = 10.0;
     const auto p_properties = CreateLinearElasticMaterialProperties(normal_stiffness, shear_stiffness);
@@ -557,4 +579,23 @@ KRATOS_TEST_CASE_IN_SUITE(InterfaceElement_CheckThrowsWhenElementIsNotInitialize
     KRATOS_EXPECT_EQ(element.Check(dummy_process_info), 0);
 }
 
+KRATOS_TEST_CASE_IN_SUITE(InterfaceElement_CheckDoesNotThrowWhenElementIsNotActive, KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    constexpr auto normal_stiffness = 20.0;
+    constexpr auto shear_stiffness  = 10.0;
+    const auto p_properties = CreateLinearElasticMaterialProperties(normal_stiffness, shear_stiffness);
+
+    Model model;
+    auto element = CreateHorizontalUnitLength3Plus3NodedLineInterfaceElementWithDisplacementDoF(model, p_properties);
+
+    // In the integrated workflow, the elements are not initialized, when they are not active.
+    // However, the Check method is always called on all elements, even if they are not active.
+    // Therefore, the Check method should not throw an exception in this case, even though the
+    // constitutive laws are not initialized.
+    element.Set(ACTIVE, false);
+
+    const auto dummy_process_info = ProcessInfo{};
+    KRATOS_EXPECT_EQ(element.Check(dummy_process_info), 0);
+}
+
 } // namespace Kratos::Testing