[GeoMechanicsApplication] Add 1D functionality to existing 2D/3D Pw element

applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.cpp

@@ -0,0 +1,27 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   John van Esch
+//                   Gennady Markelov
+//
+
+#include "custom_elements/transient_Pw_line_element.h"
+
+namespace Kratos
+{
+
+template class TransientPwLineElement<2, 2>;
+template class TransientPwLineElement<2, 3>;
+template class TransientPwLineElement<2, 4>;
+template class TransientPwLineElement<2, 5>;
+template class TransientPwLineElement<3, 2>;
+template class TransientPwLineElement<3, 3>;
+
+} // Namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp

@@ -44,7 +44,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     template< unsigned int TDim, unsigned int TNumNodes >
     static inline void InterpolateVariableWithComponents(array_1d<double, TDim>& rVector,
                                                          const Matrix& NContainer,
@@ -61,7 +60,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     template< unsigned int TDof, unsigned int TNumNodes >
     static inline void InterpolateVariableWithComponents(Vector& rVector,
                                                          const Matrix& NContainer,
@@ -81,7 +79,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void FillArray1dOutput(array_1d<double,3>& rOutputValue,
                                          const array_1d<double,2>& ComputedValue)
     {
@@ -90,7 +87,7 @@ class GeoElementUtilities
         rOutputValue[2] = 0.0;
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
     static inline void FillArray1dOutput(array_1d<double,3>& rOutputValue,
                                          const array_1d<double,3>& ComputedValue)
     {
@@ -99,7 +96,6 @@ class GeoElementUtilities
         rOutputValue[2] = ComputedValue[2];
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     template< unsigned int TDim, unsigned int TNumNodes >
     static inline void GetNodalVariableVector(array_1d<double,TDim*TNumNodes>& rNodalVariableVector,
                                               const Element::GeometryType& Geom,
@@ -116,7 +112,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------
     template< unsigned int TDim, unsigned int TNumNodes >
     static inline void GetNodalVariableMatrix(Matrix& rNodalVariableMatrix,
                                               const Element::GeometryType& rGeom,
@@ -133,7 +128,14 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    static inline void FillPermeabilityMatrix(BoundedMatrix<double, 1, 1>&   rPermeabilityMatrix,
+                                              const Element::PropertiesType& Prop)
+    {
+        // 1D
+        rPermeabilityMatrix(0, 0) = Prop[PERMEABILITY_XX];
+    }
+
+
     static inline void FillPermeabilityMatrix(BoundedMatrix<double,2,2>& rPermeabilityMatrix,
                                               const Element::PropertiesType& Prop)
     {
@@ -145,7 +147,6 @@ class GeoElementUtilities
         rPermeabilityMatrix(1,0) = rPermeabilityMatrix(0,1);
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void FillPermeabilityMatrix(BoundedMatrix<double,3,3>& rPermeabilityMatrix,
                                               const Element::PropertiesType& Prop)
     {
@@ -164,8 +165,6 @@ class GeoElementUtilities
         rPermeabilityMatrix(0,2) = rPermeabilityMatrix(2,0);
     }
 
-
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void InvertMatrix2(BoundedMatrix<double,2,2>& rInvertedMatrix,
                                      const BoundedMatrix<double,2,2>& InputMatrix,
                                      double &InputMatrixDet)
@@ -188,7 +187,6 @@ class GeoElementUtilities
         KRATOS_CATCH("")
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void InvertMatrix2(BoundedMatrix<double,2,2>& rInvertedMatrix,
                                      const BoundedMatrix<double,2,2>& InputMatrix)
     {
@@ -201,7 +199,6 @@ class GeoElementUtilities
         KRATOS_CATCH("")
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     template< unsigned int TDim>
     static inline void AssembleDensityMatrix(BoundedMatrix<double,TDim+1, TDim+1> &DensityMatrix,
                                              double Density)
@@ -269,7 +266,6 @@ class GeoElementUtilities
         AddVectorAtPosition(rPBlockVector, rRightHandSideVector, offset);
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateNewtonCotesLocalShapeFunctionsGradients(BoundedMatrix<double,2,2>& DN_DeContainer)
     {
         //Line 2-noded
@@ -286,7 +282,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateNewtonCotesLocalShapeFunctionsGradients(BoundedMatrix<double,3,3>& DN_DeContainer)
     {
         //Line 3-noded
@@ -304,7 +299,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateNewtonCotesShapeFunctions(BoundedMatrix<double,3,3>& NContainer)
     {
         //Line 3-noded
@@ -319,7 +313,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateEquallyDistributedPointsLineShapeFunctions3N(Matrix& NContainer)
     {
         //Line 3-noded
@@ -339,7 +332,6 @@ class GeoElementUtilities
         }
     }
 
-//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateEquallyDistributedPointsLineGradientShapeFunctions3N(GeometryData::ShapeFunctionsGradientsType& DN_DeContainer)
     {
         //Line 3-noded
@@ -388,8 +380,6 @@ class GeoElementUtilities
         return Circumference;
     }
 
-
-    //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateExtrapolationMatrixTriangle(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
     {
         /// The matrix contains the shape functions at each GP evaluated at each node.
@@ -427,7 +417,6 @@ class GeoElementUtilities
 
     }
 
-    //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateExtrapolationMatrixQuad(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
     {
         if (rIntegrationMethod == GeometryData::IntegrationMethod::GI_GAUSS_1)
@@ -462,7 +451,6 @@ class GeoElementUtilities
         }
     }
 
-
     static inline void CalculateExtrapolationMatrixTetra(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
     {
         if (rIntegrationMethod == GeometryData::IntegrationMethod::GI_GAUSS_1)
@@ -498,8 +486,6 @@ class GeoElementUtilities
 
     }
 
-
-    //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
     static inline void CalculateExtrapolationMatrixHexa(Matrix& rExtrapolationMatrix, const GeometryData::IntegrationMethod& rIntegrationMethod)
     {
         if (rIntegrationMethod == GeometryData::IntegrationMethod::GI_GAUSS_1)
@@ -557,9 +543,6 @@ class GeoElementUtilities
         }
      }
 
-
-    //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-
     static Vector CalculateNodalHydraulicHeadFromWaterPressures(const GeometryType& rGeom, const Properties& rProp)
     {
         const auto NumericalLimit = std::numeric_limits<double>::epsilon();

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -49,6 +49,13 @@ void KratosGeoMechanicsApplication::Register() {
     KRATOS_REGISTER_ELEMENT("TransientPwElement3D20N", mTransientPwElement3D20N)
     KRATOS_REGISTER_ELEMENT("TransientPwElement3D27N", mTransientPwElement3D27N)
 
+    KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D2N", mTransientPwLineElement2D2N)
+    KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D3N", mTransientPwLineElement2D3N)
+    KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D4N", mTransientPwLineElement2D4N)
+    KRATOS_REGISTER_ELEMENT("TransientPwLineElement2D5N", mTransientPwLineElement2D5N)
+    KRATOS_REGISTER_ELEMENT("TransientPwLineElement3D2N", mTransientPwLineElement3D2N)
+    KRATOS_REGISTER_ELEMENT("TransientPwLineElement3D3N", mTransientPwLineElement3D3N)
+
     KRATOS_REGISTER_ELEMENT("TransientPwInterfaceElement2D4N", mTransientPwInterfaceElement2D4N)
     KRATOS_REGISTER_ELEMENT("TransientPwInterfaceElement3D6N", mTransientPwInterfaceElement3D6N)
     KRATOS_REGISTER_ELEMENT("TransientPwInterfaceElement3D8N", mTransientPwInterfaceElement3D8N)

applications/GeoMechanicsApplication/geo_mechanics_application.h

@@ -94,6 +94,7 @@
 #include "custom_elements/transient_Pw_interface_element.hpp"
 #include "custom_elements/undrained_U_Pw_small_strain_element.hpp"
 #include "custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp"
+#include "custom_elements/transient_Pw_line_element.h"
 
 // Element policies
 #include "custom_elements/axisymmetric_stress_state.h"
@@ -288,6 +289,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoMechanicsApplication : publ
     const TransientPwElement<3,20> mTransientPwElement3D20N{ 0, Kratos::make_shared< Hexahedra3D20    <NodeType> >(Element::GeometryType::PointsArrayType(20)), std::make_unique<ThreeDimensionalStressState>() };
     const TransientPwElement<3,27> mTransientPwElement3D27N{ 0, Kratos::make_shared< Hexahedra3D27    <NodeType> >(Element::GeometryType::PointsArrayType(27)), std::make_unique<ThreeDimensionalStressState>() };
 
+    const TransientPwLineElement<2, 2> mTransientPwLineElement2D2N{ 0, Kratos::make_shared<Line2D2<NodeType>>(Element::GeometryType::PointsArrayType(2)) };
+    const TransientPwLineElement<2, 3> mTransientPwLineElement2D3N{ 0, Kratos::make_shared<Line2D3<NodeType>>(Element::GeometryType::PointsArrayType(3)) };
+    const TransientPwLineElement<2, 4> mTransientPwLineElement2D4N{ 0, Kratos::make_shared<Line2D4<NodeType>>(Element::GeometryType::PointsArrayType(4)) };
+    const TransientPwLineElement<2, 5> mTransientPwLineElement2D5N{ 0, Kratos::make_shared<Line2D5<NodeType>>(Element::GeometryType::PointsArrayType(5)) };
+    const TransientPwLineElement<3, 2> mTransientPwLineElement3D2N{ 0, Kratos::make_shared<Line3D2<NodeType>>(Element::GeometryType::PointsArrayType(2))};
+    const TransientPwLineElement<3, 3> mTransientPwLineElement3D3N{ 0, Kratos::make_shared<Line3D3<NodeType>>(Element::GeometryType::PointsArrayType(3))};
+
     const TransientPwInterfaceElement<2,4> mTransientPwInterfaceElement2D4N{ 0, Kratos::make_shared< QuadrilateralInterface2D4 <NodeType> >(Element::GeometryType::PointsArrayType(4)), std::make_unique<PlaneStrainStressState>() };
     const TransientPwInterfaceElement<3,6> mTransientPwInterfaceElement3D6N{ 0, Kratos::make_shared< PrismInterface3D6         <NodeType> >(Element::GeometryType::PointsArrayType(6)), std::make_unique<ThreeDimensionalStressState>() };
     const TransientPwInterfaceElement<3,8> mTransientPwInterfaceElement3D8N{ 0, Kratos::make_shared< HexahedraInterface3D8     <NodeType> >(Element::GeometryType::PointsArrayType(8)), std::make_unique<ThreeDimensionalStressState>() };

applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py

@@ -46,6 +46,7 @@
 from test_prescribed_derivatives import KratosGeoMechanicsPrescribedDerivatives
 from test_dirichlet_u import KratosGeoMechanicsDirichletUTests
 from test_normal_load_on_hexa_element import KratosGeoMechanicsNormalLoadHexaTests
+from test_pressure_line_element import KratosGeoMechanicsTransientPressureLineElementTests
 
 def AssembleTestSuites():
     ''' Populates the test suites to run.
@@ -104,7 +105,8 @@ def AssembleTestSuites():
                         TestSellmeijersRule,
                         TestElementaryGroundWaterFlow,
                         KratosGeoMechanicsTransientThermalTests,
-                        KratosGeoMechanicsTimeIntegrationTests
+                        KratosGeoMechanicsTimeIntegrationTests,
+                        KratosGeoMechanicsTransientPressureLineElementTests
                         ]
     night_test_cases.extend(small_test_cases)
 

applications/GeoMechanicsApplication/tests/test_pressure_line_element.py

@@ -0,0 +1,63 @@
+import os
+
+import KratosMultiphysics.KratosUnittest as KratosUnittest
+import test_helper
+
+class KratosGeoMechanicsTransientPressureLineElementTests(KratosUnittest.TestCase):
+    """
+    This class contains benchmark tests which are checked with the regression on a previously obtained value.
+    """
+    etalon_value1 = -20000.0
+
+    def setUp(self):
+        # Code here will be placed BEFORE every test in this TestCase.
+        pass
+
+    def tearDown(self):
+        # Code here will be placed AFTER every test in this TestCase.
+        pass
+
+    def check_water_pressure(self, test_name, etalon_value):
+        file_path = test_helper.get_file_path(os.path.join('test_pressure_line_element', test_name))
+        simulation = test_helper.run_kratos(file_path)
+        pressure = test_helper.get_water_pressure(simulation)
+        self.assertAlmostEqual(etalon_value, pressure[2])
+
+    def test_oblique_line_element2D2N(self):
+        self.check_water_pressure("test_oblique_line_element2D2N", self.etalon_value1)
+
+    def test_oblique_line_element2D3N(self):
+        self.check_water_pressure("test_oblique_line_element2D3N", self.etalon_value1)
+
+    def test_oblique_line_element2D4N(self):
+        self.check_water_pressure("test_oblique_line_element2D4N", self.etalon_value1)
+
+    def test_oblique_line_element2D5N(self):
+        self.check_water_pressure("test_oblique_line_element2D5N", self.etalon_value1)
+
+    def test_vertical_line_element2D2N(self):
+        self.check_water_pressure("test_vertical_line_element2D2N", self.etalon_value1)
+
+    def test_vertical_line_element2D3N(self):
+        self.check_water_pressure("test_vertical_line_element2D3N", self.etalon_value1)
+
+    def test_vertical_line_element2D4N(self):
+        self.check_water_pressure("test_vertical_line_element2D4N", self.etalon_value1)
+
+    def test_vertical_line_element2D5N(self):
+        self.check_water_pressure("test_vertical_line_element2D5N", self.etalon_value1)
+
+    def test_oblique_line_element3D2N(self):
+        self.check_water_pressure("test_oblique_line_element3D2N", self.etalon_value1)
+
+    def test_oblique_line_element3D3N(self):
+        self.check_water_pressure("test_oblique_line_element3D3N", self.etalon_value1)
+
+    def test_vertical_line_element3D2N(self):
+        self.check_water_pressure("test_vertical_line_element3D2N", self.etalon_value1)
+
+    def test_vertical_line_element3D3N(self):
+        self.check_water_pressure("test_vertical_line_element3D3N", self.etalon_value1)
+
+if __name__ == '__main__':
+    KratosUnittest.main()

applications/GeoMechanicsApplication/tests/test_pressure_line_element/README.md

@@ -0,0 +1,26 @@
+# Test Cases for Water Pore Pressure Line Element
+
+**Author:** [Mohamed Nabi](https://github.com/mnabideltares)
+
+**Source files:** [Water pore pressure line element](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_pressure_line_element)
+
+## Case Specification
+In this water-pressure test case, a 3 [m] deep soil is considered, with everywhere set to 10 $\mathrm{[Pa]}$ as initial condition. Then the top boundary is set to 0 $\mathrm{[Pa]}$ and the bottom boundary is left free (it automatically becomes Neumann boundary with zero flux). The simulation spans 50 hours to allow for a transition from an exponential to a linear pressure profile along the depth. This test is conducted for various configurations, including 2D2N, 2D3N, 2D4N, 2D5N, 3D2N and 3D3N line elements. The pressure distribution along the depth is then evaluated with its own result.
+
+As the water pressure is influenced by gravity force in the vertical direction (here Y-direction), the gravity then needs to be projected in the direction of the element. Therefor here we tested two configurations, namely a case with vertically-oriented elements, and a case with elements with a slope of 45 degrees. The gravity is considered to be 10 $\mathrm{[m/s^2]}$ 
+
+The boundary conditions are shown below:
+
+<img src="documentation_data/test_pressure_line_element.svg" alt="Visualization of the Boundary conditions" title="Visualization of the Boundary conditions" width="600">
+
+## Results
+
+The picture below illustrates the pressure profile resulting from the simulation (as an example the 2D3N test is shown below).
+
+<img src="documentation_data/test_pressure_line_element_2d3n_result.png" alt="Pressure along depth at the last time step" title="Pressure along the depth at the last time step" width="600">
+
+These results are associated with the final time step after the solution reaches a steady state. The results for both test configurations are identical. The analytical solution is:
+
+$P = 10000 y$
+
+In this test case, the result at node number 3 at location $y = -2 \mathrm{[m]}$ is compared with the analytical solution. The value of the pressure at node 3 is -20,000 $\mathrm{[Pa]}$