The loop I can replace by numpy operations:

>>> v_array
array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
>>> n_array
array([[ 0.1, 0.2, 0.3],
[ 0.4, 0.5, 0.6]])
>>> f_list
array([[0, 1, 2, 0],
[2, 1, 0, 1]])

Retrieving the v1 vectors:

>>> v1s = v_array[f_list[:, 0]]
>>> v1s
array([[1, 2, 3],
[7, 8, 9]])

Retrieving the normal vectors:

>>> ns = n_array[f_list[:, 3]]
>>> ns
array([[ 0.1, 0.2, 0.3],
[ 0.4, 0.5, 0.6]])

I don't think you understand quite what I'm looking for here. Every vec4 in f_list describes a triangle. The first, second, and third are indices of vertices in v_array. The fourth is an index of n_array. From your code, I've learned that I can do this, which is more what I want:
v1s = v_array[f_list[:,0:3]]
ns = n_array[f_list[:,3]]
Obviously, this changes the arrays quite a bit. With changes, the code doesn't actually crash until the corners = line. Do the following lines that do the dot product and comparison still behave correctly? They should be finding, for each triangle, whether or not the associated normal is less than 90 degrees. The triangle's edges should then be added to an array. See end.

Now how to calculate the pairwise dot product (I suppress the
difference of v1 to some_point for now):

>>> inner = numpy.inner(ns, v1s)
>>> inner
array([[ 1.4, 5. ],
[ 3.2, 12.2]])

This calculates *all* pairwise dot products, we have to select the
diagonal of this square ndarray:

>>> dotprods = inner[[numpy.arange(0, 2), numpy.arange(0, 2)]]
>>> dotprods
array([ 1.4, 12.2])

Now we can create a boolean array saying where the dotprod is > 0
(i.e, angle < 90°), and select those triangles:

>>> select = dotprods > 0
>>> select
array([ True, True], dtype=bool)
>>> selected = f_list[select]
>>> selected
array([[0, 1, 2, 0],
[2, 1, 0, 1]])

This seems like a clever idea.

In this case it's the full list. Now build the triangles corners array:

>>> corners = v_array[selected[:, :3]]
>>> corners
array([[[1, 2, 3],
[4, 5, 6],
[7, 8, 9]],

[[7, 8, 9],
[4, 5, 6],
[1, 2, 3]]])
>>>

This has indices [triangle, vertex number (0, 1, or 2), xyz].
And compute the edges (I think you can make use of them):

>>> edges_dist = numpy.asarray([corners[:, 1] - corners[:, 0], corners[:, 2] - corners[:, 0], corners[:, 2] - corners[:, 1]])
>>> edges_dist
array([[[ 3, 3, 3],
[-3, -3, -3]],

[[ 6, 6, 6],
[-6, -6, -6]],

[[ 3, 3, 3],
[-3, -3, -3]]])

This has indices [corner number, triangle, xyz].
I think it's easier to compare then "reversed" edges, because then
edge[i, j] == -edge[k, l]?

But of course:

>>> edges = numpy.asarray([[corners[:, 0], corners[:, 1]], [corners[:, 1], corners[:, 2]], [corners[:, 2], corners[:, 0]]])
>>> edges
array([[[[1, 2, 3],
[7, 8, 9]],

[[4, 5, 6],
[4, 5, 6]]],

[[[4, 5, 6],

[4, 5, 6]],

[[7, 8, 9],

[1, 2, 3]]],

[[[7, 8, 9],
[1, 2, 3]],

[[1, 2, 3],
[7, 8, 9]]]])

This has indices [edge number (0, 1, or 2), corner number in edge (0
or 1), triangle].
But this may not be what you want (not flattened in triangle number).
Therefore:

>>> edges0 = numpy.asarray([corners[:, 0], corners[:, 1]])
>>> edges1 = numpy.asarray([corners[:, 1], corners[:, 2]])
>>> edges2 = numpy.asarray([corners[:, 2], corners[:, 0]])
>>> edges0
array([[[1, 2, 3],
[7, 8, 9]],

[[4, 5, 6],
[4, 5, 6]]])
>>> edges1
array([[[4, 5, 6],

[4, 5, 6]],

[[7, 8, 9],

[1, 2, 3]]])
>>> edges2
array([[[7, 8, 9],
[1, 2, 3]],

[[1, 2, 3],
[7, 8, 9]]])

>>> edges = numpy.concatenate((edges0, edges1, edges2), axis = 0)
>>> edges
array([[[1, 2, 3],
[7, 8, 9]],

[[4, 5, 6],
[4, 5, 6]],

[[4, 5, 6],

[4, 5, 6]],

[[7, 8, 9],

[1, 2, 3]],

[[7, 8, 9],
[1, 2, 3]],

[[1, 2, 3],
[7, 8, 9]]])

This should be as intended.
The indices are [flat edge number, edge side (left or right), xyz].

Now I guess you have to iterate over all pairs of them, don't know a
numpy accelerated method. Maybe it's even faster to draw the edges
twice than to invest O(N_edges ** 2) complexity for comparing?

Unfortunately, no. The whole point of the algorithm is to extrude back-facing triangles (those with normals facing away from the light) backward, leaving polygons behind in a light volume. Although a shadow volume tutorial can explain this in a more detailed way, basically, for every triangle, if it is back-facing, add its edges to a list. Remove the duplicate edges. So, the edge between two back-facing triangles-that-meet-along-an-edge is not kept. However, if a back-facing triangle and a front-facing triangle meet, only the back-facing triangle's edge is present in the list, and so it is not removed. Thus, only the border edges between the front-facing triangles and the back facing triangles remain in the list. As front-facing triangles face the light and back-facing triangles don't, a silhouette edge is built up. When these edges are extruded, they're extremely useful. The following image http://www.gamedev.net/reference/articles/1873/image010.gif shows the process. The four triangles are all back facing. If duplicate edges are removed, only edges I0-I2, I2-I4, I4-I3, and I3-I0 remain--the silhouette edges.

Still to do, remove the duplicate edges (actually where a good deal of the optimization lies too).

So, for every back-facing triangle [v1,v2,v3,n], (where vn is a vertex index), the edges [v1,v2], [v2,v3], and [v3,v1] need to be added to a list. I.e., f_list needs to be converted into a list of edges in this way. Then, duplicate edge pairs need to be removed, noting that [v1,v2] and [v2,v1] are still a pair (in my Python code, I simply sorted the edges before removing duplicates: [123,4] -> [4,123] and [56,968] -> [56,968]). The final edge list then should be converted back into actual vertices by indexing it into v_array (I think I understand how to do this now!): [ [1,2], [4,6], ... ] -> [ [[0.3,1.6,4.5],[9.1,4.7,7.7]], [[0.4,5.5,8.3],[9.6,8.1,0.3]], ... ]

It may seem a bit complicated, but I hope this impression is mainly
because of the many outputs ...

I double-checked everything, *hope* everything is correct.

So far from me,

Friedrich

Once again, thanks so much,
Ian