Skip to content

Commit ed79031

Browse files
loumalouomegaloumalouomega
committed
Refactor CalculateReactions method to improve Epetra handling and streamline index management
1 parent c8543aa commit ed79031

1 file changed

Lines changed: 117 additions & 145 deletions

File tree

‎applications/TrilinosApplication/custom_strategies/builder_and_solvers/trilinos_block_builder_and_solver.h‎

Lines changed: 117 additions & 145 deletions
Original file line numberDiff line numberDiff line change
@@ -812,66 +812,66 @@ class TrilinosBlockBuilderAndSolver
812812
BuildRHS(pScheme, rModelPart, rb);
813813

814814
if constexpr (TSparseSpace::LinearAlgebraLibrary() == TrilinosLinearAlgebraLibrary::EPETRA) {
815-
// Initialize the Epetra importer
816-
// TODO: this part of the code has been pasted until a better solution
817-
// is found
818-
int system_size = TSparseSpace::Size(rb);
819-
int number_of_dofs = BaseType::mDofSet.size();
820-
std::vector<int> index_array(number_of_dofs);
821-
822-
// Filling the array with the global ids
823-
int counter = 0;
824-
int id = 0;
825-
for (const auto& dof : BaseType::mDofSet) {
826-
id = dof.EquationId();
827-
if (id < system_size) {
828-
index_array[counter++] = id;
815+
// Initialize the Epetra importer
816+
// TODO: this part of the code has been pasted until a better solution
817+
// is found
818+
int system_size = TSparseSpace::Size(rb);
819+
int number_of_dofs = BaseType::mDofSet.size();
820+
std::vector<int> index_array(number_of_dofs);
821+
822+
// Filling the array with the global ids
823+
int counter = 0;
824+
int id = 0;
825+
for (const auto& dof : BaseType::mDofSet) {
826+
id = dof.EquationId();
827+
if (id < system_size) {
828+
index_array[counter++] = id;
829+
}
829830
}
830-
}
831831

832-
std::sort(index_array.begin(), index_array.end());
833-
std::vector<int>::iterator NewEnd = std::unique(index_array.begin(), index_array.end());
834-
index_array.resize(NewEnd - index_array.begin());
832+
std::sort(index_array.begin(), index_array.end());
833+
std::vector<int>::iterator NewEnd = std::unique(index_array.begin(), index_array.end());
834+
index_array.resize(NewEnd - index_array.begin());
835835

836-
int check_size = -1;
837-
int tot_update_dofs = index_array.size();
838-
rb.Comm().SumAll(&tot_update_dofs, &check_size, 1);
839-
KRATOS_ERROR_IF(check_size < system_size)
840-
<< "Dof count is not correct. There are less dofs than expected.\n"
841-
<< "Expected number of active dofs = " << system_size
842-
<< " dofs found = " << check_size << std::endl;
836+
int check_size = -1;
837+
int tot_update_dofs = index_array.size();
838+
rb.Comm().SumAll(&tot_update_dofs, &check_size, 1);
839+
KRATOS_ERROR_IF(check_size < system_size)
840+
<< "Dof count is not correct. There are less dofs than expected.\n"
841+
<< "Expected number of active dofs = " << system_size
842+
<< " dofs found = " << check_size << std::endl;
843843

844-
// Defining a map as needed
845-
Epetra_Map dof_update_map(-1, index_array.size(),
846-
&(*(index_array.begin())), 0, rb.Comm());
844+
// Defining a map as needed
845+
Epetra_Map dof_update_map(-1, index_array.size(),
846+
&(*(index_array.begin())), 0, rb.Comm());
847847

848-
// Defining the importer class
849-
Kratos::shared_ptr<Epetra_Import> pDofImporter = Kratos::make_shared<Epetra_Import>(dof_update_map, rb.Map());
848+
// Defining the importer class
849+
Kratos::shared_ptr<Epetra_Import> pDofImporter = Kratos::make_shared<Epetra_Import>(dof_update_map, rb.Map());
850850

851-
// Defining a temporary vector to gather all of the values needed
852-
Epetra_Vector temp_RHS(pDofImporter->TargetMap());
851+
// Defining a temporary vector to gather all of the values needed
852+
Epetra_Vector temp_RHS(pDofImporter->TargetMap());
853853

854-
// Importing in the new temp_RHS vector the values
855-
int ierr = temp_RHS.Import(rb, *pDofImporter, Insert);
856-
KRATOS_ERROR_IF(ierr != 0) << "Epetra failure found - error code: " << ierr << std::endl;
854+
// Importing in the new temp_RHS vector the values
855+
int ierr = temp_RHS.Import(rb, *pDofImporter, Insert);
856+
KRATOS_ERROR_IF(ierr != 0) << "Epetra failure found - error code: " << ierr << std::endl;
857857

858-
double* temp_RHS_values; // DO NOT make delete of this one!!
859-
temp_RHS.ExtractView(&temp_RHS_values);
858+
double* temp_RHS_values; // DO NOT make delete of this one!!
859+
temp_RHS.ExtractView(&temp_RHS_values);
860860

861-
rb.Comm().Barrier();
861+
rb.Comm().Barrier();
862862

863-
const int ndofs = static_cast<int>(BaseType::mDofSet.size());
863+
const int ndofs = static_cast<int>(BaseType::mDofSet.size());
864864

865-
// Store the RHS values in the reaction variable
866-
// NOTE: dofs are assumed to be numbered consecutively in the
867-
// BlockBuilderAndSolver
868-
for (int k = 0; k < ndofs; k++) {
869-
auto dof_iterator = BaseType::mDofSet.begin() + k;
865+
// Store the RHS values in the reaction variable
866+
// NOTE: dofs are assumed to be numbered consecutively in the
867+
// BlockBuilderAndSolver
868+
for (int k = 0; k < ndofs; k++) {
869+
auto dof_iterator = BaseType::mDofSet.begin() + k;
870870

871-
const int i = (dof_iterator)->EquationId();
872-
// (dof_iterator)->GetSolutionStepReactionValue() = -(*b[i]);
873-
const double react_val = temp_RHS[pDofImporter->TargetMap().LID(i)];
874-
(dof_iterator->GetSolutionStepReactionValue()) = -react_val;
871+
const int i = (dof_iterator)->EquationId();
872+
// (dof_iterator)->GetSolutionStepReactionValue() = -(*b[i]);
873+
const double react_val = temp_RHS[pDofImporter->TargetMap().LID(i)];
874+
(dof_iterator->GetSolutionStepReactionValue()) = -react_val;
875875
}
876876
} else {
877877
KRATOS_ERROR << "CalculateReactions only implemented for Epetra" << std::endl;
@@ -914,49 +914,49 @@ class TrilinosBlockBuilderAndSolver
914914

915915
if constexpr (TSparseSpace::LinearAlgebraLibrary() == TrilinosLinearAlgebraLibrary::EPETRA) {
916916
// Here we construct and fill a vector "fixed local" which cont
917-
Epetra_Map localmap(-1, global_ids.size(), global_ids.data(), 0, rA.Comm());
918-
Epetra_IntVector fixed_local(Copy, localmap, is_dirichlet.data());
917+
Epetra_Map localmap(-1, global_ids.size(), global_ids.data(), 0, rA.Comm());
918+
Epetra_IntVector fixed_local(Copy, localmap, is_dirichlet.data());
919919

920-
Epetra_Import dirichlet_importer(rA.ColMap(), fixed_local.Map());
920+
Epetra_Import dirichlet_importer(rA.ColMap(), fixed_local.Map());
921921

922-
// defining a temporary vector to gather all of the values needed
923-
Epetra_IntVector fixed(rA.ColMap());
922+
// defining a temporary vector to gather all of the values needed
923+
Epetra_IntVector fixed(rA.ColMap());
924924

925-
// Detect if there is a line of all zeros and set the diagonal to a certain number (1 if not scale, some norms values otherwise) if this happens
926-
const auto& r_process_info = rModelPart.GetProcessInfo();
927-
mScaleFactor = TSparseSpace::CheckAndCorrectZeroDiagonalValues(r_process_info, rA, rb, mScalingDiagonal);
925+
// Detect if there is a line of all zeros and set the diagonal to a certain number (1 if not scale, some norms values otherwise) if this happens
926+
const auto& r_process_info = rModelPart.GetProcessInfo();
927+
mScaleFactor = TSparseSpace::CheckAndCorrectZeroDiagonalValues(r_process_info, rA, rb, mScalingDiagonal);
928928

929-
// Importing in the new temp vector the values
930-
int ierr = fixed.Import(fixed_local, dirichlet_importer, Insert);
931-
KRATOS_ERROR_IF(ierr != 0) << "Epetra failure found" << std::endl;
929+
// Importing in the new temp vector the values
930+
int ierr = fixed.Import(fixed_local, dirichlet_importer, Insert);
931+
KRATOS_ERROR_IF(ierr != 0) << "Epetra failure found" << std::endl;
932932

933-
for (int i = 0; i < rA.NumMyRows(); i++) {
934-
int numEntries; // Number of non-zero entries
935-
double* vals; // Row non-zero values
936-
int* cols; // Column indices of row non-zero values
937-
rA.ExtractMyRowView(i, numEntries, vals, cols);
933+
for (int i = 0; i < rA.NumMyRows(); i++) {
934+
int numEntries; // Number of non-zero entries
935+
double* vals; // Row non-zero values
936+
int* cols; // Column indices of row non-zero values
937+
rA.ExtractMyRowView(i, numEntries, vals, cols);
938938

939-
const int row_gid = rA.RowMap().GID(i);
940-
const int row_lid = localmap.LID(row_gid);
939+
const int row_gid = rA.RowMap().GID(i);
940+
const int row_lid = localmap.LID(row_gid);
941941

942-
if (fixed_local[row_lid] == 0) { // Not a dirichlet row
943-
for (int j = 0; j < numEntries; j++) {
944-
if (fixed[cols[j]] == true)
945-
vals[j] = 0.0;
946-
}
947-
} else { // This IS a dirichlet row
948-
// Set to zero the rhs
949-
rb[0][i] = 0.0; // note that the index of i is expected to be
950-
// coherent with the rows of A
951-
952-
// Set to zero the whole row
953-
for (int j = 0; j < numEntries; j++) {
954-
int col_gid = rA.ColMap().GID(cols[j]);
955-
if (col_gid != row_gid)
956-
vals[j] = 0.0;
942+
if (fixed_local[row_lid] == 0) { // Not a dirichlet row
943+
for (int j = 0; j < numEntries; j++) {
944+
if (fixed[cols[j]] == true)
945+
vals[j] = 0.0;
946+
}
947+
} else { // This IS a dirichlet row
948+
// Set to zero the rhs
949+
rb[0][i] = 0.0; // note that the index of i is expected to be
950+
// coherent with the rows of A
951+
952+
// Set to zero the whole row
953+
for (int j = 0; j < numEntries; j++) {
954+
int col_gid = rA.ColMap().GID(cols[j]);
955+
if (col_gid != row_gid)
956+
vals[j] = 0.0;
957+
}
957958
}
958959
}
959-
}
960960
} else {
961961
KRATOS_ERROR << "Only Epetra_MpiComm is supported for now" << std::endl;
962962
}
@@ -1027,9 +1027,9 @@ class TrilinosBlockBuilderAndSolver
10271027

10281028
// Compute T' A T
10291029
if constexpr (TSparseSpace::LinearAlgebraLibrary() == TrilinosLinearAlgebraLibrary::EPETRA) {
1030-
const TSystemMatrixType copy_A(rA);
1031-
TSparseSpace::BtDBProductOperation(rA, copy_A, r_T, true, true);
1032-
} else {
1030+
const TSystemMatrixType copy_A(rA);
1031+
TSparseSpace::BtDBProductOperation(rA, copy_A, r_T, true, true);
1032+
} else {
10331033
KRATOS_ERROR << "Only EPETRA is supported for now" << std::endl;
10341034
}
10351035

@@ -1264,82 +1264,54 @@ class TrilinosBlockBuilderAndSolver
12641264
// Generate indices database
12651265
const IndexType number_of_local_rows = mLocalSystemSize;
12661266

1267-
// Generate map - use the "temp" array here
1268-
const int temp_size = (number_of_local_rows < 1000) ? 1000 : number_of_local_rows;
1269-
std::vector<int> temp_primary(temp_size, 0);
1270-
std::vector<int> temp_secondary(temp_size, 0);
1271-
for (IndexType i = 0; i != number_of_local_rows; i++) {
1272-
temp_primary[i] = mFirstMyId + i;
1273-
}
1274-
auto& r_map = GetMap();
1275-
std::fill(temp_primary.begin(), temp_primary.begin() + number_of_local_rows, 0);
1276-
1277-
// The T graph
1278-
Epetra_FECrsGraph Tgraph(Copy, r_map, mGuessRowSize);
1279-
1280-
// Adding diagonal values
1281-
int ierr;
1282-
int index_diagonal[1] = {0};
1283-
for (IndexType i = 0; i < number_of_local_rows; ++i) {
1284-
index_diagonal[0] = mFirstMyId + i;
1285-
ierr = Tgraph.InsertGlobalIndices(1, index_diagonal, 1, index_diagonal);
1286-
KRATOS_ERROR_IF(ierr < 0) << ": Epetra failure in Graph.InsertGlobalIndices. Error code: " << ierr << std::endl;
1287-
}
1267+
// Ensure map is initialized
1268+
GetMap();
12881269

12891270
// Vector containing indices belonging to slave DoFs, not used for graph, but for master/slave index identification
1290-
std::unordered_set<std::size_t> indices;
1271+
std::unordered_set<IndexType> indices;
12911272

12921273
// TODO: Check if these should be local constraints
12931274
auto& r_constraints_array = rModelPart.MasterSlaveConstraints();
12941275

1276+
// Auxiliary vectors to construct the relation matrix structure
1277+
std::vector<std::vector<int>> all_slave_equation_ids;
1278+
std::vector<std::vector<int>> all_master_equation_ids;
1279+
all_slave_equation_ids.reserve(r_constraints_array.size());
1280+
all_master_equation_ids.reserve(r_constraints_array.size());
1281+
12951282
// Assemble all constraints
12961283
Element::EquationIdVectorType slave_equation_ids_vector, master_equation_ids_vector;
12971284
for (auto& r_const : r_constraints_array) {
12981285
r_const.EquationIdVector(slave_equation_ids_vector, master_equation_ids_vector, r_current_process_info);
12991286

1300-
// Filling the list of active global indices (non fixed)
1301-
IndexType num_active_slave_indices = 0;
1302-
for (auto& r_slave_id : slave_equation_ids_vector) {
1303-
temp_primary[num_active_slave_indices] = r_slave_id;
1304-
++num_active_slave_indices;
1305-
}
1306-
IndexType num_active_master_indices = 0;
1307-
for (auto& r_master_id : master_equation_ids_vector) {
1308-
temp_secondary[num_active_master_indices] = r_master_id;
1309-
++num_active_master_indices;
1287+
std::vector<int> slave_equation_ids;
1288+
slave_equation_ids.reserve(slave_equation_ids_vector.size());
1289+
for (const auto r_slave_id : slave_equation_ids_vector) {
1290+
indices.insert(r_slave_id);
1291+
slave_equation_ids.push_back(static_cast<int>(r_slave_id));
13101292
}
13111293

1312-
// Adding cross master-slave dofs
1313-
if (num_active_slave_indices > 0 && num_active_master_indices > 0) {
1314-
int slave_index[1] = {0};
1315-
for (IndexType i = 0; i < num_active_slave_indices; ++i) {
1316-
slave_index[0] = temp_primary[i];
1317-
indices.insert(temp_primary[i]);
1318-
ierr = Tgraph.InsertGlobalIndices(1, slave_index, num_active_master_indices, temp_secondary.data());
1319-
KRATOS_ERROR_IF(ierr < 0) << ": Epetra failure in Graph.InsertGlobalIndices. Error code: " << ierr << std::endl;
1320-
}
1294+
std::vector<int> master_equation_ids;
1295+
master_equation_ids.reserve(master_equation_ids_vector.size());
1296+
for (const auto r_master_id : master_equation_ids_vector) {
1297+
master_equation_ids.push_back(static_cast<int>(r_master_id));
13211298
}
1322-
std::fill(temp_primary.begin(), temp_primary.begin() + num_active_slave_indices, 0);
1323-
std::fill(temp_secondary.begin(), temp_secondary.begin() + num_active_master_indices, 0);
1324-
}
13251299

1326-
// Finalizing graph construction
1327-
ierr = Tgraph.GlobalAssemble();
1328-
KRATOS_ERROR_IF(ierr < 0) << ": Epetra failure in Epetra_FECrsGraph.GlobalAssemble. Error code: " << ierr << std::endl;
1329-
ierr = Tgraph.FillComplete();
1330-
KRATOS_ERROR_IF(ierr < 0) << ": Epetra failure in Epetra_FECrsGraph.FillComplete. Error code: " << ierr << std::endl;
1331-
ierr = Tgraph.OptimizeStorage();
1332-
KRATOS_ERROR_IF(ierr < 0) << ": Epetra failure in Epetra_FECrsGraph.OptimizeStorage. Error code: " << ierr << std::endl;
1333-
1334-
// Generate a new matrix pointer according to this non-zero values
1335-
TSystemMatrixPointerType p_new_T = TSystemMatrixPointerType(new TSystemMatrixType(Copy, Tgraph));
1336-
1337-
// Swap matrix
1338-
mpT.swap(p_new_T);
1300+
all_slave_equation_ids.push_back(std::move(slave_equation_ids));
1301+
all_master_equation_ids.push_back(std::move(master_equation_ids));
1302+
}
13391303

1340-
// Generate the constant vector equivalent
1341-
TSystemVectorPointerType p_new_constant_vector = TSystemVectorPointerType(new TSystemVectorType(r_map));
1342-
mpConstantVector.swap(p_new_constant_vector);
1304+
TSparseSpace::BuildConstraintsStructure(
1305+
mrComm,
1306+
number_of_local_rows,
1307+
mFirstMyId,
1308+
mGuessRowSize,
1309+
all_slave_equation_ids,
1310+
all_master_equation_ids,
1311+
mpT,
1312+
mpConstantVector,
1313+
mpMap
1314+
);
13431315

13441316
/* Fill ids for master/slave */
13451317

0 commit comments

Comments
 (0)