Let's first address the side note. "3D length" means the actual length of the path on the earth's surface represented by the polyline, accounting for the additional length contributed by its motion up and down. "Geodesic length" usually refers to the length of a path within a Riemannian manifold. There are at least two useful and reasonable ways in which the earth is viewed as such a manifold. The more conventional one uses the ellipsoidal model. In this model, elevation is ignored, "straight" means "moving along the ellipsoid without any sideways acceleration," and distance is computed by dividing the path into a sequence of very small straight segments. From this point of view 3D distance and geodesic distance usually differ.
A less conventional viewpoint is that the earth's topographic surface is the manifold. Now the 3D distance and geodesic distance almost coincide. A tiny quibble is that lines that appear to be straight on the map (or even on the ellipsoid) usually will not be straight on the earth's surface. This leads to such a small error, though, that for most purposes you can safely consider the 3D length to be geodesic in this sense.
Now, about computing that 3D length. When you account for elevation, each segment of the polyline is really (assumed to be) the hypotenuse of a triangle. The length on the map, dx, is the length of the triangle's base; the difference in elevations between endpoints, dz, is the length of the triangle's height (its "rise"); and the 3D length is the length of the hypotenuse, l. The Pythagorean Theorem tells us that
(1) l = Sqrt( dx^2 + dz^2 ).
Therefore, by breaking any polyline into its constituent segments and computing l for each segment, you can obtain the 3D length by summing the results.
dx is just the length (on the map) of the polyline segment; there are many ways to compute it in ArcGIS. Pick your favorite.
Computing dz accurately can be done in several ways. The best perhaps is to find the coordinates of the endpoints of each segment and use them to extract the DEM's elevation. This would require a bit of scripting in ArcGIS. A more straightforward method (for non self intersecting polylines) is to obtain a zonal range of the DEM using the segments as the zones. In principle, the segments should be so short that when lifted to the earth's surface they still are fairly straight. Thus, the zonal range should be the difference in elevations at the endpoints. There's a bit of bias in this approach, though, because any two successive segments share a common node. The elevation of that node will be assigned to the range of one or the other, but not both, of the segments attached to that node. This systematically shrinks the estimated ranges. You could correct for this by multiplying the zonal range (let's call it r) by (l + c)/l where c is the cellsize; that is,
(2) dz = r (l + c) / l.
This extrapolates the zonal range to the full range of each segment: its rise.
The mechanics of this operation in ArcGIS are:
Break the polyline into individual segments, each with a unique ID.
Obtain the zonal range of the DEM based on the segments layer.
Join the ranges to the segment layer attribute table. Let's call this field [r].
Calculate the length [l] of each segment.
Calculate the rise via equation (2) above.
Calculate the 3D length via equation (1) above.
Sum the 3D lengths.
Perhaps this is a little less "cumbersome and complex" than the model currently used.
I should mention there are other approaches. One of them involves computing the bearing of the polyline at each segment. This can be combined with the slope and aspect grids (derived directly from the DEM) to compute the 3D length. What is noteworthy about this strategy is that you can avoid dissecting the polyline into short segments and that all the operations are performed entirely on grids. That would make it a good solution in raster-only GIS settings.