Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Core] Adding debug check in CheckAndCorrectZeroDiagonalValues #12782

29 changes: 23 additions & 6 deletions kratos/spaces/ublas_space.h
Original file line number Diff line number Diff line change
Expand Up @@ -708,10 +708,16 @@ class UblasSpace
const SCALING_DIAGONAL ScalingDiagonal = SCALING_DIAGONAL::NO_SCALING
)
{
// The system size
const std::size_t system_size = rA.size1();

const auto& Avalues = rA.value_data();
const auto& Arow_indices = rA.index1_data();
// The matrix data
auto& r_Avalues = rA.value_data();
const auto& r_Arow_indices = rA.index1_data();
const auto& r_Acol_indices = rA.index2_data();

// Define the iterators
const auto& it_Acol_indices_begin = r_Acol_indices.begin();
loumalouomega marked this conversation as resolved.
Show resolved Hide resolved

// Define zero value tolerance
const double zero_tolerance = std::numeric_limits<double>::epsilon();
Expand All @@ -723,18 +729,29 @@ class UblasSpace
IndexPartition(system_size).for_each([&](std::size_t Index){
bool empty = true;

const std::size_t col_begin = Arow_indices[Index];
const std::size_t col_end = Arow_indices[Index + 1];
const std::size_t col_begin = r_Arow_indices[Index];
const std::size_t col_end = r_Arow_indices[Index + 1];

for (std::size_t j = col_begin; j < col_end; ++j) {
if(std::abs(Avalues[j]) > zero_tolerance) {
if(std::abs(r_Avalues[j]) > zero_tolerance) {
empty = false;
break;
}
}

if(empty) {
rA(Index, Index) = scale_factor;
const auto& it_Acol_indices_row_begin = r_Acol_indices.begin() + col_begin;
const auto& it_Acol_indices_row_end = r_Acol_indices.begin() + col_end;
loumalouomega marked this conversation as resolved.
Show resolved Hide resolved
const auto pos = std::lower_bound(it_Acol_indices_row_begin, it_Acol_indices_row_end, Index);
if (pos != it_Acol_indices_row_end) { // It was found
loumalouomega marked this conversation as resolved.
Show resolved Hide resolved
r_Avalues[pos - it_Acol_indices_begin] = scale_factor;
#ifdef KRATOS_DEBUG
} else {
KRATOS_ERROR << "Diagonal term (" << Index << ", " << Index << ") is not defined in the system matrix" << std::endl;
}
#else
}
loumalouomega marked this conversation as resolved.
Show resolved Hide resolved
#endif
rb[Index] = 0.0;
}
});
Expand Down
Loading