def compute_global_dof(num_elem, elem, row, col):
    """
    Given the two LOCAL indices row and col, the function return two arrays containing 
    row GLOBAL indices and col GLOBAL indices for all elem. In other words it map 
    the local dof indices to global dof indices.
    
    Parameters
    ----------
    num_elem : int
        number of elem in mesh
    elem : (num_elem, nodes_per_element) numpy.ndarray
        elem map (connectivity map), n-th row contains tags of nodes in n-th element 
    row : int
        row index of the local stiffness matrix entry
    col : int
        column index of the local stiffness matrix entry
    
    Returns
    -------
    row_ind: (num_elem,) numpy.ndarray
        array of global row indices related to a certain local entry
    col_ind: (num_elem,) numpy.ndarray
        array of global column indices related to a certain local entry
    """

    if (row % 2 == 0):
        row_ind = elem[:, row // 2] * 2
    else:
        row_ind = elem[:, row // 2] * 2 + 1
    
    if (col % 2 == 0):
        col_ind = elem[:, col // 2] * 2
    else:
        col_ind = elem[:, col // 2] * 2 + 1

    return row_ind, col_ind