Clean up piping strategy and piping interface element (#12785)

Cleaned up the piping strategy as well as the piping interface element. This is to prepare for the transition from piping interface elements to piping line elements. For some time (while the piping line element is developed further and validated), these two element types will co-exist.

Author: avdg81

applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp

@@ -13,7 +13,9 @@
 // Application includes
 #include "custom_elements/steady_state_Pw_piping_element.hpp"
 #include "custom_utilities/element_utilities.hpp"
+#include "custom_utilities/transport_equation_utilities.hpp"
 #include "utilities/math_utils.h"
+
 #include <cmath>
 
 namespace Kratos
@@ -272,23 +274,6 @@ double SteadyStatePwPipingElement<3, 8>::CalculateHeadGradient(const PropertiesT
                  << std::endl;
 }
 
-/// <summary>
-///  Calculate the particle diameter for the particles in the pipe. The particle diameter equals d70,
-/// when the unmodified sellmeijer piping rule is used.
-/// </summary>
-/// <param name="Prop"></param>
-/// <param name="Geom"></param>
-/// <returns></returns>
-template <unsigned int TDim, unsigned int TNumNodes>
-double SteadyStatePwPipingElement<TDim, TNumNodes>::CalculateParticleDiameter(const PropertiesType& Prop)
-{
-    double diameter;
-
-    if (Prop[PIPE_MODIFIED_D]) diameter = 2.08e-4 * pow((Prop[PIPE_D_70] / 2.08e-4), 0.4);
-    else diameter = Prop[PIPE_D_70];
-    return diameter;
-}
-
 /// <summary>
 /// Calculates the equilibrium pipe height of a piping element according to Sellmeijers rule
 /// </summary>
@@ -310,7 +295,7 @@ double SteadyStatePwPipingElement<TDim, TNumNodes>::CalculateEquilibriumPipeHeig
     double dhdx = CalculateHeadGradient(Prop, Geom, pipe_length);
 
     // calculate particle diameter
-    double particle_d = CalculateParticleDiameter(Prop);
+    double particle_d = GeoTransportEquationUtilities::CalculateParticleDiameter(Prop);
 
     // todo calculate slope of pipe, currently pipe is assumed to be horizontal
     const double pipeSlope = 0;

applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp

@@ -119,8 +119,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwPipingElement
                                       std::vector<double>&    rValues,
                                       const ProcessInfo&      rCurrentProcessInfo) override;
 
-    double CalculateParticleDiameter(const PropertiesType& Prop);
-
     double pipe_initialised = false;
 
 private:

applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp

@@ -22,6 +22,7 @@
 #include "boost/range/adaptor/filtered.hpp"
 #include "custom_elements/steady_state_Pw_piping_element.hpp"
 #include "custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp"
+#include "custom_utilities/transport_equation_utilities.hpp"
 #include "geo_mechanics_application_variables.h"
 
 #include <chrono>
@@ -45,11 +46,9 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
     using DofsArrayType            = typename BaseType::DofsArrayType;
     using TSystemMatrixType        = typename BaseType::TSystemMatrixType;
     using TSystemVectorType        = typename BaseType::TSystemVectorType;
-    using filtered_elements =
-        typename boost::range_detail::filtered_range<std::function<bool(Element*)>, std::vector<Element*>>;
-    using PropertiesType = Properties;
-    using NodeType       = Node;
-    using GeometryType   = Geometry<NodeType>;
+    using PropertiesType           = Properties;
+    using NodeType                 = Node;
+    using GeometryType             = Geometry<NodeType>;
 
     GeoMechanicsNewtonRaphsonErosionProcessStrategy(ModelPart&                    model_part,
                                                     typename TSchemeType::Pointer pScheme,
@@ -76,7 +75,17 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
 
         // get piping elements
         std::vector<Element*> PipeElements = GetPipingElements();
-        unsigned int          n_el         = PipeElements.size(); // number of piping elements
+
+        auto piping_interface_elements = std::vector<SteadyStatePwPipingElement<2, 4>*>{};
+        std::transform(
+            PipeElements.begin(), PipeElements.end(), std::back_inserter(piping_interface_elements),
+            [](auto p_element) { return dynamic_cast<SteadyStatePwPipingElement<2, 4>*>(p_element); });
+        KRATOS_DEBUG_ERROR_IF(std::any_of(piping_interface_elements.begin(),
+                                          piping_interface_elements.end(), [](auto p_element) {
+            return p_element == nullptr;
+        })) << "Not all open piping elements could be downcast to SteadyStatePwPipingElement<2, 4>*\n";
+
+        unsigned int n_el = PipeElements.size(); // number of piping elements
 
         // get initially open pipe elements
         unsigned int openPipeElements = this->InitialiseNumActivePipeElements(PipeElements);
@@ -104,7 +113,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
             std::function<bool(Element*)> filter = [](Element* i) {
                 return i->Has(PIPE_ACTIVE) && i->GetValue(PIPE_ACTIVE);
             };
-            filtered_elements OpenPipeElements = PipeElements | boost::adaptors::filtered(filter);
+            auto OpenPipeElements = piping_interface_elements | boost::adaptors::filtered(filter);
             KRATOS_INFO_IF("PipingLoop", this->GetEchoLevel() > 0 && rank == 0)
                 << "Number of Open Pipe Elements: " << boost::size(OpenPipeElements) << std::endl;
 
@@ -242,21 +251,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
         return nOpenElements;
     }
 
-    /// <summary>
-    /// Calculates the pipe particle diameter according to the modified Sellmeijer rule if selected.
-    /// Else the pipe particle diameter is equal to the d70.
-    /// </summary>
-    /// <param name="Prop"></param>
-    /// <returns></returns>
-    double CalculateParticleDiameter(const PropertiesType& Prop)
-    {
-        double diameter;
-
-        if (Prop[PIPE_MODIFIED_D]) diameter = 2.08e-4 * pow((Prop[PIPE_D_70] / 2.08e-4), 0.4);
-        else diameter = Prop[PIPE_D_70];
-        return diameter;
-    }
-
     /// <summary>
     /// Calculates the maximum pipe height which the algorithm will allow
     /// </summary>
@@ -270,8 +264,8 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
         // loop over all elements
         for (Element* pipe_element : pipe_elements) {
             // calculate pipe particle diameter of pipe element
-            PropertiesType prop              = pipe_element->GetProperties();
-            double         particle_diameter = this->CalculateParticleDiameter(prop);
+            PropertiesType prop = pipe_element->GetProperties();
+            double particle_diameter = GeoTransportEquationUtilities::CalculateParticleDiameter(prop);
 
             // get maximum pipe particle diameter of all pipe elements
             if (particle_diameter > max_diameter) {
@@ -310,7 +304,8 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
         GeoMechanicsNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>::FinalizeSolutionStep();
     }
 
-    bool check_pipe_equilibrium(filtered_elements open_pipe_elements, double amax, unsigned int mPipingIterations)
+    template <typename FilteredElementType>
+    bool check_pipe_equilibrium(const FilteredElementType& open_pipe_elements, double amax, unsigned int mPipingIterations)
     {
         bool         equilibrium = false;
         bool         converged   = true;
@@ -338,13 +333,12 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
                 // Update depth of open piping Elements
                 equilibrium = true;
                 for (auto OpenPipeElement : open_pipe_elements) {
-                    auto pElement = static_cast<SteadyStatePwPipingElement<2, 4>*>(OpenPipeElement);
                     // get open pipe element geometry and properties
                     auto& Geom = OpenPipeElement->GetGeometry();
                     auto& prop = OpenPipeElement->GetProperties();
 
                     // calculate equilibrium pipe height and get current pipe height
-                    double eq_height = pElement->CalculateEquilibriumPipeHeight(
+                    double eq_height = OpenPipeElement->CalculateEquilibriumPipeHeight(
                         prop, Geom, OpenPipeElement->GetValue(PIPE_ELEMENT_LENGTH));
                     double current_height = OpenPipeElement->GetValue(PIPE_HEIGHT);
 
@@ -420,13 +414,14 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy
     /// <param name="open_pipe_elements"> open pipe elements</param>
     /// <param name="grow"> boolean to check if pipe grows</param>
     /// <returns></returns>
-    void save_or_reset_pipe_heights(filtered_elements open_pipe_elements, bool grow)
+    template <typename FilteredElementType>
+    void save_or_reset_pipe_heights(const FilteredElementType& open_pipe_elements, bool grow)
     {
-        for (Element* OpenPipeElement : open_pipe_elements) {
+        for (auto p_element : open_pipe_elements) {
             if (grow) {
-                OpenPipeElement->SetValue(PREV_PIPE_HEIGHT, OpenPipeElement->GetValue(PIPE_HEIGHT));
+                p_element->SetValue(PREV_PIPE_HEIGHT, p_element->GetValue(PIPE_HEIGHT));
             } else {
-                OpenPipeElement->SetValue(PIPE_HEIGHT, OpenPipeElement->GetValue(PREV_PIPE_HEIGHT));
+                p_element->SetValue(PIPE_HEIGHT, p_element->GetValue(PREV_PIPE_HEIGHT));
             }
         }
     }

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.cpp

@@ -0,0 +1,25 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Anne van de Graaf
+//
+
+#include "transport_equation_utilities.hpp"
+
+namespace Kratos
+{
+
+double GeoTransportEquationUtilities::CalculateParticleDiameter(const Properties& rProperties)
+{
+    return rProperties.Has(PIPE_MODIFIED_D) && rProperties[PIPE_MODIFIED_D]
+               ? 2.08e-4 * std::pow((rProperties[PIPE_D_70] / 2.08e-4), 0.4)
+               : rProperties[PIPE_D_70];
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp

@@ -18,12 +18,13 @@
 #include "custom_retention/retention_law.h"
 #include "custom_utilities/stress_strain_utilities.h"
 #include "geo_mechanics_application_variables.h"
+#include "includes/kratos_export_api.h"
 #include "includes/variables.h"
 
 namespace Kratos
 {
 
-class GeoTransportEquationUtilities
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoTransportEquationUtilities
 {
 public:
     template <unsigned int TDim, unsigned int TNumNodes>
@@ -156,6 +157,8 @@ class GeoTransportEquationUtilities
         return result;
     }
 
+    [[nodiscard]] static double CalculateParticleDiameter(const Properties& rProperties);
+
 private:
     [[nodiscard]] static double CalculateBiotCoefficient(const Matrix&     rConstitutiveMatrix,
                                                          const Properties& rProperties)