Refactor BuildSystemStructure to support row/column block connectivity and improve error handling

Author: loumalouomega

applications/TrilinosApplication/tests/cpp_tests/test_trilinos_space.cpp

@@ -845,7 +845,7 @@ KRATOS_TEST_CASE_IN_SUITE(TrilinosBuildSystemStructure, KratosTrilinosApplicatio
     pMap = Kratos::make_shared<Epetra_Map>(-1, local_size, local_ids.data(), 0, epetra_comm);
 
     TrilinosSparseSpaceType::BuildSystemStructure(
-        epetra_comm, local_size, first_my_id, 5, all_equation_ids,
+        epetra_comm, local_size, first_my_id, 5, all_equation_ids, all_equation_ids,
         pA, pb, pDx, pReactions, system_size, pMap
     );
 
@@ -856,6 +856,43 @@ KRATOS_TEST_CASE_IN_SUITE(TrilinosBuildSystemStructure, KratosTrilinosApplicatio
     KRATOS_EXPECT_EQ(static_cast<std::size_t>(system_size), TrilinosSparseSpaceType::Size1(*pA));
 }
 
+KRATOS_TEST_CASE_IN_SUITE(TrilinosBuildSystemStructureRowColumnBlocks, KratosTrilinosApplicationMPITestSuite)
+{
+    const auto& r_comm = Testing::GetDefaultDataCommunicator();
+    auto raw_mpi_comm = MPIDataCommunicator::GetMPICommunicator(r_comm);
+    Epetra_MpiComm epetra_comm(raw_mpi_comm);
+
+    const int local_size = 2;
+    const int first_my_id = r_comm.Rank() * local_size;
+    const int system_size = 2 * r_comm.Size();
+
+    // Rectangular row/column blocks: one diagonal and one cross entry per rank
+    std::vector<std::vector<int>> all_row_equation_ids = {{first_my_id}, {first_my_id}};
+    std::vector<std::vector<int>> all_col_equation_ids = {{first_my_id}, {first_my_id + 1}};
+
+    TrilinosSparseSpaceType::MatrixPointerType pA;
+    TrilinosSparseSpaceType::VectorPointerType pb;
+    TrilinosSparseSpaceType::VectorPointerType pDx;
+    TrilinosSparseSpaceType::VectorPointerType pReactions;
+
+    TrilinosSparseSpaceType::MapPointerType pMap;
+    std::vector<int> local_ids(local_size);
+    for (int i = 0; i < local_size; i++) local_ids[i] = first_my_id + i;
+    pMap = Kratos::make_shared<Epetra_Map>(-1, local_size, local_ids.data(), 0, epetra_comm);
+
+    TrilinosSparseSpaceType::BuildSystemStructure(
+        epetra_comm, local_size, first_my_id, 5, all_row_equation_ids, all_col_equation_ids,
+        pA, pb, pDx, pReactions, system_size, pMap
+    );
+
+    KRATOS_EXPECT_FALSE(TrilinosSparseSpaceType::IsNull(pA));
+    KRATOS_EXPECT_FALSE(TrilinosSparseSpaceType::IsNull(pb));
+    KRATOS_EXPECT_FALSE(TrilinosSparseSpaceType::IsNull(pDx));
+    KRATOS_EXPECT_FALSE(TrilinosSparseSpaceType::IsNull(pReactions));
+    KRATOS_EXPECT_EQ(static_cast<std::size_t>(system_size), TrilinosSparseSpaceType::Size1(*pA));
+    KRATOS_EXPECT_EQ(2 * r_comm.Size(), pA->NumGlobalNonzeros());
+}
+
 KRATOS_TEST_CASE_IN_SUITE(TrilinosBuildConstraintsStructure, KratosTrilinosApplicationMPITestSuite)
 {
     const auto& r_comm = Testing::GetDefaultDataCommunicator();

applications/TrilinosApplication/trilinos_space.h

@@ -907,12 +907,13 @@ class TrilinosSpace
     }
 
     /**
-     * @brief Build Epetra FECrsGraph and create new system matrix + vectors.
+     * @brief Build Epetra FECrsGraph and create new system matrix + vectors from row/column block connectivity.
      * @param rComm The communicator considered
      * @param LocalSize The local size of the system
      * @param FirstMyId The first global id owned by this rank
      * @param GuessRowSize The guess row size for the graph construction
-     * @param rAllEquationIds The list of lists of equation ids for each local row
+     * @param rAllRowEquationIds The list of row equation-id blocks
+     * @param rAllColEquationIds The list of column equation-id blocks
      * @param rpA The pointer to the matrix to be created
      * @param rpb The pointer to the right-hand side vector to be created
      * @param rpDx The pointer to the solution vector to be created
@@ -925,7 +926,8 @@ class TrilinosSpace
         const IndexType LocalSize,
         const int FirstMyId,
         const int GuessRowSize,
-        const std::vector<std::vector<int>>& rAllEquationIds,
+        const std::vector<std::vector<int>>& rAllRowEquationIds,
+        const std::vector<std::vector<int>>& rAllColEquationIds,
         MatrixPointerType& rpA,
         VectorPointerType& rpb,
         VectorPointerType& rpDx,
@@ -934,11 +936,16 @@ class TrilinosSpace
         MapPointerType pMap
         )
     {
+        KRATOS_ERROR_IF(rAllRowEquationIds.size() != rAllColEquationIds.size())
+            << "BuildSystemStructure: row and column block lists must have the same size" << std::endl;
+
         // Create graph
         Epetra_FECrsGraph graph(::Copy, *pMap, GuessRowSize);
-        for (const auto& r_ids : rAllEquationIds) {
-            if (r_ids.size() != 0) {
-                graph.InsertGlobalIndices(r_ids.size(), r_ids.data(), r_ids.size(), r_ids.data());
+        for (std::size_t i_block = 0; i_block < rAllRowEquationIds.size(); ++i_block) {
+            const auto& r_row_ids = rAllRowEquationIds[i_block];
+            const auto& r_col_ids = rAllColEquationIds[i_block];
+            if (!r_row_ids.empty() && !r_col_ids.empty()) {
+                graph.InsertGlobalIndices(r_row_ids.size(), r_row_ids.data(), r_col_ids.size(), r_col_ids.data());
             }
         }
         graph.GlobalAssemble();