[GeoMechanicsApplication] Create 3d line piping element 3d2n

applications/GeoMechanicsApplication/custom_elements/geo_steady_state_Pw_piping_element.cpp

@@ -16,5 +16,6 @@ namespace Kratos
 {
 
 template class GeoSteadyStatePwPipingElement<2, 2>;
+template class GeoSteadyStatePwPipingElement<3, 2>;
 
 } // Namespace Kratos

applications/GeoMechanicsApplication/custom_elements/geo_steady_state_Pw_piping_element.h

@@ -104,7 +104,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
         CheckHasDofsFor(WATER_PRESSURE);
         CheckProperties();
         // conditional on model dimension
-        CheckForNonZeroZCoordinateIn2D();
+        if constexpr (TDim == 2) {
+            CheckForNonZeroZCoordinate();
+        }
 
         KRATOS_CATCH("")
 
@@ -216,6 +218,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
 
     using Element::CalculateOnIntegrationPoints;
 
+    std::string Info() const override { return "GeoSteadyStatePwPipingElement"; }
+
 private:
     void CheckDomainSize() const
     {
@@ -249,6 +253,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
         CheckProperty(DENSITY_WATER);
         CheckProperty(DYNAMIC_VISCOSITY);
         CheckProperty(PIPE_HEIGHT);
+        if constexpr (TDim == 3) {
+            CheckProperty(PIPE_WIDTH_FACTOR);
+        }
     }
 
     void CheckProperty(const Kratos::Variable<double>& rVariable) const
@@ -260,15 +267,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
             << ") is not in the range [0,-> at element " << Id() << std::endl;
     }
 
-    void CheckForNonZeroZCoordinateIn2D() const
+    void CheckForNonZeroZCoordinate() const
     {
-        if constexpr (TDim == 2) {
-            const auto& r_geometry = GetGeometry();
-            auto        pos        = std::find_if(r_geometry.begin(), r_geometry.end(),
-                                                  [](const auto& node) { return node.Z() != 0.0; });
-            KRATOS_ERROR_IF_NOT(pos == r_geometry.end())
-                << "Node with non-zero Z coordinate found. Id: " << pos->Id() << std::endl;
-        }
+        const auto& r_geometry = GetGeometry();
+        auto        pos        = std::find_if(r_geometry.begin(), r_geometry.end(),
+                                              [](const auto& node) { return node.Z() != 0.0; });
+        KRATOS_ERROR_IF_NOT(pos == r_geometry.end())
+            << "Node with non-zero Z coordinate found. Id: " << pos->Id() << std::endl;
     }
 
     static double CalculateLength(const GeometryType& Geom)
@@ -317,9 +322,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoSteadyStatePwPipingElement : publ
         return result;
     }
 
-    static Matrix FillPermeabilityMatrix(double pipe_height)
+    Matrix FillPermeabilityMatrix(double pipe_height) const
     {
-        return ScalarMatrix{1, 1, std::pow(pipe_height, 3) / 12.0};
+        if constexpr (TDim == 2) {
+            return ScalarMatrix{1, 1, std::pow(pipe_height, 3) / 12.0};
+        }
+        return ScalarMatrix{1, 1, this->GetProperties()[PIPE_WIDTH_FACTOR] * std::pow(pipe_height, 4) / 12.0};
     }
 
     BoundedMatrix<double, TNumNodes, TNumNodes> CalculatePermeabilityMatrix(

applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp

@@ -81,6 +81,12 @@ int SteadyStatePwPipingElement<TDim, TNumNodes>::Check(const ProcessInfo& rCurre
     return ierr;
 }
 
+template <unsigned int TDim, unsigned int TNumNodes>
+std::string SteadyStatePwPipingElement<TDim, TNumNodes>::Info() const
+{
+    return "SteadyStatePwPipingElement";
+}
+
 template <unsigned int TDim, unsigned int TNumNodes>
 void SteadyStatePwPipingElement<TDim, TNumNodes>::Initialize(const ProcessInfo& rCurrentProcessInfo)
 {

applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp

@@ -95,6 +95,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwPipingElement
 
     int Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
+    std::string Info() const override;
+
     bool InEquilibrium(const PropertiesType& Prop, const GeometryType& Geom);
 
     double CalculateHeadGradient(const PropertiesType& Prop, const GeometryType& Geom, double pipe_length);

applications/GeoMechanicsApplication/custom_python/geo_mechanics_python_application.cpp

@@ -138,6 +138,7 @@ PYBIND11_MODULE(KratosGeoMechanicsApplication, m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PIPE_IN_EQUILIBRIUM)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PIPE_MODIFIED_D)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PIPE_MODEL_FACTOR)
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PIPE_WIDTH_FACTOR)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PIPE_HEIGHT)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, PIPE_ACTIVE)
 

applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp

@@ -67,6 +67,24 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
         mPipingIterations = rParameters["max_piping_iterations"].GetInt();
     }
 
+    template <typename PipingElementType>
+    std::optional<std::vector<PipingElementType*>> TryDownCastToPipingElement(const std::vector<Element*>& rPipeElements)
+    {
+        std::vector<PipingElementType*> result;
+        std::transform(rPipeElements.begin(), rPipeElements.end(), std::back_inserter(result),
+                       [](auto p_element) { return dynamic_cast<PipingElementType*>(p_element); });
+
+        const auto number_of_piping_elements = static_cast<std::size_t>(std::count_if(
+            result.begin(), result.end(), [](auto p_element) { return p_element != nullptr; }));
+        if (number_of_piping_elements == 0) return std::nullopt;
+
+        KRATOS_ERROR_IF(number_of_piping_elements != rPipeElements.size())
+            << "Unexpected number of piping elements of type " << rPipeElements[0]->Info() << ": expected "
+            << rPipeElements.size() << " but got " << number_of_piping_elements << std::endl;
+
+        return result;
+    }
+
     void FinalizeSolutionStep() override
     {
         KRATOS_INFO_IF("PipingLoop", this->GetEchoLevel() > 0 && rank == 0)
@@ -80,42 +98,31 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
             return;
         }
 
-        auto piping_interface_elements = std::vector<SteadyStatePwPipingElement<2, 4>*>{};
-        std::transform(piping_elements.begin(), piping_elements.end(),
-                       std::back_inserter(piping_interface_elements), [](auto p_element) {
-            return dynamic_cast<SteadyStatePwPipingElement<2, 4>*>(p_element);
-        });
-
-        const auto number_of_piping_interface_elements = static_cast<std::size_t>(
-            std::count_if(piping_interface_elements.begin(), piping_interface_elements.end(),
-                          [](auto p_element) { return p_element != nullptr; }));
-        KRATOS_ERROR_IF(number_of_piping_interface_elements != 0 &&
-                        number_of_piping_interface_elements != piping_interface_elements.size())
-            << "Unexpected number of piping interface elements: expected either 0 or "
-            << piping_interface_elements.size() << " but got "
-            << number_of_piping_interface_elements << std::endl;
-
-        if (number_of_piping_interface_elements > 0) {
-            this->DetermineOpenPipingElements(piping_interface_elements);
+        auto piping_interface_elements =
+            TryDownCastToPipingElement<SteadyStatePwPipingElement<2, 4>>(piping_elements);
+        if (piping_interface_elements) {
+            this->DetermineOpenPipingElements(piping_interface_elements.value());
             this->BaseClassFinalizeSolutionStep();
             return;
         }
 
-        auto piping_line_elements = std::vector<GeoSteadyStatePwPipingElement<2, 2>*>{};
-        std::transform(piping_elements.begin(), piping_elements.end(),
-                       std::back_inserter(piping_line_elements), [](auto p_element) {
-            return dynamic_cast<GeoSteadyStatePwPipingElement<2, 2>*>(p_element);
-        });
+        auto piping_2D_line_elements =
+            TryDownCastToPipingElement<GeoSteadyStatePwPipingElement<2, 2>>(piping_elements);
+        if (piping_2D_line_elements) {
+            this->DetermineOpenPipingElements(piping_2D_line_elements.value());
+            this->BaseClassFinalizeSolutionStep();
+            return;
+        }
 
-        const auto number_of_piping_line_elements = static_cast<std::size_t>(
-            std::count_if(piping_line_elements.begin(), piping_line_elements.end(),
-                          [](auto p_element) { return p_element != nullptr; }));
-        KRATOS_ERROR_IF(number_of_piping_line_elements != piping_line_elements.size())
-            << "Unexpected number of piping line elements: expected " << piping_line_elements.size()
-            << " but got " << number_of_piping_line_elements << std::endl;
+        auto piping_3D_line_elements =
+            TryDownCastToPipingElement<GeoSteadyStatePwPipingElement<3, 2>>(piping_elements);
+        if (piping_3D_line_elements) {
+            this->DetermineOpenPipingElements(piping_3D_line_elements.value());
+            this->BaseClassFinalizeSolutionStep();
+            return;
+        }
 
-        this->DetermineOpenPipingElements(piping_line_elements);
-        this->BaseClassFinalizeSolutionStep();
+        KRATOS_ERROR << "Unexpected type of piping element\n";
     }
 
     /**
@@ -323,6 +330,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
             }
             ++piping_iteration;
         }
+
         return equilibrium;
     }
 

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -85,6 +85,7 @@ void KratosGeoMechanicsApplication::Register()
     KRATOS_REGISTER_ELEMENT("SteadyStatePwPipingElement3D6N", mSteadyStatePwPipingElement3D6N)
     KRATOS_REGISTER_ELEMENT("SteadyStatePwPipingElement3D8N", mSteadyStatePwPipingElement3D8N)
     KRATOS_REGISTER_ELEMENT("GeoSteadyStatePwPipingElement2D2N", mGeoSteadyStatePwPipingElement2D2N)
+    KRATOS_REGISTER_ELEMENT("GeoSteadyStatePwPipingElement3D2N", mGeoSteadyStatePwPipingElement3D2N)
 
     // Small strain elements
     KRATOS_REGISTER_ELEMENT("UPwSmallStrainElement2D3N", mUPwSmallStrainElement2D3N)
@@ -522,6 +523,7 @@ void KratosGeoMechanicsApplication::Register()
     KRATOS_REGISTER_VARIABLE(PIPE_IN_EQUILIBRIUM)
     KRATOS_REGISTER_VARIABLE(PIPE_MODIFIED_D)
     KRATOS_REGISTER_VARIABLE(PIPE_MODEL_FACTOR)
+    KRATOS_REGISTER_VARIABLE(PIPE_WIDTH_FACTOR)
     KRATOS_REGISTER_VARIABLE(PIPE_HEIGHT)
     KRATOS_REGISTER_VARIABLE(PREV_PIPE_HEIGHT)
     KRATOS_REGISTER_VARIABLE(DIFF_PIPE_HEIGHT)

applications/GeoMechanicsApplication/geo_mechanics_application.h

@@ -403,6 +403,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoMechanicsApplication : publ
         std::make_unique<ThreeDimensionalStressState>()};
     const GeoSteadyStatePwPipingElement<2, 2> mGeoSteadyStatePwPipingElement2D2N{
         0, Kratos::make_shared<Line2D2<NodeType>>(Element::GeometryType::PointsArrayType(2))};
+    const GeoSteadyStatePwPipingElement<3, 2> mGeoSteadyStatePwPipingElement3D2N{
+        0, Kratos::make_shared<Line3D2<NodeType>>(Element::GeometryType::PointsArrayType(2))};
 
     // small strain elements:
     const UPwSmallStrainElement<2, 3> mUPwSmallStrainElement2D3N{

applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp

@@ -165,6 +165,7 @@ KRATOS_CREATE_VARIABLE(int, PIPE_START_ELEMENT)
 KRATOS_CREATE_VARIABLE(double, PIPE_ELEMENT_LENGTH)
 KRATOS_CREATE_VARIABLE(bool, PIPE_MODIFIED_D)
 KRATOS_CREATE_VARIABLE(double, PIPE_MODEL_FACTOR)
+KRATOS_CREATE_VARIABLE(double, PIPE_WIDTH_FACTOR)
 KRATOS_CREATE_VARIABLE(double, PIPE_HEIGHT)
 KRATOS_CREATE_VARIABLE(double, PREV_PIPE_HEIGHT)
 KRATOS_CREATE_VARIABLE(double, DIFF_PIPE_HEIGHT)

applications/GeoMechanicsApplication/geo_mechanics_application_variables.h

@@ -183,6 +183,7 @@ KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, PIPE_ELEME
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, bool, PIPE_IN_EQUILIBRIUM)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, bool, PIPE_MODIFIED_D)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, PIPE_MODEL_FACTOR)
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, PIPE_WIDTH_FACTOR)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, PIPE_HEIGHT)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, PREV_PIPE_HEIGHT)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, DIFF_PIPE_HEIGHT)

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

@@ -69,6 +69,18 @@ intrusive_ptr<GeoSteadyStatePwPipingElement<2, 2>> CreateGeoSteadyStatePwPipingE
     return p_result;
 }
 
+intrusive_ptr<GeoSteadyStatePwPipingElement<3, 2>> CreateGeoSteadyStatePwPipingElement3D2NWithPWDofs(
+    const ModelPart& rModelPart, const Properties::Pointer& rProperties, const Geometry<Node>::Pointer& rGeometry)
+{
+    auto p_result = make_intrusive<GeoSteadyStatePwPipingElement<3, 2>>(
+        NextElementNumber(rModelPart), rGeometry, rProperties);
+    for (auto& node : p_result->GetGeometry()) {
+        node.AddDof(WATER_PRESSURE);
+    }
+
+    return p_result;
+}
+
 intrusive_ptr<GeoSteadyStatePwPipingElement<2, 2>> CreateHorizontalUnitLengthGeoSteadyStatePwPipingElementWithPWDofs(
     ModelPart& rModelPart, const Properties::Pointer& rProperties)
 {
@@ -79,6 +91,16 @@ intrusive_ptr<GeoSteadyStatePwPipingElement<2, 2>> CreateHorizontalUnitLengthGeo
     return p_element;
 }
 
+intrusive_ptr<GeoSteadyStatePwPipingElement<3, 2>> CreateHorizontalUnitLengthGeoSteadyStatePwPipingElement3D2NWithPWDofs(
+    ModelPart& rModelPart, const Properties::Pointer& rProperties)
+{
+    const auto p_geometry = std::make_shared<Line3D2<Node>>(CreateNodesOnModelPart(rModelPart));
+    auto p_element = CreateGeoSteadyStatePwPipingElement3D2NWithPWDofs(rModelPart, rProperties, p_geometry);
+
+    rModelPart.AddElement(p_element);
+    return p_element;
+}
+
 intrusive_ptr<GeoSteadyStatePwPipingElement<2, 2>> CreateHorizontalUnitLengthGeoSteadyStatePwPipingElementWithoutPWDofs(
     ModelPart& rModelPart, const Properties::Pointer& rProperties)
 {
@@ -304,6 +326,45 @@ KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementReturnsTheExpectedLeftHan
     KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(actual_right_hand_side, expected_right_hand_side, Defaults::relative_tolerance)
 }
 
+KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElement3D2NReturnsTheExpectedLeftHandSideAndRightHandSide,
+                          KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    // Arrange
+    const auto dummy_process_info = ProcessInfo{};
+    const auto p_properties       = std::make_shared<Properties>();
+
+    Model model;
+    auto& r_model_part = CreateModelPartWithWaterPressureVariableAndVolumeAcceleration(model);
+    auto  p_element =
+        CreateHorizontalUnitLengthGeoSteadyStatePwPipingElement3D2NWithPWDofs(r_model_part, p_properties);
+    p_element->GetProperties().SetValue(DENSITY_WATER, 1.0E3);
+    p_element->GetProperties().SetValue(DYNAMIC_VISCOSITY, 1.0E-2);
+    p_element->GetProperties().SetValue(PIPE_WIDTH_FACTOR, 1.5);
+    p_element->SetValue(PIPE_HEIGHT, 1.0E-1);
+    // Set gravity perpendicular to the line ( so no fluid body flow vector from this )
+    const auto gravity_acceleration = array_1d<double, 3>{0.0, -10.0, 0.0};
+    p_element->GetGeometry()[0].FastGetSolutionStepValue(VOLUME_ACCELERATION) = gravity_acceleration;
+    p_element->GetGeometry()[1].FastGetSolutionStepValue(VOLUME_ACCELERATION) = gravity_acceleration;
+    // Create a head gradient of -10.
+    p_element->GetGeometry()[0].FastGetSolutionStepValue(WATER_PRESSURE) = 10.0;
+    p_element->GetGeometry()[1].FastGetSolutionStepValue(WATER_PRESSURE) = 0.0;
+
+    // Act
+    Vector actual_right_hand_side;
+    Matrix actual_left_hand_side;
+    p_element->CalculateLocalSystem(actual_left_hand_side, actual_right_hand_side, dummy_process_info);
+
+    // Assert
+    auto expected_left_hand_side  = Matrix{ScalarMatrix{2, 2, 0.015 / 12.}};
+    expected_left_hand_side(0, 0) = -expected_left_hand_side(0, 0);
+    expected_left_hand_side(1, 1) = -expected_left_hand_side(1, 1);
+    KRATOS_EXPECT_MATRIX_RELATIVE_NEAR(actual_left_hand_side, expected_left_hand_side, Defaults::relative_tolerance)
+
+    auto expected_right_hand_side = Vector{ScalarVector{2, 0.15 / 12.}};
+    expected_right_hand_side(1)   = -expected_right_hand_side(1);
+    KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(actual_right_hand_side, expected_right_hand_side, Defaults::relative_tolerance)
+}
+
 KRATOS_TEST_CASE_IN_SUITE(GeoSteadyStatePwPipingElementHasValuesAfterInitialize, KratosGeoMechanicsFastSuiteWithoutKernel)
 {
     // Arrange

applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py

@@ -55,6 +55,7 @@
 from test_integration_node_extrapolation import KratosGeoMechanicsExtrapolationTests
 from test_truss_backbone_mat import KratosGeoMechanicsTrussBackboneMaterialTests
 from test_line_interface_elements import KratosGeoMechanicsInterfaceElementTests
+from test_three_dimensional_piping_validation import KratosGeoMechanicsThreeDimensionalPipingValidation
 
 def AssembleTestSuites():
     ''' Populates the test suites to run.
@@ -134,7 +135,8 @@ def AssembleTestSuites():
                         KratosGeoMechanicsBenchmarkSet2,
                         KratosGeoMechanicsTransientGroundWaterFlowTests,
                         TestSellmeijersRuleValidation,
-                        KratosGeoMechanicsDynamicsLongTests
+                        KratosGeoMechanicsDynamicsLongTests,
+                        KratosGeoMechanicsThreeDimensionalPipingValidation
                         ]
 
     # Create an array that contains all the tests from every testCase

applications/GeoMechanicsApplication/tests/test_three_dimensional_piping_validation.py

@@ -0,0 +1,32 @@
+import os
+
+import KratosMultiphysics.KratosUnittest as KratosUnittest
+import test_helper
+
+
+class KratosGeoMechanicsThreeDimensionalPipingValidation(KratosUnittest.TestCase):
+    def test_three_dimensional_piping(self):
+        file_path = test_helper.get_file_path(os.path.join('three_dimensional_piping'))
+        simulation = test_helper.run_kratos(file_path)
+        result_file_name = os.path.join(file_path, 'Cube_10sS_1.post.res')
+
+        reader = test_helper.GiDOutputFileReader()
+        actual_data = reader.read_output_from(result_file_name)
+
+        pipe_elements = [24001, 24002, 24003, 24004, 24005, 24006, 24007, 24008, 24009, 24010, 24011, 24012, 24013,
+                         24014, 24015, 24016, 24017, 24018, 24019]
+        pipe_active_before_critical_head = \
+            reader.element_integration_point_values_at_time("PIPE_ACTIVE", 0.5, actual_data, pipe_elements, [1])
+        for active, element in zip(pipe_active_before_critical_head, pipe_elements):
+            if element > 24004:
+                self.assertEqual(active[0], 1)
+            else:
+                self.assertEqual(active[0], 0)
+
+        pipe_active_after_critical_head = \
+            reader.element_integration_point_values_at_time("PIPE_ACTIVE", 1.0, actual_data, pipe_elements, [1])
+        for active in pipe_active_after_critical_head:
+            self.assertEqual(active[0], 1)
+
+if __name__ == '__main__':
+    KratosUnittest.main()