[GeoMechanicsApplication] Line interface geometry can now generate edges

applications/GeoMechanicsApplication/custom_geometries/README.md

@@ -15,4 +15,6 @@ Similarly, the 3+3 line interface geometry has the following node numbering for
 
 ![3Plus3NodedGeometry](3Plus3NodedLineGeometry.svg)
 
-One thing to note, is that this line interface geometry does not implement functions from the Geometry base class which are related to the integration scheme. That is because most of the time, interface geometries are used with a Lobatto integration scheme, which is not supported by the Geometry base class.
+One thing to note, is that this line interface geometry does not implement functions from the Geometry base class which are related to the integration scheme. That is because most of the time, interface geometries are used with a Lobatto integration scheme, which is not supported by the Geometry base class.
+
+Any line interface geometry has two edges. The first edge coincides with the first side (i.e. bottom side in the above figures) of the geometry. The first edge's nodes are identical to the nodes at the first edge. For the 3+3 line interface geometry this means the list of node IDs equals [1, 2, 3]. The second edge references the nodes of the second side (i.e. top side in the above figures) of the line interface geometry. However, this edge traverses the second side in opposite direction. This implies that the end nodes will be swapped and any high-order nodes are reversed. For instance, for the 3+3 line interface geometry the second edge's nodes are ordered as follows: [5, 4, 6];

applications/GeoMechanicsApplication/custom_geometries/line_interface_geometry.h

@@ -208,6 +208,26 @@ class LineInterfaceGeometry : public Geometry<Node>
         KRATOS_ERROR << IntegrationSchemeFunctionalityNotImplementedMessage();
     }
 
+    GeometriesArrayType GenerateEdges() const override
+    {
+        const auto points = this->Points();
+
+        // The first edge coincides with the first side of the element
+        auto begin_of_second_side = points.ptr_begin() + (points.size() / 2);
+        const auto nodes_of_first_edge = PointerVector<Node>{points.ptr_begin(), begin_of_second_side};
+
+        // The second edge coincides with the second side of the element. However, the nodes must be
+        // traversed in opposite direction.
+        auto nodes_of_second_edge = PointerVector<Node>{begin_of_second_side, points.ptr_end()};
+        std::swap(nodes_of_second_edge(0), nodes_of_second_edge(1)); // end nodes
+        std::reverse(nodes_of_second_edge.ptr_begin() + 2, nodes_of_second_edge.ptr_end()); // any high-order nodes
+
+        auto result = GeometriesArrayType{};
+        result.push_back(std::make_shared<MidGeometryType>(nodes_of_first_edge));
+        result.push_back(std::make_shared<MidGeometryType>(nodes_of_second_edge));
+        return result;
+    }
+
     void CreateIntegrationPoints(IntegrationPointsArrayType& rIntegrationPoints,
                                  IntegrationInfo&            rIntegrationInfo) const override
     {

applications/GeoMechanicsApplication/tests/cpp_tests/custom_geometries/test_interface_geometry.cpp

@@ -46,6 +46,17 @@ auto CreateThreePlusThreeNoded2DLineInterfaceGeometry()
     return LineInterfaceGeometry<Line2D3<Node>>{1, nodes};
 }
 
+void AssertNodeIdsOfGeometry(const Geometry<Node>::Pointer&  rGeometryPtr,
+                             const std::vector<std::size_t>& rExpectedNodeIds)
+{
+    KRATOS_EXPECT_NE(rGeometryPtr, nullptr);
+
+    auto node_ids = std::vector<std::size_t>{};
+    std::transform(rGeometryPtr->begin(), rGeometryPtr->end(), std::back_inserter(node_ids),
+                   [](const auto& rNode) { return rNode.Id(); });
+    KRATOS_EXPECT_VECTOR_EQ(node_ids, rExpectedNodeIds)
+}
+
 } // namespace
 
 using namespace Kratos;
@@ -363,6 +374,30 @@ KRATOS_TEST_CASE_IN_SUITE(GetLocalCoordinatesOfAllNodesOfThreePlusThreeNodedLine
     KRATOS_EXPECT_MATRIX_NEAR(result, expected_result, 1e-6)
 }
 
+KRATOS_TEST_CASE_IN_SUITE(TwoPlusTwoLineInterfaceGeometryHasTwoEdgesWithOppositeOrientations,
+                          KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    const auto geometry = CreateTwoPlusTwoNoded2DLineInterfaceGeometry();
+
+    const auto edges = geometry.GenerateEdges();
+
+    KRATOS_EXPECT_EQ(edges.size(), 2);
+    AssertNodeIdsOfGeometry(edges(0), {1, 2});
+    AssertNodeIdsOfGeometry(edges(1), {4, 3});
+}
+
+KRATOS_TEST_CASE_IN_SUITE(ThreePlusThreeLineInterfaceGeometryHasTwoEdgesWithOppositeOrientations,
+                          KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    const auto geometry = CreateThreePlusThreeNoded2DLineInterfaceGeometry();
+
+    const auto edges = geometry.GenerateEdges();
+
+    KRATOS_EXPECT_EQ(edges.size(), 2);
+    AssertNodeIdsOfGeometry(edges(0), {1, 2, 3});
+    AssertNodeIdsOfGeometry(edges(1), {5, 4, 6});
+}
+
 KRATOS_TEST_CASE_IN_SUITE(InterfaceGeometry_Throws_WhenCallingArea, KratosGeoMechanicsFastSuiteWithoutKernel)
 {
     const auto geometry = CreateThreePlusThreeNoded2DLineInterfaceGeometry();