def get_decouplings_for_line(grid, line, settings, fids=None):
    """Compute the element pairs for regulazisation decoupling, given a grid
    object and a line, denoted by start and end coordinates, in the grid.

    Parameters
    ----------
    grid: a crtomo.grid.crt_grid object
    line: np.array/list with 4 entries: x1, y1, x2, x2, denoting start and end
          point
    fids: None, or: list with file ids for
        debug_nodes



    Returns
    -------
    neighbors: Mx2 array with the numbers of adjoining elements along the given
               line
    """
    key = 'presort'
    # x/z coordinates of 'macro'-line
    x = [line[0], line[2]]
    y = [line[1], line[3]]

    # x/z coordinates of all nodes in the grid
    nx = grid.nodes[key][:, 1]
    ny = grid.nodes[key][:, 2]
    nxy = np.vstack((nx, ny)).T

    # shortest distance of all nodes to this line
    dist = np.array(distances(x, y, nx, ny))

    # set the epsilon environment
    if settings.get('eps', None) is not None:
        eps = settings['eps']
        print('Taking user supplied eps value: {0}'.format(eps))
    else:
        import scipy.spatial.distance as spdist
        pair_distances = spdist.pdist(nxy)
        print('automatic eps determination')
        print('minimal distance: {0}'.format(pair_distances.min()))
        eps = np.sqrt(
            pair_distances.min() ** 2 - (pair_distances.min() ** 2) / 4
        )
        print('estimated smallest distance to next node: {0}'.format(eps))
        eps *= 0.5
        print('final eps (0.5 * estimate): {0}'.format(eps))
        if eps == 0:
            raise Exception('eps == 0')

    # only consider nodes that lie in the eps environment
    indices = np.where(dist < eps)[0]

    # extract these nodes
    nodes = grid.nodes[key][indices, 0]
    nodes = indices

    nodes_full = np.hstack((
        grid.nodes[key][indices, :].squeeze(),
        np.atleast_2d(dist[indices]).T,
    ))
    if fids is not None:
        np.savetxt(
            fids[0],
            nodes_full.squeeze()
        )
        np.savetxt(
            fids[1],
            np.atleast_2d(dist).T,
        )

    elm_indices = []
    elm_nodes = []
    # find elements with two nodes in it
    for elmnr, element in enumerate(grid.elements):
        if len(np.intersect1d(element, nodes)) == 2:
            elm_indices.append(elmnr)
            elm_nodes.append(element)

    elms = np.array(elm_nodes)

    # find neighboring elements
    neighbors = []
    eta = settings['eta']
    # for each element on the line
    for index, elm in enumerate(elms):
        found_it = False
        for index1, elm1 in enumerate(elms):
            ints = np.intersect1d(elm, elm1)
            # only two nodes can be on a line
            if len(ints) == 2:
                # this check ensures that we do not identify the
                # boundary between two adjacent cells on the line
                # as a decoupling line.
                if np.all(dist[ints] < eps):
                    found_it = True
                    break
                else:
                    pass
                    # print(
                    #     'distances of neighbor not correct:', dist[ints], eps
                    # )
            # probably the element itself
            elif len(ints) > 2:
                # print('found more than two common nodes!')
                # print('element nodes:', elm)
                # print('element1 nodes:', elm1)
                pass
        if found_it:
            nb = (elm_indices[index] + 1, elm_indices[index1] + 1, eta)
            # reversed neighbor
            # nb_rev = (nb[1], nb[0], nb[2])
            # if nb_rev not in neighbors:
            if True:
                neighbors.append(
                    nb
                )
        else:
            pass
            print('No neighbors found. Strange...')

    return np.array(neighbors)
def get_decouplings_for_line ( grid , line , settings , fids = None ) : key = 'presort' x = [ line [ 0 ] , line [ 2 ] ] y = [ line [ 1 ] , line [ 3 ] ] nx = grid . nodes [ key ] [ : , 1 ] ny = grid . nodes [ key ] [ : , 2 ] nxy = np . vstack ( ( nx , ny ) ) . T dist = np . array ( distances ( x , y , nx , ny ) ) if settings . get ( 'eps' , None ) is not None : eps = settings [ 'eps' ] print ( 'Taking user supplied eps value: {0}' . format ( eps ) ) else : import scipy . spatial . distance as spdist pair_distances = spdist . pdist ( nxy ) print ( 'automatic eps determination' ) print ( 'minimal distance: {0}' . format ( pair_distances . min ( ) ) ) eps = np . sqrt ( pair_distances . min ( ) ** 2 - ( pair_distances . min ( ) ** 2 ) / 4 ) print ( 'estimated smallest distance to next node: {0}' . format ( eps ) ) eps *= 0.5 print ( 'final eps (0.5 * estimate): {0}' . format ( eps ) ) if eps == 0 : raise Exception ( 'eps == 0' ) indices = np . where ( dist < eps ) [ 0 ] nodes = grid . nodes [ key ] [ indices , 0 ] nodes = indices nodes_full = np . hstack ( ( grid . nodes [ key ] [ indices , : ] . squeeze ( ) , np . atleast_2d ( dist [ indices ] ) . T , ) ) if fids is not None : np . savetxt ( fids [ 0 ] , nodes_full . squeeze ( ) ) np . savetxt ( fids [ 1 ] , np . atleast_2d ( dist ) . T , ) elm_indices = [ ] elm_nodes = [ ] for elmnr , element in enumerate ( grid . elements ) : if len ( np . intersect1d ( element , nodes ) ) == 2 : elm_indices . append ( elmnr ) elm_nodes . append ( element ) elms = np . array ( elm_nodes ) neighbors = [ ] eta = settings [ 'eta' ] for index , elm in enumerate ( elms ) : found_it = False for index1 , elm1 in enumerate ( elms ) : ints = np . intersect1d ( elm , elm1 ) if len ( ints ) == 2 : if np . all ( dist [ ints ] < eps ) : found_it = True break else : pass elif len ( ints ) > 2 : pass if found_it : nb = ( elm_indices [ index ] + 1 , elm_indices [ index1 ] + 1 , eta ) if True : neighbors . append ( nb ) else : pass print ( 'No neighbors found. Strange...' ) return np . array ( neighbors )
Compute the element pairs for regulazisation decoupling, given a grid
    object and a line, denoted by start and end coordinates, in the grid.

    Parameters
    ----------
    grid: a crtomo.grid.crt_grid object
    line: np.array/list with 4 entries: x1, y1, x2, x2, denoting start and end
          point
    fids: None, or: list with file ids for
        debug_nodes



    Returns
    -------
    neighbors: Mx2 array with the numbers of adjoining elements along the given
               line