This is a very natural question, but unfortunately at the moment the answer is a bit complicated. This answer assumes sage 9.1+:

As this is a combinatorial property, you can go into `CombinatorialPolyhedron`

, which is now a method of a polyhedron. There you can obtain the answer via the face iterator. The face iterator iterates through all faces of the polyhedron with a depth-first search. Of course you could just iterate until you what you are looking for, but the iterator can also be manipulated. For this I need to explain a bit, how it works. Assume for now that the iterator is in non-dual mode, it generates all faces from the facets.

The iterator yields first all facets. Then it inductively yields all faces of the first facet (treating it as it's own polyhedron) and then marks it as being visited. Then it yields all faces of the second facet, but the faces contained in the first facet. After that it marks the second facet as being visited and so on.

Now in our situation it might happen, that the first facet does not contain the goal face. Then we want to mark it immediately as being visited. Likewise we proceed with all facets. The below code will proceed like this for the ridges and so on, which might be easier to adapt, when you need different properties.

You might observe that in this particular situation it suffices to mark the correct facets as visited.

Now, if your polyhedron contains many facets but few vertices, generating all faces from the facets might be a bad idea. This is why the iterator can work in dual mode as well. In this case you can decide to skip supfaces instead of subfaces.

As mentioned before, this is a natural question and eventually I think there should be a method like find-join-of-vertices or find-meet-of-facets.

The following should work (it assumes that the face really does exist, otherwise you need to alter the termination condition of the while loop to fit your need):

```
from sage.geometry.polyhedron.face import combinatorial_face_to_polyhedral_face
def is_subset(a,b):
return all(x in b for x in a)
def obtain_face(P, *indices):
indices = tuple(indices)
C = P.combinatorial_polyhedron()
it = C.face_iter()
# Let f be the first facet in non-dual mode or the first vertex in dual mode.
f = next(it)
while f.ambient_V_indices() != indices:
if not it.dual:
# Skipping faces unless our goal face is a subface.
if not is_subset(indices, f.ambient_V_indices()):
it.ignore_subfaces()
else: # the face iterator starts with the vertices
# Skipping faces unless our goal face is a supface.
if not is_subset(f.ambient_V_indices(), indices):
it.ignore_supfaces()
f = next(it)
# f is now the combinatorial face we where looking for (assuming it exists).
# Obtain a polyhedral face from it.
return combinatorial_face_to_polyhedral_face(P, f)
```

Then you get:

```
sage: P = polytopes.cross_polytope(4)
sage: f = obtain_face(P, 2,4,6)
sage: f.ambient_V_indices()
(2, 4, 6)
sage: f
A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 3 vertices
sage: P = polytopes.hypercube(3)
sage: f = obtain_face(P, 0,1,5,6)
sage: f.ambient_V_indices()
(0, 1, 5, 6)
sage: f
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices
```

This should be pretty fast. Note, that no options of the face iterator are exposed in the polyhedron and one has to go into combinatorial polyhedron (skipping subfaces or supfaces is crucial, because we only want one depth search and find our face immediately).

The the subset check function is not optimized deliberately.

Please provide a polyhedron

`p`

so that other users have a starting point to try and answer your question.@slelievre Thanks for the suggestion. I've updated my question with an example.