Line interface geometry can now generate edges (#12684)

The new line interface geometry now supports generating its edges. Note the that first edge's points are identical to the points of the first side. The second edge's points reference the points at the second side. However, since this edge has opposite orientation, the end points will be swapped and any intermediate points will be reversed, too.

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 ordering of the first edge's nodes is identical to the one of the first side. 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 has opposite orientation compared to the side, and therefore a node numbering according to the opposite direction. With respect to the node numbering of the side, 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();