def bc_robin(matrix_lhs, matrix_rhs, positions, grad_n_bc, a_val, b_val, c_val
    ):
    """
    Function to enforce the Boundary Condition in the Robin formalism:

                    a * (df / dn) + b * f + c = O

    Parameters
    ----------
    matrix_lhs: Left Hand Side square matrix (shape (size_i_w * size_j_ns)^2)
    matrix_rhs: Right Hand Side vector (shape (size_i_w * size_j_ns))
    positions: indices of lines corresponding to the patch being processed
    grad_n_bc: normal gradient to boundary nodes
    a_val: value of the "a" Robin parameter
    b_val: value of the "b" Robin parameter
    c_val: value of the "a" Robin parameter

    Returns
    -------
    Altered left-hand side matrix and right-hand side vector
    """
    for row in positions:
        if grad_n_bc.count_nonzero():
            matrix_lhs[(row), :] = a_val * grad_n_bc[(row), :]
        matrix_lhs[row, row] += b_val
    matrix_rhs[positions] = -c_val
    return matrix_lhs, matrix_rhs
