def _ResolutionOrder(self, variables_to_solve):
        """
        return a list of lists of tuples (block,output,ndof) to be solved

        """
#    Gp=nx.DiGraph()
#
#    for i in range(nvar):
#        Gp.add_node('v'+str(i),bipartite=0)
#
#    for i in range(neq):
#        Gp.add_node('e'+str(i),bipartite=1)
#        for j in range(nvar):
#            if Mo[i,j]==1:
#                Gp.add_edge('e'+str(i),'v'+str(j))

        Gp = nx.DiGraph()
        for variable in self.variables:
            Gp.add_node(variable, bipartite=0)
        for block in self.blocks:
            for iov, output_variable in enumerate(block.outputs):
                Gp.add_node((block, iov), bipartite=1)
                Gp.add_edge((block, iov), output_variable)
                Gp.add_edge(output_variable, (block, iov))
                for input_variable in block.inputs:
                    if not isinstance(input_variable, Signal):
                        Gp.add_edge(input_variable, (block, iov))

    #    for n1,n2 in M.items():
    #        Gp.add_edge(n1,n2)

        sinks = []
        sources = []
        for node in Gp.nodes():
            if Gp.out_degree(node) == 0:
                sinks.append(node)
            elif Gp.in_degree(node) == 0:
                sources.append(node)

        G2 = sources[:]
        for node in sources:
            for node2 in nx.descendants(Gp, node):
                if node2 not in G2:
                    G2.append(node2)

        if G2 != []:
            print(G2)
            raise ModelError('Overconstrained variables')

        G3 = sinks[:]
        for node in sinks:
            for node2 in nx.ancestors(Gp, node):
                if node2 not in G3:
                    G3.append(node2)

        if G3 != []:
            raise ModelError('Underconstrained variables')

#        vars_resolvables=[]
#        for var in vars_resoudre:
#            if not 'v'+str(var) in G2+G3:
#                vars_resolvables.append(var)


#        G1=Gp.copy()
#        G1.remove_nodes_from(G2+G3)
#
#        M1=nx.bipartite.maximum_matching(G1)
#        G1p=nx.DiGraph()
#
#        G1p.add_nodes_from(G1.nodes())
#        for e in G1.edges():
#            # equation vers variable
#            if e[0][0]=='v':
#                G1p.add_edge(e[0],e[1])
#            else:
#                G1p.add_edge(e[1],e[0])
#    #    print(len(M))
#        for n1,n2 in M1.items():
#    #        print(n1,n2)
#            if n1[0]=='e':
#                G1p.add_edge(n1,n2)
#            else:
#                G1p.add_edge(n2,n1)

        scc = list(nx.strongly_connected_components(Gp))
    #    pos=nx.spring_layout(G1p)
    #    plt.figure()
    #    nx.draw(G1p,pos)
    #    nx.draw_networkx_labels(G1p,pos)
#        print(scc)
        if scc != []:
            C = nx.condensation(Gp, scc)
            isc_vars = []
            for isc, sc in enumerate(scc):
                for var in variables_to_solve:
                    if var in sc:
                        isc_vars.append(isc)
                        break
            ancestors_vars = isc_vars[:]

            for isc_var in isc_vars:
                for ancetre in nx.ancestors(C, isc_var):
                    if ancetre not in ancestors_vars:
                        ancestors_vars.append(ancetre)

            order_sc = [sc for sc in nx.topological_sort(
                C) if sc in ancestors_vars]
            order_ev = []
            for isc in order_sc:
                # liste d'équations et de variables triées pour être séparées
                evs = list(scc[isc])
#                print(evs)
#                levs=int(len(evs)/2)
                eqs = []
                var = []
                for element in evs:
                    if type(element) == tuple:
                        eqs.append(element)
                    else:
                        var.append(element)
                order_ev.append((len(eqs), eqs, var))

            return order_ev

        raise ModelError
return a list of lists of tuples (block,output,ndof) to be solved