[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("")
 
@@ -249,6 +251,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 +265,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 +320,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_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: 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 line element\n";
     }
 
     /**

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,46 @@ 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 )
+    p_element->GetGeometry()[0].FastGetSolutionStepValue(VOLUME_ACCELERATION) =
+        array_1d<double, 3>{0.0, -10.0, 0.0};
+    p_element->GetGeometry()[1].FastGetSolutionStepValue(VOLUME_ACCELERATION) =
+        array_1d<double, 3>{0.0, -10.0, 0.0};
+    // 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