def collision_response_self(pair, liver, tool):
    new_vtx_delta = None
    move_tris = {}
    nv_aves = {}
    new_vtx_deltas = {}
    for (key, value) in pair.items():
        new_vtx_delta = np.zeros(liver.x.shape)
        (i_tet, p_index) = (int(key), np.array(value))
        p_n = p_index.shape[0]
        col_tri_index = np.argwhere((liver.tri_tet[:, 0] == i_tet)).flatten()
        if (col_tri_index.size == 0):
            raise ValueError('Update time step too big, vertices skip the surface tetrahedron elements')
        col_tri_n = col_tri_index.shape[0]
        col_tri_nv = liver.tri_normal_vec[col_tri_index]
        col_tri_p = liver.x[liver.tri_elements[col_tri_index].T[0]]
        nv_ave = tool.vtx_normal_vec[p_index].sum(axis=0)
        nv_ave = (nv_ave / np.linalg.norm(nv_ave))
        nv_aves[key] = nv_ave
        dn = nv_ave.dot(col_tri_nv.T)
        ap = (col_tri_p[np.newaxis, :] - tool.x[(p_index, np.newaxis)])
        dotn = np.tile(col_tri_nv, p_n).reshape((- 1), p_n, 3).transpose(1, 0, 2)
        apn = (ap * dotn).sum(axis=(- 1))
        ts = (apn * (1 / dn))
        int_col_p = ((ts[:, :, np.newaxis] * nv_ave[np.newaxis, np.newaxis, :]) + tool.vertices[p_index][:, np.newaxis, :])
        tri_vertices = liver.x[liver.tri_elements[col_tri_index]]
        v1 = (tri_vertices[:, 1] - tri_vertices[:, 0])
        v2 = (tri_vertices[:, 2] - tri_vertices[:, 0])
        tri_areax2 = np.linalg.norm(np.cross(v1, v2, axis=(- 1)), axis=(- 1))
        bc_temp = np.zeros((p_n, col_tri_n, 3, 3, 3))
        bc_temp[:] = np.tile(tri_vertices, 3).reshape((- 1), 3, 3, 3).transpose(0, 2, 1, 3)
        for itemp in range(p_n):
            bc_temp[itemp, :, [0, 1, 2], [0, 1, 2]] = int_col_p[itemp]
        v1 = (bc_temp[:, :, :, 1] - bc_temp[:, :, :, 0])
        v2 = (bc_temp[:, :, :, 2] - bc_temp[:, :, :, 0])
        areax2 = np.linalg.norm(np.cross(v1, v2, axis=(- 1)), axis=(- 1))
        bc_co = (areax2 * (1.0 / tri_areax2)[np.newaxis, :, np.newaxis])
        check1 = np.argwhere((abs((bc_co.sum(axis=(- 1)) - 1)) < 0.001))
        check2 = np.argwhere((dn < 0))
        inter_tri_index = np.intersect1d(check1[:, 1], check2)
        if (inter_tri_index.size == 0):
            the_best_tri = dn.argmin()
            move_tri = liver.tri_elements[col_tri_index[the_best_tri]]
            tri_nv = liver.tri_normal_vec[col_tri_index[the_best_tri]].flatten()
            tri_vtx = liver.x[move_tri].reshape(3, 3)
            v = (nv_ave - tri_nv)
            v = (v / np.linalg.norm(v))
            dn_t = v.dot(tri_nv)
            ap_t = (tri_vtx[0] - tool.x[p_index])
            t_t = (ap_t.dot(tri_nv) / dn_t)
            move_t = t_t.min()
            new_vtx_delta[move_tri] += ((- move_t) * v)
            new_vtx_deltas.setdefault(key, []).append(new_vtx_delta)
            move_tris.setdefault(key, []).append(move_tri.flatten())
            print(' None ', end='')
        elif (len(inter_tri_index) > 1):
            temp_delta = np.zeros((liver.x.shape[0], len(inter_tri_index)))
            itemp = 0
            for inter_tri_i in inter_tri_index:
                part_p_index = check1[((check1[:, 1] == inter_tri_i), 0)]
                move_t = ts[(part_p_index, inter_tri_i)].min()
                move_tri = liver.tri_elements[col_tri_index[inter_tri_i]]
                temp_delta[(move_tri, itemp)] = (- move_t)
                move_tris.setdefault(key, []).append(move_tri.flatten())
                itemp += 1
            new_vtx_delta += (temp_delta.max(axis=(- 1))[:, np.newaxis] * nv_ave[np.newaxis, :])
            new_vtx_deltas.setdefault(key, []).append(new_vtx_delta)
            print(' Multi ', end='')
        else:
            move_t = ts[:, inter_tri_index].min()
            move_tri = liver.tri_elements[col_tri_index[inter_tri_index]]
            new_vtx_delta[move_tri] += ((- move_t) * nv_ave)
            new_vtx_deltas.setdefault(key, []).append(new_vtx_delta)
            move_tris.setdefault(key, []).append(move_tri.flatten())
            print(' Single ', end='')
    return (new_vtx_delta, move_tris, nv_aves, new_vtx_deltas)
def collision_response_self(pair, liver, tool):
    new_vtx_delta = None
    move_tris = {}
    nv_aves = {}
    new_vtx_deltas = {}
    for (key, value) in pair.items():
        new_vtx_delta = np.zeros(liver.x.shape)
        (i_tet, p_index) = (int(key), np.array(value))
        p_n = p_index.shape[0]
        col_tri_index = np.argwhere((liver.tri_tet[:, 0] == i_tet)).flatten()
        if (col_tri_index.size == 0):
            raise ValueError('Update time step too big, vertices skip the surface tetrahedron elements')
        col_tri_n = col_tri_index.shape[0]
        col_tri_nv = liver.tri_normal_vec[col_tri_index]
        col_tri_p = liver.x[liver.tri_elements[col_tri_index].T[0]]
        nv_ave = tool.vtx_normal_vec[p_index].sum(axis=0)
        nv_ave = (nv_ave / np.linalg.norm(nv_ave))
        nv_aves[key] = nv_ave
        dn = nv_ave.dot(col_tri_nv.T)
        ap = (col_tri_p[np.newaxis, :] - tool.x[(p_index, np.newaxis)])
        dotn = np.tile(col_tri_nv, p_n).reshape((- 1), p_n, 3).transpose(1, 0, 2)
        apn = (ap * dotn).sum(axis=(- 1))
        ts = (apn * (1 / dn))
        int_col_p = ((ts[:, :, np.newaxis] * nv_ave[np.newaxis, np.newaxis, :]) + tool.vertices[p_index][:, np.newaxis, :])
        tri_vertices = liver.x[liver.tri_elements[col_tri_index]]
        v1 = (tri_vertices[:, 1] - tri_vertices[:, 0])
        v2 = (tri_vertices[:, 2] - tri_vertices[:, 0])
        tri_areax2 = np.linalg.norm(np.cross(v1, v2, axis=(- 1)), axis=(- 1))
        bc_temp = np.zeros((p_n, col_tri_n, 3, 3, 3))
        bc_temp[:] = np.tile(tri_vertices, 3).reshape((- 1), 3, 3, 3).transpose(0, 2, 1, 3)
        for itemp in range(p_n):
            bc_temp[itemp, :, [0, 1, 2], [0, 1, 2]] = int_col_p[itemp]
        v1 = (bc_temp[:, :, :, 1] - bc_temp[:, :, :, 0])
        v2 = (bc_temp[:, :, :, 2] - bc_temp[:, :, :, 0])
        areax2 = np.linalg.norm(np.cross(v1, v2, axis=(- 1)), axis=(- 1))
        bc_co = (areax2 * (1.0 / tri_areax2)[np.newaxis, :, np.newaxis])
        check1 = np.argwhere((abs((bc_co.sum(axis=(- 1)) - 1)) < 0.001))
        check2 = np.argwhere((dn < 0))
        inter_tri_index = np.intersect1d(check1[:, 1], check2)
        if (inter_tri_index.size == 0):
            the_best_tri = dn.argmin()
            move_tri = liver.tri_elements[col_tri_index[the_best_tri]]
            tri_nv = liver.tri_normal_vec[col_tri_index[the_best_tri]].flatten()
            tri_vtx = liver.x[move_tri].reshape(3, 3)
            v = (nv_ave - tri_nv)
            v = (v / np.linalg.norm(v))
            dn_t = v.dot(tri_nv)
            ap_t = (tri_vtx[0] - tool.x[p_index])
            t_t = (ap_t.dot(tri_nv) / dn_t)
            move_t = t_t.min()
            new_vtx_delta[move_tri] += ((- move_t) * v)
            new_vtx_deltas.setdefault(key, []).append(new_vtx_delta)
            move_tris.setdefault(key, []).append(move_tri.flatten())
            print(' None ', end=)
        elif (len(inter_tri_index) > 1):
            temp_delta = np.zeros((liver.x.shape[0], len(inter_tri_index)))
            itemp = 0
            for inter_tri_i in inter_tri_index:
                part_p_index = check1[((check1[:, 1] == inter_tri_i), 0)]
                move_t = ts[(part_p_index, inter_tri_i)].min()
                move_tri = liver.tri_elements[col_tri_index[inter_tri_i]]
                temp_delta[(move_tri, itemp)] = (- move_t)
                move_tris.setdefault(key, []).append(move_tri.flatten())
                itemp += 1
            new_vtx_delta += (temp_delta.max(axis=(- 1))[:, np.newaxis] * nv_ave[np.newaxis, :])
            new_vtx_deltas.setdefault(key, []).append(new_vtx_delta)
            print(' Multi ', end=)
        else:
            move_t = ts[:, inter_tri_index].min()
            move_tri = liver.tri_elements[col_tri_index[inter_tri_index]]
            new_vtx_delta[move_tri] += ((- move_t) * nv_ave)
            new_vtx_deltas.setdefault(key, []).append(new_vtx_delta)
            move_tris.setdefault(key, []).append(move_tri.flatten())
            print(' Single ', end=)
    return (new_vtx_delta, move_tris, nv_aves, new_vtx_deltas)