I am also interested in knowing if there are any good reasons to pick one ordering of polygons over some other ordering.
I'll post what little I know, but I'm really hoping someone else will post something much better.
I have heard of 3 "standards" for ordering these polygons, but as far as I can tell they were arbitrarily picked:
Your icosahedron-based method,
the quadrilateralized spherical cube, and
Tegmark's icosahedron-based method for pixelizing the celestial sphere.
All three are examples of geodesic grids.
Ordering polygons on a polyhedron
icosahedron-based method for pixelizing the celestial sphere
"What is the best way to pixelize a sphere?" by Max Tegmark
gives a format that appears extremely similar to your proposal:
Tegmark starts with an icosahedron, just as you do,
and recursively subdivides it, just as you do.
So I find it hard to believe that his method is "computationally and geometrically far more complex than what [you] are proposing" or that it has significantly "greater distortion in either area or geometry".
Given some specified resolution (the same as your depth of recursion)
and some point on a sphere, Tegmark's system
assigns that point to a "pixel number" (the number indicating what major icosahedral face it is on, and which particular minor sub-triangle on that face that point is in).
quadrilateralized spherical cube
Wikipedia: "quadrilateralized spherical cube"
The all-sky, skyward looking, unfolded cube is reassembled by
arranging the faces in a sideways T with the bar on the right side:
0
4 3 2 1
5
where square face 0 is centered on the North pole,
square face 5 is centered on the South pole,
the vernal equinox -- at latitude=0, longitude=0 --
lies at the center of square face 1.
Ecliptic longitude increases from face 1 to face 4.
other polyhedra
Tegmark mentions that "it is clearly desirable to use as small faces as possible".
Subdividing each face of a (30-faced) rhombic triacontahedron into two nearly-equilateral isosceles triangles gives the (60-faced) pentakis dodecahedron, as you mentioned,
which one might expect to give less error than starting with a (20-faced) icosahedron.
Further subdividing each of those faces into two irregular triangles gives the (120-faced) disdyakis triacontahedron (dual to the truncated icosidodecahedron).
I've been told that the strictly convex polyhedron with the maximum number of identical faces (exactly equal-shaped and exactly equal-area) is that disdyakis triacontahedron.
However, apparently none of these shapes gives any improvement over your icosahedron proposal.
Tegmark's code immediately breaks up each equilateral triangle of the 20-sided icosahedron into 6 identical irregular triangles.
The resulting 20*6 = 120 identical triangles are the same triangles as the faces of the disdyakis triacontahedron, right?
Ordering polygons on a recursively subdivided polygon
The quadrilateralized spherical cube recursively subdivides each square -- each division of 4 by area adding 2 more bits of localization information -- using a Morton code (Z-order curve).
The Maidenhead Locator System uses a system similar to the Morton code.
You might also want to look at the military grid reference system.
Other ways of ordering such subdivided square (i.e., drawing a path that hits the center of each square exactly once) include the
Sierpiński curve, the Hilbert curve,
the Peano curve,
the Moore curve, etc.
There exist other space-filling curves based on recursively subdividing a triangle -- each division of 4 by area adding 2 more bits of localization area.
(Do any of them have popular names?)
There are 2 popular ways of subdividing squares and triangles:
( a, b )
"Class I"/"alternate" method:
each of the original triangles (or squares) is replaced with n^2 triangles (or squares) that exactly fit inside the original triangle (or square).
(n=1 gives the original base shape, n=2 gives a shape with 4 times as many faces).
The small triangles have edges (very close to) parallel to the original big triangles.
"Class II"/"triacon" method, employed on almost every major geodesic dome project in the 1950s ( a p. 30 ),
each of the original triangles is replaced with 3*n^2 triangles.
The small triangles have edges (very close to) right angles to the original big triangles.
(Unreliable sources tell me the triacon method can be applied to the cube or rhombic solids,
each of the original squares or rhombuses is replaced with 2*n^2 rhombuses).
(n=1 converts the icosahedron to the pentakis dodecahedron I mentioned earlier; n=8 converts the icosahedron to an approximation of the Epcot "Spaceship Earth"; n=1 converts the cube into the rhombic dodecahedron, etc.).
EDIT:
The most common application of numbering the sides of a platonic solid is while making dice.
One particular way of putting numbers of the sides of an icosahedron is shown on DiceCollector's page, but I don't know if that is particularly "standard".
Have you considered asking
"What is the proper way to number the sides of a 1d20 die?" or
"What is the proper way to number the sides of a 1d120 die?",
or both,
on the https://rpg.stackexchange.com/ ?
I find your use of the term "orientation" slightly confusing, since it seems to oscillate between an absolute sense (roughly synonymous with "direction") and a relative sense (roughly synonymous with "alignment"). Let me paraphrase how I understand your question:
You have triangles for more than one cylinder in a file, and you don't know which triangle belongs to which cylinder. The normals of the triangles are all roughly orthogonal to the axis of their cylinder. For a given direction, you want to know how well the axes of the cylinders are aligned with this direction. In an extreme case, the axes of all cylinders might have the same direction; in that case, the test with that direction should yield "fully aligned" and a test with any direction orthogonal to it should yield "fully non-aligned". The problem is that though the cylinders may be fully aligned, the normals of the triangles form different angles with any given direction.
It seems to me the "incorrect way" that you're using isn't quite so bad. I don't fully understand how you arrive at $0$ and $\pi/4$, respectively, but an indicator that ranges between $0$ for aligned and $\pi/4$ for non-aligned seems useful. I gather that what you don't like about it is that you want it to be some sort of average angle.
You could probably find a function that maps your current indicator to something like an average angle, but generally speaking, averaging trigonometric functions tends to be easier and to have nicer properties than averaging angles (or absolute values of angles, which seems to be what you're doing), so I'll do something similar to what you did, but using trigonometric functions whose averages we can easily calculate exactly, so we'll get an expression in closed form relating the values of the indicator to something like an average angle.
For this to work, the triangles for each cylinder need to be roughly equidistributed about its axis. If the normals of the triangles tend to preferentially point in one of the directions orthogonal to the axis, that will make the cylinder appear more aligned with the third direction (orthogonal to the axis and the preferred normal direction) than it actually is.
So consider a cylinder with unit vector $\vec a$ along its axis, a unit vector $\vec u$ along which we want to measure the cylinder's alignment, and unit vectors $\vec n_1$ and $\vec n_2$ orthogonal to $\vec a$, which we choose such that $\vec n_1$ is also orthogonal to $\vec u$ (which we can always do). We can write the normal of a triangle of the cylinder as $\vec n=\cos\phi\,\vec n_1+\sin\phi\,\vec n_2$, and equidistribution of the normals around the cylinder's axis amounts to equidistribution of $\phi$ over $[0,2\pi]$. We have $\vec u\cdot\vec n_1=0$ (since we chose $\vec n_1$ orthogonal to $\vec u$) and $\vec u\cdot\vec n_2=\pm\sin\theta$, were $\theta$ is the angle between $\vec a$ and $\vec u$ (you may want to draw a sketch to see why that's the case). Thus $\vec u\cdot\vec n=\pm\sin\theta\sin\phi$. Now if we average this over all the equidistributed normals, we'll just get $0$ from the average over $\phi$, independent of $\theta$; that's not what we want. However, if we square before averaging, we get the average of $\sin^2\theta\sin^2\phi$ over $\phi$, which is $\frac12\sin^2\theta$. That's good, because we shouldn't be distinguishing between negative and positive values of $\theta$ or between $\theta$ and $\pi-\theta$ anyway.
Now to turn this into something like an average angle, we just have to apply the inverse operations. That is, if $\langle(\vec u\cdot\vec n)^2\rangle$ denotes the average of the square of the scalar product of the alignment direction $\vec u$ with the triangle normals $\vec n$, the "average angle" $\hat\theta$ would be derived from $\langle(\vec u\cdot\vec n)^2\rangle=\langle\frac12\sin^2\theta\rangle$ as
$$
\hat\theta=\arcsin\sqrt{2\langle(\vec u\cdot\vec n)^2\rangle}\;.
$$
You'll need to clamp any arguments above $1$ for the arcsine to $1$; this might occur if the normals aren't perfectly equidistributed.
If all the cylinders are perfectly aligned along $\vec u$, then all the normals will be orthogonal to $\vec u$, so $\langle(\vec u\cdot\vec n)^2\rangle$ will be zero, and so will $\hat\theta$. On the other hand, if all the cylinders are perfectly aligned orthogonal to $\vec u$ and their normals are equidistributed about their axes, then $\langle(\vec u\cdot\vec n)^2\rangle$ will be $1/2$, so $\hat\theta$ will be $\pi/2$ as desired. More generally, if all the cylinder axes form an angle $\theta$ with $\vec u$, then $\hat\theta$ will be $\theta$. If the cylinders form different angles with $\vec u$, then $\hat\theta$ will be an average angle in the sense that $\sin^2\theta$ (or equivalently $\cos^2\theta$) is averaged.
Best Answer
Taking the absolute values is probably not a good idea, and the $x\to1-x$ transform at the end also doesn't look promising.
I'll assume that by "least aligned" you mean that ideally your vector would be orthogonal to all the normals (and thus in the planes of all the triangles), not that it should be opposite to the normals (and thus pointing into the planes); in the latter case you should just add up all the normals and use the direction opposite to the sum.
A first attempt might be to add up the normals and then take a vector orthogonal to the sum, but then you'd have to arbitrarily choose one of the vectors in the plane orthogonal to the sum, and any information which of these might be best would be lost.
Here's a systematic way of doing this. You want your vector $x$ to be orthogonal to the normals $n_i$, so you want the absolute value of the dot product $x\cdot n_i$ to be as small as possible. A natural ansatz would then be to minimize the sum of the squares of the dot products:
$$ \sum_i(x\cdot n_i)^2\to\min\;, $$
subject to the constraint that $x$ is a unit vector. We can write this more suggestively:
$$ \sum_ix^\top n_in_i^\top x=x^\top\left(\sum_i n_in_i^\top \right)x\to\min\;. $$
Adding a Lagrange multiplier term for the normalization constraint yields the objective function
$$ x^\top\left(\sum_i n_in_i^\top \right)x-\lambda x^\top x\;, $$
and requiring the gradient with respect to $x$ to vanish yields
$$ \left(\sum_i n_in_i^\top \right)x-\lambda x=0\;. $$
This is the eigenvalue problem for the positive semidefinite matrix $\sum_in_in_i^\top$. That's just a $3\times3$ matrix, so you can readily find its eigensystem. The direction you want is the eigenvector corresponding to its least eigenvalue.